lfe: cgsolve – R documentation

Become an expert in R — Interactive courses, Cheat Sheets, certificates and more!

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.

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

Value

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

References

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

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

lfe

Linear Group Fixed Effects

v2.8-6

Artistic-2.0

Authors

Simen Gaure [aut] (<https://orcid.org/0000-0001-7251-8747>), Grant McDermott [ctb], Karl Dunkle Werner [ctb], Matthieu Stigler [ctb, cre] (<https://orcid.org/0000-0002-6802-4290>), Daniel Lüdecke [ctb]

Initial release

2021-01-04