Algoritma Lanczos

Algoritma Lanczos adalah algoritma iteratif adaptasi dari metode daya (power method) untuk menemukan m {\displaystyle m} nilai dan vektor eigen yang "paling berguna" (umumnya yang tertinggi/terendah) dari sebuah matriks Hermite berukuran n × n {\displaystyle n\times n} , dengan m {\displaystyle m} tidak perlu jauh lebih kecil dari n {\displaystyle n} . Algoritma ini dirancang oleh Cornelius Lanczos pada tahun 1950.[1] Meskipun secara prinsip metode ini efisien dalam aspek komputasi, metode yang dirancang pada awalnya tidak berguna karena sifatnya yang tidak stabil secara numerik.

Pada tahun 1970, Ojalvo dan Newman menunjukkan cara membuat metode ini stabil secara numerik, dan mengaplikasikannya untuk menemukan mode getaran dari sebuah struktur teknik berukuran besar.[2] Hal ini dicapai dengan menggunakan sebuah teknik untuk 'memurnikan' vektor-vektor Lanczos (misal dengan mengortogonalkan secara berulang setiap vektor yang baru ditemukan dengan semua vektor yang sudah ditemukan)[2] sampai ke suatu akurasi yang diinginkan. Jika teknik tidak diterapkan, metode akan menghasilkan vektor-vektor yang sangat 'terkontaminasi' oleh vektor-vektor yang berasosiasi dengan frekuensi alami yang rendah.

Dalam karyanya, Ojalvo dan Newman juga mengusulkan cara memilih vektor awal (starting vector; misalnya dengan menggunakan pembangkit bilangan acak), dan mengusulkan metode menentukan nilai m {\displaystyle m} secara empiris (sebaiknya dipilih kurang lebih 1.5 kali dari banyak nilai eigen akurat yang diinginkan). Kontribusi lain diberikan oleh Paige, yang juga memberikan analisis galat untuk metode ini.[3][4] Pada tahun 1988, Ojalvo memberikan sejarah yang lebih akurat dari algoritma ini, dan sebuah uji galat nilai eigen yang efisien.[5]

Algoritma

Algoritma ini memerlukan sebuah matriks Hermite A {\displaystyle A} berukuran n × n {\displaystyle n\times n} , dan secara opsional sebuah bilangan m {\displaystyle m} yang menyatakan banyak iterasi yang diinginkan. Jika nilai m {\displaystyle m} tidak ditentukan, umumnya diambil m = n {\displaystyle m=n} . Sebenarnya, algoritma tidak memerlukan akses ke matriks secara eksplisit, melainkan sebuah fungsi v A v {\displaystyle v\mapsto Av} yang menghasilkan produk perkalian matriks dengan vektor. Fungsi ini dipanggil paling banyak m {\displaystyle m} kali. Pada proses iterasi, algoritma rentan terhadap ketidakstabilan numerik. Ketika dieksekusi dalam aritmatika non-eksak (non-exact arithmetic), pertimbangan-pertimbangan tambahan harus diambil untuk memastikan validitas hasil (dijelaskan di bagian selanjutnya). Di akhir iterasi, algoritma akan menghasilkan sebuah matriks ortonormal V {\displaystyle V} berukuran n × m {\displaystyle n\times m} dan sebuah matriks simetrik real tridiagonal T = V A V {\displaystyle T=V^{*}AV} berukuran m × m {\displaystyle m\times m} . Dalam kasus m = n {\displaystyle m=n} , maka V {\displaystyle V} akan berupa matriks uniter, dan matriks A = V T V {\displaystyle A=VTV^{*}} .

Pada prinsipnya ada empat cara untuk menulis prosedur iterasi. Paige dan karya-karya lainnya menunjukkan bahwa urutan operasi berikut adalah yang paling stabil secara numerik.[6][7]:

  1. Anggap v 1 C n {\displaystyle v_{1}\in \mathbb {C} ^{n}} sebagai sebarang vektor dengan norma Euklidean 1 {\displaystyle 1}
  2. Langkah awal iterasi yang disingkat:
    1. Tetapkan w 1 = A v 1 {\displaystyle w_{1}'=Av_{1}} .
    2. Tetapkan α 1 = w 1 v 1 {\displaystyle \alpha _{1}=w_{1}'^{*}v_{1}} .
    3. Tetapkan w 1 = w 1 α 1 v 1 {\displaystyle w_{1}=w_{1}'-\alpha _{1}v_{1}}
  3. Untuk nilai j = 2 , , m {\displaystyle j=2,\dots ,m} , lakukan:
    1. Tetapkan β j = w j 1 {\displaystyle \beta _{j}=\|w_{j-1}\|} (juga sebuah norma Euklidean).
    2. Jika β j 0 {\displaystyle \beta _{j}\neq 0} , tetapkan nilai v j = w j 1 / β j {\displaystyle v_{j}=w_{j-1}/\beta _{j}} . Namun jika β j = 0 {\displaystyle \beta _{j}=0} , pilih v j {\displaystyle v_{j}} sebagai sebarang vektor bernorma Euclidean 1 {\displaystyle 1} yang ortogonal dengan semua vektor v 1 , , v j 1 {\displaystyle v_{1},\dots ,v_{j-1}} .
    3. Tetapkan w j = A v j {\displaystyle w_{j}'=Av_{j}} .
    4. Tetapkan α j = w j v j {\displaystyle \alpha _{j}=w_{j}'^{*}v_{j}} .
    5. Tetapkan w j = w j α j v j β j v j 1 {\displaystyle w_{j}=w_{j}'-\alpha _{j}v_{j}-\beta _{j}v_{j-1}} .
  4. Bentuk matriks V {\displaystyle V} dengan menyusun v 1 , , v m {\displaystyle v_{1},\dots ,v_{m}} sebagai kolom-kolomnya; dan bentuk matriks
    T = ( α 1 β 2 0 β 2 α 2 β 3 β 3 α 3 β m 1 β m 1 α m 1 β m 0 β m α m ) {\displaystyle T={\begin{pmatrix}\alpha _{1}&\beta _{2}&&&&0\\\beta _{2}&\alpha _{2}&\beta _{3}&&&\\&\beta _{3}&\alpha _{3}&\ddots &&\\&&\ddots &\ddots &\beta _{m-1}&\\&&&\beta _{m-1}&\alpha _{m-1}&\beta _{m}\\0&&&&\beta _{m}&\alpha _{m}\\\end{pmatrix}}}

Catatan: A v j = w j = β j + 1 v j + 1 + α j v j + β j v j 1 {\displaystyle Av_{j}=w_{j}'=\beta _{j+1}v_{j+1}+\alpha _{j}v_{j}+\beta _{j}v_{j-1}} untuk 1 < j < m {\displaystyle 1<j<m} .

Dalam praktiknya vektor awal v 1 {\displaystyle v_{1}} dapat dianggap sebagai argumen input tambahan dari prosedur; sedangkan β j = 0 {\displaystyle \beta _{j}=0} dan indikator-indikator galat numerik dapat diikutsertakan sebagai tambahan kondisi penghentian iterasi.

Tanpa menghitung perkalian matriks-vektor, setiap iterasi melakukan O ( n ) {\displaystyle O(n)} operasi aritmetika. Perkalian matriks-vektor sendiri dapat dilakukan dalam O ( d n ) {\displaystyle O(dn)} operasi aritmetika, dengan d {\displaystyle d} menyatakan rata-rata jumlah elemen bukan nol dalam sebuah baris. Dengan demikian, algoritma Lanczos memiliki total kompleksitas sebesar O ( d m n ) {\displaystyle O(dmn)} , atau O ( d n 2 ) {\displaystyle O(dn^{2})} dalam kasus m = n {\displaystyle m=n} ; hal ini membuatnya bisa sangat cepat untuk matriks rongga. Skema-skema untuk meningkatkan stabilitas numerik biasanya dibandingkan dengan kinerja tinggi ini.

Vektor v j {\displaystyle v_{j}} disebut dengan vektor Lanczos. Vektor w j {\displaystyle w_{j}'} tidak digunakan setelah nilai w j {\displaystyle w_{j}} dihitung, dan vektor w j {\displaystyle w_{j}} tidak digunakan setelah v j + 1 {\displaystyle v_{j+1}} dihitung. Hal ini memungkinkan untuk menggunakan penyimpanan yang sama untuk ketiga vektor tersebut. Begitu juga jika kita hanya ingin mencari matriks tridiagonal T {\displaystyle T} , proses iterasi tidak memerlukan v j 1 {\displaystyle v_{j-1}} setelah menghitung w j {\displaystyle w_{j}} ; meskipun beberapa skema untuk meningkatkan stabilitas numerik akan membutuhkan vektor ini nantinya. Terkadang vektor-vektor Lanczos yang berikutnya dihitung ulang dari v 1 {\displaystyle v_{1}} jika diperlukan.

Aplikasi untuk masalah eigen

Algoritma Lanczos paling sering digunakan dalam konteks menemukan nilai eigen dan vektor eigen dari sebuah matriks. Namun berbeda dengandiagonalisasi matriks yang menghasilkan vektor-vektor eigen dan nilai-nilai eigen terlihat jelas dari pemeriksaan, tridiagonalisasi yang dilakukan algoritma Lanczos memerlukan langkah-langkah tambahan yang nontrivial bahkan untuk menghitung satu nilai atau vektor eigen. Meskipun demikian, menerapkan algoritma Lanczos sering kali merupakan langkah signifikan yang dalam menghitung dekomposisi eigen. Jika λ {\displaystyle \lambda } adalah nilai eigen dari A {\displaystyle A} dan x {\displaystyle x} adalah vektor eigen dari T {\displaystyle T} (sehingga T x = λ x {\displaystyle Tx=\lambda x} ), maka y = V x {\displaystyle y=Vx} adalah vektor eigen yang bersesuaian dari A {\displaystyle A} ; karena

A y = A V x = V T V V x = V T I x = V T x = V ( λ x ) = λ V x = λ y . {\displaystyle Ay=AVx=VTV^{*}Vx=VTIx=VTx=V(\lambda x)=\lambda Vx=\lambda y.}
Hal ini mengartikan algoritma Lanczos mengubah masalah dekomposisi eigen matriks A {\displaystyle A} menjadi masalah dekomposisi eigen matriks T {\displaystyle T} .

Terdapat sejumlah algoritma khusus untuk memroses matriks tridiagonal, sering kali dengan kompleksitas komputasi yang lebih baik daripada algoritma yang umum. Misalkan T {\displaystyle T} adalah matriks simetris tridiagonal berukuran m × m {\displaystyle m\times m} . Beberapa algoritma tersebut diantaranya:

  • Rekursi kontinu (continuant recursion) memungkinkan komputasi polinomial karakteristik dengan O ( m 2 ) {\displaystyle O(m^{2})} operasi, dan mengevaluasinya pada sebuah titik dalam O ( m ) {\displaystyle O(m)} operasi.
  • Algoritma nilai eigen bagi-dan-taklukkan (divide-and-conquer) dapat digunakan untuk menghitung seluruh dekomposisi eigen dari T {\displaystyle T} dalam O ( m 2 ) {\displaystyle O(m^{2})} operasi.
  • Metode Multipole Cepat[8] (Fast Multipole method) yang dapat menghitung semua nilai eigen hanya dengan O ( m log m ) {\displaystyle O(m\log m)} operasi.

Beberapa algoritma dekomposisi eigen yang umum, terutama algoritma QR, diketahui konvergen lebih cepat pada kasus matriks tridiagonal daripada pada matriks yang umum. Kompleksitas asimtotik QR untuk matriks tridiagonal adalah O ( m 2 ) {\displaystyle O(m^{2})} , sama seperti untuk algoritma bagi-dan-taklukkan (meskipun faktor konstantanya mungkin berbeda). Mengingat semua vektor eigen memiliki m 2 {\displaystyle m^{2}} total elemen, algoritma[butuh klarifikasi] ini optimal secara asimtotik. Bahkan untuk algoritma yang tingkat konvergensinya tidak terpengaruh oleh transformasi uniter (seperti metode daya dan iterasi invers), dapat menikmati manfaat kinerja ketika diterapkan ke matriks tridiagonal T {\displaystyle T} daripada matriks asli A {\displaystyle A} . Karena struktur T {\displaystyle T} yang sangat rongga dengan semua elemen bukan nol dalam posisi yang dapat diprediksi, hal ini memungkinkan penyimpanan yang ringkas dengan kinerja yang sangat baik jika dibandingkan dengan menggunakan tembolok (caching). Selain itu, T {\displaystyle T} berupa matriks real dengan semua vektor eigen dan nilai eigen berupa real, sedangkan matriks A {\displaystyle A} secara umum mungkin memiliki elemen kompleks dan vektor eigen; penggunaan aritmetika real cukup untuk mencari vektor eigen dan nilai eigen dari T {\displaystyle T} .

Jika n {\displaystyle n} sangat besar, nilai m {\displaystyle m} dapat dikurangi sehingga matriks T {\displaystyle T} yang dihasilkan masih memiliki ukuran yang dapat dikelola. Pengurangan ini masih memungkinkan untuk menemukan nilai eigen dan vektor eigen yang ekstrem dari A {\displaystyle A} . Dalam keadaan m n {\displaystyle m\ll n} , algoritma Lanczos dapat dianggap sebagai skema kompresi lossy untuk matriks Hermite, yang mengutamakan meyimpanan nilai-nilai eigen ekstrem.

Kombinasi kinerja yang baik untuk matriks rongga dan kemampuan untuk menghitung beberapa (tanpa menghitung semua) nilai eigen adalah alasan utama untuk memilih menggunakan algoritma Lanczos.

Aplikasi

Algoritma Lanczos sangat menarik karena perkalian dengan A {\displaystyle A\,} adalah satu-satunya operasi linier berskala besar. Karena mesin pengambilan teks jangka berbobot hanya menerapkan operasi ini, algoritma Lanczos dapat diterapkan secara efisien ke dokumen teks (lihat Pengindeksan Semantik Laten). Vektor eigen juga penting untuk metode peringkat skala besar seperti Algoritma HITS yang dikembangkan oleh Jon Kleinberg, atau PageRank yang digunakan oleh Google.

Algoritma Lanczos juga digunakan dalam Fisika Materi Terkondensasi sebagai metode untuk menyelesaikan Hamiltonians dari sistem elektron berkorelasi kuat,[9] serta dalam kode model cangkang pada fisika nuklir.[10]

Implementasi

Perpustakaan NAG berisi beberapa rutinitas[11] untuk solusi sistem linear skala besar dan masalah eigen yang menggunakan algoritma Lanczos.

MATLAB dan GNU Octave hadir dengan ARPACK bawaan. Baik matriks yang tersimpan dan implisit dapat dianalisis melalui fungsi eigs() (Matlab/Octave).

Implementasi Matlab dari algoritma Lanczos (masalah presisi catatan) tersedia sebagai bagian dari Gaussian Belief Propagation Matlab Package. The GraphLab[12] perpustakaan pemfilteran kolaboratif menggabungkan implementasi paralel skala besar dari algoritma Lanczos (dalam C ++) untuk multicore.

PRIMME perpustakaan juga menerapkan algoritma seperti Lanczos.

Referensi

  1. ^ Lanczos, C. (1950). "An iteration method for the solution of the eigenvalue problem of linear differential and integral operators" (PDF). J. Res. Nat’l Bur. Std. 45: 255–282. 
  2. ^ a b Ojalvo, I. U.; Newman, M. (1970). "Vibration modes of large structures by an automatic matrix-reduction method". AIAA Journal. 8 (7): 1234–1239. Bibcode:1970AIAAJ...8.1234N. doi:10.2514/3.5878. 
  3. ^ Paige, C. C. (1971). The computation of eigenvalues and eigenvectors of very large sparse matrices (Tesis Ph.D. thesis). U. of London. OCLC 654214109. 
  4. ^ Paige, C. C. (1972). "Computational Variants of the Lanczos Method for the Eigenproblem". J. Inst. Maths Applics. 10 (3): 373–381. doi:10.1093/imamat/10.3.373. 
  5. ^ Ojalvo, I. U. (1988). "Origins and advantages of Lanczos vectors for large dynamic systems". Proc. 6th Modal Analysis Conference (IMAC), Kissimmee, FL. hlm. 489–494. 
  6. ^ Cullum; Willoughby. Lanczos Algorithms for Large Symmetric Eigenvalue Computations. 1. ISBN 0-8176-3058-9. 
  7. ^ Yousef Saad (1992-06-22). Numerical Methods for Large Eigenvalue Problems. ISBN 0-470-21820-7. 
  8. ^ Coakley, Ed S.; Rokhlin, Vladimir (2013). "A fast divide-and-conquer algorithm for computing the spectra of real symmetric tridiagonal matrices". Applied and Computational Harmonic Analysis. 34 (3): 379–414. doi:10.1016/j.acha.2012.06.003 alt=Dapat diakses gratis. 
  9. ^ Chen, HY; Atkinson, W.A.; Wortis, R. (July 2011). "Disorder-induced zero-bias anomaly in the Anderson-Hubbard model: Numerical and analytical calculations". Physical Review B. 84 (4): 045113. arXiv:1012.1031 alt=Dapat diakses gratis. Bibcode:2011PhRvB..84d5113C. doi:10.1103/PhysRevB.84.045113. 
  10. ^ Shimizu, Noritaka (21 October 2013). "Nuclear shell-model code for massive parallel computation, "KSHELL"". arΧiv:1310.5431 [nucl-th]. 
  11. ^ The Numerical Algorithms Group. "Keyword Index: Lanczos". NAG Library Manual, Mark 23. Diakses tanggal 2012-02-09. 
  12. ^ GraphLab Diarsipkan 2011-03-14 di Wayback Machine.

Bacaan lebih lanjut

  • Golub, Gene H.; Van Loan, Charles F. (1996). "Lanczos Methods". Matrix Computations. Baltimore: Johns Hopkins University Press. hlm. 470–507. ISBN 0-8018-5414-8. 
  • Ng, Andrew Y. (2001). "Link Analysis, Eigenvectors and Stability".