Algoritma Penyeimbang Matriks

Saya telah menulis toolbox sistem kontrol dari awal dan murni di Python3 (plug shameless:) harold. Dari penelitian saya sebelumnya, saya selalu mengeluh tentang pemecah Riccati care.mkarena alasan teknis / tidak relevan.

Karenanya, saya telah menulis serangkaian rutinitas saya sendiri. Satu hal yang saya tidak dapat menemukan jalan keluar adalah untuk mendapatkan algoritma penyeimbangan kinerja tinggi, setidaknya sebaik balance.m. Sebelum Anda menyebutkannya, xGEBALkeluarga terpapar dalam Scipy dan pada dasarnya Anda dapat memanggil dari Scipy sebagai berikut, misalkan Anda memiliki array 2D tipe float A:

import scipy as sp
gebal = sp.linalg.get_lapack_funcs(('gebal'),(A,)) # this picks up DGEBAL
Ab, lo, hi, scaling , info = gebal(A, scale=1 , permute=1 , overwrite_a=0 )

Sekarang jika saya menggunakan matriks tes berikut

array([[ 6.      ,  0.      ,  0.      ,  0.      ,  0.000002],
       [ 0.      ,  8.      ,  0.      ,  0.      ,  0.      ],
       [ 2.      ,  2.      ,  6.      ,  0.      ,  0.      ],
       [ 2.      ,  2.      ,  0.      ,  8.      ,  0.      ],
       [ 0.      ,  0.      ,  0.000002,  0.      ,  2.      ]])

saya mendapat

array([[ 8.      ,  0.      ,  0.      ,  2.      ,  2.      ],
       [ 0.      ,  2.      ,  0.000002,  0.      ,  0.      ],
       [ 0.      ,  0.      ,  6.      ,  2.      ,  2.      ],
       [ 0.      ,  0.000002,  0.      ,  6.      ,  0.      ],
       [ 0.      ,  0.      ,  0.      ,  0.      ,  8.      ]])

Namun, jika saya meneruskan ini balance.m, saya mengerti

>> balance(A)

ans =

    8.0000         0         0    0.0625    2.0000
         0    2.0000    0.0001         0         0
         0         0    6.0000    0.0002    0.0078
         0    0.0003         0    6.0000         0
         0         0         0         0    8.0000

Jika Anda memeriksa pola permutasi, mereka sama namun penskalaannya tidak aktif. gebalmemberikan kerak kesatuan sedangkan Matlab memberikan kekuatan berikut 2: [-5,0,8,0,2].

Jadi ternyata, ini tidak menggunakan mesin yang sama. Saya sudah mencoba berbagai opsi seperti Lemonnier, Van Dooren skala dua sisi, Parlett-Reinsch asli dan juga beberapa metode lain yang kurang dikenal dalam literatur seperti versi padat SPBALANCE.

Satu hal yang mungkin saya tekankan adalah bahwa saya sadar akan pekerjaan Benner; khususnya Balancing Matriks Hamiltonian Symplectic khusus untuk tujuan ini. Namun, perhatikan bahwa jenis perawatan ini dilakukan di dalam gcare.m(solver Riccati umum) dan penyeimbangan dilakukan secara langsung melalui balance.m. Karena itu, saya akan sangat menghargai jika seseorang dapat mengarahkan saya ke implementasi yang sebenarnya.

Pengungkapan: Saya benar-benar tidak berusaha untuk membalikkan kode matematika insinyur: Saya benar-benar ingin menjauh darinya karena berbagai alasan termasuk motivasi dari pertanyaan ini, yaitu, saya tidak tahu apa yang dilakukannya yang membuat saya harus mengeluarkan banyak biaya waktu kembali pada hari itu. Tujuan saya adalah untuk mendapatkan algoritma balancing yang memuaskan yang memungkinkan saya untuk lulus contoh CAREX sehingga saya dapat menerapkan metode iterasi Newton di atas pemecah reguler.

Butuh waktu cukup lama untuk memikirkan ini dan seperti biasa menjadi jelas setelah Anda menemukan pelakunya.

Setelah memeriksa kasus-kasus bermasalah yang dilaporkan dalam David S. Watkins. Kasus di mana keseimbangan berbahaya. Elektron. Trans. Angka Anal, 23: 1-4, 2006 dan juga diskusi di sini (keduanya dikutip dalam arXiv: 1401.5766v1 ), ternyata matlab menggunakan penyeimbang dengan memisahkan elemen diagonal terlebih dahulu.

Pikiran awal saya adalah, sesuai dengan dokumentasi klasik terbatas pada fungsi LAPACK, GEBAL melakukan ini secara otomatis. Namun, saya kira apa yang penulis maksudkan dengan mengabaikan elemen diagonal tidak menghapusnya dari jumlah baris / kolom.

Bahkan, jika saya secara manual menghapus diagonal dari array, maka kedua hasilnya bertepatan, yaitu

import scipy as sp
gebal = sp.linalg.get_lapack_funcs(('gebal'),(A,)) # this picks up DGEBAL
Ab, lo, hi, scaling , info = gebal(A - np.diag(np.diag(A)), scale=1 , permute=1 , overwrite_a=0 )  

memberikan hasil yang sama dengan balance.m(tanpa entri diagonal tentu saja).

Jika ada pengguna Fortran-savy dapat mengkonfirmasi ini dengan memeriksa dgebal.f , saya akan berterima kasih.

EDIT: Hasil di atas tidak menyiratkan bahwa ini adalah satu-satunya perbedaan. Saya juga membuat matriks yang berbeda di mana GEBAL dan keseimbangan. Saya menghasilkan hasil yang berbeda bahkan setelah diagonal dipisahkan.

Saya cukup ingin tahu apa bedanya tetapi tampaknya tidak ada cara untuk mengetahui karena itu adalah matlab built-in dan karenanya kode tertutup.

EDIT2 : Ternyata matlab menggunakan versi LAPACK yang lebih lama (mungkin sebelum 3.5.0) dan pada 2016b mereka tampaknya ditingkatkan ke versi yang lebih baru. Sekarang hasilnya setuju sejauh yang saya bisa uji. Jadi saya pikir itu menyelesaikan masalah. Saya seharusnya mengujinya dengan versi LAPACK yang lebih lama.