Proposal: package with iterative solvers for linear systems

[Vladimír Chalupecký]

Over the last couple of months I've been occasionally working on a package for iterative solution of linear systems. We already have several direct solvers for dense systems based on matrix factorizations. With this addition we would have solvers for large sparse systems, for systems where the matrix doesn’t need to be explicitly formed, or for systems where the solution is not required to be accurate to machine precision (for example in inexact Newton methods for the optimize package).

There are certainly some corners that need to be polished but I hope that it has reached a presentable state. The code currently lives at:

https://github.com/vladimir-ch/iterative

If it is accepted here, I will then submit a PR to the gonum repository to discuss the code itself.

First of all, the package name. I settled on “iterative” very early on and so far I didn’t come up with a better name. It reads quite nicely: the main interface is iterative.Method, the main function is iterative.LinearSolve.

The overall design is similar to the optimize package. At the low level there are various solver types that implement the Method interface. It uses the reverse communication mechanism in order to avoid depending on how the system matrix is represented. Operations like (transposed) matrix times vector are delegated to the caller. In this way we can provide the iterative algorithms without having any implementation for sparse matrices available. At a higher level there is a convenience function LinearSolve which implements the main iteration loop similarly to optimize.Local. The convenience comes at the “cost” of concretising how the (transposed) matrix times vector (and other) operations are done.

The implementation of the various methods loosely follows “The Templates” http://www.netlib.org/templates/double/. Actually, I used the version in SciPy where the bugs present in the netlib version are fixed. For the start I would have been happy with providing just CG and GMRES but to better understand what operations are needed, I did also BiCG and BiCGSTAB. I plan to focus on Krylov methods, so methods like SOR don’t fit into the current design.

Regarding the API, I’d like to get some opinions mainly about how to represent the matrix operations (and the preconditioner). At the moment there is a struct MatrixOps with two fields MatVec func(dst, x []float64) and MatTransVec func(dst, x []float64). The other obvious option is to make it an interface which is probably better. I have other smaller doubts (e.g., details of the reverse communication protocol: who should compute what, what operations are necessary), and there are many TODOs sprinkled over the code but I hope these will be raised and addressed during the review.

[Brendan Tracey]

I think this is a great idea to have in gonum. These kinds of linear solvers are widely used in large-scale engineering. 

I don't think Iterative is a good name. What the package does is solve large-scale linear programs. I would think that "linear" or "linsolve" would be better. The call "linsolve.Iterative" is very invocative of what is happening. 

I think it should be an interface, and there should be an easy "dense" wrapper. I would also guess in this case the "Result" should allow for supplying the destination slice. It's pretty common that these are the inner loop, and I would guess the extra allocation could be relatively costly. 

Are we sure that []float64 is the right input and output type, as opposed to a Vector or some more general interface? I would imagine it's relatively common that a lot of the RHS variables are 0, and so a benefit may be gained from allowing a sparse input. I suspect that in general the output will be dense? One advantage of a more abstract vector is that in HPC environments it may not be possible to hold "x" or "b" in memory. Having the vector portion be more abstract could allow that use case.

[Vladimír Chalupecký]

I don't think Iterative is a good name. What the package does is solve large-scale linear programs. I would think that "linear" or "linsolve" would be better. The call "linsolve.Iterative" is very invocative of what is happening. 

I had considered those two. For me "linear" is to broad and so does not express what's in the package. "linsolve" is better but it could suggest that it has all Gonum's linear solvers which would be only partly true. In my plan, being iterative is the distinctive feature. Also, "iterative.LinearSolve" is just as invocative as "linsolve.Iterative". I'm not arguing for "iterative", I'm just comparing the options.

Perhaps what else we would like to have in the package in the future could decide the name. If it was a sparse direct solver, then "linsolve" is probably better (linsolve.Iterative, linsolve.Direct could be two main entry points).  If we wanted to have there iterative methods for computing eigenvalues, then probably "iterative" is better (iterative.LinearSolve, iterative.Eigen, ...).

I think it should be an interface, and there should be an easy "dense" wrapper.

Then, what about

type Matrix interface {

    MulVec(dst, x []float64, trans bool)

}

? Or we could keep it a struct, call it System, and add the right-hand side to it:

type System struct {

    Matrix func(dst, x []float64, trans bool)

    RHS []float64  // or B

}

I will make a mental note to add an easy dense wrapper later.

 

I would also guess in this case the "Result" should allow for supplying the destination slice. It's pretty common that these are the inner loop, and I would guess the extra allocation could be relatively costly. 

I agree. For this I could move the initial guess X0 from Settings to the signature of LinearSolve. If nil, allocate. If not nil, use as an initial guess and overwrite with the solution. In both cases return also in Result. Would this be ok?

 

Are we sure that []float64 is the right input and output type, as opposed to a Vector or some more general interface?

I would imagine it's relatively common that a lot of the RHS variables are 0, and so a benefit may be gained from allowing a sparse input. I suspect that in general the output will be dense?

Yes, the result will be dense. The right-hand side may often have many zeros but I'm not sure if the extra generality and complexity wins anything here. The rhs is used to compute the residual b-Ax once at the start and then only if the method doesn't update it iteratively (as is the case in GMRES). Compared to the matrix-vector multiply the gains in computing sparse vector minus dense vector instead of difference of two dense vectors is probably negligible.

 

One advantage of a more abstract vector is that in HPC environments it may not be possible to hold "x" or "b" in memory. Having the vector portion be more abstract could allow that use case.

Switching to a distributed computation would be more involved than just using a general interface for the vectors. The proposal considers only serial solvers, parallel solvers would need a completely separate set of interfaces and implementations (although I would be excited to have them one day). Also note that Method only cares about the type of the output solution vector but it is oblivious to what the rhs is.

[Vladimír Chalupecký]

I have moved the development to the exp repository at https://github.com/gonum/exp

The code is at the moment in the add-linsolve branch.

In the end I decided to rename the package to linsolve as suggested by Brandan. It's actually a nice package name, and describes better than "iterative" what the package does.

I also decided to introduce the System struct that holds the matrix times vector function and the right-hand side. Together with providing the initial vector and storing the solution in-place, the main entry function has signature:

func Iterative(sys System, x []float64, method Method, settings Settings) (*Result, error)

which I think is not bad.

If this is roughly acceptable, is it ok if I create a PR to merge this into exp/master? We could discuss at code level there.

[Brendan Tracey]

If you create a PR I’ll start working with you on it. I remember having a bunch of comments when I worked through it last.

--
You received this message because you are subscribed to the Google Groups "gonum-dev" group.
To unsubscribe from this group and stop receiving emails from it, send an email to gonum-dev+...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.