Bisakah solusi sistem linear persamaan diperkirakan hanya untuk beberapa variabel pertama?

15

Saya memiliki sistem persamaan linear ukuran mxm, di mana m besar. Namun, variabel yang saya minati hanyalah variabel n pertama (n kecil dibandingkan dengan m). Apakah ada cara saya dapat memperkirakan solusi untuk nilai m pertama tanpa harus menyelesaikan seluruh sistem? Jika demikian, apakah perkiraan ini lebih cepat daripada menyelesaikan sistem linear penuh?

Paul
sumber
2
Tidak, kecuali fungsi pemaksaan Anda juga terbatas pada variabel n pertama. Jika ya, Anda dapat membentuk komplemen Schur, meskipun kemungkinan padat. Jika operator asli Anda jarang, mungkin tidak sepadan.
Jack Poulson
1
Saya kira Anda bisa menggunakan eliminasi gaussian mulai dari sudut kanan bawah matriks. Ini akan ~ 2x lebih cepat daripada eliminasi gaussian biasa jika Anda hanya peduli pada beberapa elemen pertama dan berhenti di tengah jalan. Saya tidak tahu bagaimana ini akan dibandingkan dengan metode berulang.
Dan
4
@OscarB: Harap tidak. Aturan Cramer adalah kekejaman dalam aritmatika floating point. Saya belum pernah mendengarnya digunakan untuk perhitungan serius, dan dibutuhkan pemikiran yang cukup untuk menghindari kompleksitas faktorial , di mana ia masih tidak bersaing dengan eliminasi Gaussian.
Jack Poulson
1
@ Paul: Sebagian besar pengurangan model pesanan digunakan dalam konteks sistem ODE atau DAE yang besar. Kadang-kadang, metodologi reduksi dimotivasi oleh sistem ODE atau DAE yang muncul dari diskritisasi PDE. Saya belum melihat pengurangan model yang digunakan pada persamaan aljabar murni. (Jika sudah, kirimkan saya referensi, karena saya sedang melakukan tesis tentang metode pengurangan model dan akan sangat tertarik untuk melihatnya.) Jika Anda mau, saya bisa membuat sketsa seperti apa model pengurangan jika kita memperlakukan persamaan aljabar sebagai kasus degenerasi dari sistem persamaan diferensial-aljabar.
Geoff Oxberry
1
@JackPoulson - apakah Anda keberatan meringkas komentar Anda sebagai jawaban? Saya pikir itu solusi yang paling benar dan saya tidak ingin itu hilang dalam komentar.
Aron Ahmadia

Jawaban:

13

Seperti yang telah ditunjukkan orang lain, ini sulit dilakukan dengan pemecah langsung. Yang mengatakan, itu tidak sulit untuk dilakukan dengan pemecah iteratif. Untuk tujuan ini, perhatikan bahwa sebagian besar pemecah iteratif dengan satu atau lain cara meminimalkan kesalahan sehubungan dengan beberapa norma. Seringkali, norma ini disebabkan oleh matriks itu sendiri, tetapi kadang-kadang juga hanya norma vektor l2. Tetapi itu tidak harus menjadi masalah: Anda dapat memilih norma mana yang Anda inginkan untuk meminimalkan kesalahan (atau residual), dan Anda dapat, misalnya, memilih norma di mana Anda menimbang komponen yang Anda pedulikan dengan 1 dan semua yang lain dengan 1e-12, misalnya misalnya sesuatu seperti (1e-24)N i = 6 x 2 i dan produk skalar yang sesuai. Kemudian tulis semua langkah solver iteratif sehubungan dengan norma dan produk skalar ini, dan Anda mendapatkan solver iteratif yang secara signifikan lebih memperhatikan elemen vektor yang Anda pedulikan daripada yang lain.||x||2=i=15xi2+i=6Nxi2

Pertanyaannya tentu saja adalah apakah Anda membutuhkan iterasi yang lebih sedikit dibandingkan dengan produk norma / skalar yang menimbang semua komponen secara merata. Tapi itu memang harus menjadi kasus: katakanlah Anda hanya peduli pada lima elemen vektor pertama. Maka Anda harus paling banyak lima iterasi untuk mengurangi kesalahan dengan faktor 1e12 karena lima iterasi adalah apa yang diperlukan untuk sistem 5x5 yang menggambarkannya. Itu bukan bukti tetapi saya cukup yakin bahwa Anda memang harus pergi dengan jumlah iterasi yang jauh lebih kecil jika bobot dalam norma (1e-12 di atas) lebih kecil daripada toleransi yang Anda inginkan untuk menyelesaikan sistem linear secara iteratif. .

Wolfgang Bangerth
sumber
2
Hmm, poin bagus. Saya akan tertarik melihat contoh nyata, karena saya agak khawatir tentang efek dari hanya mencoba untuk menyelesaikan beberapa derajat kebebasan; walaupun residu mungkin kecil, mungkin norma kesalahan masih cukup besar (lakukan untuk secara efektif mengabaikan sebagian besar operator).
Jack Poulson
Secara intuitif, ini hanya berfungsi jika komponen-komponen sistem yang sangat kecil benar-benar mendominasi jawaban dalam L2 (atau norma yang Anda pahami kesalahan Anda untuk diukur dalam arti). Kalau tidak, saya pikir kepedulian Jack adalah valid, tetapi saya pasti akan tertarik bahkan melihat bukti numerik dari ini ...
Aron Ahmadia
Seseorang harus memastikan Anda mengambil metode yang meminimalkan kesalahan , bukan sisa. Saya pikir MinErr mungkin merupakan titik awal yang baik.
Wolfgang Bangerth
@ WolfgangBangerth: Saya tidak akrab dengan MINERR: apakah ini referensi utama?
Jack Poulson
1
Bahkan itu tidak cukup, karena Anda akan tidak akurat. Anda tidak bisa mendapatkan beberapa komponen secara akurat menggunakan bobot ini.
Matt Knepley
17

Membentuk komplemen Schur

Misalkan Anda telah mengijinkan dan mempartisi matriks Anda ke dalam formulir

A=(A11A12A21A22),

sedemikian rupa sehingga mengandung derajat kebebasan Anda dan jauh lebih kecil dari A 11 , maka seseorang dapat membentuk komplemen SchurA22A11

S22:=A22A21A111A12,

baik melalui faktorisasi LU yang tampak benar parsial atau rumus eksplisit, dan kemudian dapat dipahami dalam pengertian berikut:S22

S22x=y(A11A12A21A22)(x)=(0y),

di mana mewakili bagian 'tidak menarik' dari solusi. Dengan demikian, disediakan sisi kanan yang hanya nol dalam derajat kebebasan komplemen Schur S 22 , kita hanya perlu menyelesaikan terhadap S 22 untuk mendapatkan bagian dari solusi yang sesuai dengan derajat kebebasan tersebut.S22S22

Kompleksitas komputasi dalam kasus padat yang tidak terstruktur

Pengaturan dengan tinggi A dan n dengan tinggi A 22 , maka metode standar untuk komputasi S 22 adalah untuk Faktor pertama L 11 U 11 : = A 11 (mari kita mengabaikan berputar untuk saat ini) di sekitar 2 / 3 ( N - n ) 3 bekerja, lalu membentukNAnA22S22L11U11:=A112/3(Nn)3

S22:=A22(A21U111)(L111A12)=A22A21A111A12

menggunakan dua penyelesaian segitiga yang masing - masing membutuhkan kerja, dan kemudian melakukan pembaruan ke A 22 dalam 2 n 2 ( N - n ) bekerja.n(Nn)2A222n2(Nn)

Dengan demikian, total pekerjaan adalah sekitar . Ketika n sangat kecil, N - n N , sehingga biaya dapat dilihat secara kasar 2 / 3 N 3 , yang merupakan biaya faktorisasi penuh.2/3(Nn)3+2n(Nn)2+2n2(Nn)nNnN2/3N3

The benefit is that, if there is a very large number of right-hand sides to be solved with the same system of equations, then S22 could potentially be reused a large number of times, where each solve would only require 2n2 work (rather than 2N2 work) if S22 is factored.

Computational complexity in the (typical) sparse case

If your sparse system arose from some type of finite-difference or finite-element approximation, then sparse-direct solvers will almost certainly be able to exploit some of the structure; 2d systems can be solved with O(N3/2) work and O(NlogN) storage, while 3d systems can be solved with O(N2) work and O(N4/3) storage. The factored systems can then be solved with the same amount of work as the storage requirements.

The point of bringing up the computational complexities is that, if nN and you have a 2d system, then since the Schur complement will likely be dense, the solve complexity given the factored Schur complement will be O(n2)=O(N), which is only missing a logarithmic factor versus solving the full system! In 3d, it requires O(N) work instead of O(N4/3).

It is thus important to keep in mind that, in your case where n=N, there will only be significant savings if you're working in several dimensions and have many right-hand sides to solve.

Jack Poulson
sumber
1
This is a great summary of the schur complement method and when it is computationally efficient to use it!
Paul
6

The model reduction approach

Since Paul asked, I'll talk about what happens if you use projection-based model reduction methods on this problem. Suppose that you could come up with a projector P such that the range of P, denoted R(P), contains the solution to your linear system Ax=b, and has dimension k, where k is the number of unknowns for which you wish to solve in a linear system.

A singular value decomposition of P will yield the following partitioned matrix:

P=[V][diag(1k)000][WT].

The matrices obscured by stars matter for other things (like estimating error, etc.), but for now, we'll avoid dealing with extraneous details. It follows that

P=VWT

is a full rank decomposition of P.

Essentially, you'll solve the system

PAx=Pb

in a clever way, because V and W also have the property that WTV=I. Multiplying both sides of PAx=Pb by WT and letting y=Vx^ be an approximation for x yields

WTAx^=WTb.

Solve for x^, premultiply it by V, and you have y, your approximation for x.

Why the Schur complement approach is probably better

For starters, you have to pick P somehow. If the solution to Ax=b is in R(P), then y=x, and y isn't an approximation. Otherwise, yx, and you introduce some approximation error. This approach doesn't really leverage at all the structure you mentioned wanting to exploit. If we pick P such that its range is the standard unit basis in the coordinates of x you want to calculate, the corresponding coordinates of y will have errors in them. It's not clear how you'd want to pick P. You could use an SVD of A, for instance, and select P to be the product of the first k left singular vectors of A and the adjoint of the first k right singular vectors of A, assuming that singular vectors are arranged in decreasing order of singular value. This choice of projector would be equivalent to performing proper orthogonal decomposition on A, and it would minimize the L2-error in the approximate solution.

In addition to introducing approximation errors, this approach also introduces three extra matrix multiplies on top of the linear solve of the smaller system and the work needed to calculate V, and W. Unless you're solving the same linear system a lot, only changing the right hand side, and P is still a "good" projection matrix for all of those systems, those extra costs will probably make solving the reduced system more expensive than solving your original system.

The drawbacks are much like JackPoulson's approach, except that you're not quite leveraging the structure that you mentioned.

Geoff Oxberry
sumber
4

The long answer is...sort of.

You can re-arrange your system of equations such that the farthest right k columns are the variables which you wish to solve for.

Step 1: Perform Gaussian Elimination so that the matrix is upper triangular. Step 2: solve by back substitution for only the first (last) k variables which you are interested in

This will save you the computational complexity of having to solve for the last nk variables via back-substitution, which could be worth it if n is as large as you say. Keep in mind that a fair amount of work will still have to be done for step 1.

Also, keep in mind that restricting the order in which you are going to perform back-substituion may restrict the form of the matrix (it takes away the ability to exchange columns) which could possibly lead to an ill conditioned system, but I am not sure about that - just something to keep in mind.

drjrm3
sumber
Gaussian elimination requires O(n3) work, but backward substitution only requires O(n2). Thus, as n grows larger, the percentage of time spent in the triangle solve becomes vanishingly small.
Jack Poulson
which is why the answer is "sort of" instead of "yes" =)
drjrm3
It makes sense that it can be done this way... However, the bulk of the computation in a Gaussian Elimination is in the forward elimination phase, yielding an O(n^3) complexity despite the truncated backward substitution phase. I was hoping there was a faster method...
Paul