3
$\begingroup$

If we call LAPACK DGGEV or DGGEVX routines for two badly-conditioned matrices in a C++ code, will we get the same eigen-values & vectors as MATLAB's eig command with 'qz' algorithm like: eig(a,b,'qz')?

Is there any possibility that eigenvectors of MATLAB are different from LAPACK DGGEV or DGGEVX?

$\endgroup$
  • $\begingroup$ LAPACK calls BLAS routines to do things like matrix multiplication. Different implementations of these BLAS routines do the multiplies and adds in different order, and this can result in slightly different results from the BLAS routines given the same inputs. It's likely that the differences in results that you're seeing are do to such differences in the BLAS used by MATLAB and your C++ code. $\endgroup$ – Brian Borchers Nov 19 '14 at 2:49
  • 4
    $\begingroup$ In any case, you really don't have any reason to complain as long as the eigenvectors returned by both codes are correct. Remember that with repeated eigenvalues there are going to be infinitely many sets of eigenvectors that span the same eigenspaces, so the eigenvectors are not uniquely defined. $\endgroup$ – Brian Borchers Nov 19 '14 at 2:50
  • 1
    $\begingroup$ Have you checked that the the eigenvectors returned by the different routines are correct? That is, $Ax=\lambda x$? $\endgroup$ – Brian Borchers Nov 19 '14 at 19:13
2
$\begingroup$

A few things to note:

  • By definition A·v = λ·v, eigenvectors are not unique. You can multiply by any constant and still get another valid eigenvector.

    The convention in MATLAB is that for eig(A), the eigenvectors are scaled so that the norm of each is 1.0, and for eig(A,B), the eigenvectors are not normalized (see here for an example). Here is the relevant part in the documentation:

    V: right eigenvectors, returned as a square matrix whose columns are the right eigenvectors of A or generalized right eigenvectors of the pair (A,B). The form and normalization of V depends on the combination of input arguments:

    • [V,D] = eig(A) returns matrix V, whose columns are the right eigenvectors of A such that A*V = V*D. The eigenvectors in V are normalized so that the 2-norm of each is 1.

    • [V,D] = eig(A,'nobalance') also returns matrix V. However, the 2-norm of each eigenvector is not necessarily 1.

    • [V,D] = eig(A,B) and [V,D] = eig(A,B,algorithm) returns V as a matrix whose columns are the generalized right eigenvectors that satisfy A*V = B*V*D. The 2-norm of each eigenvector is not necessarily 1. In this case, D contains the generalized eigenvalues of the pair (A,B), along the main diagonal.

      If A is symmetric and B is symmetric positive definite, then the eigenvectors in V are normalized so that the 2-norm of each is 1.

  • In addition, eigenvalues are not sorted. You are only guaranteed that the columns of V are the corresponding right eigenvectors to the eigenvalues in D. That's not the same as svd.

    In fact, there is no total ordering of complex numbers. The convention in MATLAB is that the sort function sorts complex elements first by magnitude (i.e. abs(x)), then by phase angle on the [-pi,pi] interval (i.e. angle(x)) if magnitudes are equal.

Having considered the above, when solving the generalized eigenvalue problem, you should be getting similar results to MATLAB by using DGGEV routine in LAPACK.

$\endgroup$
  • $\begingroup$ The relevant part in the documentation says If A is symmetric and B is symmetric positive definite, then the eigenvectors in V are normalized so that the B-norm of each is 1. I wonder what they mean by B-norm is it the same as 2-norm? $\endgroup$ – user3853917 Nov 27 '14 at 12:24
  • $\begingroup$ Does B-norm mean that Z**T*B*Z = I, mentioned here in DSYGV documentation? $\endgroup$ – user3853917 Nov 27 '14 at 12:28
  • 1
    $\begingroup$ huh, good observation! I just assumed it was a typo, and as you can see corrected it as 2-norm in my quote above. I'm not so sure now... Anyway see here to understand the naming convention of LAPACK routines: xyyzzz where x is the type (single, double, complex), yy is the type of the matrix (general, symmetric, diagonal, etc..), and zzz is the computation performed (e.g: EV, GV, SVX, etc..). So if there is a version of the algorithm specialized for the matrix shape you have you should use it instead (say symmetric positive-definite) $\endgroup$ – Amro Nov 27 '14 at 13:40
  • $\begingroup$ Thanks, I asked another question here which is kind of relevant. $\endgroup$ – user3853917 Nov 27 '14 at 14:09
1
$\begingroup$

The best way to find the answer to this question was to experiment:

We calculated the Eigen values & vectors of two badly-conditioned matrices with MATLAB R2013b by using eig command like eig(a,b,'qz'). The ill-conditioned matrices were 342x342 in size.

Based on the explanations on the following link, MATLAB R2009a implemented LAPACK DGGEV for doing eig(a,b,'qz'):

https://stackoverflow.com/questions/26808510/which-algorithm-do-dggev-or-dsygv-eigen-solvers-in-lapack-implement-is-it-qz

Therefore, we implemented DGGEV (and also DGGEVX) routines from LAPACK in our C++ code by using LAPACKE which is a "C" language wrapper for LAPACK developed by INTEL-MKL team.

The conclusion was that we got exactly the same Eigen values & vectors as MATLAB eig(a,b,'qz') precisely.

Please note that if you use DGGEVX (expert version of DGGEV), you should turn the balancing off by using 'N' for BALANC argument, to get exactly the same Eigen vectors as MATLAB 'qz' algorithm.

We also checked the correctness of Eigen values & vectors by checking if they satisfy A*V-B*V*D=0, and the coclusion was that Eigen values & vectors were correct, as the residual error was close enough to zero.

$\endgroup$
  • $\begingroup$ Please move the actual answer to your question into your answer. Also, if you could condense the original question into a short summary, and move some of your learnings into your answer, that would also help with readability. Generally speaking, the practice of repeat edits where the OP just appends more and more information does two things: (1) it triggers a flag because you're editing so much, and (2) the result looks bad, because it has accreted several layers, so it's tough to follow. One or two added sections is fine; more than 10 edits adding sections is too much. $\endgroup$ – Geoff Oxberry Nov 20 '14 at 2:54
  • $\begingroup$ @GeoffOxberry Done. $\endgroup$ – user3853917 Nov 21 '14 at 2:02
  • $\begingroup$ Does anyone know how to get the 'chol' version of [V, lam]=eig(e) that matlab uses by default from DGGEV? $\endgroup$ – ggb667 Feb 17 '17 at 0:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.