Algorytm równoważenia macierzy

Pisałem przybornik systemu sterowania od zera i wyłącznie w Python3 (bezwstydna wtyczka:) harold. Na podstawie moich wcześniejszych badań zawsze narzekałem na solver Riccati care.mz powodów technicznych / nieistotnych.

Dlatego piszę własny zestaw procedur. Jednej rzeczy, której nie mogę znaleźć, jest uzyskanie algorytmu równoważenia o wysokiej wydajności, przynajmniej tak dobrego balance.m. Zanim o tym wspomnisz, xGEBALrodzina jest widoczna w Scipy i możesz w zasadzie wywoływać z Scipy w następujący sposób, załóżmy, że masz tablicę 2D typu 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 )

Teraz jeśli użyję poniższej matrycy testowej

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.      ]])

dostaję

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.      ]])

Jeśli jednak przekażę to balance.m, otrzymam

>> 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

Jeśli sprawdzisz wzorce permutacji, będą one takie same, jednak skalowanie jest wyłączone. gebalDaje skalowania jedność natomiast Matlab daje następujące uprawnienia 2: [-5,0,8,0,2].

Najwyraźniej nie używają tej samej maszyny. Próbowałem różnych opcji, takich jak Lemonnier, dwustronne skalowanie Van Doorena, oryginalny Parlett-Reinsch, a także kilka innych mniej znanych metod w literaturze, takich jak gęsta wersja SPBALANCE.

Być może chciałbym podkreślić, że znam pracę Bennera; w szczególności Symplectic Balance of Hamiltonian Matrices specjalnie do tego celu. Należy jednak pamiętać, że ten rodzaj leczenia odbywa się w ramach gcare.m(uogólnionego solvera Riccati), a równoważenie odbywa się bezpośrednio za pośrednictwem balance.m. Dlatego byłbym wdzięczny, gdyby ktoś mógł wskazać mi faktyczną implementację.

Ujawnienie: Naprawdę nie próbuję odwracać kodu matematyki: naprawdę chcę uciec od niego z różnych powodów, w tym motywacji tego pytania, to znaczy, nie wiem, co robi, co kosztowało mnie bardzo dużo czasu z powrotem w ciągu dnia. Moim zamiarem jest uzyskanie zadowalającego algorytmu równoważenia, który pozwala mi przekazywać przykłady CAREX, dzięki czemu mogę zaimplementować metody iteracyjne Newtona na zwykłym rozwiązaniu.

Zajęło mi to trochę czasu, aby to rozgryźć i jak zwykle staje się to oczywiste po znalezieniu sprawcy.

Po sprawdzeniu problematycznych przypadków zgłoszonych w David S. Watkins. Przypadek, w którym równoważenie jest szkodliwe. Elektron. Trans. Liczba Anal, 23: 1–4, 2006, a także dyskusja tutaj (oba cytowane w arXiv: 1401.5766v1 ), okazuje się, że matlab używa równoważenia, oddzielając najpierw elementy ukośne.

Moją pierwszą myślą było, zgodnie z klasyczną ograniczoną dokumentacją funkcji LAPACK, GEBAL wykonał to automatycznie. Myślę jednak, że autorzy rozumieją, ignorując elementy ukośne, nie usuwając ich z sum wierszy / kolumn.

W rzeczywistości, jeśli ręcznie usunę przekątną z tablicy, oba wyniki będą zbieżne, to znaczy

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 )  

daje taki sam wynik jak balance.m(bez oczywiście przekątnych wpisów).

Jeśli jakikolwiek użytkownik Fortran-savy może to potwierdzić, sprawdzając plik dgebal.f , byłbym wdzięczny.

EDYCJA: Powyższy wynik nie oznacza, że ​​jest to jedyna różnica. Skonstruowałem także różne macierze, w których GEBAL i balance.m dają różne wyniki nawet po rozdzieleniu przekątnych.

Jestem dość ciekawy, jaka może być różnica, ale wydaje się, że nie ma sposobu, aby wiedzieć, ponieważ jest to wbudowany matlab, a zatem zamknięty kod.

EDIT2 : Okazuje się, że Matlab używał starszej wersji LAPACK (prawdopodobnie wcześniejszej niż 3.5.0) i do 2016b wydaje się, że zostały zaktualizowane do nowszej wersji. Teraz wyniki są zgodne, o ile mogę przetestować. Myślę więc, że to rozwiązuje problem. Powinienem był przetestować go ze starszymi wersjami LAPACK.