Add an implicit residual norm criterion.

[pratikvn]

This PR adds an implicit residual norm convergence criterion, which checks the reduction in the implicit residual calculated in a solver. Some libraries (example: MFEM) use an implicit residual as their main criterion. This criterion class will allow apples to apples comparison of our solvers to theirs.

The idea is that for solvers such as CG or extensions of CG (BiCG, BiCGSTAB, CGS, FCG...) we can use scalars which are representative of the residual norm of the solution. For example, in CG, the scalar rho is always computed in each iteration and closely follows the error in the solution.

This allows us to save computing the residual norm from the residual vector, which is what we do in most cases to check for convergence in the ResidualNormReduction stopping criterion.

Of course, for some solvers, there is no implicit residual norm available, and for those solvers, this stopping criterion will throw a NOT_SUPPORTED exception.

Edited

Updated modifications:

1. The ResidualNorm class now can perform all three types of criterion checks: absolute, relative and reduction.
2. The previous classes are marked with the [[deprecated("message")]] attribute. Some byproducts of this are:
   i. Whenever this class is used a warning is thrown at compile time.
   ii. As we need to still test the older class functionality, we will always throw the warning during compilation.
   A deprecated note is added to the deprecated classes.
3. As we use ResidualNormReduction kernels everywhere, these are now replaced by ResidualNorm. The default behaviour of ResidualNorm should be the same as ResidualNormReduction, so the change in function name should not affect where ResidualNormReduction was previously called.
4. ResidualNorm and ImplicitResidualNorm use the same class structure and functions except that they call different check_impl kernels. In the future we might want to create a ResidualNormBase class and derive both from that. But I think currently that breaks interface, so I more or less duplicated the code.

Unrelated modifications:

1. fill with value functions were added for Array and Dense. I think these are really helpful to just fill the entire values arrays with a single value. This just uses the fill_array kernel that we previously had, so no new kernels are added.
2. In IDR, we were computing the residual norm in each iteration, but were not updating the stopping criterion with this information, which is now fixed.

[codecov]

Codecov Report

Merging #702 (2f74910) into develop (086cbab) will increase coverage by 0.09%.
The diff coverage is 97.21%.

@@             Coverage Diff             @@
##           develop     #702      +/-   ##
===========================================
+ Coverage    92.74%   92.83%   +0.09%     
===========================================
  Files          354      355       +1     
  Lines        25612    26046     +434     
===========================================
+ Hits         23753    24179     +426     
- Misses        1859     1867       +8     

Impacted Files	Coverage Δ
core/device_hooks/common_kernels.inc.cpp	0.00% <0.00%> (ø)
include/ginkgo/core/base/array.hpp	89.56% <ø> (ø)
include/ginkgo/core/base/types.hpp	92.59% <ø> (ø)
include/ginkgo/core/log/logger.hpp	87.50% <0.00%> (-4.61%)	⬇️
include/ginkgo/core/matrix/dense.hpp	98.13% <ø> (ø)
include/ginkgo/core/preconditioner/ic.hpp	76.34% <0.00%> (ø)
omp/components/fill_array.cpp	100.00% <ø> (ø)
reference/components/fill_array.cpp	100.00% <ø> (ø)
include/ginkgo/core/stop/residual_norm.hpp	90.12% <90.62%> (+1.98%)	⬆️
reference/test/stop/residual_norm_kernels.cpp	97.63% <97.48%> (+0.34%)	⬆️
... and 46 more

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 086cbab...2f74910. Read the comment docs.

[hartwiganzt]

I think having this is very useful! However, I am unsure whether "implicit" is the best name, or whether "recurrence residual" is better.

[pratikvn]

We always pass in the recurrent residual (the explicitly computed residual vector) to the stopping criterion. This is what we check in the relative residual norm reduction stopping criterion. Even though this is also a recurrent residual, I think "implicit residual" might be better suited because the residual norm is implicitly computed, so for example in case of preconditioned CG, the implicit residual norm would be r'*P*r, while the recurrent residual vector would be r and its norm r'*r

[upsj]

Here I would note the distinction between residual, residual norm and implicit residual norm:
(recurrent) residual is exact modulo rounding errors (every solver)
(recurrent) residual norm is exact modulo rounding errors (GMRES)
implicit residual norm is only correlated with the exact residual norm (CG, ...)

[hartwiganzt]

LGTM

[Slaedr]

LGTM, though I do have a few comments. Thanks!

[Slaedr]

Here I would note the distinction between residual, residual norm and implicit residual norm:
(recurrent) residual is exact modulo rounding errors (every solver)
(recurrent) residual norm is exact modulo rounding errors (GMRES)
implicit residual norm is only correlated with the exact residual norm (CG, ...)

I would note that, only in the case of left preconditioned GMRES, the recurrent residual norm (from the Hessenberg least-squares RHS) is the exact preconditioned residual norm. So when convergence is detected based on it, the true residual would typically still not be converged.

[nbeams]

Thanks @pratikvn, it seems this will be useful for the MFEM integration, but also more generally for all Ginkgo users in the future.

[thoasm]

Nice work!

My biggest concern is the usage of a string instead of an enum or boolean to decide between the initial guess or the right-hand side as the reference. Otherwise, I only have minor nits.

[upsj]

LGTM in general, great idea! I would prefer if we could remove the compute_absolute calls, since they allocate memory in each allocation. Another question would be: Do we want to report these implicit residuals to the loggers? Finally, in your stopping criterion, you currently represent relative residuals and residual reduction. What do you think about making this consistent with the regular residual stopping criteria and adding 3 separate stopping criteria based on the baseline value (absolute, relative, residual reduction)?

[pratikvn]

Regarding reporting the implicit residuals to the loggers is a good point and I did have a look at it, but I think that breaks our interface as we cannot add another parameter to the Updater class.

Regarding adding three separate criterion classes, I think that would be a lot of duplication and unnecessary code. But if everyone feels that they would prefer to have separate classes, I can do that.

[pratikvn]

@thoasm thank you for the updates and the fixes. I think your solution is good. I think the fact that we dont have to duplicate check_impl for the deprecated criteria is good.

I also added the strided kernels, which I had missed. I think this should be ready for reviewing again.

[pratikvn]

@tcojean , regarding the logger. Maybe I misunderstand your comment. As I see it, I need to update the criterion_check_started to also take the implicit_residual_norm from the updater, which is passed to the logger in include/.../stop/criterion.hpp line 157. I think I cannot add another GKO_REGISTER_EVENT with criterion_check_started as that sets some variables with the same name so overloading is not possible. Adding something like criterion_check_started_new is possible, I guess without breaking interface, but that again will be in the public interface and will cause naming issues in the future. I would prefer to modify the criterion_check_started and break the interface for good in 2.0.

[thoasm]

LGTM!

[tcojean]

@pratikvn What I mean is that maybe it works if you add the implicit residual at the end and use an = nullptr default case, both for the logger event and the criterion updater? Both the current and new interface would work. Of course, that's maybe not optimal in terms of positioning etc.

Changes addressed.

[upsj]

LGTM mostly, I have two ideas how we could improve and future-proof this approach:

1. logging: You can add another parameter to a logger without breaking interface as follows:

    GKO_LOGGER_REGISTER_EVENT(21, iteration_complete, ..., const LinOp *implicit_tau_sq = nullptr)
protected:                                                           \
    virtual void on_iteration_complete(const LinOp *solver, const size_type &it,
            const LinOp *r, const LinOp *x, const LinOp *tau) const {
        this->on_iteration_complete(solver, it, r, x, tau, nullptr);
    }

2. stopping criterion generation: We could avoid the with_baseline factory parameter call if we had more speaking with_tolerance functions:

Factory& with_reduction_factor(double) { baseline = initial_resnorm; ... }
Factory& with_relative_residual(double) { baseline = rhs_norm; ... }
Factory& with_absolute_residual(double) { baseline = absolute; ... }

and adding another value none to the baseline enum. When we try to generate a residual norm on that enum value, we throw an exception: "Baseline not specified" or something like that.

[upsj]

Since we don't update rho between these two calls, does it make sense to provide the same value twice?

[pratikvn]

good point

[pratikvn]

I need to always report this so that check_impl can capture it. If it is nullptr, then the check_impl for ImplicitResidualNorm will fail.

[upsj]

Do you see any downside to checking for nullptr in check_impl and, in that case, doing nothing? If the stopping status is not changed, the iteration should continue normally, right? Only then we can't stop at half-iterations with BiCGSTAB

[pratikvn]

@upsj , Regarding the logger, I guess @tcojean was also talking about iteration_complete. I completely missed that. I was thinking of criterion_check_x. I will update that.

Regarding the factory parameters, I prefer the with_baseline than separate parameters for the different baselines.

[pratikvn]

Actually, having another look at the iteration_complete and that it is a public function which is overridden in other loggers, adding parameters, actually breaks interface, I think.

[upsj]

@pratikvn That's what the second on_iteration_complete overload is for. Though I actually got it the wrong way round, this current implementation might cause iteration events to get lost on old loggers. The correct solution would be

    GKO_LOGGER_REGISTER_EVENT(21, iteration_complete, ...)
protected:                                                           \
    virtual void on_iteration_complete(const LinOp *solver, const size_type &it,
            const LinOp *r, const LinOp *x, const LinOp *tau, const LinOp *implicit_tau_sq = nullptr) const {
        this->on_iteration_complete(solver, it, r, x, tau);
    }

and overriding the new on_iteration_complete where necessary

[pratikvn]

@upsj , I think the problem is that the call of this->on_iteration_complete(solver, it, r, x, tau) inside the new on_iteration_complete is ambiguous because the compiler cannot resolve between the two when you call it with just non-specific parameters. Also I think any call to on_iteration_complete will be ambiguous unless all the parameters are specified.

[upsj]

@pratikvn Oh yeah, when you drop the default value for the last parameter, the ambiguity goes away. You only need to make sure when overriding the new interface that you override both functions, one with your actual implementation, and the other one with

void on_iteration_complete(const LinOp *solver, const size_type &it,
            const LinOp *r, const LinOp *x = nullptr, const LinOp *tau = nullptr) const {
        this->on_iteration_complete(solver, it, r, x, tau, nullptr);
    }

[pratikvn]

Maybe I am still misunderstanding what you mean here.

Consider a logging call from the gmres apply. This essentially the log function, which calls the on_x based on the Event template parameter, therefore calling on_iteration_complete with the parameters passed into the log function.

The question is as for example in the call to log in gmres, where you pass in 4 parameters in addition to the solver, the on in the logger does not know which on_iteration_complete to call, because both the on_iteration_complete have the same parameters and have only one extra in the second overload.

Unless we specify all the parameters for the iteration_complete in the log call, this will be ambiguous. And adding parameters to all log calls will break interface, because in MFEM and deal.ii, we do derive and use from the previous loggers.

[upsj]

The new overload with the additional parameter must not have any default parameters, then there is no ambiguity between them. If you specify all parameter, the new overload is used, otherwise the old one is called.

[pratikvn]

Okay, I think I finally get you mean. :)

[upsj]

LGTM! I like the unification of the interface, probably something we can finish with 2.0
Only the two overloads of iteration_complete need to be swapped, otherwise old loggers can miss new events.

[sonarcloud]

Kudos, SonarCloud Quality Gate passed!

0 Bugs
0 Vulnerabilities
0 Security Hotspots
27 Code Smells

85.9% Coverage
3.7% Duplication