Diagonalisasi Matriks Pengkondisian Padat

10

Saya mencoba mendiagonalisasi beberapa matriks yang padat dan tidak terkondisikan. Dalam presisi mesin, hasilnya tidak akurat (mengembalikan nilai eigen negatif, vektor eigen tidak memiliki simetri yang diharapkan). Saya beralih ke fungsi Mathematica Eigensystem [] untuk mengambil keuntungan dari presisi yang berubah-ubah, tetapi perhitungannya sangat lambat. Saya terbuka untuk sejumlah solusi. Apakah ada paket / algoritma yang cocok untuk masalah yang dikondisikan? Saya bukan ahli prakondisi, jadi saya tidak yakin seberapa banyak ini bisa membantu. Kalau tidak, yang bisa saya pikirkan adalah pemecah nilai eigen presisi yang diparalelkan secara acak, tetapi saya tidak terbiasa dengan apa pun di luar Mathematica, MATLAB dan C ++.

Untuk memberikan latar belakang masalah, matriksnya besar, tetapi tidak besar (paling banyak 4096x4096 hingga 32768x32768). Mereka nyata, simetris, dan nilai eigen dibatasi antara 0 dan 1 (eksklusif), dengan banyak nilai eigen sangat dekat dengan 0 dan tidak ada yang mendekati 1. Matriks pada dasarnya adalah operator konvolusi. Saya tidak perlu mendiagonalisasi semua matriks saya, tetapi semakin besar saya bisa, semakin baik. Saya memiliki akses ke cluster komputasi dengan banyak prosesor dan kemampuan komputasi terdistribusi.

Terima kasih

Leigh
sumber
2
Rutin apa yang Anda gunakan untuk mendiagonalisasi matriks simetris Anda yang sebenarnya? Dan dalam hal apa dekomposisi nilai eigen tidak akurat?
Jack Poulson
Inilah ide yang berkaitan dengan jawaban Arnold: lakukan dekomposisi Cholesky dari matriks SPD Anda, dan kemudian temukan nilai singular dari segitiga Cholesky yang baru saja Anda peroleh, mungkin menggunakan algoritma tipe dqd untuk menjaga akurasi.
JM
1
@ JM: Membentuk dekomposisi Cholesky dari matriks definitif positif tunggal numerik tidak stabil secara numerik dengan metode yang biasa, karena seseorang mungkin menemukan pivot negatif. (Misalnya, chol Matlab (A) biasanya gagal.) Seseorang harus mengaturnya ke nol dan memusnahkan baris faktor yang sesuai. Melakukan ini memberi cara untuk mendapatkan ruang nol numerik secara andal.
Arnold Neumaier
@Arnold, jika ingatanku, ada adaptasi Cholesky yang menggunakan pivot simetris untuk kasus-kasus di mana matriks positif semi- terbatas (atau hampir jadi). Mungkin itu bisa digunakan ...
JM
@ JK: Seseorang tidak perlu berputar untuk menyelesaikan kasus semidefinite; resep yang saya berikan sudah cukup. Saya hanya ingin menunjukkan bahwa seseorang tidak dapat menggunakan program kaleng standar tetapi harus memodifikasinya sendiri.
Arnold Neumaier

Jawaban:

7

Hitung SVD sebagai pengganti dekomposisi spektral. Hasilnya sama dalam aritmatika persis, karena matriks Anda pasti pasti simetris, tetapi dalam aritmatika presisi terbatas, Anda akan mendapatkan nilai eigen kecil dengan akurasi lebih banyak.

Sunting: Lihat Demmel & Kahan, Nilai Singular Akurat dari Matriks Bidiagonal, SIAM J. Sci. Stat. Komputasi. 11 (1990), 873-912.
ftp://netlib2.cs.utk.edu/lapack/lawnspdf/lawn03.pdf

Edit2; Perhatikan bahwa tidak ada metode yang dapat menyelesaikan nilai eigen yang lebih kecil dari pada norma waktu akurasi mesin yang digunakan, karena mengubah entri tunggal dengan satu ulp mungkin sudah mengubah nilai eigen kecil dengan ini. Jadi mendapatkan nol nilai eigen sebagai pengganti yang sangat kecil adalah tepat, dan tidak ada metode (kecuali bekerja dengan presisi yang lebih tinggi) akan mengurai vektor eigen yang sesuai, tetapi hanya mengembalikan dasar untuk ruang nol numerik yang umum.

Arnold Neumaier
sumber
[0,BT;B,0]
2
@JackPoulson: Intinya adalah bahwa bentuk bidiagonal menentukan nilai singular kecil jauh lebih baik. Bentuk tridiagonal simetris yang terkait memiliki nol pada diagonal, yang dipertahankan oleh reduksi bidiagonal menjadi diagonal, tetapi tidak dengan QR yang diterapkan pada tridiagonal.
Arnold Neumaier
1
Referensi? Metode Jacobi dikenal sangat akurat (meskipun lambat).
Jack Poulson
@JackPoulson: Coba dan lihat. Demmel & Kahan, Nilai Singular Akurat dari Bidiagonal Matrices, 202.38.126.65/oldmirrors/ftp.netlib.org/lapack/lawnspdf/…
Arnold Neumaier
[0,BT;B,0]
1

Terima kasih atas saran ini. Saya mencoba perintah SVD Mathematica, tetapi saya tidak mendapatkan peningkatan yang nyata (masih kehilangan simetri yang sesuai, 'nilai eigen' salah nol di mana mereka salah keluar negatif sebelumnya). Mungkin saya perlu mengimplementasikan salah satu algoritma yang Anda jelaskan di atas alih-alih fungsi bawaan? Saya mungkin ingin menghindari masalah menggunakan metode tertentu seperti ini kecuali saya yakin sebelumnya bahwa itu akan menawarkan peningkatan yang signifikan.

@ JackPoulson, saya membaca skrip tentang metode Jacobi yang Anda referensikan, dan itu terlihat menjanjikan. Bisakah Anda atau siapa pun merekomendasikan cara yang baik untuk menerapkan metode Jacobi untuk menemukan sistem eigens? Saya menduga bahwa jika saya mengkodekannya sendiri (dalam MATLAB), itu akan sangat lambat.

Leigh
sumber
Saya belum mengujinya, tetapi ada implementasi MATLAB di sini: groups.google.com/forum/?fromgroups#!msg/sci.math.num-analysis/…
Jack Poulson
Perhatikan bahwa tidak ada metode yang dapat menyelesaikan nilai eigen yang lebih kecil dari pada norma waktu akurasi mesin yang digunakan, karena mengubah entri tunggal dengan satu ulp mungkin sudah mengubah nilai eigen kecil dengan ini. Jadi mendapatkan nol nilai eigen sebagai pengganti yang sangat kecil adalah tepat, dan tidak ada metode (kecuali bekerja dengan presisi yang lebih tinggi) akan mengurai vektor eigen yang sesuai, tetapi hanya mengembalikan dasar untuk ruang nol numerik yang umum. Untuk apa Anda membutuhkan nilai eigen?
Arnold Neumaier
@ArnoldNeumaier: Saya menjalankan beberapa tes di MATLAB dengan nilai eigen di kisaran [0,1], dengan satu nilai eigen secara manual ditetapkan ke nilai-nilai seperti 6.3e-16, dan rutin SVD Octave (berdasarkan dgesvd, yang menggunakan reduksi menjadi bidiagonal dan lalu QR) mengambil nilai-nilai ini jauh lebih akurat daripada Octave eig. Kode Jacobi yang ditautkan tampaknya terlalu lambat untuk digunakan, bahkan pada matriks berukuran sedang.
Jack Poulson
@JackPoulson: Ya. Tetapi Leigh tampaknya mengeluh tentang banyak nilai eigen yang sangat kecil, dan vektor eigennya jarang yang dirancang tetapi akan bercampur dengan bebas, tidak peduli metode mana yang digunakan. Dan nilai positif positif yang sangat kecil (lebih kecil dari 1e-16) tentu saja akan ditemukan nol.
Arnold Neumaier
@ArnoldNeumaier benar bahwa saya menemukan beberapa nilai eigen yang sangat kecil, yang saya duga memperburuk masalah. Saya tidak menyadari (meskipun dalam retrospeksi itu jelas) bahwa nilai eigen kurang dari 1e-16 akan menjadi nol di floating point. Saya kira meskipun jumlahnya dapat disimpan, kesalahan pembulatan terjadi ketika menambahkannya ke nomor yang lebih besar. Vektor eigen memberi tahu saya jika masalah tertentu dapat dipecahkan. Vektor eigen memungkinkan penguraian masalah menjadi bagian yang dapat dipecahkan dan tidak dapat dipecahkan. Jika saya pada dasarnya dibatasi oleh presisi, maka dapatkah Anda merekomendasikan paket untuk solusi yang lebih cepat?
Leigh