Monday 22 August 2016

Summary

GSoC 2016 is about to end so in this post I will sum up the work we had done during this period. I will write about 3 major topics: eig, gsvd and chol2inv. The original project was about eig, and when it was done we moved on to the other ones.

eig

The goals of the project were:
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)

Calling forms:


We finished these goals before the end of GSoC, and most of the work is already merged to the main octave repository:
In my repository there are some other to be reviewed commits that reduces code duplication with the use of templates.
You can find some information about the individual tasks on this blog in previous posts, and also here:

gsvd

The function gsvd was imported from Octave-Forge linear-algebra package. However it uses deprecated LAPACK functions: *ggsvd. Moreover it is not compatible with Matlab’s gsvd. For more information check these discussions:
  1. Deprecated LAPACK routines (GSVD)
  2. new gsvd function incompatible w/Matlab
  3. bug #48807: new gsvd function is incompatible w/Matlab
Work already merged to the main Octave repository:
Other commits that are not yet merged:

chol2inv

chol2inv couldn’t handle sparse matrixes (bug # 36437), we submitted a fix for this bug and tests for chol2inv.
Commits:


GSoC 2016 with GNU Octave

This was an exeptional opportunity to gain experience and get into open source development and learn about Octave. I am really happy for this project and the time we spent working on it. I would like to thank for my mentors Carnë Draug and Nir Krakauer, without their help and guidance non of this would have been possible.

Tuesday 12 July 2016

chol/qz and matrix/vector flags

In the past weeks I worked on implementing the chol/qz and matrix/vector flags. My work can be found at my repository:

Calling forms

chol/qz flag

In case of a generalized eigenvalue problem the chol/qz flag can be used in Matlab. The aim of this task was to introduce this functionality to Octave.
The chol/qz flag only matters if your matrices are symmetric as:
Regardless of the algorithm you specify, the eig function always uses the QZ algorithm when A or B are not symmetric. [1]

If you omit the chol/qz flag it works as it worked before:
The algorithm used depends on whether there are one or two input matrices, if they are real or complex, and if they are symmetric (Hermitian if complex) or non-symmetric. [2]
Which is the same as in Matlab:
When you omit the algorithm argument, the eig function selects an algorithm based on the properties of A and B. It uses the 'chol' algorithm for symmetric (Hermitian) A and symmetric (Hermitian) positive definite B. Otherwise, it uses the 'qz' algorithm. [1]
So if it is a symmetric problem than the chol method is used otherwise the qz method is used.

matrix/vector flag

The matrix/vector flag can be used in every previous call forms. With it you can specify the format of the eigenvalues.
If you use the “vector” flag the eigenvalues will be returned in a column vector.
If you use the “matrix” flag the eigenvalues will be returned in a diagonal matrix.
If you omit this flag than the eigenvalues will be returned in the same format as before, which is the same as in Matlab:
The default behavior varies according to the number of outputs specified:
If you specify one output, such as e = eig(A), then the eigenvalues are returned as a column vector by default.
If you specify two or three outputs, such as [V,D] = eig(A), then the eigenvalues are returned as a diagonal matrix, D, by default. [3]


[1] http://www.mathworks.com/help/matlab/ref/eig.html#inputarg_algorithm

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]

Thursday 5 May 2016

Hello World!

I am Barbara Lócsi, a Software Engineering student at the Budapest University of Technology and Economics. During this summer as a GSoC student I will work on Generalized eigenvalue problem. I will blog about my progress here.
You can find more information about me and my application here: 

Here you can see my timeline
  • Community Bonding period (Until May 22)
    • I am already started working on implementing preliminary balancing I would like to finish most of this task before the coding period begins.
  • Week 1-2 (May 23 - Jun 5)
    • Finals, non-coding time
  • Week 3-4 (Jun 6 - Jun 19)
    • Finish preliminary balancing if it is not finished, start working on implementing left eigenvector calculation
  • Midterm evaluations (Jun 20 - Jun 27)
  • Week 6 (Jun 27 - Jul 3)
    • Finish the left eigenvector task (if not finished)
  • Week 7-10 (Jul 4 - Jul 31)
    • algorithm choosing for eigenvalue calculation (chol or qz)
    • creating tests
  • Week 11 (Aug 1 - Aug 7)
    • documenting
  • Week 12-13 (Aug 8 - Aug 28)
    • deciding return value format of the eigenvalues (vector or matrix)
    • testing, documenting