## Tuesday, 14 June 2016

### 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 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.” 

#### Testing

To test the preliminary balancing I used Matlab’s example for „Eigenvalues of a Matrix Whose Elements Differ Dramatically in Scale” 

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.4409e-16   2.2204e-16  -2.9055e-16  -1.6653e-16
8.8818e-16   2.7756e-16   3.1094e-17  -1.3878e-16
1.7176e-17   1.5361e-18   6.6297e-18   0.0000e+00
-6.9389e-17  -4.4409e-16   2.2204e-16   6.9389e-17

(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 , so it deals more intelligently with cases where the balancing causes loss of accuracy. 

“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.” 

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.4409e-16   2.2204e-16  -2.4947e-16  -1.6653e-16
8.8818e-16   3.8858e-16  -4.1531e-17   1.3878e-16
0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00
-8.3267e-17  -4.4409e-16   2.2204e-16  -4.1633e-17

ans =

0.0000e+00  -4.4409e-16   2.2204e-16   1.6653e-16
0.0000e+00   0.0000e+00  -2.2204e-16  -1.3878e-16
0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00
0.0000e+00   1.3878e-17   0.0000e+00   0.0000e+00

ans =

0.0000e+00  -4.4409e-16   2.2204e-16   1.6653e-16
0.0000e+00   0.0000e+00  -2.2204e-16  -1.3878e-16
0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00
0.0000e+00   1.3878e-17   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.” 

#### 2 comments:

1. Hi Barbara. Is there a way to understand which version of LAPACK Octave uses, from inside Octave?

Marco Caliari

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