I have been writing a control system toolbox from scratch and purely in Python3 (shameless plug :
harold ). 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
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. ]])
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. ]])
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:
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
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.
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 Newton iteration methods on top of the regular solver.