Visitor pattern for convergence history

[haampie]

Currently there are two or three options like plot, log and verbose for the solvers. In fact these are all different ways to show progress of the algorithm. I think it is a good design choice to separate output about progress completely from the algorithms themselves.

One way to go is to use the visitor pattern: let the user pass in a callback, and each iteration call it with the current convergence history results. Then you could do stuff like this:

p = plot(yscale = :log10)
plot_on_the_go = t::ConvergenceHistory -> push!(p, t.iterations, last(t[:resnorm]))
cg(A, b, visitor = plot_on_the_go)

Is there a drawback to this?

[haampie]

The problem with the above is that the algorithm can't know what data the user needs. Maybe a ConvergenceHistory object isn't enough:

In GMRES you might want to access the Hessenberg matrix to compute Ritz values at the end of a "macro" iteration; you could use approximate eigenvectors to apply deflation techniques in a second problem you need to solve. Or you might want to implement a multishift version of GMRES where you solve (A - tI)x = b for multiple t's, reusing the search subspace from the Ax=b problem. In BiCGStab(l) you might want to compare how the computed residual actually differs from the real residual (so you need access to intermediate solutions x). Et cetera.

All this should be possible without changing the algorithm itself or copy / pasting it into a new one.

Now @jiahao started with the iterator idea for CG, and I currently think that is the best way to make it possible. So basically you have some low-level helper function to construct an iterable:

iterable = cg_iterable(A, b, ...)

And then let the algorithm run as

for (iteration, residual) = enumerate(iterable)
    println("Iteration ", iteration, " |r| = ", residual)
end

The iterable gives the user access to x and r et cetera:

for res = iterable
  surface(reshape(iterable.x, 128, 128))
end

The original API can then just be a thin wrapper around the iterator:

function cg!(x, A, b; kwargs...)
  iterable = cg_iterable(x, A, b; kwargs...)

  for item = iterable end

  iterable.x
end

A last important benefit of having an iterable object, is that it is reusable. For instance for non-linear problems where in each Newton iterations a large linear system needs to be solved (approximately), it is very useful to pre-allocate vectors once for all Newton steps. In the current implementation this is not possible, because the methods themselves do the allocations.

[ChrisRackauckas]

Yes, that is definitely the way to do it. Note that DifferentialEquations.jl went this way with the Integrator interface actually being a (mutable) iterator, and the solve function actually just building the integrator and iterating it. So many beautiful things come out of this. For example, it makes it easier to test and debug algorithms because you can always just easily step-by-step it, even without the debugger (makes great MWEs for difficult in-package bugs). You can use this for algorithm development: between the iterations you can inject new code that basically is the new algorithm to easily try it out without having to modify the actual implementation (great for research papers). Any user can modify things as it goes along for whatever weird algorithm they are interested in.

After that change I never looked back. I would be happy to see IterativeSolvers.jl do the same.

[haampie]

Thanks Chris! Seems to make a lot of sense indeed