Repository navigation
linalg.matrix_balance gives wrong transformation matrix #7258
Description
The documentation for linalg.matrix_balance states that
B, T = matrix_balance(A, permute=True, scale=True, separate=False, overwrite_a=False)
should yield matrices B and T such that
B=T^-1 A T
However,
On my machine, the following code fails
import scipy.linalg as linalg
import numpy as np
A = np.array([[-0.5, 0., 0., 0. ],
[ 0., -1., 0. , 0. ],
[ 1., 0., -0.5, 0. ],
[ 0., 1., 0., -1. ]])
B, T = linalg.matrix_balance(A)
print 'T^-1 * A * T - B = '
print linalg.solve(T, A).dot(T) - B
# This should be all zeros, but it is not.
The output of running this code is
[[-0.5 -1. 0. 0. ]
[ 1. -0.5 0. 0. ]
[ 0. 0. 0.5 -1. ]
[ 0. 0. 1. 0.5]]
I think the issue lies in the permutation, not the scaling. Basically, I think the output matrix is defined as
T = DP
where D is the diagonal scaling matrix, and P is the permutation matrix. I think to correct the issue, we should have
T = D*P^T.
This is illustrated in the following code:
A2 = np.array([[-3., 0., 1., 0.],
[-1., -1., -0., 1.],
[-3., -0., -0., 0.],
[-1., -0., 1., -1.]])
print '------------------'
B2, T2 = linalg.matrix_balance(A2)
print 'T2^-1 * A2 * T2 - B2 = '
print linalg.solve(T2, A2).dot(T2) - B2
# This is non zero on my machine
B2, (d,p) = linalg.matrix_balance(A2, separate=True)
D = np.diag(d)
P = np.eye(4)[p]
T_fix = D.dot(P.T)
print 'T_fix^-1 * A2 * T_fix - B2 ='
print linalg.solve(T_fix, A2).dot(T_fix) - B2
# This is zero on my machine
Activity
ilayn commentedon Apr 4, 2017
Thanks for reporting it. I will take a look as soon as I have time today. It is possible that just the documentation is misleading too or there is a mistake with the generation of the outer factor T from the permutation.
ilayn commentedon Apr 4, 2017
If there are examples that work without any problems that's an indication that LAPACK is performing chain permutations and just using indices are not sufficient.
What I mean is if row 1 is swapped with row 3 and then row 2 is swapped row 3 (which is now row 1) you get a permutation that is not possible to correct with straightforward indexing.
some group theoretical explanation is possible but I thought I have tackled this issue. Apparently not :)
FIX: Corrected the transformation matrix permutation
FIX: Corrected the transformation matrix permutation