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

eigfortran

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,

Best Answer

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.
Related Question