Algorytm równoważenia macierzy

9

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.

percusse
źródło

Odpowiedzi:

7

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.

percusse
źródło