Reproducing Matlab's eig() results in Fortran for getting eigenvectors corresponding to zero eigenvalues

I've been working on speeding up some code by rewriting it in Fortran. One of the critical steps is to determine x such that Ax = 0. That is to find the nullspace of A=A'. Matlab seems to do this very well with the eig() function but I cannot reproduce the resulting nullspace eigenvectors with Lapack's DSYEV.

DSYEV has a bug where using "U" (upper diagonal) to store the matrix gives erroneous results. I'm using the lower diagonal option here.

Nonetheless, the vectors produced by eig() and DSYEV corresponding to small eigenvalues (<1e-13) are not the same. Both eig() and DSYEV produce orthonormal vectors that lead to Ax=0 for each vector. But it seems the vectors produced by eig() are much cleaner. They produce better results in the application I'm using them in.

How can I reproduce the results given by eig() so that it can be used in an optimized fortran code which runs faster than MATLAB. I'm keen on running independently of MATLAB. What is eig() doing that DSYEV is not. my matrix does not need any permutations or scaling. I've tried using DGESVD which also does not do as well as eig(). Thanks,

Answer 1

[Matt J]

If the nullity of A is greater than 1, which it seems to be for you, then the null space has an infinite continuum of possible bases. There is no way to guarantee that any eigenvector solver, be it EIG or DSYEV, will reliably make a particular choice. The output will also depend discontinuously on the data. That is, if you make a small perturbation of A, you can get a completely different looking set of eigenvectors.

The real solution to your problem is to come up with a more precise mathematical definition of what you call a "clean eigenvector" and use that definition to modify the math that your application is doing.

[Ali]

Thanks Matt. I should have gotten the hint when both methods produced orthonormal vectors and in both cases Ax=0. The thing is there was no obvious pattern in the angles between the two sets of eigenvectors as obtained by the dot product between the two sets of vectors.

I have to look for another explanation for the difference in the results because both sets of bases should be good. It may be that the code in Fortran for the actual application has a bug.