Computing min and max eigenvalues

[Michał Wichrowski]

Dear all,

I need to compute the minimum and maximum eigenvalues of the matrix-free operator (vmult is only available ). For the largest eigenvalue, I'm using the power method, but I have some problems with minimal eigenvalue. The condition number of the operator is close to 1, so I tried using the power method on inverse obtained by CG. This solution is time-consuming and for some cases I got strange results.

I have seen that CG can output eigenvalues: 

  solver_outer.connect_eigenvalues_slot( [](const std::vector<  double > & vec)->void { std::cout<< "  "<<vec.front()<< " , " <<vec.back()<< "  size: "<< vec.size()<<std::endl; } , false);

  solver_outer.solve (shat, sol, rhs,

PreconditionIdentity()

  );

but the estimate of largest eigenvalue obtained this way does not  match one from power method: 

Maximum eigenvalue from power method: 1.001

Maximum eigenvalues from CG:   0.0071333 , 1793.58  size: 249

Michał

[Bruno Turcksin]

Michal,

Have you tried to use ARPACK to get the smallest eigenvalues? That's what we are doing in one of our code and it works pretty well. We have a patch here so that you don't need to provide the inverse of the operator to use ARPACK. You still need to provide a mass matrix but it's no used.

Best,

Bruno

[Michał Wichrowski]

W dniu czwartek, 18 października 2018 17:59:40 UTC+2 użytkownik Bruno Turcksin napisał:

Michal,

Have you tried to use ARPACK to get the smallest eigenvalues? That's what we are doing in one of our code and it works pretty well. We have a patch here so that you don't need to provide the inverse of the operator to use ARPACK. You still need to provide a mass matrix but it's no used.

Best,

Bruno

Thanks a lot!

In the meantime, I found that  I was probably using the wrong solver to get inverse (the system may be not SPD) using GMRes resulted in a significant speedup. However, I still have strange results, including the smallest eigenvalue being larger than the largest one. I think using ARPACK instead of my own code will solve the problem, thank you.

Michał

[Michał Wichrowski]

Dear Bruno,

I've used ARPACK and the computation time has dropped signifficantly. But I have some strange results, I have following code:

  IterationNumberControl arpack_control(30, 1e-8);

  typename ArpackSolver::AdditionalData arpack_data_smallest(30,

                                               ArpackSolver:: smallest_magnitude,

                                               false);

  typename ArpackSolver::AdditionalData arpack_data_largest(30,

                                               ArpackSolver:: largest_magnitude,

                                               false);

    ArpackSolver arpack_smallest (arpack_control,arpack_data_smallest);

    arpack_smallest.solve(shat,mass, s_hatS_inverse, eigenvalues, eigenvectors_small, 1 );

    pcout<<"\t Aparck small "<<"coverged in "

        <<arpack_control.last_step()<<"\t "

<<" Residual: "

<< arpack_control.last_value()<< " \t eigenvalue:  "<< eigenvalues[0]<<std::endl;

    std::vector< LevelVectorType > eigenvectors_large(2, rhs);

    ArpackSolver arpack_largest (arpack_control,arpack_data_largest);

    arpack_largest.solve(shat,mass, s_hatS_inverse, eigenvalues, eigenvectors_large, 1 );

    pcout<<"\t Aparck large "<<"coverged in "

        <<arpack_control.last_step()<<"\t "

<<" Residual: "

<< arpack_control.last_value()<< " \t eigenvalue:  "<< eigenvalues[0]<<std::endl;

    typename   EigenPower<LevelVectorType>::AdditionalData data(0.);

  EigenPower<LevelVectorType> eigen(solver_control_eigen,mem,data );

  double l_max=1;

  double l_min=0;

//   rhs=1;  //starting from 1 gives same results

  rhs= eigenvectors_large[0];

  internal::PreconditionChebyshevImplementation::set_initial_guess(rhs);

  eigen.solve(l_max,shat, rhs );

  sol=1;

    pcout<<"\t Eigen power  "<<"coverged in "

        <<solver_control_eigen.last_step()<<"\t "

<<" Residual: "

<< solver_control_eigen.last_value()<<std::endl;

The output is:

Aparck small coverged in 4294967295   Residual: nan eigenvalue:  (1.17207,0) 

Aparck large  coverged in 4294967295   Residual: nan eigenvalue:  (0.714841,0)

Eigen power  coverged in 29   Residual: 8.71076e-05

  Eigen value:1.17147

1) The number of iterations and residual are obviously strange. I guess that ARPACK solver does not update SlolverControl object

2) The largest eigenvalue is computed instead of the largest one, and, I guess, the smallest eigenvalue is computed instead of the largest one. 

3) I get complex eigenvalues in some cases and I am quite sure, that all eigenvalues of the system are real.

Michał

[Michał Wichrowski]

PS. 

I noticed that in some cases the "smallest" eigenvalue does not match maximum eigenvalue obtained from power method.

[Bruno Turcksin]

Michal,

Le ven. 19 oct. 2018 à 06:33, Michał Wichrowski
<mtwich...@gmail.com> a écrit :

> 1) The number of iterations and residual are obviously strange. I guess that ARPACK solver does not update SlolverControl object

That's correct. SolverControl is used to set ARPACK's parameters but
it is not updated by ARPACK so the number of iterations and the
residual are garbage.

> 2) The largest eigenvalue is computed instead of the largest one, and, I guess, the smallest eigenvalue is computed instead of the largest one.

If you use ARPACK through deal without patching it, then the only
method available is shift-inverse. The thing is that when you know the
largest eigenvalue of the inverse you can get the smallest eigenvalue
of the operator. So yes, the names are reversed.

> 3) I get complex eigenvalues in some cases and I am quite sure, that all eigenvalues of the system are real.

I've not had this problem. Maybe it's because of numerical error?

> I noticed that in some cases the "smallest" eigenvalue does not match maximum eigenvalue obtained from power method.

Are you sure that the power method has converged? Also, try using
ARPACK's regular mode instead of shift-inverse. I think it might be
better to find the largest eigenvalue.

Best,

Bruno

[Michał Wichrowski]

Dear Bruno

> 2) The largest eigenvalue is computed instead of the largest one, and, I guess, the smallest eigenvalue is computed instead of the largest one.
If you use ARPACK through deal without patching it, then the only
method available is shift-inverse. The thing is that when you know the
largest eigenvalue of the inverse you can get the smallest eigenvalue
of the operator. So yes, the names are reversed.

Yup,  using the patch solved this issue.  At least the smallest eigenvalue is smaller than the largest.

> 3) I get complex eigenvalues in some cases and I am quite sure, that all eigenvalues of the system are real.
I've not had this problem. Maybe it's because of numerical error? 

Eigenvalue:  (1.03342,0.0069945) with arpack tolerance 1e-8. Is it possible?

> I noticed that in some cases the "smallest" eigenvalue does not match maximum eigenvalue obtained from power method.
Are you sure that the power method has converged? Also, try using
ARPACK's regular mode instead of shift-inverse. I think it might be
better to find the largest eigenvalue.

 The deal.II power method converged with residual: 9.5052e-05 I switched to regular mode and obtained 1.04192 from arpack and 1.03981 from power method. The tolerance of arpack was set to 1e-8. The difference is growing with mesh size,  so I assume the tolerance is set for eigenvector, resulting in much smaller tolerance for the eigenvalue.

Also, I got the smallest eigenvalues equal to  0 in some cases, while I'm sure that the matrix is not singular. I can solve a linear system with the same matrix using  Krylov solvers without a problem.

Michał

[Earl Fairall]

I am posing this question purely out of personal interest: what do you mean by maximum and minimum eigenvalue? I am confused by this.

Earl

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Michał Wichrowski]

Earl,

I am posing this question purely out of personal interest: what do you mean by maximum and minimum eigenvalue? I am confused by this.

 

The eigenvalues with largest and smallest magnitude. My problem is symmetric positive-definite, thus the spectrum is real and positive, so in my case "minimum eigenvalue" makes sense. 

Bruno,

I think I finally have a solution. By increasing Arnoldi space size, I obtained results that seems ok.

I think the documentation of arpack solver is not valid: 

https://www.dealii.org/developer/doxygen/deal.II/classArpackSolver.html#afdc3aa9d761c43b5b4132e731c6191a5

A	The operator for which we want to compute eigenvalues. Actually, this parameter is entirely unused.Actually, this parameter is entirely unused.

The matrix A is obviously used. 

Also, calling arpack with WhichEigenvalues::both_ends always results in a crash.

Thanks for help,

Michał

[Wolfgang Bangerth]

On 10/25/2018 05:23 AM, Michał Wichrowski wrote:
>
> I think the documentation of arpack solver is not valid:
> https://www.dealii.org/developer/doxygen/deal.II/classArpackSolver.html#afdc3aa9d761c43b5b4132e731c6191a5
> A The operator for which we want to compute eigenvalues. Actually, this

> parameter is entirely unused.*Actually, this parameter is entirely unused*.

>
> The matrix A is obviously used.
>
> Also, calling arpack with WhichEigenvalues::both_ends always results in
> a crash.

Would you mind opening github issues about both of these?

Thanks
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/