Memecahkan sistem linear dengan argumen matriks

10

Kita semua akrab dengan banyak metode komputasi untuk menyelesaikan sistem linear standar

Ax=b.
Namun, saya ingin tahu apakah ada metode komputasi "standar" untuk memecahkan sistem linier yang lebih umum (dimensi terbatas) dari formulir tersebut.

A m 1 × n 1 B m 2 × n 2 L m 1 × n 1 m 2 × n 2 A x = b

LA=B,
mana, katakanlah, adalah matriks , adalah matriks , dan adalah operator linier yang mengambil matriks ke matriks , yang tidak melibatkan vektorisasi matriks , yaitu mengkonversi semuanya ke bentuk standar . Am1×n1Bm2×n2Lm1×n1m2×n2Ax=b

Alasan saya bertanya adalah saya harus menyelesaikan persamaan berikut untuk :u

R R u f

(RR+λI)u=f
mana adalah transformasi Radon 2d, adjoinnya, dan dan adalah array 2d (gambar). Dimungkinkan untuk memvariasikan persamaan ini, tetapi ini menyebalkan, terutama jika kita beralih ke 3D.RRuf

Lebih umum, bagaimana dengan array ? Misalnya, menyelesaikan mana dan adalah array 3D (saya harus melakukan ini dengan Radon transform di beberapa titik juga).nDLA=BAB

Terima kasih sebelumnya, dan jangan ragu untuk mengirim saya ke StackExchange lain jika Anda merasa perlu.

icurays1
sumber
1
Anda mungkin dapat membangun preconditioner bertingkat yang efektif kemudian menggunakan gradien konjugat. Saya memiliki masalah serupa di mana ini cukup efektif dan sangat dapat diparalelkan. Jika Anda ingin metode langsung, pertimbangkan reduksi ke bentuk schur seperti dalam makalah ini tentang persamaan Lyapunov: cs.cornell.edu/cv/ResearchPDF/Hessenberg.Schur.Method.pdf
Nick Alger
Luar biasa, terima kasih atas referensi! Saya baru saja mendapat CG untuk bekerja secara efektif, jadi saya senang.
icurays1

Jawaban:

9

Ya, Anda sudah benar, dan itu memang akan berfungsi dengan baik ketika Anda meningkatkan ke 3-D. Bagian yang paling mudah, sebenarnya, adalah produk dalam --- lakukan saja produk titik standar pada vektor yang setara, ekuivalen . Karena Anda mungkin menyimpan data secara bersamaan, Anda dapat melakukannya di tempat. Ini bahkan bekerja dengan ruang vektor yang kompleks --- perlakukan saja nilai-nilai kompleks sebagai pasangan nilai nyata. Itu karena untuk CG Anda memerlukan produk dalam nyata .Rny,xRe(yHx)

Satu hal yang Anda harus berhati-hati ketika Anda menerapkan CG (atau pendekatan berulang serupa) dengan operator linier umum adalah menerapkan adjoint operator linier Anda dengan benar. Artinya, orang sering mendapatkan benar, tetapi melakukan kesalahan dengan mengimplementasikan .y=F(x)z=F(y)

Saya sarankan menerapkan tes sederhana yang memanfaatkan identitas berikut: untuk setiap dan sesuai , Jadi yang Anda lakukan adalah menghasilkan nilai acak dan , menjalankannya masing-masing melalui operasi forward dan adjoint Anda, dan menghitung dua produk dalam di atas. Pastikan mereka cocok dengan presisi yang masuk akal, dan ulangi beberapa kali.xy

y,F(x)=F(y),x.
xy

EDIT: apa yang Anda lakukan jika operator linier Anda seharusnya simetris? Nah, Anda juga perlu memverifikasi simetri itu. Jadi gunakan tes yang sama, hanya mencatat bahwa --- terapkan operasi yang sama untuk dan . Tentu saja, OP memiliki operator asimetris dan simetris untuk menangani ...F=Fxy

Michael Grant
sumber
Terima kasih @ChristianClason! Saya tahu dari pengalaman bagaimana kesalahan frustasi dalam perhitungan adjoint dapat terjadi. :) Dalam paket TFOCS kami menerapkan linop_test.mrutin untuk alasan ini. Paket itu juga mendukung matriks, array, dan produk Cartesian dalam ruang vektor.
Michael Grant
3

Ternyata, karena sistem saya simetris dan pasti positif (karena operator linier saya ditulis sebagai ), gradien konjugat dapat diadaptasi untuk menyelesaikan jenis persamaan ini secara berulang. Satu-satunya modifikasi muncul ketika menghitung produk dalam - yaitu perhitungan khas produk dalam di CG terlihat seperti atau . Dalam versi modifikasi, kami menggunakan produk dalam Frobenius, yang dapat dihitung dengan menjumlahkan entri dari produk Hadamard (searah). YaituRR+λIrkTrkpkTApk

A,B=i,jAijBij

Saya menduga ini akan baik-baik saja ketika saya memutakhirkan ke array 3D, meskipun saya belum melihat produk dalam Frobenius didefinisikan pada array 3D (saya akan bekerja di bawah asumsi bahwa saya dapat lagi hanya menjumlahkan produk pointwise).

Saya masih tertarik dengan metode yang lebih umum jika ada yang tahu!

icurays1
sumber