Discrepancy in Computing the Inf-Sup Constants Using Matlab and deal.II

[Najwa Alshehri]

Dear all,

Thank you for your always support and help.

I am encountering an issue while solving an eigenvalue problem related to the computation of the discrete inf-sup constant for a saddle point problem. Specifically, I am solving the following system:

Bt A^-1 B eigenvector = eigenvalue M eigenvector,

where Bt A^-1 B is clearly symmetric and at least semi-positive definite, and M is a symmetric positive definite matrix.

The problem arises when comparing the inf-sup constants obtained using Matlab and dealii (solving on the same mesh ):

1. Matlab: Using the "eig" solver, the inf-sup constant decays to zero as the mesh is refined.
2. deal.II: Using the "ArpackSolver", the inf-sup constant is bounded away from zero.

The discrepancy is puzzling, as both methods should give consistent results. To illustrate this issue, I have attached a figure comparing the two results and included in this post the two data files obtained from Matlab and deal.II.

I would greatly appreciate any insights or suggestions on what might be causing the discrepancy.

Your help is highly appreciated,

Najwa

[Wolfgang Bangerth]

On 5/30/24 09:55, Najwa Alshehri wrote:
>
> The problem arises when comparing the inf-sup constants obtained using Matlab
> and dealii (solving on the same mesh ):
>

> 1. Matlab: Using the "eig" solver, the inf-sup constant decays to zero as the
> mesh is refined.
> 2. deal.II: Using the "ArpackSolver", the inf-sup constant is bounded away

> from zero.
>
> The discrepancy is puzzling, as both methods should give consistent results.
> To illustrate this issue, I have attached a figure comparing the two results
> and included in this post the two data files obtained from Matlab and deal.II.
>
>
> I would greatly appreciate any insights or suggestions on what might be
> causing the discrepancy.

Najwa:
my first question would be if your matrices are the same. If they are not,
then you should not expect eigenvalues to be the same. So what happens if you
output the relevant matrices from both codes on, say, a 2x2 mesh? On a coarse
enough mesh, you can often even compute the matrix by hand on a piece of paper
within half an hour, so that if the matrices you output are different, you can
figure out which one is right and which one is not.

Best
W.

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

[Matthias Maier]

Do you use the same "mass matrix" in both cases, meaning do you call
eig(A, M) where A is the stiffness matrix and M is the mass matrix in
matlab?


You will need to compute the spectrum of the eigenvalue problem

A x = \lambda M x


Best,
Matthias

On Thu, May 30, 2024, at 10:55 CDT, Najwa Alshehri <najwaa...@gmail.com> wrote:

> Dear all,
>
>
> Thank you for your always support and help.
>
>
> I am encountering an issue while solving an eigenvalue problem related to
> the computation of the discrete inf-sup constant for a saddle point
> problem. Specifically, I am solving the following system:
>
>
> Bt A^-1 B eigenvector = eigenvalue M eigenvector,
>
>
> where Bt A^-1 B is clearly symmetric and at least semi-positive definite,
> and M is a symmetric positive definite matrix.
>
>
> The problem arises when comparing the inf-sup constants obtained using
> Matlab and dealii (solving on the same mesh ):
>
>

> 1. Matlab: Using the "eig" solver, the inf-sup constant decays to zero

> as the mesh is refined.

> 2. deal.II: Using the "ArpackSolver", the inf-sup constant is bounded

[Najwa Alshehri]

Dear Wolfgang:

Thank you for your answer. I have checked the matrices and I printed the case for a 2 by 2 mesh from both Matlab and deal.II. With a simple change of rows and columns, they are exactly the same ( see the files attached). Not sure yet why they result in different results of the eigenvalue problem. There is a very importent detail that I want to point out here. When I solved on the very course mesh (2 by 2) I got the following error, which is not the case when I solved on a 4 by 4 mesh for example.

An error occurred in line <791> of file </usr/local/include/deal.II/lac/arpack_solver.h> in function
    void dealii::ArpackSolver::solve(const MatrixType1&, const MatrixType2&, const INVERSE&, std::vector<std::complex<double> >&, std::vector<ValueType>&, unsigned int) [with VectorType = dealii::Vector<double>; MatrixType1 = dealii::LinearOperator<dealii::Vector<double>, dealii::Vector<double>, dealii::internal::LinearOperatorImplementation::EmptyPayload>; MatrixType2 = dealii::LinearOperator<dealii::Vector<double>, dealii::Vector<double>, dealii::internal::LinearOperatorImplementation::EmptyPayload>; INVERSE = dealii::SparseDirectUMFPACK]
The violated condition was:
    false
Additional information:
    Error with dnaupd, info -2. Check documentation of ARPACK

Dear Mathias:

Thank you for your answer. The matrices are exactly the same in both codes. The different is that I solve for the recipricals in deal.II. To be precise,

let AA= Bt A^inv B 

and let M be a mass matrix

where B, A, M are the same matrices in both Matlab and deal.II.

In Matlab: I solve eig(AA, M) (AA x = \lambda M x)  then I seeks the  square root of the smalles \lambda 

In dealii: I solve using ArpackSolver solving M x = 1/\Lambda AA x) 

so I do first 

M_inv_umfpack.initialize(M);

auto M = linear_operator(M);

auto M_inv = linear_operator(M, M_inv_umfpack);

then, I create a linear operator of the other matrices and the inverse of A and perform AA= Bt * A^inv *  B 

lastly I solve : arpack_solver.solve( M, AA, M_inv_umfpack, eigenvalues,eigenvectors);

So I am seeking the largest eigenvalue and \lambda = 1/largest eigenvalue and the infsup constant is the square root of that.

Could solving for the recipricals is the issue? if yes what should I fill in the third intery of the solving if Iam doing 

arpack_solver.solve( AA, M, ????????, eigenvalues,eigenvectors);

One more time, your questions and help are highly appreciated.

Best,

Najwa

[Najwa Alshehri]

Hello all,

One more time, thanks for your support.  I've managed to identify and resolve the issue I was facing, and I'd like to share the solution here for the benefit of others who might encounter similar challenges.

The problem stemmed from my approach of solving for the reciprocals of the eigenvalues, which I did just because the inverse of the mass matrix M was easier. Specifically,  I was solving M x = 1/\lambda AA x. The question of why this was the problem remains unsolved.

I decided to solve the exact problem directly, namely AA x = \lambda M x.  To achieve this, I computed the inverse of the matrix AA= Bt * A^inv *  B  using a Conjugate Gradient (CG) solver. Subsequently, I solved for the exact eigenvalues, and to my satisfaction, I obtained results that aligned with those obtained from MATLAB.

If anyone has insights into why solving for the reciprocal yields different results, I would greatly appreciate your input on this matter.

Best,

Najwa

[Wolfgang Bangerth]

On 6/1/24 14:49, Najwa Alshehri wrote:
>
> I decided to solve the exact problem directly, namely AA x = \lambda M x. To
> achieve this, I computed the inverse of the matrix AA= Bt * A^inv *  B using a
> Conjugate Gradient (CG) solver. Subsequently, I solved for the exact
> eigenvalues, and to my satisfaction, I obtained results that aligned with
> those obtained from MATLAB.
>
> If anyone has insights into why solving for the reciprocal yields different
> results, I would greatly appreciate your input on this matter.

Najwa,
I'm glad you figured it out because I don't really have a good idea about what
ARPACK does. Typically,

AA x = \lambda M x

is solved in exactly this way because M may be a matrix that is singular. In
many applications (not sure whether yours is one of those), M is a matrix that
can have a large null-space -- for example, if you are computing eigenvalues
of the Stokes operator, in some applications you have that M = [M_u 0 ; 0 0],
that is, it has a large null block on all pressures. This matrix is not
invertible, whereas AA is. By convention, the invertible matrix is kept on the
left side of the eigenvalue equation.

[Najwa Alshehri]

Wolfgang and all,

I have a positive definite matrix M, and its inverse can be quickly and easily found in deal.II using the following lines:

M_inv_umfpack.initialize(M);

auto M = linear_operator(M);

auto M_inv = linear_operator(M, M_inv_umfpack);

However, when solving the direct system, I need to provide the inverse of AA to the Arpack solver. This requires solving the system using the Conjugate Gradient (CG) method a number of times equal to the size of AA. As the mesh is refined, this size increases, and depending on the finite element method (FE) I am using, it could be even larger.

While I have been able to obtain at least the first three results for many FEs I am testing which is enough for the current study, I am seeking a faster solution for a general eigenvalue problem of the same type. If anyone has suggestions for a more efficient approach to finding the inverse of AA, I would greatly appreciate it.

One more time, you have always been great supporters. This group is very useful. Thank you to you Wolfgang, to Mathias for the answer, and to all members and developers of deal.ii.

Best regards,
Najwa

[Luca Heltai]

Dear Najwa, 

You don’t need to create the inverse of A. Just its action, pretty much in the exact same way you did for the mass matrix, but replacing the inverse operator done with umfpack with one done with CG. 

Luca

Il giorno 2 giu 2024, alle ore 09:42, Najwa Alshehri <najwaa...@gmail.com> ha scritto:

Wolfgang and all,

--
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.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/e448210c-be8d-4460-8b41-27a4425f7accn%40googlegroups.com.

[Luca Heltai]

To be more specific, you can use inverse_operator, to which you pass A, a solver, and a preconditioner, and it returns an operator the applies the inverse of A to a vector. Just as your mass inverse.

Luca

Il giorno 2 giu 2024, alle ore 11:25, Luca Heltai <luca....@gmail.com> ha scritto:

Dear Najwa, 

[Najwa Alshehri]

Dear Luca,

Thank you for your answer, that works perfectly. 

Best,

Najwa

[Wolfgang Bangerth]

Najwa:
on a separate note: Computing the inf-sup constant is an interesting problem.
It would be nice to have a program that shows how to do that available as the
base for others' experiments. Could I interest you in submitting your code,
once you've got it working, to the code gallery?
https://dealii.org/code-gallery.html

Best
W.

[Najwa Alshehri]

Wolfgang:

Thank you for bringing this to my attention. I would definitely love to share the code in the code gallery for the benefit of others. In this case, I will create a code for a simpler problem, such as the Stokes problem or mixed Laplace problem, to illustrate the concept of the inf-sup condition and how to compute the constant. As soon as it is ready, I will submit it there.

Best regards,

Najwa

[Wolfgang Bangerth]

On 6/3/24 02:03, Najwa Alshehri wrote:
> Thank you for bringing this to my attention. I would definitely love to share
> the code in the code gallery for the benefit of others. In this case, I will
> create a code for a simpler problem, such as the Stokes problem or mixed
> Laplace problem, to illustrate the concept of the inf-sup condition and how to
> compute the constant. As soon as it is ready, I will submit it there.

Yes, I think that would be fabulous! Let us know if you have questions or need
help!