cgsolve: Solve a symmetric linear system with the conjugate gradient...
In sgaure/lfe: Linear Group Fixed Effects

Description Usage Arguments Details Value References See Also Examples

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.

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

`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.

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|| < ε ||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.

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

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

kaczmarz

  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))