core icon indicating copy to clipboard operation
core copied to clipboard

apf::eigen not converging

Open zhang-alvin opened this issue 7 years ago • 2 comments

I'm trying to get the eigenvectors and eigenvalues for a simple matrix using the apf::eigen() function:

  apf::Matrix3x3 testMat;
  testMat[0][0] = 1.0;
  testMat[0][1] = 0.0;
  testMat[0][2] = 15.0;
  testMat[1][0] = 0.0;
  testMat[1][1] = 1.0;
  testMat[1][2] = 0.0;
  testMat[2][0] = 0.0;
  testMat[2][1] = 2.0;
  testMat[2][2] = 5.0;

  apf::Vector3 eigenVectors[3];
  double eigenValues[3];
  apf::eigen(testMat, eigenVectors, eigenValues);

For this simple case, I would run into an assertion error at https://github.com/SCOREC/core/blob/master/apf/apfMatrix.cc#L76 because the algorithm had not converged. This seems abnormal since Octave finds the eigenvalues to be integers $\lambda = {1,5,1}$.

I compiled this code with the core libraries loaded via module on the SCOREC machines.

I also noticed that the maximum number of iterations is hard-coded as maxiters = 100 which seems like something that the user should decide.

zhang-alvin avatar Nov 17 '17 21:11 zhang-alvin

I had forgotten about this issue. I remember talking to @mortezah about this and one limitation of apf::eigen is that the input matrix must be symmetric. The input might also have to be invertible. So I think what would be great is really a better error message that details such limitations.

zhang-alvin avatar Aug 17 '18 17:08 zhang-alvin

Bringing this issue back from the dead a second time. In my specific case, which does have a symmetric matrix, changing maxiters = 100 to maxiters = 1000 does not change the matrix which is supposed to converge.

The culprit appears to be the convergence criterion used: https://github.com/SCOREC/core/blob/2ca5177fd172dd310a9412ad9b549b8afe291b2c/mth/mthQR.cc#L207-L218 In lines 211-212, the limit is set to an absolute value 1e-10. Instead the value should probably be relative to some value dependent on the matrix in question, here a. I'm looking for suggestions as to what that should be.

This does not work for the specific case @zhang-alvin mentioned, but that could just be because that matrix is not symmetric.

AjinkyaDahale avatar Aug 02 '19 16:08 AjinkyaDahale