Tuesday 21 June 2016

Midterm evaluations


You can find my work here:


Project goals

Certain calling forms of the eig function are currently missing, including:
  • preliminary balancing
  • computing left eigenvectors as a third output
  • choosing among generalized eigenvalue algorithms
  • choosing among return value formats of the eigenvalues (vector or matrix) see more here.


Calling forms


The aim for this period was:

  • Finish preliminary balancing if it is not finished, start working on implementing left eigenvector calculation


The preliminary balancing task is done, I wrote a blog post about it:
I started working on implementing the left eigenvector calculation and it is mostly done, some tests needs to be added to eig.cc which will be completed until the end of this week.

Remaining tasks for the second period:

  1. choosing among generalized eigenvalue algorithms
  2. choosing among return value formats of the eigenvalues (vector or matrix) see more here.



Tuesday 14 June 2016

Preliminary balancing


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.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 [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.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.” [5]