Matrix Balancing Algorithm
I have been writing a control system toolbox from scratch and purely in Python3 (shameless plug : <PERSON> ). From my past research, I have always complaints about the Riccati solver care.m for reasons that are technical/irrelevant.
Hence, I've been writing my own set of routines. One thing I can't find a way around is to obtain a high-performance balancing algorithm, at least as good as balance.m. Before you mention it, xGEBAL family is exposed in Scipy and you can basically call from Scipy as follows, suppose you have a float type 2D array 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 )
Now if I use the following test matrix
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. ]])
I get
array([[ 8. , 0. , 0. 768
, 2. , 2. ],
[ 0. , 2. , 0.000002, 0. , 0. ],
[ 0. , 0. , 6. , 2. , 2. ],
[ 0. , 0.000002, 0. , 6. , 0. ],
[ 0. , 0. , 0. , 0. , 8. ]])
However, if I pass this to balance.m, I get
>> 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
If you check the permutation patterns, they are the same however the scaling is off. gebal gives unity scalings whereas matlab gives the following powers of 2: [-5,0,8,0,2].
So apparently, these are not using the same machinery. I've tried various options such as Lemonnier, Van Dooren two sided scaling, original Parlett-Reinsch and also some other less-known methods in the literature such as the dense version of SPBALANCE.
One point maybe I might emphasize is that I am aware of <PERSON>'s work; in particular the Symplectic Balancing of Hamiltonian Matrices specifically for this purpose. However, note that this type of treatment is done within gcare.m (generalized <PERSON> solver) and balancing is done directly via balance.m. Hence, I would appreciate if someone can point me to the actual implementation.
Disclosure: I am really not trying to reverse engineer mathworks code: I actually want to get away from it due to various reasons including the motivation of this question, that is to say, I don't know what it is doing which costed me lots of time back in the day. My intention is to get a satisfactory balancing algorithm that allows me to pass CAREX examples such that I can implement <PERSON><PHONE_NUMBER>
0 0.0003 0 6.0000 0
0 0 0 0 8.0000
If you check the permutation patterns, they are the same however the scaling is off. gebal gives unity scalings whereas matlab gives the following powers of 2: [-5,0,8,0,2].
So apparently, these are not using the same machinery. I've tried various options such as Lemonnier, Van Dooren two sided scaling, original Parlett-Reinsch and also some other less-known methods in the literature such as the dense version of SPBALANCE.
One point maybe I might emphasize is that I am aware of Benner's work; in particular the Symplectic Balancing of Hamiltonian Matrices specifically for this purpose. However, note that this type of treatment is done within gcare.m (generalized Riccati solver) and balancing is done directly via balance.m. Hence, I would appreciate if someone can point me to the actual implementation. 768