SOLUSI MASALAH SPEKTRAL MATRIX RAKSASA DENGAN ITERASI SUBRUANG KRYLOV & METODE LANCZOS M. Bunjamin *
ABSTRAK SOLUSI MASALAH SPEKTRAL MATRIX RAKSASA DENGAN ITERASI SUBRUANG KRYLOV & METODE LANCZOS. Berikut disajikan suatu algoritma komputasi untuk mencari solusi masalah spektral matrix raksasa (monster matrices), yaitu menentukan nilai pribadi (eigenvalues) dan vektor pribadi (eigenvectors) matrix real simetrik A N× N , dengan N ribuan sampai jutaan. Metode Iterasi Subruang Krylov dan algoritma Lanczos yang disajikan berikut adalah versi modern yang telah mengalami banyak pengembangan seperti yang dilaporkan oleh Lin & Gubernatis dll.
ABSTRACT THE SOLUTION OF MONSTER MATRICES PROBLEMS WITH THE KRYLOV SUBSPACE ITERATION AND THE LANCZOS METHOD. A computational algorithm to solve the spectral problems of monster matrices, i.e. finding the eigenvalues and eigenvectors of real symmetric matrices A N× N , with N in the thousands or millions, will be presented below. The Krylov subspace iteration and the Lanczos methods to be presented below are modern versions that have undergone many improvements as reported by Lin & Gubernatis and others.
SUBRUANG INVARIAN Jika matrix real simetrik 0A beroperasi di ruang vektor real R N berdimensi N, suatu konsep penting yang diperlukan dalam masalah spektral ialah subruang invarian / invariant subspace, yang didefinisikan sebagai berikut: Subruang R M ⊂ R N , M < N, disebut subruang invarian terhadap matrix A jika ∀x ∈ R M è Ax ∈ R M . Jika x ∈ R M suatu vektor pribadi dari A, maka Ax = λx juga ∈ R M . Karena itu subruang invarian berdimensi M dari A dibentang (spanned) oleh M buah vektor pribadi dari A.
*
FMIPA-UI dan FMIPA-UNAS
Jika M buah vektor {q1 , q 2 , L , q M } adalah basis dari subruang invarian
R
M
, dan mereka membentuk M kolom dari matrix Q N×M , maka M kolom dari
matrix AQ adalah M vektor ∈ R M , yaitu kombinasi linear dari {q1 , q 2 , L , q M } . Jadi dapat ditulis:
Aq 1 = t11q1 + t 21q 2 + L + t M 1q M Aq 2 = t12 q1 + t 22 q 2 + L + t M 2 q M M Aq M = t1M q1 + t 2 M q 2 + L + t MM q M t 11 t T = 21 M t M1
t12 t 22 M tM2
AQ=QT
L t 1M L t 2M . O M L t MM
Jika basis {q1 , q 2 , L , q M } ortonormal
QTQ = I
Q TM ×N A N ×N Q N ×M = TM ×M
Ortonormalisasi {q1 , q 2 ,L, q M } dibantu metode Gram-Schmidt. Dari hasil di atas dapat disimpulkan hal penting berikut. Andaikan λ dan y adalah nilai pribadi dan vektor pribadi matrix T, alias Ty = λy, maka QTy = λQy, atau A(Qy) = λ Qy). Ini berarti bahwa jika λ dan y adalah nilai pribadi dan vektor pribadi matrix T, hal ini berakibat bahwa λ dan Qy adalah nilai pribadi dan vektor pribadi matrix A. Hasil ini sangat penting karena memberi kita jalan untuk mencari sebagian dari nilai pribadi dan vektor pribadi matrix raksasa A N× N dengan bantuan matrix kecil TM ×M , sehingga lebih mudah dan dalam jangkauan komputer kita. Kunci keberhasilan metode ini antara lain terletak pada a) kemampuan kita menemukan subruang invarian dari matrix A, yang dalam praktek didekati oleh subruang Krylov (lihat di bawah), dan b) kemampuan kita membuat kolom-kolom matrix Q saling ortonormal, suatu hal yang tidak mudah karena keterbatasan presisi
komputer (= finite precision), sehingga penggunaan “double -precision” adalah suatu keharusan. Inipun sebenarnya belum cukup.
METODE LANCZOS DAN SUBRUANG KRYLOV Jika diketahui matrix real simetrik A N× N dan vektor sembarang b = b1 ≠ 0, metode Lanczos diawali dengan menentukan barisan vektor sebagai berikut: b1 ,
b 2 = Ab 1 , b 3 = A 2 b 1 , …, b M = A M −1b 1 . Jika M vektor {b1 , b 2 , L , b M } ini bebas linear (linearly independent), mereka membentang apa yang disebut subruang Krylov K berdimensi M. Karena operasi oleh A, barisan vektor b1 , b 2 , L , b M ∈ K, tetapi b M +1 = A M b ∉ K . Dari metode pangkat (lihat Bunjamin bab 5) kita tahu bahwa barisan vektor ini konvergen menuju ke vektor pribadi pasangan dari nilai pribadi dominan dari A jika M cukup besar, sehingga A M b = A(A M −1 b) ≈ λ(A M −1 b), λ nilai pribadi dominan dari A è A M b hampir proporsional dengan
A M −1 b è A M b hampir ∈ K. Jadi subruang Krylov K hampir invarian terhadap matrix A, dan aproximasi beberapa nilai pribadi & vektor pribadi dari A yang berukuran raksasa dapat diperoleh dari nilai pribadi & vektor pribadi dari matrix T yang jauh lebih kecil dan lebih mudah. Beberapa sifat dari matrix T dalam metode Lanczos adalah sebagai berikut. Matrix T = Q T AQ adalah simetrik dan tridiagonal (Stoer & Bulirsch). Suatu rumus rekurensi tiga suku tersedia untuk menghitung kolom Q (lihat di bawah). Matrix A diperlukan hanya untuk menghitung perkalian matrix dengan vektor. Konvergensi ternyata cukup cepat. Misalnya, Lin & Gubernatis melaporkan bahwa untuk N = 16 juta, cukup diambil M di antara 30 sampai 100. Jika semua unsur a ij , i ≠ j, bertanda sama, maka konvergensi terjamin. Jika diandaikan matrix T yang tridiagonal mempunyai unsur-unsur sebagai berikut :
β1 α1 T= 0 M 0
α1 β2 α2 O L
0 α2 β3 O 0
L O O O α M −1
, α M −1 β M 0 M 0
maka dapat dijabarkan beberapa rumus berikut. Dari persamaan AQ = QT serta mengingat bentuk T yang tridiagonal di atas, dapat disimpulkan berlakunya rumus rekurensi tiga-suku berikut. Aq i = α i −1q i −1 + βi q i + α i q i +1 , i = 1,…,M 1, di mana q i adalah kolom # i dari Q.
LANGKAH-LANGKAH ALGORITMA LANCZOS Algoritma Lanczos diawali dengan memilih skalar α 0 = 1, vektor q 0 = 0 dan vektor q1 sembarang asal q1 = 1. Vektor {q i } yang saling ortonormal berarti
βi = q Ti Aq i dan α i = q Ti Aq i +1 . Langkah Lanczos akan dibantu sedikit jika didefinisikan dua vektor bantu, yaitu u i = Aq i − α i −1q i −1 dan ri = u i − βi q i , sehingga βi = q Ti u i dan ri = α i q i +1 . Langkah yang lebih urut dan rinci adalah sbb. 1. a. Dengan modal α 0 = 1, q 0 = 0 dan vektor awal sembarang q1 ≠ 0 dan q1 = 1, hitung vektor u1 = Aq 1 − α 0 q 0 dan hitung skalar β1 = q1T u1 . b. Selanjutnya hitung vektor r1 = u1 − β1q 1 = α 1q 2 è hitung α 1 =
r1 dan
q 2 = r1 / α 1 = r1 / r1 . 2. a. Dengan modal q 2 , hitung u 2 = Aq 2 − α 1q1 dan β 2 = q T2 u 2 . b. Selanjutnya hitung vektor r2 = u 2 − β 2 q 2 = α 2 q 3 è hitung α 2 = r2 dan
q 3 = r2 / α 2 = r2 / r2 . 3. demikian seterusnya sampai α 1 , α 2 , …, α M −1 , β1 , β 2 , …, β M , q1 , q 2 , …,
q M semuanya diperoleh è matrix T dan Q terbentuk. Nilai pribadi & vektor pribadi matrix tridiagonal TM ×M dicari dengan cara konvensional, misalnya dengan polinom-polinom yang membentuk barisan Sturm, diikuti metode belah-dua (lihat Bunjamin bab 2 & 5). Cara lain ialah dengan metode QL atau QR (Press dkk, Stoer & Bulirsch). Nilai pribadi T mendekati nilai pribadi A, sedangkan vektor pribadi y dari T memberi vektor Qy yang mendekati vektor pribadi A. Satu hal yang masih sedikit mengganjal dalam algoritma di atas ialah, bagaimana menentukan nilai M yang tepat, mengingat adanya fakta berikut. Nilai M yang ideal ialah = N, karena dengan demikian matrix tridiagonal TM ×M = TN×N
akan menghasilkan semua nilai pribadi dari A N× N . Namun hal ini tidak mungkin dicapai mengingat besarnya N, sehingga pemilihan M << N mutlak diperlukan, karena justru di sinilah letak kekuatan metode Lanczos. Di lain pihak, jika M terlalu kecil, nilai pribadi dari TM ×M meleset semua dari nilai pribadi A N× N , sehingga algoritma gagal. Hal ini disebabkan karena subruang Krylov K makin menjauhi sifat invarian jika M makin kecil, dan makin mendekati sifat invarian dengan makin besarnya M, asal q1 , q 2 , …, q M masih tetap saling bebas linear. Karena itu disarankan strategi sebagai berikut. Jika untuk suatu M tertentu, λ 0 ( M ) adalah nilai pribadi dominan dari T, maka prosedur Lanczos perlu diulang untuk berbagai nilai M demikian hingga untuk suatu bilangan ε > 0 kecil, syarat konvergensi
λ 0 ( M ) − λ 0 ( M − 1) <ε, λ 0 (M ) dipenuhi, sejauh q1 , q 2 , …, q M masih tetap saling bebas linear. Pengalaman mengajarkan bahwa penerapan “double precision” mutlak diperlukan agar sifat ortogonalitas Q T Q = I è Q T AQ = T bukan sekedar teori di atas kertas saja. Pengalaman juga mengajarkan bahwa metode Lanczos yang diterapkan pada A akan memberi nilai pribadi dominan. Jika aplikasi menuntut ditemukannya nilai pribadi dengan nilai mutlak terkecil dari A, metode Lanczos dapat diterapkan pada
A −1 dengan menerapkan strategi semacam metode pangkat invers dengan bantuan metode LU (lihat Bunjamin bab 1). Dari sini akan didapat nilai pribadi dominan dari
A −1 , dan invers dari nilai pribadi dominan dari A −1 ini = nilai pribadi dengan nilai mutlak terkecil dari A.
HASIL KOMPUTASI METODE LANCZOS PADA PC Berhubung dengan keterbatasan sarana komputer yang saya miliki, maka laporan ini membatasi diri pada matrix real simetrik berukuran N = 80 serta M dipilih = 6. Yang penting di sini ialah, terbukti bahwa sebagian dari nilai pribadi dan vektor pribadi matrix real simetrik A80×80 dapat diperoleh dari matrix tridiagonal simetrik
T6×6 . Adapun matrix A = { a ij } dipilih dengan acak simetrik sebagai berikut: a ij = INT(500 * RND) / 100, di mana RND adalah bila ngan acak ∈ (0, 1), dan dibuat a ij = a ji agar A simetrik. Hasil yang didapat cukup meyakinkan, sebagai berikut.
Dengan metode Lanczos, didapat nilai pribadi matrix real simetrik A80×80 dengan nilai mutlak terbesar=201.2344 dan dengan nilai mutlak terkecil=0.2091635, dan hasil ini didapat dari matrix T6×6 . Sebagai perbandingan dan pengujian efektivitas algoritma Lanczos, jika hasil ini dicari dengan metode pangkat dan pangkat invers (lihat Bunjamin bab 5) langsung pada matrix asli A80×80 , diperoleh hasil yang sama. Perhitungan vektor pribadi lewat metode Lanczos atau secara langsung juga memberi hasil sesuai harapan. Dari hasil di atas dapat diduga bahwa metode Lanczos akan sangat efektif untuk membantu mencari nilai pribadi matrix raksasa, asal didukung oleh sarana komputer yang cukup besar dan cepat untuk menyimpan dan mengolah matrix raksasa.
METODE ARNOLDI DAN DAVIDSON Sebagai tambahan informasi dapat dilaporkan bahwa selain metode Lanczos, juga dikenal beberapa metode lain untuk masalah spektral sejenis. Jika metode Lanczos tersebut di atas sesuai untuk masalah spektral untuk matrix raksasa yang simetrik, maka metode Arnoldi sesuai untuk matrix taksimetrik. Akibatnya, jika metode Lanczos menghasilkan matrix T yang tridiagonal simetrik, maka metode Arnoldi menghasilkan matrix T yang berbentuk Hessenberg atas (upper-Hessenberg). Selain itu ada juga metode Davidson yang populer di bidang kimia komputasi (Davidson sendiri seorang profesor kimia, tetapi “terpaksa” juga menekuni komputasi karena terdesak oleh kebutuhan di bidang kimia). Metode Davidson lebih sesuai untuk matrix simetrik dengan diagonal dominan, dengan memanfaatkan metode RitzGalerkin untuk mendapatkan matrix T yang berukuran kecil, tetapi penuh (jadi tidak tridiagonal dan tidak Hessenberg). Semua metode solusi tersebut di atas selama dasawarsa 90-an telah mengalami pengembangan pesat karena didesak oleh kebutuhan dalam sains & teknologi modern, termasuk pengembangan metode solusi masalah spektral umum, yaitu yang memenuhi persamaan Ax = λBx, di mana λ dan x adalah pasangan pribadi yang dicari, matrix A dan B berukuran raksasa dan simetrik, dan matrix B positif definit.
DAFTAR PUSTAKA 1. WILKINSON, J. H., The Algebraic Eigenvalue Problem, Clarendon Press, Oxford, (1965) 388
2. PRESS, W. H., B. P. FLANNERY, S. A. TEUKOLSKY, W. T. VETTERLING, Numerical Recipes, Cambridge University Press, Cambridge, (1992) 469 3. STOER, J., R. BULIRSCH, Introduction to Numerical Analysis, 2nd ., SpringerVerlag, New York, (1991) 362 4. LIN, H. Q., J. E. GUBERNATIS, Exact Diagonalization Methods for Quantum Systems, Comp. Phys, 7 (4) (Jul / Aug) (1993) 5. DAVIDSON, E. R., Monster Matrices: Their Eigenvalues and Eigenvectors, Comp. Phys, 7 (5) (Sep / Oct) (1993) 6. BOOTEN, A., H.V.D. VORST, Cracking Large-scale Eigenvalue Problems, Part I & II, Comp. Phys, 10 (3) (May / Jun) (1996) 7. BOOTEN, A., H.V.D. VORST, Cracking Large-scale Eigenvalue Problems, Part I & II, Comp. Phys, 10 (4) (July / Aug) (1996) 8. SONG HE, Stabilizing the Richardson Eigenvector Algorithm by Controlling Chaos, Comp. Phys, 11 (2) (Mar / Apr) (1997) 9. DESCLOUX, J., J.-L. FATTEBERT, F. GYGI, Rayleigh Quotient Iteration, an Old Recipe for Solving Large-scale Eigenvalue Problems, Comp. Phys, 12 (1) (Jan / Feb) (1998) 10. BUNJAMIN, M., Pengantar Pemodelan & Simulasi Komputer, Diktat Kuliah Komputasi, Jakarta, (2001)
HOME
KOMPUTASI DALAM SAINS DAN TEKNOLOGI NUKLIR XIII
