Foreplay

Solusi numerik sistem persamaan linear merupakan langkah yang sangat diperlukan dalam banyak bidang matematika terapan, teknik, dan industri TI, apakah itu bekerja dengan grafik, menghitung aerodinamika pesawat terbang, atau mengoptimalkan logistik. "Mesin" yang sekarang modis juga tidak akan banyak berkembang tanpa itu. Selain itu, solusi sistem linear, sebagai aturan, memakan persentase terbesar dari semua biaya komputasi. Jelas bahwa komponen penting ini memerlukan optimasi kecepatan maksimum.

Seringkali bekerja dengan apa yang disebut matriks jarang - mereka yang nol adalah urutan besarnya lebih besar dari nilai-nilai lainnya. Ini, misalnya, tidak bisa dihindari jika Anda berurusan dengan persamaan diferensial parsial atau beberapa proses lain di mana elemen yang muncul dalam hubungan linear mendefinisikan hanya terhubung dengan "tetangga". Berikut adalah contoh yang mungkin dari matriks jarang untuk persamaan Poisson satu dimensi yang dikenal dalam fisika klasik $$$- \ phi ^ {''} = f$ pada kisi yang seragam (ya, sementara tidak ada banyak nol di dalamnya, tetapi ketika menggiling kisi itu akan sehat!):

$$

$$\ begin {pmatrix} 2 & -1 & 0 & 0 & 0 \\ -1 & 2 & -1 & 0 & 0 \\ 0 & -1 & 2 & -1 & 0 \\ 0 & 0 & -1 & 2 & -1 \\ 0 & 0 & 0 & -1 & 2 \ end {pmatrix}$$

Matriks yang berlawanan dengan mereka - matriks yang tidak memperhatikan jumlah nol dan memperhitungkan semua komponen tanpa kecuali - disebut padat.

Matriks jarang sering, karena berbagai alasan, disajikan dalam format kolom terkompresi - Compressed Sparse Column (CSC). Deskripsi ini menggunakan dua array integer dan satu floating point. Biarkan matriks memiliki segalanya $$$nnz_A$ bukan nol dan $$$N$ kolom. Elemen matriks terdaftar dalam kolom dari kiri ke kanan, tanpa pengecualian. Array pertama $$$i_A$ panjangnya $$$nnz_A$ berisi nomor baris komponen matriks bukan nol. Array kedua $$$j_A$ panjangnya $$$N + 1$ mengandung ... tidak, bukan nomor kolom, karena kemudian matriks akan ditulis dalam format koordinat (koordinat), atau ternary (triplet). Dan array kedua berisi nomor seri komponen matriks yang digunakan untuk memulai kolom, termasuk kolom dummy tambahan di bagian akhir. Akhirnya, array ketiga $$$v_A$ panjangnya $$$nnz_A$ sudah mengandung komponen itu sendiri, dalam urutan yang sesuai dengan array pertama. Di sini, misalnya, dengan asumsi bahwa penomoran baris dan kolom matriks dimulai dari nol, untuk matriks tertentu

$$

$$A = \ begin {pmatrix} 0 & 1 & 0 & 4 & -1 \\ 7 & 0 & 2.1 & 0 & 3 \\ 0 & 0 & 0 & 10 & 0 \ end {pmatrix}$$

array akan menjadi $$$i_A = \ {1, 0, 1, 0, 2, 0, 1 \}$ , $$$j_A = \ {0, 1, 2, 3, 5, 7 \}$ , $$$v_A = \ {7, 1, 2.1, 4, 10, -1, 3 \}$ .

Metode untuk memecahkan sistem linear dibagi menjadi dua kelas besar - langsung dan berulang. Garis lurus didasarkan pada kemungkinan merepresentasikan matriks sebagai produk dari dua matriks yang lebih sederhana, untuk kemudian memecah solusi menjadi dua tahap sederhana. Iteratif menggunakan sifat unik dari ruang linear dan bekerja pada yang kuno sebagai metode dunia pendekatan yang tidak diketahui untuk solusi yang diinginkan, dan dalam proses konvergensi itu sendiri, matriks digunakan, sebagai suatu peraturan, hanya untuk memperbanyaknya dengan vektor. Metode berulang jauh lebih murah daripada yang langsung, tetapi bekerja lambat pada matriks yang tidak terkondisikan. Ketika keandalan beton bertulang itu penting - gunakan garis lurus. Di sini saya di sini dan saya ingin menyentuh sedikit.

Katakanlah untuk matriks persegi $$$A$ dekomposisi formulir tersedia bagi kami $$$A = LU$ dimana $$$L$ dan $$$U$ - masing-masing, matriks segitiga bawah dan segitiga atas, masing-masing. Yang pertama berarti memiliki nol di atas diagonal, yang kedua berarti di bawah diagonal. Bagaimana tepatnya kita mendapatkan dekomposisi ini - kita tidak tertarik sekarang. Berikut ini adalah contoh dekomposisi sederhana:

$$

$$\ begin {pmatrix} 1 & -1 & -1 \\ 2 & - 1 & -0.5 \\ 4 & -2 & -1.5 \ end {pmatrix} = \ begin {pmatrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ 4 & 2 & 1 \ end {pmatrix} \ begin {pmatrix} 1 & -1 & -1 \\ 0 & 1 & 1,5 \\ 0 & 0 & -0 & \0,5 \ end {pmatrix}$$

Bagaimana dalam hal ini menyelesaikan sistem persamaan $$$Ax = f$ misalnya dengan sisi kanan $$$f = \ begin {pmatrix} 4 \\ 2 \\ 3 \ end {pmatrix}$ ? Tahap pertama adalah langkah langsung (maju memecahkan = substitusi maju). Kami menunjukkan $$$y: = Ux$ dan bekerja dengan sistem $$$Ly = f$ . Karena $$$L$ - segitiga bawah, berturut-turut dalam siklus kami menemukan semua komponen $$$y$ atas ke bawah:

$$

$$\ begin {pmatrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ 4 & 2 & 1 \ end {pmatrix} \ begin {pmatrix} y_1 \\ y_2 \\ y_3 \ end {pmatrix} = \ begin {pmatrix} 4 \\ 2 \\ 3 \ end {pmatrix} \ menyiratkan y = \ mulai {pmatrix} 4 \\ -6 \\ -1 \ end {pmatrix}$$

Gagasan utamanya adalah, setelah menemukan $$$i$ komponen vektor $$$y$ itu dikalikan dengan kolom dengan nomor matriks yang sama $$$L$ yang kemudian dikurangi dari sisi kanan. Matriks itu sendiri tampaknya runtuh dari kiri ke kanan, berkurang dalam ukuran karena semakin banyak komponen vektor ditemukan $$$y$ . Proses ini disebut eliminasi kolom.

Tahap kedua adalah langkah mundur (mundur memecahkan = substitusi mundur). Ditemukan vektor $$$y$ kami yang memutuskan $$$Ux = y$ . Di sini kita sudah membahas komponen dari bawah ke atas, tetapi idenya tetap sama: $$$i$ Kolom th dikalikan dengan komponen yang baru saja ditemukan $$$x_i$ dan bergerak ke kanan, dan matriks itu runtuh dari kanan ke kiri. Seluruh proses diilustrasikan untuk matriks yang disebutkan dalam contoh dalam gambar di bawah ini:

$$

$$\ small \ begin {pmatrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ 4 & 2 & 1 \ end {pmatrix} \ begin {pmatrix} y_1 \\ y_2 \\ y_3 \ end {pmatrix} = \ begin {pmatrix} 4 \\ 2 \\ 3 \ end {pmatrix} \ menyiratkan \ begin {pmatrix} 1 & 0 \\ 2 & 1 \ end {pmatrix} \ begin {pmatrix} y_2 \\ y_3 \ end {pmatrix } = \ begin {pmatrix} -6 \\ -13 \ end {pmatrix} \ menyiratkan \ begin {pmatrix} 1 \ end {pmatrix} \ begin {pmatrix} y_3 \ end {pmatrix} = \ begin {pmatrix} -1 \ end {pmatrix}$$

$$

$$\ small \ begin {pmatrix} 1 & -1 & -1 \\ 0 & 1 & 1.5 \\ 0 & 0 & -0.5 \ end {pmatrix} \ begin {pmatrix} x_1 \\ x_2 \\ x_3 \ end { pmatrix} = \ begin {pmatrix} 4 \\ -6 \\ -1 \ end {pmatrix} \ Rightarrow \ begin {pmatrix} 1 & -1 \\ 0 & 1 \ end {pmatrix} \ begin {pmatrix} x_1 \ \ x_2 \ end {pmatrix} = \ begin {pmatrix} 6 \\ -9 \ end {pmatrix} \ menyiratkan \ begin {pmatrix} 1 \ end {pmatrix} \ begin {pmatrix} x_1 \ end {pmatrix} = \ begin {pmatrix} -3 \ end {pmatrix}$$

dan keputusan kita akan menjadi $$$x = \ begin {pmatrix} -3 \\ - 9 \\ 2 \ end {pmatrix}$ .

Jika matriks padat, artinya, ia sepenuhnya direpresentasikan dalam bentuk semacam array, satu dimensi atau dua dimensi, dan akses ke elemen tertentu di dalamnya berlangsung dari waktu ke waktu. $$$O (1)$ , maka prosedur solusi serupa dengan dekomposisi yang ada tidak sulit dan mudah dikodekan, jadi kami bahkan tidak akan membuang waktu untuk itu. Bagaimana jika matriksnya jarang? Di sini, pada prinsipnya, juga tidak sulit. Berikut adalah kode C ++ untuk langkah maju di mana solusinya $$$x$ ditulis untuk tempat di sisi kanan, tanpa memeriksa kebenaran data input (array CSC sesuai dengan matriks segitiga bawah):

Algoritma 1:

void forward_solve (const std::vector<size_t>& iA, const std::vector<size_t>& jA, const std::vector<double>& vA, std::vector<double>& x) { size_t j, p, n = x.size(); for (j = 0; j < n; j++) //    { x[j] /= vA[jA[j]]; for (p = jA[j]+1; p < jA[j+1]; p++) x[iA[p]] -= vA[p] * x[j] ; } } 

Semua diskusi lebih lanjut hanya akan membahas langkah maju untuk menyelesaikan sistem segitiga bawah $$$Lx = f$ .

Dasi

Tetapi bagaimana jika sisi kanan, mis. vektor di sebelah kanan tanda sama dengan $$$Lx = f$ Apakah itu sendiri memiliki jumlah nol yang besar? Maka masuk akal untuk melewati perhitungan yang terkait dengan posisi nol. Perubahan kode dalam kasus ini minimal:

Algoritma 2:

 void fake_sparse_forward_solve (const std::vector<size_t>& iA, const std::vector<size_t>& jA, const std::vector<double>& vA, std::vector<double>& x) { size_t j, p, n = x.size(); for (j = 0; j < n; j++) //    { if(x[j]) { x[j] /= vA[jA [j]]; for (p = jA[j]+1; p < jA[j+1]; p++) x[iA[p]] -= vA[p] * x[j]; } } } 

Satu-satunya hal yang kami tambahkan adalah if , yang tujuannya adalah untuk mengurangi jumlah operasi aritmatika ke angka aktualnya. Jika, misalnya, seluruh sisi kanan terdiri dari nol, maka Anda tidak perlu menghitung apa pun: solusinya adalah sisi kanan.

Semuanya tampak hebat dan tentu saja itu akan berhasil, tetapi di sini, setelah foreplay yang lama, masalahnya akhirnya terlihat - kinerja pemecah ini secara asimptot dalam kasus sistem besar. Dan semua karena fakta memiliki for loop. Apa masalahnya? Bahkan jika kondisi di dalam if ternyata benar sangat jarang, tidak ada jalan keluar dari siklus itu sendiri, dan ini menciptakan kompleksitas algoritma $$$O (N)$ dimana $$$N$ Apakah ukuran matriks. Tidak peduli bagaimana loop dioptimalkan oleh kompiler modern, kompleksitas ini akan membuatnya terasa besar $$$N$ . Terutama ofensif ketika seluruh vektor $$$f$ seluruhnya terdiri dari nol, karena, seperti yang kami katakan, maka jangan lakukan apa pun dalam kasus ini! Bukan operasi aritmatika tunggal! Apa-apaan ini $$$O (N)$ aksi?

Baiklah, katakan saja. Meski begitu, mengapa Anda tidak bisa hanya bertahan for menjalankan, karena perhitungan nyata dengan bilangan real, mis. mereka yang jatuh if masih sedikit? Tetapi kenyataannya adalah bahwa prosedur aliran langsung ini dengan sisi kanan yang jarang digunakan sendiri sering digunakan dalam siklus eksternal dan mendasari dekomposisi Cholesky. $$$A = LL ^ T$ dan dekomposisi LU yang tampak kiri. Ya, ya, salah satu dari ekspansi itu, tanpa kemampuan untuk melakukan semua gerakan langsung dan mundur dalam menyelesaikan sistem linear, kehilangan semua makna praktis.

Teorema Jika matriks adalah symmetric positive definite (SPD), maka dapat direpresentasikan sebagai $$$A = LL ^ T$ satu-satunya jalan kemana $$$L$ - matriks segitiga bawah dengan elemen positif di diagonal.

Untuk matriks SPD yang sangat jarang, digunakan dekomposisi Cholesky. Secara skematis merepresentasikan dekomposisi dalam bentuk matriks-blok

$$

$$\ begin {pmatrix} \ mathbf {L} _ {11} & \ mathbf {0} \\ \ mathbf {l} _ {12} ^ T & l_ {22} \ end {pmatrix} \ begin {pmatrix} \ mathbf {L} _ {11} ^ T & \ mathbf {l} _ {12} \\ \ mathbf {0} & l_ {22} \ end {pmatrix} = \ begin {pmatrix} \ mathbf {A} _ { 11} & \ mathbf {a} _ {12} \\ \ mathbf {a} _ {12} ^ T & a_ {22} \ end {pmatrix},$$

keseluruhan proses faktorisasi dapat secara logis dibagi menjadi hanya tiga langkah.

Algoritma 3:

1. metode atas dekomposisi Cholesky dimensi yang lebih kecil $$$\ mathbf {L} _ {11} \ mathbf {L} _ {11} ^ T = \ mathbf {A} _ {11}$ (rekursi!)
2. sistem segitiga bawah dengan sisi kanan jarang $$$\ mathbf {L} _ {11} \ mathbf {l} _ {12} = \ mathbf {a} _ {12}$ (ini dia! )
3. perhitungan $$$l_ {22} = \ sqrt {a_ {22} - \ mathbf {l} _ {12} ^ T \ mathbf {l} _ {12}}$ (operasi sepele)

Dalam praktiknya, ini dilaksanakan sedemikian rupa sehingga langkah 3 dan 2 dilakukan dalam satu siklus besar, dan dalam urutan ini. Dengan demikian, matriks dijalankan secara diagonal dari atas ke bawah, meningkatkan matriks $$$L$ baris demi baris pada setiap iterasi loop.

Jika dalam agoritme seperti itu kompleksitas ke depan pada langkah 2 akan $$$O (N)$ dimana $$$N$ Adalah ukuran matriks segitiga bawah $$$\ mathbf {L} _ {11}$ pada iterasi sembarang dari siklus besar, maka kompleksitas seluruh dekomposisi akan setidaknya $$$O (N ^ 2)$ ! Oh, aku tidak akan menyukainya!

Canggih

Banyak algoritma yang berdasarkan pada meniru tindakan manusia dalam menyelesaikan masalah. Jika Anda memberi seseorang sistem linear segitiga bawah dengan sisi kanan, di mana hanya 1-2 yang tidak nol, maka tentu saja ia pertama-tama menjalankan vektor sisi kanan dengan matanya dari atas ke bawah (siklus kerumitan terkutuk itu) $$$O (N)$ ! ) untuk menemukan nonzeros ini. Maka dia akan menggunakan hanya mereka, tanpa membuang waktu pada komponen nol, karena yang terakhir tidak akan mempengaruhi solusi: tidak masuk akal untuk membagi nol menjadi komponen diagonal dari matriks, serta memindahkan kolom dikalikan dengan nol ke kanan. Ini adalah algoritma di atas 2. Tidak ada keajaiban. Tetapi bagaimana jika seseorang segera diberi nomor komponen bukan nol dari beberapa sumber lain? Misalnya, jika sisi kanan adalah kolom dari beberapa matriks lain, seperti halnya dengan dekomposisi Cholesky, maka di sana kami memiliki akses instan ke komponen bukan nolnya, jika diminta secara berurutan:

 //     j,     : for (size_t p = jA[j]; p < jA[j+1]; p++) //    vA[p]     iA[p]   j 

Kompleksitas akses ini adalah $$$O (nnz_j)$ dimana $$$nnz_j$ Apakah jumlah komponen bukan nol dalam kolom j. Terima kasih Tuhan untuk format CSC! Seperti yang Anda lihat, itu tidak hanya digunakan untuk menghemat memori.

Dalam hal ini, kita dapat menangkap esensi dari apa yang terjadi ketika menyelesaikan dengan metode kursus langsung $$$Lx = f$ untuk matriks jarang $$$L$ dan sisi kanan $$$f$ . Tahan napas: kami mengambil komponen yang bukan nol $$$f_i$ di sisi kanan, kami menemukan variabel yang sesuai $$$x_i$ membaginya dengan $$$L_ {ii}$ dan kemudian, mengalikan seluruh kolom ke-i dengan variabel yang ditemukan ini, kami memperkenalkan non-nol tambahan di sisi kanan dengan mengurangi kolom ini dari sisi kanan! Proses ini dijelaskan dengan baik dalam bahasa ... grafik. Apalagi berorientasi dan non-siklik.

Untuk matriks segitiga bawah yang tidak memiliki nol pada diagonal, kami mendefinisikan grafik yang terhubung. Kami menganggap bahwa penomoran baris dan kolom dimulai dari nol.
Definisi Grafik konektivitas untuk matriks segitiga bawah $$$L$ ukurannya $$$N$ , yang tidak memiliki nol pada diagonal, disebut himpunan node $$$V = \ {0, ..., N-1 \}$ dan tepi berorientasi $$$E = \ {(j, i) | L_ {ij} \ ne 0 \}$ .

Di sini, misalnya, adalah seperti apa grafik koneksi untuk matriks segitiga yang lebih rendah.

di mana angka-angka hanya sesuai dengan nomor urut elemen diagonal, dan titik-titik menunjukkan komponen tidak nol di bawah diagonal:



Definisi Grafik berorientasi reachability $$$G$ pada beberapa indeks $$$W \ subset V$ panggil begitu banyak indeks $$$R_G (W) \ subset V$ bahwa dalam indeks apa pun $$$z \ dalam R_G (W)$ Anda bisa mendapatkan dengan melewati grafik dari beberapa indeks $$$w (z) \ dalam W$ .

Contoh: untuk grafik dari gambar $$$R_G (\ {0, 4 \}) = \ {0, 4, 5, 6 \}$ . Jelas bahwa itu selalu berlaku $$$W \ subset R_G (W)$ .

Jika kami mewakili setiap simpul grafik sebagai nomor kolom dari matriks yang menghasilkannya, maka simpul tetangga yang ujung-ujungnya diarahkan sesuai dengan nomor baris komponen bukan nol dalam kolom ini.

Biarkan $$$nnz (x) = \ {j | x_j \ ne 0 \}$ menunjukkan himpunan indeks yang sesuai dengan posisi bukan nol dalam vektor x.

Teorema Hilbert (tidak, tidak dengan nama yang diberi nama spasi) $$$nnz (x)$ dimana $$$x$ ada vektor solusi dari sistem segitiga bawah yang jarang $$$Lx = f$ dengan sisi kanan jarang $$$f$ bertepatan dengan $$$R_G (nnz (f))$ .

Tambahan kita sendiri: dalam teorema kita tidak memperhitungkan kemungkinan yang tidak mungkin bahwa angka bukan nol di sisi kanan, ketika menghancurkan kolom, dapat dikurangi menjadi nol, misalnya, 3 - 3 = 0. Efek ini disebut pembatalan numerik. Memperhatikan nol yang terjadi secara spontan seperti itu adalah buang-buang waktu, dan mereka dianggap sebagai semua angka lain di posisi nol.

Metode yang efektif untuk melakukan perpindahan langsung dengan sisi kanan yang jarang diberikan, dengan asumsi bahwa kita memiliki akses langsung ke komponen-komponennya yang bukan nol menurut indeks, adalah sebagai berikut.

1. Kami menjalankan melalui grafik menggunakan metode "pencarian pertama kedalaman", berurutan dimulai dari indeks $$$i \ in nnz (f)$ setiap komponen bukan nol di sisi kanan. Menulis node yang ditemukan ke array $$$R_G (nnz (f))$ saat melakukan ini dalam urutan di mana kita kembali ke grafik! Dengan analogi dengan pasukan penjajah: kami menduduki negara tanpa kekejaman, tetapi ketika kami mulai diusir kembali, kami, marah, menghancurkan segala sesuatu di jalannya.
   Perlu dicatat bahwa sama sekali tidak perlu bahwa daftar indeks $$$nnz (f)$ diurutkan dengan menaik ketika diumpankan ke input dari algoritma "pencarian pertama-kedalaman". Anda dapat memulai dalam urutan apa pun di set $$$nnz (f)$ . Urutan berbeda milik set $$$nnz (f)$ indeks tidak mempengaruhi solusi akhir, seperti yang akan kita lihat pada contoh.
   
   Langkah ini tidak memerlukan pengetahuan tentang array nyata. $$$v_A$ ! Selamat datang di dunia analisis simbolik ketika bekerja dengan pemecah jarang langsung!
   
2. Kami melanjutkan ke solusi itu sendiri, yang siap membantu kami $$$R_G (nnz (f))$ dari langkah sebelumnya. Kami menghancurkan kolom dalam urutan terbalik dari catatan array ini . Voila!

Contoh

Pertimbangkan contoh di mana kedua langkah diperlihatkan secara rinci. Misalkan kita memiliki matriks 12 x 12 dengan bentuk berikut:

Grafik konektivitas yang sesuai memiliki bentuk:

Biarkan bukan nol di bagian kanan hanya pada posisi 3 dan 5, yaitu $$$nnz (f) = \ {3, 5 \}$ . Mari kita jalankan melalui grafik, mulai dari indeks ini dalam urutan tertulis. Maka metode "pencarian dalam" akan terlihat sebagai berikut. Pertama kita akan mengunjungi indeks secara berurutan $$$3 \ hingga 8 \ hingga 11$ , tidak lupa untuk menandai indeks yang dikunjungi. Pada gambar di bawah mereka diarsir dengan warna biru:

Ketika kembali, kami memasukkan indeks ke dalam array kami indeks komponen solusi nol $$$nnz (x): = \ {\ color {blue} {11}, \ color {blue} 8, \ color {blue} 3 \}$ . Selanjutnya, coba jalankan $$$5 \ to8 \ hingga ...$ , tapi kami menemukan simpul 8 yang sudah ditandai, jadi kami tidak menyentuh rute ini dan pergi ke cabang $$$5 \ to9 \ hingga ...$ . Hasil dari proses ini adalah $$$5 \ to9 \ to10$ . Kami tidak dapat mengunjungi Node 11, karena sudah diberi label. Jadi, kembali dan lengkapi array $$$nnz (x)$ tangkapan baru yang ditandai dengan warna hijau: $$$nnz (x): = \ {\ color {blue} {11}, \ color {blue} 8, \ color {blue} 3, \ color {green} {10}, \ color {green} 9, \ color {green} 5 \}$ . Dan inilah gambarnya:

Node hijau berlumpur 8 dan 11 adalah yang ingin kami kunjungi selama run kedua, tetapi tidak bisa, karena sudah dikunjungi selama yang pertama.

Jadi susunannya $$$nnz (x) = R_G (nnz (f))$ terbentuk. Kami lolos ke langkah kedua. Langkah sederhana: kita jalankan melalui array $$$nnz (x)$ dalam urutan terbalik (dari kanan ke kiri), menemukan komponen yang sesuai dari vektor solusi $$$x$ pembagian oleh komponen diagonal dari matriks $$$L$ dan memindahkan kolom ke kanan. Komponen lainnya $$$x$ karena mereka nol, mereka tetap demikian. Secara skematis, ini ditunjukkan di bawah ini, di mana panah bawah menunjukkan urutan penghancuran kolom:

Catatan: dalam urutan penghancuran kolom ini, angka 3 akan ditemui paling lambat dari angka 5, 9, dan 10! Kolom tidak dihancurkan dalam urutan menaik, yang akan menjadi kesalahan bagi yang padat, mis. matriks yang tidak lengkap. Tetapi untuk matriks yang jarang, ini adalah hal biasa. Seperti dapat dilihat dari struktur matriks non-nol dalam contoh ini, menggunakan kolom 5, 9, dan 10 ke kolom 3 tidak akan mengubah komponen dalam respons $$$x_5$ , $$$x_9$ , $$$x_ {10}$ tidak sama sekali, mereka memiliki c $$$x_3$ hanya "tidak ada persimpangan." Metode kami memperhitungkan ini! Demikian pula, menggunakan kolom 8 setelah kolom 9 tidak akan merusak komponen $$$x_9$ . Algoritma hebat, bukan?

Jika kita berkeliling grafik pertama dari indeks 5, dan kemudian dari indeks 3, maka array kita akan menjadi $$$nnz (x): = \ {\ color {blue} {11}, \ color {blue} 8, \ color {blue} {10}, \ color {blue} {9}, \ color {blue} 5, \ color {green} 3 \}$ . Menghancurkan kolom dalam urutan terbalik dari array ini akan memberikan solusi yang persis sama seperti pada kasus pertama!

Kompleksitas semua operasi kami diskalakan dengan jumlah operasi aritmatika aktual dan jumlah komponen bukan nol di sisi kanan, tetapi tidak dengan ukuran matriks! Kami telah mencapai tujuan kami.

Kritik

NAMUN! Pembaca kritis akan melihat bahwa asumsi di awal bahwa kita memiliki "akses langsung ke komponen bukan nol dari sisi kanan dengan indeks" sudah berarti bahwa sekali sebelum kita menjalankan sisi kanan dari atas ke bawah untuk menemukan indeks ini dan mengaturnya dalam array $$$nnz (f)$ , yaitu, sudah menghabiskan $$$O (N)$ aksi! Selain itu, grafik berjalan dengan sendirinya mengharuskan kita terlebih dahulu mengalokasikan memori dari panjang maksimum yang mungkin (di tempat lain Anda perlu menuliskan indeks yang ditemukan dengan mencari secara mendalam!) Agar tidak menderita realokasi lebih lanjut karena array bertambah $$$nnz (x)$ , dan ini juga membutuhkan $$$O (N)$ operasi! Lalu, mengapa, kata mereka, semua ribut?

Tetapi memang, untuk solusi satu kali dari sistem segitiga bawah yang jarang dengan sisi kanan jarang, awalnya didefinisikan sebagai vektor padat, tidak masuk akal untuk membuang waktu pengembang pada semua algoritma yang disebutkan di atas. Mereka mungkin bahkan lebih lambat daripada metode dahi yang disajikan oleh algoritma 2 di atas. Tapi, seperti yang telah disebutkan, alat ini sangat diperlukan untuk faktorisasi Cholesky, jadi Anda tidak harus melemparkan tomat ke arahku. Memang, sebelum menjalankan algoritma 3, semua memori yang diperlukan dengan panjang maksimum dialokasikan segera dan membutuhkan $$$O (N)$ waktu. Dalam siklus kolom selanjutnya $$$A$ semua data hanya ditimpa dalam array dengan panjang tetap $$$N$ , dan hanya pada posisi-posisi di mana penulisan ulang ini relevan, berkat akses langsung ke elemen yang tidak nol. Dan justru karena inilah efisiensi muncul!

Implementasi C ++

Grafik itu sendiri sebagai struktur data dalam kode tidak perlu dibangun. Sudah cukup untuk menggunakannya secara implisit ketika bekerja dengan matriks secara langsung. Semua konektivitas yang diperlukan akan diperhitungkan oleh algoritma. Pencarian mendalam tentu saja mudah diimplementasikan menggunakan rekursi dangkal. Di sini, misalnya, tampak seperti pencarian kedalaman rekursif berdasarkan indeks mulai tunggal:

 /* j -   iA, jA -    ,    CSC top -     nnz(x);     0 result -   N,    nnz(x)   0  top-1  marked -    N;      */ void DepthFirstSearch(size_t j, const std::vector<size_t>& iA, const std::vector<size_t>& jA, size_t& top, std::vector<size_t>& result, std::vector<int>& marked) { size_t p; marked[j] = 1; //   j   for (p = jA[j]; p < jA[j+1]; p++) //       j { if (!marked[iA[p]]) //  iA[p]   { DepthFirstSearch(iA[p], iA, jA, top, result, marked); //      iA[p] } } result[top++] = j; //  j   nnz(x) } 

Jika Anda melewatkan variabel top sama dengan nol pada panggilan pertama ke DepthFirstSearch , maka setelah penyelesaian seluruh prosedur rekursif, variabel top akan sama dengan jumlah indeks yang ditemukan dalam array result . Pekerjaan rumah: tulis fungsi lain yang mengambil array indeks komponen bukan nol di sisi kanan dan meneruskannya secara berurutan ke fungsi DepthFirstSearch . Tanpa ini, algoritme tidak lengkap. Harap dicatat: array bilangan real $$$vA$ kami tidak meneruskannya ke fungsi sama sekali, karena itu tidak diperlukan dalam proses analisis simbolik.

Terlepas dari kesederhanaannya, cacat dalam implementasinya jelas: untuk sistem yang besar, stack overflow baru saja tiba. Kalau begitu ada pilihan berdasarkan loop bukan rekursi. Ini lebih sulit untuk dipahami, tetapi sudah cocok untuk semua kesempatan:

 /* j -   N -   iA, jA -    ,    CSC top -     nnz(x) result -   N,     nnz(x)   top  ret-1 ,  ret -    marked -    N work_stack -     N */ size_t DepthFirstSearch(size_t j, size_t N, const std::vector<size_t>& iA, const std::vector<size_t>& jA, size_t top, std::vector<size_t>& result, std::vector<int>& marked, std::vector<size_t>& work_stack) { size_t p, head = N - 1; int done; result[N - 1] = j; //    while (head < N) { j = result[head]; //  j     if (!marked[j]) { marked[j] = 1; //   j   work_stack[head] = jA[j]; } done = 1; //    j      for (p = work_stack[head] + 1; p < jA[j+1]; p++) //    j { //   c   iA[p] if (marked[iA[p]]) continue; //  iA[p]  ,   work_stack[head] = p; //       j result[--head] = iA[p]; //       iA[p] done = 0; //   j    break; } if (done) //      j  { head++; //  j    result[top++] = j; //  j    } } return (top); } 

Dan inilah yang terlihat dari generator struktur vektor bukan nol $$$nnz (x)$ :

 /* iA, jA -   ,    CSC iF -        f nnzf -      f result -   N,    nnz(x)   0  ret-1 ,  ret -    marked -    N,    work_stack -     N */ size_t NnzX(const std::vector<size_t>& iA, const std::vector<size_t>& jA, const std::vector<size_t>& iF, size_t nnzf, std::vector<size_t>& result, std::vector<int>& marked, std::vector<size_t>& work_stack) { size_t p, N, top; N = jA.size() - 1; top = 0; for (p = 0; p < nnzf; p++) if (!marked[iF[p]]) //        top = DepthFirstSearch(iF[p], N, iA, jA, vA, top, result, marked, work_stack); for (p = 0; p < top; p++) marked[result[p]] = 0; //   return (top); } 

Akhirnya, menggabungkan semuanya menjadi satu, kami menulis solver segitiga bawah untuk kasus sisi kanan yang dijernihkan:

 /* iA, jA, vA -  ,    CSC iF -        f nnzf -      f vF -       f result -   N,    nnz(x)   0  ret-1 ,  ret -    marked -    N,    work_stack -     N x -    N,    ;     */ size_t lower_triangular_solve (const std::vector<size_t>& iA, const std::vector<size_t>& jA, const std::vector<double>& vA, const std::vector<size_t>& iF, size_t nnzf, const std::vector<double>& vF, std::vector<size_t>& result, std::vector<int>& marked, std::vector<size_t>& work_stack, std::vector<double>& x) { size_t top, p, j; ptrdiff_t px; top = NnzX(iA, jA, iF, nnzf, result, marked, work_stack); for (p = 0; p < nnzf; p++) x[iF[p]] = vF[p]; //    for (px = top; px > -1; px--) //     { j = result[px]; // x(j)   x [j] /= vA[jA[j]]; //   x(j) for (p = jA[j]+1; p < jA[j+1]; p++) { x[iA[p]] -= vA[p]*x[j]; //  j-  } } return (top) ; } 

Kami melihat bahwa siklus kami hanya berjalan melalui indeks array $$$nnz (x)$ , tidak seluruh set $$$0, 1, ..., N - 1$ . Selesai!

Ada implementasi yang tidak menggunakan array tag markeduntuk menghemat memori. Sebagai gantinya, satu set indeks tambahan digunakan.$$$V_1$ tidak berpotongan dengan set $$$V = \{0, ..., N - 1\}$di mana ada pemetaan satu-ke-satu dengan operasi aljabar sederhana sebagai prosedur menandai node. Namun, di era memori murah kita, simpan dalam satu larik panjang$$$N$ tampaknya sepenuhnya berlebihan.

Kesimpulannya

Proses memecahkan sistem persamaan linear yang jarang dengan metode langsung, sebagai suatu peraturan, dibagi menjadi tiga tahap:

1. Analisis karakter
2. Faktorisasi numerik berdasarkan data analisis sivol
3. Solusi dari sistem segitiga yang diperoleh dengan sisi kanan yang padat

Langkah kedua, faktorisasi numerik, adalah bagian yang paling memakan sumber daya dan memakan paling banyak (> 90%) dari perkiraan waktu. Tujuan dari langkah pertama adalah untuk mengurangi biaya tinggi yang kedua. Contoh analisis simbolik disajikan dalam posting ini. Namun, ini adalah langkah pertama yang membutuhkan waktu pengembangan paling lama dan biaya mental maksimum dari pihak pengembang. Analisis simbolis yang baik membutuhkan pengetahuan tentang teori grafik dan pohon serta memiliki "naluri algoritmik". Langkah kedua jauh lebih mudah diimplementasikan.
Nah, langkah ketiga, baik dalam implementasi dan dalam perkiraan waktu, dalam banyak kasus adalah yang paling bersahaja.

Pengantar yang baik untuk metode langsung untuk sistem jarang ada di Tim Davis , seorang profesor di Departemen Ilmu Komputer di Texas A&M University , berjudul "Metode Langsung untuk Sistem Linier Jarang ". Algoritme yang dijelaskan di atas dijelaskan di sana. Saya sendiri tidak mengenal Davis, walaupun saya sendiri dulu bekerja di universitas yang sama (walaupun di fakultas yang berbeda). Jika saya tidak salah, Davis sendiri berpartisipasi dalam pengembangan solver (!). digunakan dalam Matlab selain itu, ia terlibat bahkan untuk gambar pembangkit listrik di Google Street View (OLS) dengan cara, di fakultas yang sama yang tercantum tidak lain dirinya Bjarne Stroustrup - .. pencipta C ++

Mungkin suatu hari nanti Saya akan menulis sesuatu yang lain tentang matriks yang jarang atau metode numerik di umum. Baik untuk semua!