Matrix Balancing Algorithm

Question

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

Answer 1

Took me quite a while to figure this out and as usual it becomes obvious after you find the culprit.

After checking the problematic cases reported in David S. Watkins. A case where balancing is harmful. Electron. Trans. Numer. Anal, 23:1–4, 2006 and also the discussion here (both being cited in arXiv:1401.5766v1), it turns out that matlab uses the balancing by separating the diagonal elements first.

My initial thought was, as per the classical limited documentation on LAPACK functions, GEBAL performed this automatically. However, I guess what the authors mean by ignoring the diagonal elements is not removing them from the row/column sums.

In fact, if I manually remove the diagonal from the array, then both results coincide, that is

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 )  

gives the same result as balance.m (without the diagonal entries of course).

If any Fortran-savy user can confirm this by checking dgebal.f, I would be grateful.

EDIT: Above result does not imply that this is the only difference. I have constructed also different matrices where GEBAL and balance.m produces different results even after the diagonals are separated.

I am quite curious what the difference might be but it seems there is no way to know since it is a matlab built-in and hence closed code.

EDIT2 : It turns out matlab was using an older version of LAPACK (probably pre 3.5.0) and by 2016b they seem to be upgraded to the newer version. Now results agree as far as I can test. So I think that settles the issue. I should have tested it with older LAPACK versions.