lfe (version 2.5-1998)

cgsolve: Solve a symmetric linear system with the conjugate gradient method

Description

cgsolve uses a conjugate gradient algorithm to solve the linear system $A x = b$ where $A$ is a symmetric matrix. cgsolve is used internally in lfe in the routines fevcov and bccorr, but has been made public because it might be useful for other purposes as well.

Usage

cgsolve(A, b, eps = 0.001, init = NULL, symmtest = FALSE, name = "")

Arguments

A

matrix, Matrix or function.

b

vector or matrix of columns vectors.

eps

numeric. Tolerance.

init

numeric. Initial guess.

symmtest

logical. Should the matrix be tested for symmetry?

name

character. Arbitrary name used in progress reports.

Value

A solution $x$ of the linear system $A x = b$ is returned.

Details

The argument A can be a symmetric matrix or a symmetric sparse matrix inheriting from "Matrix" of the package Matrix. It can also be a function which performs the matrix-vector product. If so, the function must be able to take as input a matrix of column vectors.

If the matrix A is rank deficient, some solution is returned. If there is no solution, a vector is returned which may or may not be close to a solution. If symmtest is FALSE, no check is performed that A is symmetric. If not symmetric, cgsolve is likely to raise an error about divergence.

The tolerance eps is a relative tolerance, i.e. $||x - x_0|| < \epsilon ||x_0||$ where $x_0$ is the true solution and $x$ is the solution returned by cgsolve. Use a negative eps for absolute tolerance. The termination criterion for cgsolve is the one from Kaasschieter (1988), Algorithm 3.

Preconditioning is currently not supported.

If A is a function, the test for symmetry is performed by drawing two random vectors x,y, and testing whether $|(Ax, y) - (x, Ay)| < 10^{-6} sqrt((||Ax||^2 + ||Ay||^2)/N)$, where $N$ is the vector length. Thus, the test is neither deterministic nor perfect.

References

Kaasschieter, E. (1988) A practical termination criterion for the conjugate gradient method, BIT Numerical Mathematics, 28(2):308-322. http://link.springer.com/article/10.1007%2FBF01934094

See Also

kaczmarz

Examples

  N <- 100000
# create some factors
  f1 <- factor(sample(34000,N,replace=TRUE))
  f2 <- factor(sample(25000,N,replace=TRUE))
# a matrix of dummies, which probably is rank deficient
  B <- makeDmatrix(list(f1,f2))
  dim(B)
# create a right hand side
  b <- as.matrix(B %*% rnorm(ncol(B)))
# solve B' B x = B' b
  sol <- cgsolve(crossprod(B), crossprod(B, b), eps=-1e-2)
  #verify solution
  sqrt(sum((B %*% sol - b)^2))