Calling IterativeSolvers.cg! for N iterations in one call produces a different per-iteration residual sequence than performing the same total number of iterations via multiple cg! calls.
A single call with 6 iterations produces a different result from two calls of thee iterations, or six calls of one iteration. Continuing with the same statevars or reseting them, does not seem to have an effect.
Is this normal behavior in the cg! implementation?
It looks like this issue has been raised five years ago
and the solution outlined:
Resuming iterative solvers requires not only the last approximate solution xₙ, but also other state variables. For example in CG, these state variables are he last residual vector rₙ and the last search direction vector pₙ. Luckily, these state variables are already stored in CGStateVariables, but currently the code is not written to resume the CG algorithm from the state stored in the passed statevars::CGStateVariables.
CG is a Krylov subspace method which tries to update the solution in orthogonal directions, which means that behind the scenes, a orthogonalization process runs which delivers orthogonal residuals and A-ortohogonal search directions. A restart without the the context accumulated so far loses the orthogonal basis built so far.
CG for symmetric positive definite A has the advantage that the new orthogonal direction only depends on the residual and the previous direction which occurs as “state vector”, and does not need to carry the residual of all the previous iterations (unlike GMRES). If the state vector is lost, the orthogonality information is discarded.
With restart(s), you thus you should observe a slowdown of the convergence.
You are absolutely right. In CG, the new orthogonal search direction only depends on the residual rₙ and the previous search direction pₙ. In the cg! implementation this information is maintained in the statevars::CGStateVariables.
The genius in our group figured out how to use the iterator constructed and used by cg! to achieve resumption of iterations from the latest xₙ, rₙ and pₙ. But the question still remains. Why not expose this functionality directly?
If the CG algorithm uses short recurrences to compute iterands as @j-fu and the theory suggests, it should suffice to provide the lastest iterand x_n at restart. Thus only limited advantage to keep track of r_n and p_n.
The situation does change in case of long term recurrence methods like GMRES and friends or in case of solving linear systems with same matrix and various right-hand hands.
BTW there are two more iterative solver packages: Krylov.jl and KrylovKit.jl. I is possible that one of the provides an API which allows to keep the context with restarts.
For some solvers like CG, it should not be too hard to have an option to resume the solver from the current workspace but it may be way harder for some solvers where we update the solution with a lag like MINRES-QLP or only at the end like GMRES.
I like to compare what you request as a reverse-communication API, where you have the full control on the solve and can do one by one the iteration and check the status / information from the solver.
To be able to do what you request in Krylov.jl or IterativeSolvers.jl, we need to store a few additional scalar variables in the workspace, they are the key missing part to resume the solver from a current state.
For example, \alpha and \beta in CG.
In practice we only store vectors in the workspace and keep the local scalar variable on the stack.
Please open an issue / discussion an Krylov.jl if you are interested by this feature.
