J
Saat ini saya melakukan ini dengan perpustakaan armadillo tetapi ternyata terlalu lambat. Masalahnya adalah bahwa saya perlu melakukan ini untuk satu triliun matriks dan ternyata menghitung dua penentu adalah hambatan dari program saya. Karena itu saya punya dua pertanyaan
Apakah ada trik yang bisa saya gunakan untuk menghitung penentu lebih cepat mengingat saya tahu ukurannya? Apakah mungkin ekspansi berantakan untuk matriks yang dapat bekerja dalam kasus ini?
Apakah ada cara lain yang efisien untuk menguji kesetaraan
Edit. Untuk menjawab komentar. Saya perlu menghitung semua graf saling melengkapi yang tidak terhubung sendiri dengan urutan sehingga dan memiliki jumlah pohon rentang yang sama. Motivasi untuk ini dapat ditemukan di pos mathoverflow ini . Sedangkan untuk mesin saya menjalankan ini pada mesin 8 core 3.4GHh secara paralel.13 G ¯ G
Edit. Saya dapat mengurangi waktu berjalan yang diharapkan sebesar 50% dengan membuat program C untuk secara spesifik menghitung determinan matriks . Saran masih diterima.
sumber
Jawaban:
Karena Anda sudah menggunakan C ++ dan matriks Anda pasti positif simetris, saya akan melakukan faktorisasi tidak dipvotivasikan dari dan juga . Di sini saya berasumsi bahwa juga pasti positif, jika tidak, akan memerlukan pivoting untuk stabilitas numerik (juga mungkin meskipun itu tidak pasti positif, pivoting tidak diperlukan, tetapi Anda harus mencobanya). Q 12 I - Q - J 12 I - Q - J L D L TLDLT Q 12I−Q−J 12I−Q−J LDLT
Ini lebih cepat daripada faktorisasi LU, dan juga lebih cepat daripada Cholesky karena akar kuadrat dihindari. Penentu hanyalah produk dari elemen-elemen matriks diagonal . Kode untuk melakukan faktorisasi LDL sangat sederhana sehingga Anda dapat menulisnya dalam waktu kurang dari 50 baris C. Halaman Wikipedia di dalamnya menjelaskan algoritme, dan saya memiliki beberapa kode templated sederhana untuk melakukan Cholesky di sini . Anda dapat menyederhanakannya dan memodifikasinya untuk menghindari akar kuadrat untuk mengimplementasikan faktorisasi .L D L TD LDLT
Karena Anda juga dapat mengontrol format penyimpanan, Anda dapat lebih mengoptimalkan rutin untuk menyimpan hanya setengah matriks dan mengemasnya dalam array linier untuk memaksimalkan memori lokalitas. Saya juga akan menulis produk kustom titik sederhana dan rutinitas pembaruan peringkat-1 karena ukuran masalahnya sangat kecil, Anda harus membiarkan kompiler menyatukan rutinitas untuk mengurangi overhead panggilan. Karena itu adalah loop ukuran tetap, kompiler harus dapat secara otomatis inline dan membuka gulungan barang bila perlu.
Saya akan menghindari mencoba memainkan trik untuk mengambil keuntungan dari kenyataan bahwa berisi di dalam ekspresi. Sangat mungkin bahwa untuk ukuran masalah sekecil itu, trik-trik ini pada akhirnya menjadi lebih lambat daripada hanya melakukan dua perhitungan penentu terpisah. Tentu saja, satu-satunya cara untuk memverifikasi klaim ini adalah dengan mencobanya.Q12I−Q−J Q
sumber
Tanpa beberapa informasi tentang konstruksi matriks simetrik nyata pasti positif ini, saran-saran yang akan dibuat adalah kebutuhan yang cukup terbatas.12×12
Saya mengunduh paket Armadillo dari Sourceforge dan melihat dokumentasi. Cobalah untuk meningkatkan kinerja komputasi secara terpisah dan , di mana adalah matriks peringkat satu dari semua yang ada, dengan menetapkan misalnya . Dokumentasi mencatat bahwa ini adalah default untuk matriks hingga ukuran , jadi karena kelalaian saya menganggap opsi adalah default untuk kasus .det(Q) det(12I−Q−J) J 4×4 12×12
det(Q,slow=false)
slow=true
Apa yangQ 12×12 penentu akan dilakukan karena perhitungan akan "dilempar ke atas tembok" ke LAPACK (atau ATLAS) pada saat itu tanpa indikasi bahwa pivot tidak diperlukan; lihat
slow=true
mungkin dilakukan adalah berputar sebagian atau penuh dalam mendapatkan bentuk eselon baris, dari mana determinan mudah ditemukan. Namun Anda tahu sebelumnya bahwa matriks adalah positif pasti, jadi pivot tidak diperlukan untuk stabilitas (setidaknya secara dugaan untuk sebagian besar perhitungan Anda. Tidak jelas apakah paket Armadillo melempar pengecualian jika pivot menjadi terlalu kecil, tetapi ini harus menjadi fitur paket aljabar linear numerik yang masuk akal. EDIT: Saya menemukan kode Armadillo yang mengimplementasikan dalam file header , menggunakan template C ++ untuk fungsionalitas substansial. Pengaturan ini tampaknya tidak mempengaruhi bagaimanadet
include\armadillo_bits\auxlib_meat.hpp
slow=false
det_lapack
dan doanya dalam file itu.Poin lainnya adalah mengikuti rekomendasi mereka untuk membangun paket Armadillo yang menghubungkan ke penggantian kecepatan tinggi untuk BLAS dan LAPACK, jika Anda memang menggunakannya; lihat Sec. 5 file Armadillo README.TXT untuk detailnya. [Penggunaan BLAS atau LAPACK versi 64-bit khusus juga disarankan untuk kecepatan pada mesin 64-bit saat ini.]
Reduksi baris ke bentuk eselon pada dasarnya adalah eliminasi Gaussian, dan memiliki kompleksitas aritmatika . Untuk kedua matriks ini, berarti dua kali lipat yang berfungsi, atau . Operasi-operasi ini mungkin menjadi "penghambat" dalam pemrosesan Anda, tetapi ada sedikit harapan bahwa tanpa struktur khusus di (atau beberapa hubungan yang diketahui di antara triliun kasus uji yang memungkinkan amortisasi) pekerjaan dapat dikurangi menjadi .23n3+O(n2) 43n3+O(n2) Q O(n2)
Sebagai perbandingan, ekspansi oleh kofaktor umum matriks melibatkanoperasi perkalian (dan kira-kira sebanyak penambahan / pengurangan), jadi untuk perbandingan ( vs ) jelas mendukung penghapusan lebih dari kofaktor.n×n n! n=12 12!=479001600 23n3=1152
Pendekatan lain yang membutuhkan bekerja akan mengurangi ke bentuk tridiagonal dengan transformasi Householder, yang juga menempatkan ke dalam bentuk tridiagonal. Komputasi dan selanjutnya dapat dilakukan dalam operasi . [Efek dari pembaruan peringkat satu dalam penentu kedua dapat dinyatakan sebagai faktor skalar yang diberikan dengan menyelesaikan satu sistem tridiagonal.]43n3+O(n2) Q 12I−Q det(Q) det(12I−Q−J) O(n) −J
Menerapkan perhitungan independen semacam itu mungkin bermanfaat sebagai pemeriksaan atas hasil panggilan yang berhasil (atau gagal) ke
det
fungsi Armadillo .Kasus Khusus: Seperti yang disarankan oleh Komentar Jernej, anggaplah bahwa mana seperti sebelumnya adalah matriks (peringkat 1) dari semua yang ada dan adalah matriks diagonal nonsingular (positif). Memang untuk aplikasi yang diusulkan dalam teori grafik ini akan menjadi matriks bilangan bulat. Maka rumus eksplisit untuk adalah:Q=D−J J D=diag(d1,…,dn) det(Q)
Sebuah sketsa pembuktiannya memberikan peluang untuk menggambarkan penerapan yang lebih luas, yaitu setiap kali memiliki determinan yang diketahui dan sistem dengan cepat diselesaikan. Mulailah dengan memperhitungkan:D Dv=(1…1)T
Sekarang kembali menjadi peringkat 1, yaitu . Perhatikan bahwa penentu kedua sederhana:D−1J (d−11…d−1n)T(1…1)
di mana adalah polinomial karakteristik dari . Sebagai matriks peringkat 1, harus memiliki (setidaknya) faktor untuk memperhitungkan ruang nolnya. Nilai eigen "hilang" adalah , seperti yang dapat dilihat dari perhitungan:f(x) D−1J f(x) n−1 x ∑d−1i
Karena itu karakteristik polinomial , dan adalah seperti yang ditunjukkan di atas untuk , .f(x)=xn−1(x−∑d−1i) f(1) det(I−D−1J) 1−∑d−1i
Perhatikan juga bahwa jika , maka , sebuah matriks diagonal yang determinannya hanyalah produk dari entri diagonalnya.Q=D−J 12I−Q−J=12I−D+J−J=12I−D
sumber
Jika Anda memiliki cara terstruktur untuk menghitung grafik yang ingin Anda hitung determinannya, mungkin Anda bisa menemukan pembaruan peringkat rendah yang memindahkan Anda dari satu grafik ke grafik lainnya.
Jika demikian, maka Anda dapat menggunakan lemma penentu matriks untuk menghitung dengan murah determinan grafik berikutnya yang akan dihitung menggunakan pengetahuan Anda tentang determinan grafik saat ini.
Yaitu, untuk matriks dan vektor : Ini dapat digeneralisasi jika U dan V adalah matriks dan adalah :A u,v det(A+uvT)=(1+vTA−1u)det(A) n×m A n×n det(A+UVT)=det(Im+VTA−1U)det(A)
Untuk menghitung invers secara efisien, Anda dapat menggunakan rumus Sherman-Morrison untuk mendapatkan invers dari matriks berikutnya dari matriks saat ini:(A+uvT)−1=A−1−A−1uvTA−11+vTA−1u
sumber