Mengapa sistem linear yang dikondisikan dengan buruk dapat diselesaikan dengan tepat?

13

Menurut jawaban di sini , angka kondisi besar (untuk penyelesaian sistem linier) mengurangi jumlah digit yang benar dalam solusi floating point. Matriks diferensiasi orde tinggi dalam metode pseudospectral biasanya sangat dikondisikan. Mengapa saat itu mereka masih merupakan metode yang sangat akurat?

Saya memahami bahwa presisi rendah yang berasal dari matriks yang dikondisikan buruk hanya merupakan nilai yang dijamin , tetapi tetap saja membuat saya bertanya-tanya mengapa matriks yang dikondisikan secara akurat diselesaikan dengan metode langsung dalam praktiknya - misalnya, LCOLkolom Tabel 3.1 pada halaman 11 dari Wang et al., METODE KOLOKASI YANG BAIK DENGAN KONDISI MENGGUNAKAN MATRIX INTEGRASI PSEUDOSPECTRAL , SIAM J. Sci. Komputasi, 36 (3) .

Zoltán Csáti
sumber
2
Intuisi saya adalah bahwa solvabilitas / akurasi sistem Ax = b terkait dengan vektor pemaksaan b, bukan hanya matriks A. Mungkin jika b tidak "menyelidiki" atau "membangkitkan" mode A yang dikondisikan dengan buruk, maka solusi yang akurat tetap mungkin. Sebagai contoh pembatas, A bisa saja singular (angka kondisi tak terbatas), namun Ax = b mungkin masih memiliki solusi, yang dapat dihitung secara akurat, jika data pemaksaan b berada di kisaran A. Saya akui ini cukup bagus -wavy, itulah sebabnya saya hanya berkomentar alih-alih menjawab.
rchilton1980
@ rchilton1980 "namun Ax = b mungkin masih memiliki solusi" Tapi solusi itu tidak unik. Dan contoh yang saya maksudkan memiliki solusi unik.
Zoltán Csáti
Itu tandingan yang adil - mungkin artefak memilih nomor kondisi tak terbatas (nilai eigen yang persis nol). Namun saya pikir Anda dapat mengganti nilai eigen nol dengan epsilon mesin dan poin saya masih berdiri. (Yaitu, sistem memiliki angka kondisi yang sangat besar, sistem tidak berbentuk dengan solusi unik, yang dapat kita hitung dengan sangat akurat asalkan tidak memiliki komponen di sepanjang eigenpair kecil itu).
rchilton1980
1
Untuk lebih spesifik, eksperimen pemikiran saya di sini adalah sesuatu seperti A = diag ([1 1 1 1 1 eps]), b = [b1 b2 b3 b4 b5 0]. Ini dibuat-buat, tapi saya pikir cukup untuk membenarkan klaim aslinya: "kadang-kadang kondisi buruk A dapat diselesaikan secara akurat untuk pilihan b"
rchilton1980
1
Hanya memberikan contoh lain dari Moler ini blog blogs.mathworks.com/cleve/2015/02/16/...
percusse

Jawaban:

7

Ditambahkan setelah jawaban awal saya:

Tampak bagi saya sekarang bahwa penulis makalah yang direferensikan memberikan angka kondisi (angka kondisi 2-norma tetapi angka kondisi norma tak terhingga) dalam tabel sambil memberikan kesalahan absolut maksimum daripada kesalahan relatif normal atau kesalahan relatif elemenwise maksimum ( ini semua adalah ukuran yang berbeda.) Perhatikan bahwa kesalahan relatif elemen-elemen maksimum tidak sama dengan kesalahan relatif norma-tak terhingga. Selain itu, kesalahan dalam tabel relatif terhadap solusi yang tepat untuk masalah nilai batas persamaan diferensial asli daripada sistem linear persamaan diskritisasi. Dengan demikian informasi yang disediakan dalam makalah benar-benar tidak sesuai untuk digunakan dengan batas kesalahan berdasarkan nomor kondisi.

Namun, dalam replikasi perhitungan saya, saya memang melihat situasi di mana kesalahan norma tak terhingga relatif (atau kesalahan relatif dua norma) secara substansial lebih kecil daripada terikat yang ditetapkan oleh nomor kondisi tak terhingga-norma (masing-masing nomor kondisi 2-norma.) Terkadang Anda beruntung.

Saya menggunakan paket DMSUITE MATLAB dan memecahkan contoh masalah dari makalah ini menggunakan metode pseudospectral dengan polinomial Chebyshev. Nomor kondisi saya dan kesalahan absolut maksimum mirip dengan yang dilaporkan di koran.

Saya juga melihat kesalahan relatif norma yang agak lebih baik dari yang diperkirakan berdasarkan nomor kondisi. Sebagai contoh, pada contoh masalah dengan , menggunakan N = 1024 , saya dapatkanϵ=0,01N=1024

cond (A, 2) = 7.9e + 8

cond (A, inf) = 7.8e + 8

norm (u-uexact, 2) / norm (uexact, 2) = 3.1e-12

norm (u-uexact, inf) / norm (uexact, inf) = 2.7e-12

Tampaknya solusinya bagus untuk sekitar 11-12 digit, sedangkan nomor kondisinya ada di urutan 1e8.

Namun, situasi dengan kesalahan elemen-elemen lebih menarik.

maks (abs (u-uexact)) = 2.7e-12

Itu masih terlihat bagus.

maks (abs ((u-uexact) ./ uexact) = 6.1e + 9

Wow- ada kesalahan relatif yang sangat besar dalam setidaknya satu komponen solusi.

Apa yang terjadi? Solusi yang tepat dari persamaan ini memiliki komponen yang kecil (mis. 1.9e-22) sedangkan solusi perkiraan keluar pada nilai yang jauh lebih besar dari 9e-14. Ini disembunyikan oleh pengukuran kesalahan relatif norma (apakah itu 2-norma atau tak terhingga-norma) dan hanya menjadi terlihat ketika Anda melihat kesalahan relatif elemen-elemen dan mengambil maksimum.

Jawaban asli saya di bawah ini menjelaskan mengapa Anda bisa mendapatkan kesalahan relatif normal dalam solusi yang kurang dari batas yang diberikan oleh nomor kondisi.


Seperti yang telah Anda catat dalam pertanyaan, nomor kondisi, , dari matriks non-singular memberikan kesalahan relatif kasus terburuk yang terikat untuk solusi untuk sistem persamaan yang terganggu. Yaitu, jika kita menyelesaikan A ( x + Δ x ) = b + Δ b dengan tepat dan menyelesaikan A x = b dengan tepat, makaκ(SEBUAH)SEBUAH(x+Δx)=b+ΔbSEBUAHx=b

Δxxκ(A)Δbb

Nomor kondisi dapat dihitung sehubungan dengan berbagai norma, tetapi nomor kondisi dua norma sering digunakan, dan itulah nomor kondisi yang digunakan dalam kertas yang Anda rujuk.

Kesalahan kasus terburuk terjadi ketika adalah vektor singular kiri A yang sesuai dengan nilai singular terkecil dari A . Kasus terbaik terjadi ketika Δ b adalah vektor singular kiri A yang sesuai dengan nilai singular A terbesar . Ketika Δ b adalah acak, maka Anda harus melihat proyeksi Δ b pada semua vektor singular kiri A dan nilai singular yang sesuai. Bergantung pada spektrum A , segalanya mungkin berjalan sangat buruk atau sangat baik. ΔbAAΔbAAΔbΔbAA

Pertimbangkan dua matriks , keduanya dengan nomor kondisi 2-norma 1,0 × 10 10 . Matriks pertama memiliki nilai singular 1 , 1 × 10 - 10 , , 1 × 10 - 10 . Matriks kedua memiliki nilai singular 1 , 1 , , 1 , 1 × 10 - 10 . A1.0×101011×10101×10101111×1010

Dalam kasus pertama, gangguan acak tidak mungkin berada dalam arah vektor singular kiri pertama, dan lebih cenderung dekat dengan salah satu vektor singular dengan nilai singular . Dengan demikian perubahan relatif dalam solusi cenderung sangat besar. Dalam kasus kedua, hampir semua gangguan akan dekat ke arah vektor singular dengan nilai singular 1 , dan perubahan relatif dalam solusi akan kecil. 1×10101

PS (ditambahkan kemudian setelah saya kembali dari kelas yoga ...)

Rumus untuk solusi untuk adalahAΔx=Δb

Δx=VΣ1UTΔb=i=1nUiTΔbσiVi

Oleh teorema Pythagoras,

Δx22=saya=1n(UsayaTΔbσsaya)2

Jika kita tetap , maka jumlah ini dimaksimalkan ketika Δ b = U n dan diminimalkan ketika Δ b = U 1 .Δb2=1Δb=UnΔb=U1

Dalam situasi yang dipertimbangkan di sini, adalah hasil dari kesalahan pembulatan acak, sehingga nilai U T i Δ b semuanya harus besarnya kira-kira sama. Istilah dengan nilai yang lebih kecil dari σ i akan berkontribusi banyak pada kesalahan, sementara istilah dengan nilai yang lebih besar dari σ saya tidak akan berkontribusi banyak. Bergantung pada spektrumnya, ini bisa dengan mudah jauh lebih kecil dari batas kasus terburuk. ΔbUsayaTΔbσsayaσsaya

Brian Borchers
sumber
Bukankah argumen ini menyiratkan bahwa adalah mungkin (bahkan jika tidak mungkin) untuk mencapai batas kasus terburuk dari untuk matriks dalam contoh ini? AFAIU, berdasarkan jawaban saya dan berdasarkan pada dokumentasi ini seharusnya tidak mungkin. κ(A)?getrs
Kirill
@BrianBorchers Bisakah Anda menjelaskan mengapa "Kesalahan kasus terburuk terjadi ketika adalah vektor singular kiri A yang sesuai dengan nilai singular terkecil dari A. Kasus terbaik terjadi ketika Δ b adalah vektor singular kiri A yang sesuai dengan yang terbesar nilai singular dari A. " memegang? Dari contoh di bawah ini adalah logis, tetapi saya perlu beberapa rumus. Biarkan SVD dari A menjadi A = U Σ V T . Dalam kasus pertama, A = Δ b σ 1 v T 1 +ΔbAAΔbAAAA=UΣVT . Bagaimana untuk melanjutkan? A=Δbσ1v1T+i=2NuiσiviT
Zoltán Csáti
Saya belum pernah membahas kesalahan pembulatan dalam matriks , tetapi efek umumnya serupa - kecuali jika Anda benar-benar beruntung dalam kesalahan pembulatan, Anda biasanya melakukan sedikit lebih baik daripada ikatan kasus terburuk pesimistis. A
Brian Borchers
(-1) Diskusi kesalahan relatif komponen-bijaksana dalam output sangat menyesatkan.
Kirill
1

tl; dr Mereka melaporkan sebuah nomor kondisi, belum tentu tepat jumlah kondisi matriks, karena ada perbedaan.

Ini khusus untuk matriks dan vektor sisi kanan. Jika Anda melihat dokumentasi untuk*getrs , dikatakan batas kesalahan maju adalah Di sinicond(A,x)tidak cukup seperti bilangan kondisiκ(A), melainkan cond(A,x)=| A - 1

x-x0xcHaind(SEBUAH,x)kamucHaind(SEBUAH)kamu.
cHaind(SEBUAH,x)κ(SEBUAH) (Di sini di dalam norma ini adalah nilai absolut komponen-bijaksana.) Lihat, misalnya,perbaikan berulang untuk sistem linier dan LAPACKoleh Higham, atauAkurasi dan Stabilitas Algoritma NumerikHigham(7.2).
cond(A,x)=|A1||A||x|x,cond(A)=|A1||A|.

Sebagai contoh Anda, saya mengambil operator diferensial pseudospectral untuk masalah yang sama dengan , dan sebenarnya ada perbedaan besar antara | A - 1 | | A | Dan κ ( A ) , saya menghitung 7 × 10 3 dan 2,6 × 10 7n=128|A1||A|κ(A)7×1032.6×107, yang cukup untuk menjelaskan pengamatan bahwa ini terjadi untuk semua sisi kanan, karena urutan besarnya kira-kira sesuai dengan yang terlihat pada Tabel 3.1 (3-4 urutan kesalahan yang lebih baik). Ini tidak bekerja ketika saya mencoba sama untuk hanya random matrix sakit-AC, sehingga harus menjadi milik .A

Contoh eksplisit di mana kedua nomor kondisi tidak cocok, yang saya ambil dari Higham (7.17, hal.124), karena Kahan adalah Contoh lain yang saya temukan hanyalah matriks Vandermonde polosdenganbacak. Saya melewatidan beberapa matriks berkondisi buruk lainnya juga menghasilkan jenis hasil ini, sepertidan.

(2111ϵϵ1ϵϵ),(2+2ϵϵϵ).
[1:10]bMatrixDepot.jltriwmoler

Pada dasarnya, apa yang terjadi adalah ketika Anda menganalisis stabilitas penyelesaian sistem linier sehubungan dengan gangguan, Anda harus terlebih dahulu menentukan gangguan mana yang Anda pertimbangkan. Saat memecahkan sistem linier dengan LAPACK, batas kesalahan ini mempertimbangkan gangguan komponen-bijaksana dalam , tetapi tidak ada gangguan di b . Jadi ini berbeda dari yang biasa κ ( A ) = A - 1A , yang menganggap gangguan normwise di kedua A dan b .Abκ(A)=A1AAb

Pertimbangkan (sebagai contoh tandingan) juga apa yang akan terjadi jika Anda tidak membuat perbedaan. Kita tahu bahwa menggunakan perbaikan berulang dengan presisi ganda (lihat tautan di atas) kita bisa mendapatkan kesalahan relatif ke depan yang paling baik dari untuk matriks dengan κ ( A ) 1 / u . Jadi jika kita mempertimbangkan gagasan bahwa sistem linier tidak dapat diselesaikan dengan akurasi lebih baik daripada κ ( A ) u , bagaimana mungkin solusi pemurnian bekerja?O(u)κ(A)1/uκ(A)u

PS Itu penting yang ?getrsmengatakan solusi dihitung adalah solusi yang benar (A + E)x = bdengan perturbasi di A , tetapi tidak ada gangguan di b . Keadaan akan berbeda jika gangguan diizinkan di b .EAbb

Mengedit Untuk menunjukkan ini bekerja lebih langsung, dalam kode, bahwa ini bukan kebetulan atau masalah keberuntungan, melainkan (biasa) konsekuensi dari dua angka kondisi yang sangat berbeda untuk beberapa matriks tertentu, yaitu,

cond(A,x)cond(A)κ(A).
function main2(m=128)
    A = matrixdepot("chebspec", m)^2
    A[1,:] = A[end,:] = 0
    A[1,1] = A[end,end] = 1
    best, worst = Inf, -Inf
    for k=1:2^5
        b = randn(m)
        x = A \ b
        x_exact = Float64.(big.(A) \ big.(b))
        err = norm(x - x_exact, Inf) / norm(x_exact, Inf)
        best, worst = min(best, err), max(worst, err)
    end
    @printf "Best relative error:       %.3e\n" best
    @printf "Worst relative error:      %.3e\n" worst
    @printf "Predicted error κ(A)*ε:    %.3e\n" cond(A, Inf)*eps()
    @printf "Predicted error cond(A)*ε: %.3e\n" norm(abs.(inv(A))*abs.(A), Inf)*eps()
end

julia> main2()
Best relative error:       2.156e-14
Worst relative error:      2.414e-12
Predicted error κ(A)*ε:    8.780e-09
Predicted error cond(A)*ε: 2.482e-12

Sunting 2 Berikut adalah contoh lain dari fenomena yang sama di mana kondisi yang berbeda jumlahnya berbeda-beda. Kali ini, Di sini A adalah matriks Vandermonde 10 × 10 pada 1 : 10 , dan ketika x dipilih secara acak, c o n d ( A , x ) secara nyata lebih kecil dari κ

cond(A,x)cond(A)κ(A).
A1:10xcond(A,x) , dan kasus terburuk x diberikan oleh x i = i a untuk beberapa a .κ(A)xxi=iaa
function main4(m=10)
    A = matrixdepot("vand", m)
    lu = lufact(A)
    lu_big = lufact(big.(A))
    AA = abs.(inv(A))*abs.(A)
    for k=1:12
        # b = randn(m) # good case
        b = (1:m).^(k-1) # worst case
        x, x_exact = lu \ b, lu_big \ big.(b)
        err = norm(x - x_exact, Inf) / norm(x_exact, Inf)
        predicted = norm(AA*abs.(x), Inf)/norm(x, Inf)*eps()
        @printf "relative error[%2d]    = %.3e (predicted cond(A,x)*ε = %.3e)\n" k err predicted
    end
    @printf "predicted κ(A)*ε      = %.3e\n" cond(A)*eps()
    @printf "predicted cond(A)*ε   = %.3e\n" norm(AA, Inf)*eps()
end

Kasus rata-rata (hampir 9 urutan kesalahan lebih besar lebih baik):

julia> T.main4()
relative error[1]     = 6.690e-11 (predicted cond(A,x)*ε = 2.213e-10)
relative error[2]     = 6.202e-11 (predicted cond(A,x)*ε = 2.081e-10)
relative error[3]     = 2.975e-11 (predicted cond(A,x)*ε = 1.113e-10)
relative error[4]     = 1.245e-11 (predicted cond(A,x)*ε = 6.126e-11)
relative error[5]     = 4.820e-12 (predicted cond(A,x)*ε = 3.489e-11)
relative error[6]     = 1.537e-12 (predicted cond(A,x)*ε = 1.729e-11)
relative error[7]     = 4.885e-13 (predicted cond(A,x)*ε = 8.696e-12)
relative error[8]     = 1.565e-13 (predicted cond(A,x)*ε = 4.446e-12)
predicted κ(A)*ε      = 4.677e-04
predicted cond(A)*ε   = 1.483e-05

Kasus terburuk ( ):a=1,,12

julia> T.main4()
relative error[ 1]    = 0.000e+00 (predicted cond(A,x)*ε = 6.608e-13)
relative error[ 2]    = 1.265e-13 (predicted cond(A,x)*ε = 3.382e-12)
relative error[ 3]    = 5.647e-13 (predicted cond(A,x)*ε = 1.887e-11)
relative error[ 4]    = 8.895e-74 (predicted cond(A,x)*ε = 1.127e-10)
relative error[ 5]    = 4.199e-10 (predicted cond(A,x)*ε = 7.111e-10)
relative error[ 6]    = 7.815e-10 (predicted cond(A,x)*ε = 4.703e-09)
relative error[ 7]    = 8.358e-09 (predicted cond(A,x)*ε = 3.239e-08)
relative error[ 8]    = 1.174e-07 (predicted cond(A,x)*ε = 2.310e-07)
relative error[ 9]    = 3.083e-06 (predicted cond(A,x)*ε = 1.700e-06)
relative error[10]    = 1.287e-05 (predicted cond(A,x)*ε = 1.286e-05)
relative error[11]    = 3.760e-10 (predicted cond(A,x)*ε = 1.580e-09)
relative error[12]    = 3.903e-10 (predicted cond(A,x)*ε = 1.406e-09)
predicted κ(A)*ε      = 4.677e-04
predicted cond(A)*ε   = 1.483e-05

A=(010000100001ϵ000).
A=1A1=ϵ1κ(A)=ϵ1|A1|=A1=|A|1cond(A)=1Ax=bκ(A)

cond(A)κ(A)

A = matrixdepot("kahan", 48)
κ, c = cond(A, Inf), norm(abs.(inv(A))*abs.(A), Inf)
@printf "κ=%.3e c=%.3e ratio=%g\n" κ c (c/κ)

κ=8.504e+08 c=4.099e+06 ratio=0.00482027
Kirill
sumber
Nomor kondisi dalam makalah yang disebut oleh OP adalah angka kondisi dua norma. Jika Anda kembali ke referensi [17] oleh ElBarbary Anda akan melihat bahwa pada makalah sebelumnya ini adalah nomor kondisi dua norma. Juga, saya menyiapkan contoh dari makalah ini menggunakan DMsuite, dan mendapatkan angka kondisi 2-norma yang hampir sama persis seperti yang dilaporkan dalam makalah ini.
Brian Borchers
Nomor norma kondisi norma tak terbatas untuk contoh-contoh ini yang saya dapatkan menggunakan interpolasi dmsuite dan Chebyshev serupa dalam besarnya dengan angka kondisi dua norma. Saya tidak berpikir bahwa perbedaan antara 2-norma dalam angka kondisi infinity adalah yang penting untuk contoh khusus ini.
Brian Borchers
Saya percaya bahwa kesalahan yang dilaporkan dalam makalah ini absolut daripada kesalahan relatif (tidak membuat banyak perbedaan kecuali untuk ϵ=0,01, di mana solusinya turun mendekati 0.
Brian Borchers
Untuk ϵ=0,01 dan N=1024, kesalahan relatif untuk bagian-bagian dari solusi yang mendekati 0 besar, tetapi kesalahan absolut kecil. Saya setuju bahwa makalah itu sangat kabur tentang apa nomor kondisi yang digunakan dan tentang apa "kesalahan" itu sebenarnya (kesalahan relatif atau absolut.)
Brian Borchers
@BrianBorchers Saya tidak yakin apa yang Anda maksud: ini bukan perbedaan antara angka kondisi 2-norm dan infty-norm, melainkan nomor kondisi bijaksana dan komponen-bijaksana (gangguan relatif komponen-bijaksana dalam input, bukan komponen Kesalahan relatif-bijaksana dalam output seperti dalam jawaban Anda).
Kirill