Repository
I have created a repository, you can take a
look at my work here:
Preliminary balancing
Currently in octave eig uses preliminary
balancing and it can’t be turned off, while in Matlab[1] it can be.
§ e = eig(A)
§ [V,D] =
eig(A)
§ [V,D,W] =
eig(A)
§ e =
eig(A,B)
§ [V,D] =
eig(A,B)
§ [V,D,W] =
eig(A,B)
§ [___] = eig(A,balanceOption)
§ [___] =
eig(A,B,algorithm)
§ [___] =
eig(___,eigvalOption)

§ lambda =
eig (A)
§ lambda =
eig (A, B)
§ [V,
lambda] = eig (A)
§ [V,
lambda] = eig (A, B)

The aim of this task was to change the *geev LAPACK calls to the extended *geevx which allows us to turn off the balancing. The
ability to do this is important becuase the balancing can be harmful:
“Balancing sometimes seriously
degrades accuracy. In particular, one should not balance a matrix after it has
been transformed to Hessenberg form. However, we must emphasize that balancing
is usually not harmful and often very beneficial. When in doubt, balance.” [2]
Testing
To test the preliminary balancing I used Matlab’s
example for „Eigenvalues of a Matrix Whose Elements Differ Dramatically in
Scale” [3]
A = [ 3.0
2.0 0.9 2*eps;
2.0 4.0 1.0 eps;
eps/4 eps/2 1.0 0;
0.5 0.5 0.1 1.0];
[VN,DN] = eig(A,'balance');
A*VN  VN*DN
But regardless of the balancing is on or
off the results were the same.
4.4409e16 2.2204e16
2.9055e16 1.6653e16
8.8818e16 2.7756e16
3.1094e17 1.3878e16
1.7176e17 1.5361e18
6.6297e18 0.0000e+00
6.9389e17 4.4409e16
2.2204e16 6.9389e17
(Which is the
more accurate one, the one we want)
It seems like the matrix is not balanced.
Moreover balance(A) is the same as A:
balance(A)A
ans =
0
0 0 0
0
0 0 0
0
0 0 0
0
0 0 0
The reason for this is that the stopping
criteria for balancing was changed in LAPACK 3.5.0 [4], so it deals more intelligently with
cases where the balancing causes loss of accuracy. [5][6][7]
“However, for the case where A is
dense and poorly scaled, the new algorithm will still balance the matrix and
improve the eigenvalue condition number. If accurate eigenvectors are desired,
then one should consider not balancing the matrix.” [5]
Other matrix
A = [3 2 0.9 0; 2 4 1 0; 0 0 1 0; 0.5 0.5 0.1 1]; % drop the eps terms for
now
Balancing this matrix is not harmful, so LAPACK
will balance as we expect it to do.
A = [3 2 0.9 0; 2 4 1 0; 0 0 1 0; 0.5 0.5 0.1 1]; % drop the eps terms for
now
[VN,DN] = eig(A,'nobalance');
A*VN  VN*DN
[VN,DN] = eig(A,'balance');
A*VN  VN*DN
[VN,DN] = eig(A);
A*VN  VN*DN
ans =
4.4409e16 2.2204e16
2.4947e16 1.6653e16
8.8818e16 3.8858e16
4.1531e17 1.3878e16
0.0000e+00 0.0000e+00
0.0000e+00 0.0000e+00
8.3267e17 4.4409e16
2.2204e16 4.1633e17
ans =
0.0000e+00 4.4409e16
2.2204e16 1.6653e16
0.0000e+00 0.0000e+00
2.2204e16 1.3878e16
0.0000e+00 0.0000e+00
0.0000e+00 0.0000e+00
0.0000e+00 1.3878e17
0.0000e+00 0.0000e+00
ans =
0.0000e+00 4.4409e16
2.2204e16 1.6653e16
0.0000e+00 0.0000e+00
2.2204e16 1.3878e16
0.0000e+00 0.0000e+00
0.0000e+00 0.0000e+00
0.0000e+00 1.3878e17
0.0000e+00 0.0000e+00
On this example we can see that the eig(A) is
working as the eig(A,'balance') which is a behaviour we want.
So the ability to turn off the balancing is not
as important as it was before LAPACK 3.5.0. but there are some cases when the balancing is
bad, so the ability to turn off the balancing is still important.
“However, for the case where A is
dense and poorly scaled, the new algorithm will still balance the matrix and
improve the eigenvalue condition number. If accurate eigenvectors are desired,
then one should consider not balancing the matrix.” [5]
Hi Barbara. Is there a way to understand which version of LAPACK Octave uses, from inside Octave?
ReplyDeleteMarco Caliari
You can try "ldd" to see the libraries linked to the Octave executable, but in general it is difficult. There was discussion about it in the bug tracker: https://savannah.gnu.org/bugs/?func=detailitem&item_id=45659
ReplyDelete