BerandaComputers and TechnologyMenerapkan LU dan Dekomposisi Cholesky di D

Menerapkan LU dan Dekomposisi Cholesky di D

Penulis: Dr Chibisi Chima-Okereke Dibuat: 26 November 2020 00:24:47 GMT Diterbitkan: 26 November 2020 00:24:47 GMT

Pengantar

Artikel ini berfokus pada implementasi Dekomposisi LU dan Cholesky menggunakan bahasa pemrograman D. Tujuan utama faktorisasi matriks ini adalah untuk mengubah persamaan linier ke dalam bentuk (Ax=b ), di mana (A in mathbb {R} ^ {n times n} ), (x in mathbb {R} ^ {n} ), dan (b in mathbb {R} ^ {n} ) ke dalam sistem segitiga yang dapat diselesaikan dengan substitusi maju atau mundur. Dalam kasus dekomposisi LU solver (Eliminasi Gaussian ), satu-satunya kendala adalah bahwa matriks (A ) memiliki peringkat penuh, meskipun dekomposisi itu sendiri dapat diterapkan ke matriks peringkat penuh dimensi (A in mathbb {R} ^ {m times n} ). Dekomposisi Cholesky dibatasi hanya untuk faktor matriks simetris persegi.

Referensi:

  • Algoritma untuk Eliminasi Gaussian dan Dekomposisi LU dari Algoritma 3 p65 dari Aljabar Linear Numerik dan Pengantar, oleh Holger Wendland.
  • PA=Dekomposisi LU, p51-9 Buku Pegangan Aljabar Linear, Edisi ke-2, Leslie Hogden.
  • Dekomposisi Cholesky , hal.78 Aljabar Linear Numerik Pengantar, oleh Holger Wendland.

Eliminasi Gaussian

Dekomposisi LU adalah formalisasi matriks Eliminasi Gaussian, yang beroperasi pada matriks yang diperbesar, penggabungan berdasarkan kolom (A | b ). Ini membuat segitiga atas di bagian (A ) dari matriks yang diperbesar dengan berulang kali mengalikan baris tertentu dengan faktor yang sesuai dan menguranginya dari baris lain. Ini menghilangkan koefisien dalam (A ) sedemikian rupa sehingga berakhir memiliki struktur segitiga atas yang sesuai. Pada akhir proses eliminasi, persamaan linier menjadi (Ux=c ), dimana (U ) merupakan matriks segitiga atas (dibentuk dari (A )) dan (c ) merupakan hasil dari (b ) telah mengambil bagian dalam proses eliminasi.

Sistem persamaan linier dapat dibuat dipersembahkan oleh:

$$ begin {matrix} a_ {11} x_ {1} & + & a_ {12} x_ {2} & + ldots & + & a_ {1n} x_ {n} &=& b_ {1} \
a_ {21} x_ {1} & + & a_ {22} x_ {2} & + ldots & + & a_ {2n} x_ {n} &=& b_ {2} \
vdots & + & vdots & + ldots & + & vdots &=& vdots \
a_ {n1} x_ {1} & + & a_ {n2} x_ {2} & + ldots & + & a_ {nn} x_ {n} &=& b_ {n} \
end {matrix} $$

Sejak (a_ {1, *} ) akan berakhir sebagai satu-satunya baris dengan koefisien lengkap, baris ini digunakan untuk menghilangkan koefisien pertama di baris lainnya. Iterasi baris (k=2, ldots, n ), kita mengalikan baris (1 ) dengan ( frac {a_ {k, 1}} {a_ {1,1}} ), ini membuat koefisien pertama dari baris pertama yang diubah sama dengan (a_ {k, 1} ). Baris yang diubah ini dikurangi dari baris (k ^ {th} ) yang menginduksi nol di posisi pertama baris itu. Kami kemudian mengambil baris asli (1 ) dan mengalikannya dengan ( frac {a_ {k + 1, 1}} {a_ {1, 1}} ) dan menguranginya dari baris (k + 1 ), terus melakukan ini sepanjang baris hingga baris (n ). Baris (2 ) kemudian menjadi baris faktor dan kita mengalikannya dengan ( frac {a_ {k, 2} ^ {2}} {a_ {2,2} ^ {2}} ), dan melaksanakan pengurangan seperti sebelumnya untuk (k=3, ldots, n ), proses ini diulang terus menerus sampai kita mendapatkan matriks segitiga atas (U ) dan vektor yang diubah terkait (c ). Yang dapat diselesaikan dengan substitusi belakang, topik dan implementasinya tercakup dalam artikel lain . Bagian (L ) dari dekomposisi LU terletak di subdiagonal matriks (A ^ {(n)} ) (di akhir algoritma). Eliminasi Gaussian adalah pemecah dan dekomposisi terkait adalah LU: (A=LU ).

Dalam artikel ini kami menggunakan alternatif nama implementasi backsubstitution tersebut di atas. Ada juga fungsi untuk pergantian pemain depan. Dalam kode di bawah ini, hanya tanda tangan yang ditampilkan di sini karena penerapan substitusi maju hanyalah kebalikan dari substitusi belakang. Alih-alih menamai pemecah kembali dan maju substitusi, penamaan telah didasarkan pada penyelesaian matriks segitiga atas / bawah masing-masing.

   enum  UPLO  {Atas, 

Bawah } alias Atas = UPLO . Atas ; alias Lebih Rendah = UPLO . Menurunkan ; T [] pengganti ( UPLO uplo , T

) (Matriks! ( T ) U , T [] y ) jika ( uplo == Atas ) { // … } T [] pengganti (

UPLO uplo , T

) (Matriks! ( T ) L , T [] y ) jika ( uplo == Bawah ) { // … }

Berputar sebagian

Di (k ^ { th} ) langkah Eliminasi Gaussian, persamaan (k ) harus dibagi dengan (a_ {k, k} ^ {(k)} ), dan jika (a_ {k, k} ^ {(k )} ) nol, Eliminasi Gaussian akan gagal dan jika mendekati nol akan menghasilkan kesalahan numerik dan menyebabkan algoritma menjadi tidak stabil. Untuk mengurangi ini, pada awal setiap set eliminasi untuk kolom, kita menukar baris (k ^ {th} ) saat ini dengan baris (l ) sehingga (l in {k, k + 1, ldots, n} ) dengan (| a_ {l, k} |=max | a_ {i, k} ^ {(k)} | ), artinya pilih baris yang memiliki nilai absolut terbesar di kolom (k ). Baris (l ) ditukar dengan baris (k ), yang secara nyata mengubah struktur matriks (A ^ {(k)} ). Untuk menjelaskan hal ini, matriks permutasi (P ) digunakan untuk merekam pertukaran baris. Ini dimulai sebagai matriks identitas, tetapi setiap baris di (A ^ {(k)} ) ditukar, baris terkait di (P ) juga ditukar. Dalam praktiknya, permutasi swap dicatat dalam vektor indeks yang dapat digunakan untuk mengubah (A ) sesuai kebutuhan. Dekomposisi LU kemudian menjadi (PA=LU ).

Dekomposisi Cholesky

Dekomposisi cholesky terkait dengan dekomposisi LU . Sekarang (A=LU ), dan (U ) adalah matriks segitiga atas yang mencakup diagonal ( (L ) adalah subdiagonal). Untuk matriks simetris (A ), (U=DL ^ {T} ), dan (A=LDL ^ {T} ), di mana (D ) adalah matriks diagonal. Dekomposisi Cholesky kemudian adalah (A=(LD ^ {1/2}) (LD ^ {1/2}) ^ {T}= tilde {L} tilde {L} ^ {T} ). Lebih formal (A=LL ^ {T} ) (jelas tidak sama (L ) seperti di (LU )). Dekomposisi Cholesky standar membutuhkan penggunaan akar kuadrat untuk (D ^ {1/2} ), yang berarti bahwa semua elemennya harus positif – matriks (A ) harus pasti positif. Dekomposisi alternatif adalah (A=LDL ) yang tidak memerlukan penggunaan akar kuadrat tetapi berarti bahwa semua elemen (D ) harus bukan nol karena mereka membagi elemen.

Implementasi dari faktorisasi matriks dijelaskan di atas sekarang mengikuti. Seperti biasa, algoritme ini untuk tujuan studi bukan untuk digunakan dalam aplikasi dunia nyata.

Menerapkan Dekomposisi LU

Penerapan dekomposisi LU sebagai pemecah diberikan di bawah ini. Ini memberi pengguna opsi untuk menimpa matriks (A ) dan (b ) dengan (U ) dan (c ).

  mobil  luSolve  

( bool menimpa =

Salah , T ) (Matriks! ( T ) SEBUAH, T [] b ) { tetap mobil n = SEBUAH. ncol ; // Opsi timpa statis jika ( menimpa ) { Matriks! ( T ) U = SEBUAH; T [] c = b ; } lain

{ Matriks! ( T ) U = SEBUAH. dup ; T [] c = b . dup ; } // Menghitung U dan c untuk setiap

( k ;

0 . . ( n 1 )) { T d yang tidak dapat diubah = 1 / U [k, k]; untuk setiap (saya ; ( k + 1 ) .. n ) { U [i, k] = U [i, k] d ; c [i] = c [i] c [k] U [i, k]; untuk setiap ( j ; ( k + 1 ) .. n ) { U [i, j] -= U [i, k] U [k, j]; } } } // Mendapatkan L mobil

L = mata ! ( T ) ( n ); untuk ( panjang saya = 0 ; i n ;

++ saya) { untuk ( panjang j = i ; j n ; ++ j ) { jika ( i == j ) terus ; L [j, i] = U [j, i]; // Menyetel subdiagonal ke nol U [j, i] = 0

; } } “U: n” . writeln ( U

); “L: n” . writeln ( L

); // Pemecah substitusi kembali kembali pengganti ! (Atas) ( U , c ); }

Dekomposisi LU dengan pivot parsial diberikan di bawah ini:

  mobil  lu  

( Metode LUAlgoritma , bool

menimpa = Salah , T ) ( ref Matriks ! ( T ) Ap ) { Matriks! ( T ) SEBUAH; statis jika

( menimpa ) { SEBUAH = Ap ; } lain

{ SEBUAH = Ap . dup ; } tetap mobil m = SEBUAH. sekarang ; tetap mobil n = SEBUAH. ncol ; mobil L = mata ! ( T ) ( m

); mobil U = nol ! ( T ) ( m , n ); / Dalam algoritma , poros disimpan di sebuah vektor / mobil p = seq ! ( size_t ) ( 0 , m 1 , 1 ); panjang r = 0 ; untuk ( panjang k = 0 ; k n ; ++ k ) { / Mendapatkan indeks untuk baris yang sesuai untuk diputar . Yang mana

fungsi mengembalikan indeks file       item menggunakan perbandingan berpasangan      /      mobil mu = yang ! ((Sebuah,  b 

)

=> abs (Sebuah)

> abs ( b )) ( SEBUAH [k] [r..$]); mu += r ; jika (SEBUAH [mu, k] != 0 ) { // Pertukaran baris mobil

tmp = SEBUAH [mu, k..$]; SEBUAH [mu, k..$] = SEBUAH [r, k..$]; SEBUAH [r, k..$] = tmp ; jika ( r > 0 ) { tmp = L [mu, 0..r]; L [mu, 0..r] = L [r, 0..r]; L [r, 0..r]

= tmp ; } // Tukar dengan vektor pivot mobil swp = p [mu]; p [mu] = p [r]; p [r]

= swp ; // Perhitungan untuk L dan U L [(r + 1)..$, r] = SEBUAH [(

r + 1)..$, k] / SEBUAH [r, k]; U [r, k..$] = SEBUAH [r, k..$]; untuk ( panjang saya = r + 1 ; saya m ; ++ saya) { untuk ( panjang j = k + 1 ; j n ; ++ j ) { SEBUAH [i, j] -=

L [i, r] U [r, j]; } } r += 1 ; } } / Memulihkan matriks permutasi dari vektor / mobil P = nol ! ( T ) ( m , m ); untuk ( panjang saya = 0 ; i m ; ++ saya) { mobil baris = baru T [m]; baris [] = 0 ; baris [ p [i]] = 1 ; P [i, 0..$] = matriks (baris, 1 , m ); } kembali LU ! ( T ) ( P , L , U ); }

Objek LU dalam kode di atas adalah:

   struct   LU  

( T ) jika ( isFloatingPoint ! (

T )) { Matriks! ( T ) P ; Matriks! ( T ) L ; Matriks! ( T ) U ; }

Implementasi Dekomposisi Cholesky

Di bawah ini adalah enumerasi yang menunjukkan dua varian dekomposisi Cholesky yang diimplementasikan:

   enum  Cholesky  {Vanila, 

LDL } alias Vanila = Cholesky . Vanila ; alias LDL = Cholesky . LDL ;

Dasar (Vanila) dekomposisi cholesky diberikan di bawah ini. Itu ditulis sedemikian rupa sehingga (L ) ditulis di segitiga bawah matriks.

  mobil  cholesky  

( Metode Cholesky , bool

menimpa = Salah , T ) (Matriks! ( T ) Ap ) jika (metode == Vanila) { tetap mobil n = Ap . ncol ; statis jika

(

menimpa ) { Matriks! ( T ) SEBUAH = Ap ; } lain

{ Matriks! ( T ) SEBUAH = Ap . dup ; } untuk ( panjang j = 0 ; j n ; ++ j ) { untuk ( panjang k = 0 ; k j ; ++ k ) { SEBUAH [j, j] -= SEBUAH [j, k] SEBUAH [j, k]; } // Akar kuadrat berpotensi menyebabkan masalah // seperti yang dibahas di atas SEBUAH [j, j] = sqrt (SEBUAH [j, j]); untuk ( panjang saya = ( j + 1 ); saya n ; ++ saya

{ untuk ( panjang k = 0 ; k j ; ++ k ) { SEBUAH [i, j] -= SEBUAH [i, k] SEBUAH [j, k]; } SEBUAH [i, j] /= SEBUAH [j, j]; } } // Mengumpulkan L untuk hasil Matriks ! ( T ) L = mata ! ( T ) ( n

; untuk setiap (saya ; 0. . n ) { untuk setiap ( j ; saya .. n ) { L [j, i] = SEBUAH [j, i]; } } kembali L ; }

Di bawah ini adalah LDL varian – diadaptasi dari kode Eliminasi Gaussian:

  mobil  cholesky  

( Metode Cholesky , bool

menimpa = Salah , T ) (Matriks! ( T ) SEBUAH) jika (metode == LDL ) { tetap mobil n = SEBUAH. ncol ; statis jika

(

menimpa ) { Matriks! ( T ) U = SEBUAH; } lain

{ Matriks! ( T ) U = SEBUAH. dup ; } T [] d = baru T [n]; // Mirip dengan kode di luSolve untuk setiap

( k ;

0 . . ( n 1 )) { // Kumpulkan diagonal d [k] = U [k, k]; untuk setiap (saya ; ( k + 1 ) .. n ) { U [i, k] = U [i, k] / d [k]; untuk setiap ( j ; ( k + 1 ) .. n ) { U [i, j] -= U [i, k] U [k, j]; } } } d [$ 1] = U [$ 1, $ 1]; // Kumpulkan elemen untuk L mobil

L = mata ! ( T ) ( n ); untuk ( panjang saya = 0 ; i n ; ++ saya) { untuk ( panjang j = i ; j n ; ++ j ) { jika ( i == j ) terus ; L [j, i] = U [j, i]; } } / Ini adalah Sebuah “tipe voldermort” di D . Itu adalah tidak tersedia untuk konstruksi di luar

ini

fungsi . / struct Cholesky

(

T ) { Matriks! ( T ) L ; T [] d ; } kembali Cholesky ! ( T ) ( L , d

); }

Beberapa keluaran

Di sini beberapa nomor dimasukkan melalui algoritma yang ditulis untuk tujuan demonstrasi. Data yang dihasilkan secara acak:

  mobil SEBUAH = createRandomMatrix ! ( dua kali lipat ) (

7 , 7 ); mobil x = createRandomMatrix ! ( dua kali lipat ) ( 7 , 1 ); mobil b = gemm (SEBUAH, x ). Himpunan ; “Sebuah” . writeln (SEBUAH

); “x: n” . writeln ( x

); “b: n” . writeln ( b

);

keluaran:

 SEBUAH: Matriks (7 x 7) 0.78760895 0.16791495 0.51998612 0.90592727 0.35697508 0.89517028 0.63818117 0,94644681 0,65230649 0,70400738 0,74370928 0,93061186 0,95699091 0,51571648 0.7791607 0.021011612 0.95554588 0.88994508 0.54347486 0.20869079 0.67902032 0.30326116 0.76147071 0.22517496 0.36843291 0.9923756 0.13921349 0.30278154 0.85278273 0.15511374 0.053926307 0.6661747 0.94426026 0.8441975 0.68269395 0,082958874 0,17310216 0,39019568 0,82435547 0,93347347 0,094591039 0,95515174 0,25554523 0,69263844 0,52331036 0,32446879 0,45709359 0,32742983 0,49883607  x: Matriks (7 x 1) 0.40713286 0.45848681 0.79254078 0.84800791 0.17773008 0,59358732 0,060270934  b: [2.21126, 2.63757, 2.10024, 1.24075, 1.73605, 1.40107, 1.41717]  

Selesaikan selanjutnya menggunakan Gaussian Elimination:

  mobil xOut = luSolve  (

SEBUAH, b ); “x: n” . writeln ( x . Himpunan ); “Gaussian Elimination Solver: n” . writeln ( xOut );

keluaran:

  U: Matriks (7 x 7) 0.7876 0.1679 0.52 0.9059 0.357 0.8952 0.6382 0 0,4505 0,07916 -0,3449 0,5016 -0,1187 -0,2512 0 0 0.4666 -0.1174 0.3519 -0.7151 -0.03321 0 0 0 0,5286 0,1525 -0,1712 0,4386 0 0 0 0 1,101 -1,055 0,3243 0 0 0 0 0 1.027 0.1629 0 0 0 0 0 0 0,3834  L: Matriks (7 x 7) 1 0 0 0 0 0 0 1.202 1 0 0 0 0 0 0,9893 -0,3221 1 0 0 0 0 0,385 1,547 -0,2089 1 0 0 0 1.083 -0.05925 -1.081 -0.8741 1 0 0 0,1053 0,345 0,6603 1,751 0,2029 1 0 0,3245 ​​1,416 0,5196 1,097 -0,6534 0,07343 1  x: [0.407133, 0.458487, 0.792541, 0.848008, 0.17773, 0.593587, 0.0602709] Pemecah Eliminasi Gaussian: [0.407133, 0.458487, 0.792541, 0.848008, 0.17773, 0.593587, 0.0602709]  

(PA=LU ) dekomposisi:

  mobil plu = SEBUAH. 

lu ! ( Produk Luar ); “PLU decomp:” . writeln ( plu ); impor std.algorithm : mengurangi; // Membandingkan elemen PA dengan LU “Kesalahan PLU:” . writeln ( mapply ! (( a1 , a2 ) => pow ( a1 a2 , 2.0

)) (

gemm ( plu . L , plu . U ). Himpunan , gemm ( plu . P , SEBUAH). Himpunan ) . mengurangi ! (( a1 , a2 ) => a1 + a2 ) . sqrt ) ;

keluaran:

  PLU decomp: LU! Double (Matrix (7 x 7) 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 , Matriks (7 x 7) 1 0 0 0 0 0 0 0,3204 1 0 0 0 0 0 0,901 -0,7831 1 0 0 0 0 0,08765 0,2098 -0,5658 1 0 0 0 0.8322 -0.6786 0.1139 0.4626 1 0 0 0,8232 -0,934 -0,6468 0,5876 -0,4156 1 0 0,27 0,9349 -0,5745 0,07403 0,2895 -0.1709 1 , Matriks (7 x 7) 0,9464 0,6523 0,704 0,7437 0,9306 0,957 0,5157 0 0.5525 -0.0004036 0.1301 0.6942 -0.1674 0.1375 0 0 -0.5807 0.09797 0.6494 -0.1492 0.3257 0 0 0 0.7873 1.074 -0.03858 1.065 0 0 0 0 -0,517 0,02 -0,2276 0 0 0 0 0 -0,801 -0,127 0 0 0 0 0 0 0,3834 ) Kesalahan PLU: 3.55513e-16  

Keluaran untuk Cholesky, pertama pemecah “vanilla”:

  SEBUAH = createRandomMatrix ! ( dua kali lipat ) (  7 

, 7 ); // Buat matriks acak menjadi simetris SEBUAH = gemm ! ( CblasTrans , CblasNoTrans ) (SEBUAH, SEBUAH); “Sebuah” . writeln (SEBUAH

); mobil L = cholesky ! (Vanila) (SEBUAH

); “Cholesky: n” . writeln ( L

); mobil calcA = gemm ! ( CblasNoTrans , CblasTrans ) ( L , L ); “LL ^ T: n” . writeln ( calcA

; // Membandingkan elemen A dengan LL ^ T “Kesalahan Cholesky (Vanilla):” . writeln ( peta! (( x1 , x2 ) => abs ( x1 x2 )) (SEBUAH. Himpunan , calcA . Himpunan ) . mengurangi ! (( x1 , x2 ) => x1 + x2 ) );

keluaran:

  Cholesky: Matriks (7 x 7) 1,73 0 0 0 0 0 0 1,327 0,3089 0 0 0 0 0 1,239 -0,1535 0,5974 0 0 0 0 1.804 0.4317 0.1312 0.6613 0 0 0 1.1 0.1185 0.4336 -0.3596 0.3032 0 0 1.529 0.2632 0.1535 0.5298 0.1423 0.3863 0 1.501 0.4592 -0.1566 0.5664 0.4249 -0.3999 0.1552  LL ^ T: Matriks (7 x 7) 2.9932443 2.2959181 2.1435358 3.120934 1.9031914 2.6449977 2.5977037 2.2959181 1.8564814 1.5967375 2.5272337 1.4964281 2.1101046 2.1343782 2.1435358 1.5967375 1.9155117 2.2470581 1.6037649 1.9454282 1.696221 3.120934 2.5272337 2.2470581 3.8950476 1.8545822 3.2419384 3.260778 1.9031914 1.4964281 1.6037649 1.8545822 1.6334361 1.6321453 1.5633707 2.6449977 2.1101046 1.9454282 3.2419384 1.6321453 2.880207 2.598331 2.5977037 2.1343782 1.696221 3.260778 1.5633707 2.598331 3.1751308  Kesalahan Cholesky (Vanilla): 2.44249e-15  

lalu LDL

pemecah:

   // Kami mendapatkan matriks A yang sama dari penghitungan LDL ^ T     "LDL:  n" 

. writeln (SEBUAH. cholesky ! ( LDL ));

keluaran:

  A=LDL ^ T: Matriks (7 x 7) 2.9932443 2.2959181 2.1435358 3.120934 1.9031914 2.6449977 2.5977037 2.2959181 1.8564814 1.5967375 2.5272337 1.4964281 2.1101046 2.1343782 2.1435358 1.5967375 1.9155117 2.2470581 1.6037649 1.9454282 1.696221 3.120934 2.5272337 2.2470581 3.8950476 1.8545822 3.2419384 3.260778 1.9031914 1.4964281 1.6037649 1.8545822 1.6334361 1.6321453 1.5633707 2.6449977 2.1101046 1.9454282 3.2419384 1.6321453 2.880207 2.598331 2.5977037 2.1343782 1.696221 3.260778 1.5633707 2.598331 3.1751308  // Ini adalah objek LDL yang memegang matriks bawah dan diagonal d LDL: Cholesky! Dobel (Matriks (7 x 7) 1 0 0 0 0 0 0 0,767 1 0 0 0 0 0 0.7161 -0.4969 1 0 0 0 0 1.043 1.398 0.2196 1 0 0 0 0,6358 0,3837 0,7258 -0,5438 1 0 0 0.8837 0.8519 0.2569 0.801 0.4695 1 0 0.8679 1.486 -0.2622 0.8564 1.401 -1.035 1 , [2.99324, 0.0954357, 0.356905, 0.437381, 0.0919369, 0.149217, 0.0240929])  

Itu dia. Terima kasih telah membaca!

Read More

RELATED ARTICLES

LEAVE A REPLY

Please enter your comment!
Please enter your name here

Most Popular

Recent Comments