Bagaimana menangani kompleksitas dalam kode numerik, misalnya, ketika berhadapan dengan matriks Jacobian yang besar?

Saya memecahkan sistem persamaan linear berganda, dan telah menghitung Jacobian dari sistem diskritisasi. Hasilnya sangat rumit, di bawah ini (hanya!) 3 kolom pertama dari matriks $3\times 9$ ,

(Kerumitan muncul, sebagian, karena skema numerik memerlukan pemasangan eksponensial untuk stabilitas.)

Saya punya pertanyaan umum tentang implementasi kode numerik menggunakan Jacobian.

Saya dapat melanjutkan dan mengimplementasikan matriks ini dalam kode. Tetapi intuisi saya mengatakan kepada saya untuk mengharapkan debugging yang membosankan beberapa hari (mungkin beberapa minggu!) Karena kerumitannya dan tidak terhindarkannya kesalahan yang muncul. Bagaimana seseorang mengatasi kompleksitas seperti ini dalam kode numerik, tampaknya tidak dapat dihindari ?! Apakah Anda menggunakan pembuatan kode otomatis dari paket simbolis (kemudian men-tweak kode dengan tangan)?

Pertama saya berencana untuk men-debug Jacobian analitik dengan perkiraan perbedaan yang terbatas, haruskah saya mengetahui adanya jebakan? Bagaimana Anda menangani masalah serupa dalam kode Anda?

Memperbarui

Saya mengkodekan ini dengan Python dan menggunakan sympy untuk menghasilkan Jacobian. Mungkin saya bisa menggunakan fitur pembuatan kode ?

Satu kata: Modularitas .

Ada banyak ekspresi berulang dalam Jacobian Anda yang dapat ditulis sebagai fungsi mereka sendiri. Tidak ada alasan bagi Anda untuk menulis operasi yang sama lebih dari sekali, dan itu akan membuat proses debug lebih mudah; jika Anda hanya menulis sekali hanya ada satu tempat untuk kesalahan (secara teori).

Kode modular juga akan membuat pengujian lebih mudah; Anda dapat menulis tes untuk setiap komponen Jacobian Anda sebagai lawan mencoba untuk menguji seluruh matriks. Misalnya, jika Anda menulis fungsi saya () dengan cara modular, Anda dapat dengan mudah menulis tes kewarasan untuk itu, periksa apakah Anda membedakannya dengan benar, dll.

Saran lain adalah untuk melihat perpustakaan diferensiasi otomatis untuk merakit Jacobian. Tidak ada jaminan mereka bebas dari kesalahan, tetapi mungkin akan ada lebih sedikit kesalahan debugging / lebih sedikit daripada menulis sendiri. Berikut ini beberapa yang mungkin ingin Anda lihat:

• Sacado (Sandia Labs)
• ADIC (Argonne)

Maaf, baru saja melihat bahwa Anda menggunakan python. ScientificPython memiliki dukungan untuk AD.

• ScientificPython

Biarkan saya menimbang di sini dengan beberapa kata hati-hati, diawali dengan sebuah cerita. Dulu, saya bekerja dengan seorang rekan ketika saya baru memulai. Dia memiliki masalah optimisasi untuk dipecahkan, dengan tujuan yang agak berantakan. Solusinya adalah menghasilkan turunan analitik untuk optimasi.

Masalah yang saya lihat adalah turunan ini buruk. Dihasilkan menggunakan Macsyma, dikonversi ke kode fortran, mereka masing-masing lusinan pernyataan lanjutan. Bahkan, kompiler Fortran menjadi marah pada saat itu, karena melebihi jumlah maksimum pernyataan lanjutan. Sementara kami menemukan bendera yang memungkinkan kami untuk mengatasi masalah itu, ada masalah lain.

• Dalam ekspresi panjang, seperti yang umumnya dihasilkan oleh sistem CA, ada risiko pembatalan subtraktif masif. Hitung banyak angka besar, hanya untuk menemukan mereka semua membatalkan satu sama lain untuk menghasilkan angka kecil.

• Seringkali derivatif yang dihasilkan secara analitik sebenarnya lebih mahal untuk dievaluasi daripada derivatif yang dihasilkan secara numerik menggunakan perbedaan hingga. Gradien untuk n variabel dapat memakan waktu lebih dari n kali biaya mengevaluasi fungsi objektif Anda. (Anda mungkin dapat menghemat waktu karena banyak istilah dapat digunakan kembali di berbagai turunan, tetapi itu juga akan memaksa Anda untuk melakukan pengkodean tangan secara hati-hati, alih-alih menggunakan ekspresi yang dihasilkan komputer. Dan setiap kali Anda membuat kode matematika tidak menyenangkan ekspresi, probabilitas kesalahan tidak sepele. Pastikan Anda memverifikasi turunan ini untuk akurasi.)

Inti cerita saya adalah ekspresi yang dihasilkan CA ini memiliki masalah sendiri. Yang lucu adalah kolega saya sebenarnya bangga dengan kerumitan masalahnya, bahwa ia jelas memecahkan masalah yang sangat sulit karena aljabarnya begitu jahat. Apa yang saya pikir tidak dia pikirkan adalah apakah aljabar itu benar-benar menghitung hal yang benar, apakah itu melakukannya dengan akurat, dan apakah itu melakukannya dengan efisien.

Seandainya saya menjadi orang senior saat itu dalam proyek ini, saya akan membacakan kepadanya aksi kerusuhan. Kebanggaannya menyebabkan dia menggunakan solusi yang mungkin tidak perlu rumit, bahkan tanpa memeriksa bahwa gradien berbasis perbedaan yang terbatas sudah memadai. Saya berani bertaruh kami telah menghabiskan mungkin satu minggu minggu waktu untuk menjalankan optimasi ini. Paling tidak, saya akan menasihatinya untuk hati-hati menguji gradien yang dihasilkan. Apakah itu akurat? Seberapa akurat itu, dibandingkan dengan turunan beda hingga? Bahkan, ada alat di sekitar hari ini yang juga akan mengembalikan perkiraan kesalahan dalam prediksi turunannya. Ini tentu benar untuk kode diferensiasi adaptif, (turunan) yang saya tulis dalam MATLAB.

Uji kodenya. Verifikasi turunannya.

Tetapi sebelum Anda melakukan APAPUN ini, pertimbangkan apakah skema optimasi lain yang lebih baik adalah sebuah pilihan. Misalnya, jika Anda melakukan pemasangan eksponensial, maka ada peluang yang sangat baik Anda dapat menggunakan kuadrat terkecil non-linier yang dipartisi (kadang-kadang disebut kuadrat terkecil yang dapat dipisah. Saya pikir itulah istilah yang digunakan oleh Seber dan Wild dalam buku mereka.) Idenya adalah untuk memecah set parameter menjadi set intrinsik linier dan intrinsik nonlinier. Gunakan pengoptimalan yang hanya berfungsi pada parameter nonlinear. Mengingat parameter tersebut "diketahui", maka parameter linear intrinsik dapat diperkirakan dengan menggunakan kuadrat terkecil linier sederhana. Skema ini akan mengurangi ruang parameter dalam optimasi. Itu membuat masalah lebih kuat, karena Anda tidak perlu menemukan nilai awal untuk parameter linier. Ini mengurangi dimensi ruang pencarian Anda, sehingga membuat masalah berjalan lebih cepat. Sekali lagi saya berikanalat untuk tujuan ini , tetapi hanya di MATLAB.

Jika Anda menggunakan turunan analitik, beri kode untuk menggunakan kembali istilah. Ini bisa menjadi penghematan waktu yang serius, dan sebenarnya dapat mengurangi bug, menghemat waktu Anda sendiri. Tapi periksa angka-angka itu!

Ada beberapa strategi untuk dipertimbangkan:

1. Temukan turunan dalam bentuk simbolis menggunakan CAS, lalu ekspor kode untuk menghitung turunannya.

2. Gunakan alat diferensiasi otomatis (AD) untuk menghasilkan kode yang menghitung turunan dari kode untuk menghitung fungsi.

3. Gunakan perkiraan perbedaan hingga untuk mendekati Jacobian.

Diferensiasi otomatis dapat menghasilkan kode yang lebih efisien untuk menghitung seluruh Jacobian kemudian menggunakan perhitungan simbolik untuk menghasilkan formula untuk setiap entri dalam matriks. Perbedaan terbatas adalah cara yang baik untuk memeriksa turunan Anda.

Berikut adalah contoh di mana kami menggunakan diferensiasi otomatis menggunakan Sacado dalam satu kode: http://www.dealii.org/developer/doxygen/deal.II/step_33.html

Selain saran bagus BrianBorcher, pendekatan lain yang mungkin dilakukan untuk fungsi bernilai riil adalah dengan menggunakan pendekatan turunan langkah-kompleks (lihat artikel ini (paywalled) dan artikel ini ). Dalam beberapa kasus, pendekatan ini menghasilkan turunan numerik yang lebih akurat dengan biaya mengubah nilai variabel dalam fungsi Anda dari nyata menjadi kompleks. Artikel kedua berisi daftar beberapa kasus di mana pendekatan fungsi langkah kompleks mungkin rusak.