Bagaimana operator backslash MATLAB memecahkan

Saya membandingkan beberapa kode saya dengan "stok" kode MATLAB. Saya terkejut dengan hasilnya.

Saya menjalankan kode sampel (Matriks Jarang)

n = 5000;
a = diag(rand(n,1));
b = rand(n,1);
disp('For a\b');
tic;a\b;toc;
disp('For LU');
tic;LULU;toc;
disp('For Conj Grad');
tic;conjgrad(a,b,1e-8);toc;
disp('Inv(A)*B');
tic;inv(a)*b;toc;

Hasil:

    For a\b
    Elapsed time is 0.052838 seconds.

    For LU
    Elapsed time is 7.441331 seconds.

    For Conj Grad
    Elapsed time is 3.819182 seconds.

    Inv(A)*B
    Elapsed time is 38.511110 seconds.

Untuk Matriks Padat:

n = 2000;
a = rand(n,n);
b = rand(n,1);
disp('For a\b');
tic;a\b;toc;
disp('For LU');
tic;LULU;toc;
disp('For Conj Grad');
tic;conjgrad(a,b,1e-8);toc;
disp('For INV(A)*B');
tic;inv(a)*b;toc;

Hasil:

For a\b
Elapsed time is 0.575926 seconds.

For LU
Elapsed time is 0.654287 seconds.

For Conj Grad
Elapsed time is 9.875896 seconds.

Inv(A)*B
Elapsed time is 1.648074 seconds.

Bagaimana sih a \ b begitu mengagumkan?

Dalam Matlab, perintah '\' memanggil algoritma yang tergantung pada struktur matriks A dan mencakup pemeriksaan (overhead kecil) pada properti A.

1. Jika A jarang dan berpita, gunakan pemecah berpita.
2. Jika A adalah matriks segitiga atas atau bawah, gunakan algoritma substitusi mundur.
3. Jika A simetris dan memiliki elemen diagonal positif nyata, cobalah faktorisasi Cholesky. Jika A jarang, gunakan pemesanan ulang terlebih dahulu untuk meminimalkan pengisian.
4. Jika tidak ada kriteria di atas yang dipenuhi, lakukan faktorisasi segitiga umum menggunakan eliminasi Gaussian dengan pivoting parsial.
5. Jika A jarang, maka gunakan perpustakaan UMFPACK.
6. Jika A tidak kuadrat, gunakan algoritma berdasarkan faktorisasi QR untuk sistem yang tidak ditentukan.

Untuk mengurangi overhead, dimungkinkan untuk menggunakan perintah linsolve di Matlab dan pilih sendiri solver di antara opsi-opsi ini.

Jika Anda ingin melihat apa yang a\bdilakukan untuk matriks khusus Anda, Anda dapat mengatur spparms('spumoni',1)dan mencari algoritma apa yang membuat Anda terkesan. Sebagai contoh:

spparms('spumoni',1);
A = delsq(numgrid('B',256));
b = rand(size(A,2),1);
mldivide(A,b);  % another way to write A\b

akan menampilkan

sp\: bandwidth = 254+1+254.
sp\: is A diagonal? no.
sp\: is band density (0.01) > bandden (0.50) to try banded solver? no.
sp\: is A triangular? no.
sp\: is A morally triangular? no.
sp\: is A a candidate for Cholesky (symmetric, real positive diagonal)? yes.
sp\: is CHOLMOD's symbolic Cholesky factorization (with automatic reordering) successful? yes.
sp\: is CHOLMOD's numeric Cholesky factorization successful? yes.
sp\: is CHOLMOD's triangular solve successful? yes.

jadi saya bisa melihat bahwa "\" akhirnya menggunakan "CHOLMOD" dalam kasus ini.

Untuk matriks jarang, Matlab menggunakan UMFPACK untuk \operasi " ", yang, dalam contoh Anda, pada dasarnya dijalankan melalui nilai - nilai a, membalikkannya, dan mengalikannya dengan nilai - nilai dari b. Namun, untuk contoh ini, Anda harus menggunakan b./diag(a), yang jauh lebih cepat.

Untuk sistem yang padat, operator backslash sedikit lebih rumit. Deskripsi singkat tentang apa yang dilakukan saat diberikan di sini . Menurut deskripsi itu, dalam contoh Anda, Matlab akan menyelesaikan a\bmenggunakan substitusi mundur. Untuk matriks kuadrat umum, dekomposisi LU digunakan.

Sebagai aturan praktis, jika Anda memiliki matriks tipis kompleksitas yang masuk akal (yaitu, itu tidak harus menjadi stensil 5 poin tetapi sebenarnya bisa menjadi diskritisasi dari persamaan Stokes yang jumlah nonzeros per baris adalah jauh lebih besar dari 5), maka pemecah langsung yang jarang seperti UMFPACK biasanya mengalahkan pemecah Krylov yang berulang jika masalahnya tidak lebih besar dari sekitar mungkin 100.000 yang tidak diketahui.

Dengan kata lain, untuk sebagian besar matriks jarang yang dihasilkan dari diskresi 2d, pemecah langsung adalah alternatif tercepat. Hanya untuk masalah 3d di mana Anda dengan cepat mendapatkan di atas 100.000 yang tidak diketahui apakah penting untuk menggunakan pemecah berulang.