Add chebyshev Iteration

[yhmtsai]

It adds the Chebshev iteration in https://en.wikipedia.org/wiki/Chebyshev_iteration
The second-order richardson uses a similar formula but the scalars are constant in all iteration. Chebyshev Iteration update the scalar from the previous one (the scalar are the same if upper/lower eigval does not change)

[MarcelKoch]

Do you intend adding an eigenvalue estimation? I think that would be very helpful, because most times users don't have that available. I think PETSc is also doing that.
I briefly searched for that and here are some references that might be interesting for that:

• https://doi.org/10.1137/0907057
• https://doi.org/10.1137/0613050
• https://doi.org/10.1007/BF01397475

I think this is the GMRES version used by PETSC to compute the estimate: https://petsc.org/release/docs/manualpages/KSP/KSPAGMRES/

[yhmtsai]

@MarcelKoch thanks for providing these reference.

I do not think I will put them in this pull request.

From https://doi.org/10.1137/0907057, they introduce the algorithm
adaptive chebyshev iteration
Perform m-step GMRES and use the information from the Hessenberg matrix to get the estimation of eigenvalue (I may be wrong),
and then perform https://link.springer.com/article/10.1007/BF01389971 to get the parameter for chebchev
Repeat these two steps until converge
This will be more like a standalone solver not a smoother for multigrid because it will require several steps on GMRES, which is too expensive IMO.
We can introduce by with_adaptive(GMRES LinOpFactory)

From https://petsc.org/release/docs/manualpages/KSP/KSPCHEBYSHEV/#kspchebyshev,
More precisely, https://petsc.org/release/docs/manualpages/KSP/KSPChebyshevEstEigSet/
They use Lanczo or GMRES to get the maximum estimation to decide the parameter by following
min-eig-boundary = a * min-eigest + b * max-eigest = 0.1 max-eigest (default), and
max-eig-boundary = c * min-eigest + d * max-eigest = 1.1 max-eigest (default)
(because min_eigest is usually inaccurate)
By this, user can use the information by GMRES and generate min/max eig boundary for Chebyshev, so we do not need additional parameter?
It only does once in generation, so this is more usable for Multigrid

For GMRES, we need to get the Hessenberg out and compute eigenvalue on it.
GMRES eigenvalue and Lanczo are related to #135

[sonarcloud]

Kudos, SonarCloud Quality Gate passed!   

0 Bugs
0 Vulnerabilities
0 Security Hotspots
8 Code Smells

44.3% Coverage
7.1% Duplication

[MarcelKoch]

Besides the comments left below, I want to mention the num_keep and num_generated mechanic. That needs some explaination, because right now I don't understand what that is for. It seems like some sort of restart mechanic, but I might be wrong there. Also, exposing this to the user seems confusing to me.

[MarcelKoch]

It is not clear what this parameter is for.

[MarcelKoch]

Perhaps this means, construct Chebyshev polynomials until degree num_keep, and after that just do IR with these polynomial?

[yhmtsai]

no, num_keep is to keep the generated scalar in the storage.
alpha and beta are changed by iteration, but they are fixed by the given bound.
It's used for fixed iteration runs because we do not need to refill these scalars for different apply

[MarcelKoch]

But shouldn't we then just keep all scalars?

[yhmtsai]

when using it as normal solver, we may not have the maximum iteration information.
allocating one big dense matrix and then moving them to another one when full may not be a good approach.
we can also add the ability of workspace such that it can handle the std::vector then it should be more flexible for uncertain size.

[yhmtsai]

I can not assume that there's an iteration criterion. It can contain the residual norm or time criterion.
If we assume that, we should provide the standalone iteration parameter.
Users do not need to know the implementation details.
The statement is to keep how many scalars for chebyshev to avoid the refill overhead.
I can introduce this usage later when we have the use case.
I think it should help the situation when kernel launch overhead is noticeable.

[MarcelKoch]

You are right, this is very useful to negate the overhead. Still, I think we can assume that there is an iteration criterion somewhere in the combined one and use that. If not, we would just throw.

[MarcelKoch]

You could set the default value to unspecified, which you define before. When the solver is generated, you can check if the parameter is unspecified and then try to extract an iteration criterion from the criteria. If that doesn't work, you throw with a message to either pass an iteration criterion or set this parameter to something else.
Or alternatively, just replace the criteria parameter with a iterations parameter and move the class into the preconditioner namespace. I think that is also a fine option, since that would be the main use case anyway.

[yhmtsai]

I have moved it to based on the given iteration from stopping criterion.
It is only increased after creating object.
I also think whether staying a fixed number is enough or not.

[MarcelKoch]

Thanks, I think it should be fine the way it is now.

[MarcelKoch]

Suggested change

	* residual. It has another term for the difference of solution. Moreover, this
	* residual. The solution is then updated using the Chebyshev polynomials. Moreover, this

Is that what the sentence is trying to say? Anyway, I think it should be mentioned somewhere that this uses these polynomials.

[yhmtsai]

the solution x_i is also based on the $x_{i-1} - x_{i-2}$

[MarcelKoch]

TBH I don't see the relevance of that. Has that any effect for the user?

[yhmtsai]

no, I only try to explain the algorithm's difference from IR. IR uses $x += \alpha M^{-1}r$, but chebyshev uses $x += \alpha_i M^{-1} r + \beta_i (x_{i-1} - x_{i-2})$ $(x_{i-1} - x_{i-2})$ is the additional term against IR

[MarcelKoch]

I still think this has to be rephrased

[yhmtsai]

I try to rephrase it again. Could you take a look?

[vasilisge0]

I am fine with merging this PR when corrections in parts of the documentation that are pointed out (comments) take place. I could understand what num_keep variable was used for but perhaps a more descriptive name would be better.

[sonarcloud]

Kudos, SonarCloud Quality Gate passed!   

0 Bugs
0 Vulnerabilities
0 Security Hotspots
21 Code Smells

42.6% Coverage
4.4% Duplication

The version of Java (11.0.3) you have used to run this analysis is deprecated and we will stop accepting it soon. Please update to at least Java 17.
Read more here

[MarcelKoch]

LGTM, but I would like to resolve some of my earlier issues. Other than that only minor nits.

[MarcelKoch]

I still think this has to be rephrased

[sonarcloud]

Kudos, SonarCloud Quality Gate passed!   

0 Bugs
0 Vulnerabilities
1 Security Hotspot
27 Code Smells

59.4% Coverage
5.2% Duplication

The version of Java (11.0.3) you have used to run this analysis is deprecated and we will stop accepting it soon. Please update to at least Java 17.
Read more here

[sonarcloud]

Quality Gate passed

The SonarCloud Quality Gate passed, but some issues were introduced.

27 New issues
1 Security Hotspot
59.4% Coverage on New Code
5.7% Duplication on New Code

See analysis details on SonarCloud

[sonarcloud]

Quality Gate passed

Issues
29 New issues
0 Accepted issues

Measures
1 Security Hotspot
55.6% Coverage on New Code
6.4% Duplication on New Code

See analysis details on SonarCloud

[upsj]

There are a few questions around data recycling and mutability I would like to discuss before merging this

[upsj]

I would rearrange this slightly (the fact that we need to know the spectrum is important, so we should mention it early), and try to be consistent with capitalization of iteration.

Suggested change

	* Chebyshev iteration is an iterative method that can solving nonsymeetric
	* problems. Chebyshev Iterations avoids the inner products for computation
	* which may be the bottleneck for distributed system. Chebyshev Iteration is
	* developed based on the Chebyshev polynomials of the first kind. Moreover,
	* this method requires knowledge about the spectrum of the preconditioned
	* matrix. This implementation follows the algorithm in "Templates for the
	* Solution of Linear Systems: Building Blocks for Iterative Methods, 2nd
	* Edition".
	* Chebyshev iteration is an iterative method for solving nonsymmeetric
	* problems based on some knowledge of the spectrum of the (preconditioned) system matrix. It avoids the computation of inner products,
	* which may be a performance bottleneck for distributed system. Chebyshev iteration is
	* developed based on Chebyshev polynomials of the first kind.
	* This implementation follows the algorithm in "Templates for the
	* Solution of Linear Systems: Building Blocks for Iterative Methods, 2nd
	* Edition".

[upsj]

I'm a bit wary of mutable members beyond well-controlled cases like the workspace. They might make future work on thread-safety more complicated. Not that I have a good alternative for it, but I want to bring it up.

[yhmtsai]

Isn't wokrspace handled in similar way?

[upsj]

We don't have any mutability in our LinOps right now beyond vector caches and workspaces, most of the mutable uses are in loggers etc. So the mutability is constrained into a single class. That's not the case here.

In general, if I write code, I expect a const object to not change internally, and the only reasons we don't follow this requirement fully right now seems to be the need for cache data structures to avoid reallocation, and the fact that loggers are const and can't change themselves without mutable.

[upsj]

this entire file can use deferred factories

[upsj]

deferred?

[upsj]

The naming solver is a bit of a stretch here - maybe even more so than IR. It makes sense for us that solvers are iterative algorithms that use a preconditioner, but from the user perspective, it's confusing. Maybe we need a concept in-between, like smoother, but that could also include classical preconditioners like Gauß-Seidel

[upsj]

Is this actually mixed precision? Doesn't look like it to me on a first glance

[upsj]

Is this mixed precision?

[upsj]

Here solver == this, so I would suggest removing the capture

[upsj]

I think this should be a separate parameter, not derived from the stopping criterion. The stopping criteria are not really intended to be queried this way.

[yhmtsai]

It was. but @MarcelKoch suggested that getting it from stopping criterion to reduce the parameter, I think.

[MarcelKoch]

If we put it as a separate parameter, then we should not allow for a stopping criteria parameter. This would be fine with me, but I think it causes other issues, as this can't be considered a solver anymore.

[upsj]

The smoother doesn't have much purpose with norm-based stopping criteria, because that would negate the benefit of avoiding reductions, right? And can it run correctly when it uses more iterations than the number of stored scalar values? If yes, then we could have two separate parameters iterations and subspace_dimension (or similar), if no, we could have just a single parameter iterations, and internally create an Iteration stopping criterion

[upsj]

Am I understanding correctly that this us potentially re-using alpha and beta values from previous iterations? If so, that seems like it is using the workspace it was not originally intended for - it's intended for avoiding reallocation, not recycling values from previous iterations. That also relates to the mutable members I would prefer to avoid.

[yhmtsai]

@upsj To avoid the mutable and reusing value from workspace, I think moving them into generatio to pre-fill all alpha and beta before apply. I think it will also give better performance because we do not need to do some tiny fill before the limitation during apply.
Also, refill it when the foci or criterion change (or specfic factor for that).

[upsj]

If those are just static coefficients independent of the matrix, that sounds very reasonable. I tend to agree with what I believe @MarcelKoch said, which is that we shouldn't consider changeable stopping criteria here. We also currently don't seem to have a way to propagate information from other parts of the solver hierarchy to compute appropriate focus values.