feat (linalg_iterative) : add gmres solver to the iterative methods library

[jalvesz]

This PR proposes adding the generalized minimal residual method (GMRES) to the linalg_iterative module.

It was drafted with the help of CoPilot and received the help of @AlexanderGSC to add the Modified Gram Schmidt scheme for increased robustness.

• Public core kernel implementation
• Standard interface for dense and CSR matrix types
• Unit test
• Examples
• Documentation

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 68.20%. Comparing base (a2ef8aa) to head (c4734a0).
⚠️ Report is 1 commits behind head on master.

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1186      +/-   ##
==========================================
- Coverage   68.46%   68.20%   -0.26%     
==========================================
  Files         412       19     -393     
  Lines       13816     2378   -11438     
  Branches     1559        0    -1559     
==========================================
- Hits         9459     1622    -7837     
+ Misses       4357      756    -3601     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[AlexanderGSC]

Thanks @jalvesz!

If more tests or documentation need to be added, or if the maintainers think it’s necessary to make any changes, I’d be happy to help.

[jvdp1]

Thank you for this addition. Here are some suggestions.

[jvdp1]

is it on purpose the same as stdlib_size_wksp_cg?

[jalvesz]

It is not on purpose the as stdlib_size_wksp_cg. The workspace for the GMRES algorithm is bit more tricky than for the other gradient descent methods as not only are there reference buffer vectors with the size of the problem but also there are arrays depending on the kdim (number of internal iterations per restart cycle).

The current workspace design :
r(:) → residual (1 vector)
w(:) → Arnoldi work vector (1 vector)
v(:,1:kdim+1) → Krylov basis (kdim+1 vectors)
z(:,1:kdim) → preconditioned basis (kdim vectors)

Actually the last vector z could be reduced to z(:) as the preconditionner is applied per each j inner iteration.

I'll review this and try to propose a more lean version in terms of internal storage size. This won't change the fact that the workspace depends not only on n (problem size) but also on the selected number of internal iterations kdim.

Maybe @AlexanderGSC you have some suggestions here?

[AlexanderGSC]

Hi @jvdp1, @jalvesz!

@jalvesz is right; the size of the workspace depends on the size of the restart (kdim), and there isn’t much that can be done about it. Storing only the last vector from the preconditioner saves a lot of memory – well spotted!

I can think of a few ideas for keeping the interface unchanged, each with its own pros and cons:

• move the allocation of v inside the kernel:
  pros: keeps the current interface unchanged.
  cons: I don’t know if it’s allowed to reserve memory inside the kernel. At the moment there are variables such as the Hessemberg matrix, the rotation vectors, etc. that are allocated inside the kernel, but they are small, depending only on kdim, which by default is 30, meaning very little memory is used.

• Explicit workspace handling:
  Pros: no memory allocation within the kernel.
  Cons: we would need to document that the workspace size is 3+(kdim+1). And we would need to add a check in the code to ensure this is the case. Furthermore, in the interface, the value would not reflect the workspace size.

• Implicit kdim handling: the workspace size is checked, and the kdim size is deduced implicitly. If a value of 33 is set in the interface, gmres would run with a default kdim=30.
  Pros: allows running ‘by default’ or with ‘fine-tuning’ without allocating memory in the kernel.
  Cons: the way the workspace is managed becomes a little odd and less straightforward. Restart is a parameter that every version of gmres has.

In my opinion, if it isn’t a problem, I would move the v memory allocation into the kernel. If that is a problem, then we could explore other options.

[jalvesz]

I gave another round to the proposal. I added a compact option to allow switching between a memory efficient or a speed oriented implementation versions. Let me know your thoughts on this.

[AlexanderGSC]

Hello @jalvesz,

What you’re suggesting is much better; it’s a definitive solution that would allow for much greater flexibility and provide a clean interface for any method requiring a variable number of vectors.

Is the Hilbert test causing problems again in single precision? If so, I think the best thing is to remove it, or if necessary, I’ll modify it so that the test is skipped entirely in single precision. Adding the interface is much more important. There’ll be time to fix the test later on.

Great work!

[jalvesz]

Yes, the Hilbert test in single precision is causing problems with two of the intel jobs, don't know what is best, if just remove it or try with a jacobi preconditioner ?

[AlexanderGSC]

I don’t think using a preconditioner will make any difference.

The fact that it fails with Intel compilers suggests to me that they may be using aggressive instruction fusion optimisations that cause some numerical drift. The condition number of the Hilbert's matrix must be of the order of $10^{12}$, which is impossible for SP.

The simplest solution: reduce the size of the Hilbert matrix to something manageable for SP. With $n=4$ it should work.

If you agree, I can push the change to update the test to $n=4$ directly so we can see if the Intel CI jobs finally turn green.

[jalvesz]

Indeed! I'm thinking that another option would be to check the condition number estimate and pick a size such that the condition is below some critical threshold ? Using
$$\kappa \approx \frac{(1+\sqrt{2})^{4*N}}{\sqrt{N}}$$
a limit size could be precomputed at compile time

... or simply just do the test for kind >sp with fixed size.

[jalvesz]

@loiseaujc thanks for your input! Indeed you are right on both comments, I added the apply_givens_rotation following your suggestion, including a small comment to LightKrylov. I also changed the triangular solve to rely on trtrs, we have to decide how to properly exit in case of a singular triangular matrix. For the moment I left an exit path if (info/=0) exit

@AlexanderGSC I pushed a change bypassing the Hilbert test for sp following your suggestion.

[loiseaujc]

As far as I know, the triangular matrix should never be singular. The cheap residual check norm_sq = g(j+1) * g(j+1) is precisely here for that. As soon as g(j+1) is small enough (indicating an almost invariant subspace of A), the triangular matrix will be of size j x j and invertible. It comes for the on-the-fly QR factorization of the upper Hessenberg matrix.

[jalvesz]

Thanks @loiseaujc for the clarification. So in that case the exit path is useless. Would you say it is ok to remove it or what would you advise here?

[loiseaujc]

Well, I guess we could have the standard linalg_state_type stuff just to be on the safe side but, as far as I know, we should never encounter an error.

[jalvesz]

Right, though introducing that would imply changing a bit the API. If needed we could include that in a future PR. So I guess we can leave it as is for the moment.

[AlexanderGSC]

Hello @jalvesz @loiseaujc ,

@loiseaujc is absolutely right: the cheap residual check (norm_sq < tolsq) should trigger a breakdown exit ($h_{j+1,j} \approx 0$) well before we ever feed a singular system to trtrs.

What I'm not sure will work from a practical standpoint in floating-point computation is using quadratic tolerances (tolsq) and a hard-coded epsilon test within the Arnoldi loop. For example, if $rtol=10^{-6}$ is set in single precision, for an epsilon that will be on the order of $10^{-7}$, it drives the stopping threshold tolsq down to $10^{-12}$ (assuming a normalized system). In single precision, $10^{-12}$ is indistinguishable from pure round-off noise.

Then, the cheap quadratic check (norm_sq < tolsq) might never be mathematically satisfied. The algorithm will keep iterating, the Krylov basis will lose orthogonality, and the Hessenberg matrix will become highly ill-conditioned, and we end up calling trtrs on a nearly-singular matrix where the diagonal elements $H(j,j)$ have collapsed to near-zero. This can lead to division-by-zero, extreme scaling issues, or NaNs.

Keeping the if (info /= 0) exit path after trtrs is actually a very cheap and robust safety net to prevent runtime crashes, or we could avoid the quadratic problem keeping the tolerances and the residual estimate on a linear scale.

These are minor considerations; for now, we can leave it at that and keep refining it step by step. What do you think?

[jalvesz]

Thanks for the feedback @AlexanderGSC so I think we agree and the PR is ready!
Any comments @loiseaujc @jvdp1 ?

[jvdp1]

Nice addition to stdlib. Thank you! I only have some minor suggestions/questions. Feel free to ignore them.