v0.2.0

- changed interface, now we use f_resid() and f_BAx()
- unit test and manual are rewritte in consequence
- added attribute 'converged'

Signed-off-by: Serguei Sokol <sokol@insa-toulouse.fr>

DESCRIPTION

 Package: gmresls
 Title: Solve Least Squares with GMRES(k)
-Version: 0.1.1
+Version: 0.2
 Authors@R: 
     person("Serguei", "Sokol", , "sokol@insa-toulouse.fr", role = c("aut", "cre"),
            comment = c(ORCID = "0000-0002-5674-3327"))

NEWS.md

-## v0.1.1 2024-??-??
+## v1.0.0 2024-??-??
 
  - initial R implemetation

R/gmresls.R

 #' Solve a Least Squares System with a Preconditioner.
 #'
-#' Solve a least squares system Ax=b (dim(A)=(m,n) with m >= n) with a
+#' 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.
 #'
-#' @param f_Ax A function calculating matrix-vector product Ax for a given x
-#' @param f_Br A function calculating matrix-vector product Br for a given r
-#' @param b A vector of right hand side
+#' @param 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.
+#' @param f_BAx A function f_BAx(x, ...) calculating matrix-matrix-vector product BAx for a given x.
 #' @param x0 A vector or NULL (which means 0), initial approximation for Ax=b
 #' @param 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.
-#' @param maxit A maximal iteration number. Default (0) means length(x).
-#' @param tol A tolerance for solution, estimated as ||B(Ax-b)||/||Bb||, default 1.e-7
-#' @param ... Parameters passed through to f_Ax and f_Br
-#' @returns The solution x, having the structure of Bb.
+#' @param maxit A maximal iteration number. Here, itereation number continues to increment even after a possible GMRES restart. Default (0) means length(x).
+#' @param tol A tolerance for solution x, estimated as ||B(Ax-b)||/||Bb||, default 1.e-7
+#' @param ... Parameters passed through to f_BAx and f_resid
+#' @returns The solution x, having the structure of f_resid(x,...).
 #' @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
@@ -22,33 +21,32 @@
 #' 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.
 #'
 #' @examples
-#' # prepare a 4x3 toy problem Ax=b
+#' # 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
-#' f_Ax=function(x,...) list(...)$A%*%x
-#' f_Br=function(r,...) t(list(...)$A)%*%r
-#' x=gmresls(f_Ax, f_Br, b, A=A)
+#' 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))))
 #' @export
-gmresls=function(f_Ax, f_Br, b, x0=NULL, k=0, maxit=0, tol=1.e-7, ...) {
-  f_ba=function(x, ...) f_Br(f_Ax(x, ...), ...)
-  bb=f_Br(b, ...)
-  nbb=sqrt(sum(bb*bb)) # useful for convergence detection
-  nx=length(bb)
-  if (maxit == 0L)
-    maxit=nx
+gmresls=function(f_resid, f_BAx, x0=NULL, k=0, maxit=0, tol=1.e-7, ...) {
   it=0L
   warn=""
   repeat { # restart till convergence or maxit reached
     it=it+1L
-    r0=bb
-    nr0=nbb
-    if (!is.null(x0)) {
-      r0=r0-f_ba(x0, ...)
-      nr0=sqrt(sum(r0*r0))
+    r0=f_resid(x0, ...)
+    nr0=sqrt(sum(r0*r0))
+    if (it == 1L) {
+      nx=length(r0)
+      if (maxit == 0L)
+        maxit=nx
+      nbb=if (is.null(x0)) nr0 else sqrt(crossprod(c(f_resid(NULL, ...))))
     }
     
-    v=as.matrix(r0/nr0)
+    v=as.matrix(c(r0))*(1./nr0)
     h=matrix(0., 1L, 0L) # Heisenberg matrix, will be growing
     re1=nr0
     if (k == 0)
@@ -57,7 +55,7 @@ gmresls=function(f_Ax, f_Br, b, x0=NULL, k=0, maxit=0, tol=1.e-7, ...) {
     for (i in seq_len(k)) {
       if (converged)
         return(if (is.null(x0)) double(nx) else x0)
-      bavi=f_ba(v[,i], ...)
+      bavi=f_BAx(v[,i], ...)
       h=cbind(h, crossprod(v, bavi))
       v=cbind(v, bavi-v%*%h[,i])
       h=rbind(h, 0.)
@@ -88,8 +86,7 @@ gmresls=function(f_Ax, f_Br, b, x0=NULL, k=0, maxit=0, tol=1.e-7, ...) {
     warn=paste0(warn, "Maximal iteration number (", maxit, ") is reached without convergence")
   if (nchar(warn))
     warning(attr(x, "warning") <- warn)
-
-  attr(x, "resid")=nri
+  attr(x, "norm_res")=nri
   attr(x, "iteration")=it
   attr(x, "converged")=converged
   x

README.md

@@ -21,12 +21,12 @@ devtools::install_git(git@forgemia.inra.fr:mathscell/gmresls.git)
 This is a basic example which shows you how to solve a common problem:
 
 ``` r
-library(gmresls)
 # prepare a 4x3 toy problem Ax=b
 A=rbind(diag(1:3)+matrix(1, 3,3), rep(1, 3))
 xsol=1:3
-b=A%*%xsol
-f_Ax=function(x,...) list(...)$A%*%x
-f_Br=function(r,...) t(list(...)$A)%*%r
-x=gmresls(f_Ax, f_Br, b, A=A)
+b=A%*%xsol+rnorm(4, 0., 0.1)
+f_resid=function(x,...) with(list(...), if (length(x) == 0) 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))))
 ```

inst/unitTests/runit.gmresls.R

@@ -2,10 +2,9 @@
 # prepare a 4x3 toy problem Ax=b
 A=rbind(diag(1:3)+matrix(1, 3,3), rep(1, 3))
 xsol=1:3
-b=A%*%xsol
-f_Ax=function(x,...) list(...)$A%*%x
-f_Br=function(r,...) t(list(...)$A)%*%r
-x=gmresls(f_Ax, f_Br, b, A=A)
-test.doc <- function() {
-  checkEqualsNumeric(x, xsol)
-}
+b=A%*%xsol+rnorm(4, 0., 0.1)
+f_resid=function(x,...) with(list(...), if (length(x) == 0) 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)
+test.doc <- function()
+  checkEqualsNumeric(x, qr.solve(A, b))

man/gmresls.Rd

@@ -4,30 +4,28 @@
 \alias{gmresls}
 \title{Solve a Least Squares System with a Preconditioner.}
 \usage{
-gmresls(f_Ax, f_Br, b, x0 = NULL, k = 0, maxit = 0, tol = 1e-07, ...)
+gmresls(f_resid, f_BAx, x0 = NULL, k = 0, maxit = 0, tol = 1e-07, ...)
 }
 \arguments{
-\item{f_Ax}{A function calculating matrix-vector product Ax for a given x}
+\item{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.}
 
-\item{f_Br}{A function calculating matrix-vector product Br for a given r}
-
-\item{b}{A vector of right hand side}
+\item{f_BAx}{A function f_BAx(x, ...) calculating matrix-matrix-vector product BAx for a given x.}
 
 \item{x0}{A vector or NULL (which means 0), initial approximation for Ax=b}
 
 \item{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.}
 
-\item{maxit}{A maximal iteration number. Default (0) means length(x).}
+\item{maxit}{A maximal iteration number. Here, itereation number continues to increment even after a possible GMRES restart. Default (0) means length(x).}
 
-\item{tol}{A tolerance for solution, estimated as ||B(Ax-b)||/||Bb||, default 1.e-7}
+\item{tol}{A tolerance for solution x, estimated as ||B(Ax-b)||/||Bb||, default 1.e-7}
 
-\item{...}{Parameters passed through to f_Ax and f_Br}
+\item{...}{Parameters passed through to f_BAx and f_resid}
 }
 \value{
-The solution x, having the structure of Bb.
+The solution x, having the structure of f_resid(x,...).
 }
 \description{
-Solve a least squares system Ax=b (dim(A)=(m,n) with m >= n) with a
+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.
 }
@@ -40,11 +38,14 @@ In case of non convergence after maxit iterations, attr(x) will contain a field
 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.
 }
 \examples{
-# prepare a 4x3 toy problem Ax=b
+# 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
-f_Ax=function(x,...) list(...)$A\%*\%x
-f_Br=function(r,...) t(list(...)$A)\%*\%r
-x=gmresls(f_Ax, f_Br, b, A=A)
+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))))
 }