Incorrect eigenvalues for a symmetric positive-definite matrix

Summary

It appears I am getting incorrect eigenvalues for a symmetric positive-definite matrix. It results in incorrect eigenvectors which end up being linearly dependant (as if the matrix was defective). This impacts the turbulentInletDFSEM as eddy objects lose a dimension.

Steps to reproduce

Here's a minimal example, I concocted:

#include "fvCFD.H"

int main(int argc, char *argv[])
{
    Info<< "\nCalculating eigenvalues of a symmetric positive-definite matrix\n";

    symmTensor A = symmTensor(
        0.005,
        -0.0003,
        0.001,
        0.002,
        0.0001,
        0.002);
    vector lambda(eigenValues(A));
    Info<< "Eigenvalues (OpenFOAM): " << lambda << endl;
    Info<< "Eigenvalues (numpy.linalg): (0.00532289, 0.00160802, 0.0020691)" << endl;
    return 0;
}

What is the current bug behaviour?

Output I am getting:

Eigenvalues (OpenFOAM): (0.00183856 0.00183856 0.00532289)
Eigenvalues (numpy.linalg): (0.00532289, 0.00160802, 0.0020691)

What is the expected correct behaviour?

See above. I compared it with numpy.linalg.

Environment information

• OpenFOAM version : v2006
• Operating system : Arch Linux
• Hardware info : Intel(R) Core(TM) i7-4810MQ CPU @ 2.80GHz
• Compiler : gcc :

Possible fixes

Use alternatives eigensolvers. I currently went with Eigen but I am looking for another OpenFOAM native solution.

changed title from Incorrect eigenvalues for symmetric positive-definite matrix to Incorrect eigenvalues for a symmetric positive-definite matrix

assigned to @kuti

I found OpenFOAM native solution with EigenMatrix:

#include "fvCFD.H"
#include "EigenMatrix.H"

int main(int argc, char *argv[])
{
    Info<< "\nCalculating eigenvalues of a symmetric positive-definite matrix\n";

    symmTensor A = symmTensor(
        0.005,
        -0.0003,
        0.001,
        0.002,
        0.0001,
        0.002);
    vector lambda(eigenValues(A));
    SquareMatrix<scalar> As(3);
    As(0, 0) = 0.005;
    As(0, 1) = -0.0003;
    As(0, 2) = 0.001;
    As(1, 0) = -0.0003;
    As(1, 1) = 0.002;
    As(1, 2) = 0.0001;
    As(2, 0) = 0.001;
    As(2, 1) = 0.0001;
    As(2, 2) = 0.002;
    EigenMatrix<scalar> D(As, true);
    Info<< "Eigenvalues w. OpenFOAM eigenValue: " << lambda << endl;
    Info<< "Eigenvalues w. OpenFOAM EigenMatrix: " << D.EValsRe() << endl;
    Info<< "Eigenvalues w. numpy.linalg: (0.00532289, 0.00160802, 0.0020691)" << endl;
    return 0;
}

The eigenvalues are not sorted, but the numerical values coincide with numpy.linalg.eig.

I haven't read into details of how this cubicEqn is solved but this substitution makes me highly suspicious that this is not general enough. I am not sure if this in some near term plans but it might be a good idea to reimplement eigenValues with EigenMatrix class. I am going to do it and test on my local version for now.

Many thanks for your report, suggestions and contribution, @robertsawko.

• In the first test of yours, the eigenvalues were computed by an implementation of the closed-form solution of the eigenvalue problem - a direct solution. Such errors are expected, I'm afraid, and such errors will always occur from that algorithm without any discernible pattern. The reason is the difference between the closed-form equations involving real numbers and the implementation of the closed-form equations involving floating-point numbers (Wikipedia.) The advantage of this method is its relatively low computational cost.
• In the second test of yours, the eigenvalues were computed by an iterative method (i.e. EigenMatrix). numpy (through LAPACK) does compute eigen decompositions only in an iterative manner. More robust, still cheap, but new in OpenFOAM (6-month old whereas the first impl is there at least for 10 years), therefore did not propagate across the entire code.

May be it is a good time to refactor things.

Yet the first error will keep appearing in one way or another. No generic cure for that. Though, I am happy to be proven incorrect :).

Thanks for a quick reply.

Yes. I realised that that eigenValues uses a direct method and the other approaches are iterative. What I didn't understand though is how this direct approach solves the characteristic polynomial. From the quick glance at it, I don't think it is using Cardano formulas but rather uses some substitution which makes the equation quadratic. That's what made it a little suspcious to my mind.

I confess, I don't remember Cardano formulas so well now, but I am happy to revise. If solving that quadratic equation is essentially the same as Cardano, then the problem must be related to computing in finite precision?

yes, it is due to the floating-point arithmetic.

having said that we had a tiny change in the routines - if possible, try to test the case with the develop branch.

many thanks.

Ok. Will test the dev branch today then and let you know.

By the way, note that the equation in the documentation don't render quite right in doxygen(?). The line break in two equations is omitted and they get amalgamated into a single line.

Hi @robertsawko,

Could you please update us whenever you complete your test, if possible?

Many thanks.

Sorry. I am stuck with recompilation of the branch. I've already been there. I think it's the parallel build which freezes machine. Even when I set the WM_NCOMPPROCS=1. I remember having to manually update some makefiles to fix it.

Oh, I think I may know the reason.

At the time of writing, develop branch (HEAD:7f17a71f) could not be compiled. If you can wait a day or two (or even the new version v2012), it will be OK.

Sorry for forgetting to notify you.

Dear @robertsawko, the develop branch (HEAD:9ee0023b, at the time of writing) can be compiled, if need be.

It crashes for me at:

g++ -std=c++11 -m64 -pthread -DOPENFOAM=2012 -DWM_DP -DWM_LABEL_SIZE=32 -Wall -Wextra -Wold-style-cast -Wnon-virtual-dtor -Wno-unused-parameter -Wno-invalid-offsetof -Wno-attributes -Wno-unknown-pragmas  -O3  -DNoRepository -ftemplate-depth-100 -I/home/rsa/OpenFOAM/openfoam/build/linux64GccDPInt32Opt/src/OpenFOAM -DHAVE_LIBZ -iquote. -IlnInclude -I/home/rsa/OpenFOAM/openfoam/src/OpenFOAM/lnInclude -I/home/rsa/OpenFOAM/openfoam/src/OSspecific/POSIX/lnInclude   -fPIC -c primitives/functions/Function1/Function1/function1Base.C -o /home/rsa/OpenFOAM/openfoam/build/linux64GccDPInt32Opt/src/OpenFOAM/primitives/functions/Function1/Function1/function1Base.o
make: *** No rule to make target 'lnInclude/SquareI.H', needed by '/home/rsa/OpenFOAM/openfoam/build/linux64GccDPInt32Opt/src/OpenFOAM/primitives/functions/Function1/makeFunction1s.C.dep'.  Stop

Hi @robertsawko - this often happens when a file which was previously required (in this case SquareI.H) has been removed entirely. The old dependency check still tries to see if it changed, but since it is gone make figures it must be an error. In various places we try to provide a dummy file to help with that sort of thing. In the case of Function1, there were several changes with refactoring so essentially have to rebuild that bit anyhow.

Can use what Kuti mentioned (remove the offending .dep file), or use the command wrmdep SquareI.H to remove depend files mentioning SquareI.H. Sometimes I use one approach, sometimes the other.

may the following help?

rm -rf /home/rsa/OpenFOAM/openfoam/build/linux64GccDPInt32Opt/src/OpenFOAM/primitives/functions/Function1`

foam

wmakeLnIncludeAll -u

./Allwmake -j -l

There were a few more files I had to treat this way, but it worked. And what's more my eigentest works as well.

Calculating eigenvalues of a symmetric positive-definite matrix
Eigenvalues w. OpenFOAM eigenValue: (0.00160802 0.0020691 0.00532289)
Eigenvalues w. OpenFOAM EigenMatrix: 3(0.00160802 0.0020691 0.00532289)
Eigenvalues w. numpy.linalg: (0.00532289, 0.00160802, 0.0020691)

Could you please elaborate a little what has changed? Just for my own education. The diff on the equation shows the only change in a discriminant check:

<     const scalar discr = (mag(numDiscr) > VSMALL) ? numDiscr : 0;
---
>     const scalar discr = (mag(numDiscr) > sqr(SMALL)) ? numDiscr

closed

Hi @kuti

Since OpenFOAM tensor is fixed 3x3, there's no need to use Eigen as the call overhead is considerably large. While you are fixing the eigenVectors() bug, I'm also looking for external diagonalization solvers and found this:

https://www.mpi-hd.mpg.de/personalhomes/globes/3x3/index.html

It works for both real symm matrices and Hermitian matrices.(OpenFOAM only need the real symm part.) The main idea is to use analytical calculation and falls back to iterative method if the matrix is ill-conditioned.

I've compiled it as an external C library and tested it in OpenFOAM. It's matching MATLAB result under extreme conditions while OpenFOAM eigenValues() can't. Please consider porting this C code in OpenFOAM. Thanks.

Many thanks @xuegy.

Unfortunately, we are aware of the publication, the associated code, and its benefits from the outset. The problem is the lack of time necessary to perform the required tasks for a seamless operation, i.e. integration, testing etc.

Maybe you can integrate and test the module locally - then I believe all the team members are more than happy to merge your contribution into the main code. :)

@kuti Yes I'm willing to do that. It will benefit my current project as well.

The problem is my code will be in "C styled" C++. So even after I port the code, your team still have a lot of work to do.