Package 'gmresls'

Solve a Least Squares System with a Preconditioner.

Description

Solve a least squares system Ax~=b (dim(A)=(m,n) with m >= n) with a preconditioner B: BAx=Bb (dim(B)=(n,m)). Implemented method uses GMRES(k) with callback functions, i.e. no explicit A or B are required. GMRES can be restarted after k iterations.

Usage

gmresls(f_resid, f_BAx, x0 = NULL, k = 0, maxit = 0, tol = 1e-07, ...)

Arguments

f_resid	A function f_resid(x, ...) calculating B(b-Ax) for a given x. If x is of length 0 (e.g. NULL), it must be considered as 0.
f_BAx	A function f_BAx(x, ...) calculating matrix-matrix-vector product BAx for a given x.
x0	A vector or NULL (which means 0), initial approximation for Ax=b
k	An integer, parameter for restarting GMRES. Value 0 (default) means no restart, i.e. at most length(x) basis vectors will be constructed and used.
maxit	A maximal iteration number. Here, itereation number continues to increment even after a possible GMRES restart. Default (0) means length(x).
tol	A tolerance for solution x, estimated as ||B(Ax-b)||/||Bb||, default 1.e-7
...	Parameters passed through to f_BAx and f_resid

Details

Implemented method is equivalent to a classical GMRES(k) method with restart after constructing k basis vectors and applied to a square system BAx=Bb. Dense matrices constructed and stored by this method are of size (length(x), k) and (k+1, k) where k is GMRES current basis vector number. If maxit > k, GMRES will be restarted after each k iterations Particularity of this implementation that matrices A and B have no to be stored explicitly. User provides just callback function mimicking their multiplication by adequate vectors. In case of non convergence after maxit iterations, attr(x) will contain a field 'warning' with the message which will be also issued with warning() If the operator BA is not of full rank, iterations will be stopped before reaching convergence or maxit. A warning will be emitted in this case.

Value

The solution x, having the structure of f_resid(x,...).

Examples

# prepare a 4x3 toy least squares (LS) problem Ax=b
A=rbind(diag(1:3)+matrix(1, 3,3), rep(1, 3))
xsol=1:3
b=A%*%xsol+rnorm(4, 0., 0.1) # add some noise as it is often the case in LS
f_resid=function(x,...)
   with(list(...), if (length(x) == 0L) crossprod(A, b) else crossprod(A, b-A%*%x))
f_BAx=function(x,...)
   with(list(...), crossprod(A, A%*%x))
x=gmresls(f_resid, f_BAx, A=A, b=b)
stopifnot(all.equal(c(x), c(qr.solve(A,b))))

---