Apa itu algoritma sederhana untuk menghitung SVD matriks?
Idealnya, saya ingin algoritma yang kuat secara numerik, tetapi saya ingin melihat implementasi yang sederhana dan yang tidak terlalu sederhana. Kode C diterima.
Adakah referensi untuk makalah atau kode?
Jawaban:
Lihat /math/861674/decompose-a-2d-arbitrary-transform-into-only-scaling-and-rotation (maaf, saya akan menuliskannya di komentar tapi saya sudah mendaftar hanya untuk memposting ini jadi saya belum bisa memposting komentar).
Tetapi karena saya menulisnya sebagai jawaban, saya juga akan menulis metode:
Itu menguraikan matriks sebagai berikut:
Satu-satunya hal yang harus dijaga dengan metode ini adalah bahwa atau untuk atan2.G=F=0 H=E=0
Saya ragu itu bisa lebih kuat dari itu( Perbarui: lihat jawaban Alex Eftimiades!).Referensi tersebut adalah: http://dx.doi.org/10.1109/38.486688 (diberikan oleh Rahul di sana) yang berasal dari bagian bawah posting blog ini: http://metamerist.blogspot.com/2006/10/linear-algebra -for-graphics-geeks-svd.html
Pembaruan: Sebagaimana dicatat oleh @VictorLiu dalam komentar, mungkin negatif. Itu terjadi jika dan hanya jika penentu matriks input negatif juga. Jika itu masalahnya dan Anda ingin nilai singular positif, ambil saja nilai absolut dari .sy sy
sumber
@Pedro Gimeno
"Aku ragu itu bisa lebih kuat dari itu."
Tantangan diterima.
Saya perhatikan pendekatan yang biasa digunakan adalah menggunakan fungsi trigonometri seperti atan2. Secara intuitif, seharusnya tidak perlu menggunakan fungsi trigonometri. Memang, semua hasil berakhir sebagai sinus dan cosinus arctans - yang dapat disederhanakan menjadi fungsi aljabar. Butuh waktu cukup lama, tetapi saya berhasil menyederhanakan algoritma Pedro untuk hanya menggunakan fungsi aljabar.
Kode python berikut melakukan triknya.
sumber
y1
= 0,x1
= 0,h1
= 0, dant1
= 0/0 =NaN
.The GSL memiliki SVD 2-by-2 solver yang mendasari bagian QR dekomposisi algoritma SVD utama untuk
gsl_linalg_SV_decomp
. Lihatsvdstep.c
file dan carisvd2
fungsinya. Fungsi ini memiliki beberapa kasus khusus, tidak sepele, dan terlihat melakukan beberapa hal untuk berhati-hati secara numerik (misalnya, menggunakanhypot
untuk menghindari luapan).sumber
ChangeLog
file jika Anda mengunduh GSL. Dan Anda dapat melihatsvd.c
detail dari keseluruhan algoritma. Satu-satunya dokumentasi yang benar tampaknya untuk fungsi-fungsi yang dapat dipanggil pengguna tingkat tinggi, misalnyagsl_linalg_SV_decomp
,.Ketika kita mengatakan "kuat secara numerik" biasanya kita maksudkan algoritma di mana kita melakukan hal-hal seperti berputar untuk menghindari penyebaran kesalahan. Namun, untuk matriks 2x2, Anda dapat menuliskan hasilnya dalam bentuk rumus eksplisit - yaitu, menuliskan rumus untuk elemen SVD yang menyatakan hasil hanya dalam hal input , bukan dalam hal nilai menengah yang sebelumnya dihitung . Itu berarti bahwa Anda mungkin memiliki pembatalan tetapi tidak ada propagasi kesalahan.
Intinya adalah bahwa untuk sistem 2x2, tidak perlu khawatir tentang ketahanan.
sumber
Kode ini didasarkan pada kertas Blinn ini , kertas Ellis , SVD kuliah , dan tambahan perhitungan. Algoritma cocok untuk matriks nyata reguler dan tunggal. Semua versi sebelumnya berfungsi 100% dan juga yang ini.
sumber
Saya membutuhkan algoritma yang dimiliki
Kami ingin menghitung dan sebagai berikut:c1,s1,c2,s2,σ1 σ2
Gagasan utamanya adalah menemukan matriks rotasi yang mendiagonalisasi , yaitu adalah diagonal.V ATA VATAVT=D
Ingat itu
Mengalikan kedua sisi dengan kita dapatkanS−1
Karena adalah diagonal, pengaturan ke akan memberi kita , artinya adalah matriks rotasi, adalah matriks diagonal, adalah matriks rotasi dan , seperti yang kita cari untuk.D S D−−√ UTU=Identity U S V USV=A
Menghitung rotasi diagonalisasi dapat dilakukan dengan menyelesaikan persamaan berikut:
dimana
dan adalah tangen dari sudut . Ini dapat diturunkan dengan memperluas dan membuat elemen off-diagonal sama dengan nol (mereka sama satu sama lain).t2 V VATAVT
Masalah dengan metode ini adalah bahwa ia kehilangan presisi floating point yang signifikan ketika menghitung dan untuk matriks tertentu, karena pengurangan dalam perhitungan. Solusi untuk ini adalah untuk melakukan dekomposisi RQ ( , segitiga atas dan orthogonal) pertama, kemudian menggunakan algoritma untuk pd . Hal ini memberikan . Perhatikan bagaimana pengaturan to 0 (seperti pada ) menghilangkan beberapa penambahan / pengurangan. (Dekomposisi RQ cukup sepele dari perluasan produk matriks).β−α γ A=RQ R Q USV′=R USV=USV′Q=RQ=A d R
Algoritme yang diimplementasikan secara naif dengan cara ini memiliki beberapa anomali numerik dan logis (misalnya atau ), yang saya perbaiki dalam kode di bawah ini.S +D−−√ −D−−√
Saya melemparkan sekitar 2000 juta matriks acak pada kode, dan kesalahan numerik terbesar yang dihasilkan adalah sekitar (dengan floats 32 bit, ). Algoritma ini berjalan di sekitar 340 siklus clock (MSVC 19, Ivy Bridge).6⋅10−7 error=||USV−M||/||M||
Gagasan dari:
http://www.cs.utexas.edu/users/inderjit/public_papers/HLA_SVD.pdf
http://www.math.pitt.edu/~sussmanm/2071Spring08/lab09/index.html
http: // www.lucidarme.me/singular-value-decomposition-of-a-2x2-matrix/
sumber
Saya telah menggunakan deskripsi di http://www.lucidarme.me/?p=4624 untuk membuat kode C ++ ini. Matriksnya adalah dari perpustakaan Eigen, tetapi Anda dapat dengan mudah membuat struktur data Anda sendiri dari contoh ini:
Dengan fungsi tanda standar
Ini menghasilkan nilai yang persis sama dengan
Eigen::JacobiSVD
(lihat https://eigen.tuxfamily.org/dox-devel/classEigen_1_1JacobiSVD.html ).sumber
S2 = hypot( a*a + b*b - c*c - d*d, 2*(a*c + b*d))
I have pure C code for the 2x2 real SVD here. See line 559. It essentially computes the eigenvalues ofATA by solving a quadratic, so it's not necessarily the most robust, but it seems to work well in practice for not-too-pathological cases. It's relatively simple.
sumber
For my personal need, I tried to isolate the minimum computation for a 2x2 svd. I guess it is probably one of the simplest and fastest solution. You can find details on my personal blog : http://lucidarme.me/?p=4624.
Advantages : simple, fast and you can only calculate one or two of the three matrices (S, U or D) if you don't need the three matrices.
Drawback it uses atan2, which may be inacurate and may require an external library (typ. math.h).
sumber
Here is an implementation of a 2x2 SVD solve. I based it off of Victor Liu's code. His code was not working for some matrices. I used these two documents as mathematical reference for the solve: pdf1 and pdf2.
The matrix
setData
method is in row-major order. Internally, I represent the matrix data as a 2D array given bydata[col][row]
.sumber