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

In lfe: Linear Group Fixed Effects

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

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

lfe documentation built on Feb. 16, 2023, 7:32 p.m.