From e9ac86a1711d137a17fdcff084bd1cbbd2fc549b Mon Sep 17 00:00:00 2001
From: Serguei Sokol <sokol@insa-toulouse.fr>
Date: Mon, 14 Oct 2024 16:08:11 +0200
Subject: [PATCH] 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 | 2 +-
NEWS.md | 2 +-
R/gmresls.R | 55 ++++++++++++++++------------------
README.md | 10 +++----
inst/unitTests/runit.gmresls.R | 13 ++++----
man/gmresls.Rd | 31 +++++++++----------
6 files changed, 55 insertions(+), 58 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
index 0835ef2..89c0f41 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
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"))
diff --git a/NEWS.md b/NEWS.md
index a88c804..4ec817b 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,3 @@
-## v0.1.1 2024-??-??
+## v1.0.0 2024-??-??
- initial R implemetation
diff --git a/R/gmresls.R b/R/gmresls.R
index 5daa3ef..5429c86 100644
--- a/R/gmresls.R
+++ b/R/gmresls.R
@@ -1,18 +1,17 @@
#' 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
diff --git a/README.md b/README.md
index 01e694a..e166583 100644
--- a/README.md
+++ b/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))))
```
diff --git a/inst/unitTests/runit.gmresls.R b/inst/unitTests/runit.gmresls.R
index 0289ff1..cc5fbeb 100644
--- a/inst/unitTests/runit.gmresls.R
+++ 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))
diff --git a/man/gmresls.Rd b/man/gmresls.Rd
index f61b9c2..5076bda 100644
--- a/man/gmresls.Rd
+++ 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))))
}
--
GitLab