lm.fit: Fitter for Linear Models

Description

Fits a real linear model via QR with a "limited pivoting strategy", as in R's DQRDC2 (fortran).

Usage

## S3 method for class 'ddmatrix,ddmatrix':
lm.fit(x, y, tol = 1e-07, singular.ok = TRUE)

Arguments

x, y

numeric distributed matrices

tol

tolerance for numerical rank estimation in QR decomposition.

singular.ok

logical. If FALSE then a singular model (rank-deficient x) produces an error.

Value

Returns a list of values similar to R's lm.fit(). Namely, the list contains: coefficients - distributed matrix, residuals - distributed matrix, effects - distributed matrix, rank - global numeric, fitted.values - distributed matrix, assign - NULL if lm.fit() is called directly, qr - list, same as return from qr(), df.residual - global numeric. The return values are, respectively: (1) a solution x to the linear least squares problem, (2) the difference in the numerical fit A %*% x and the observed b, (3) t(Q) %*% b, where Q is the orthogonal matrix from a QR-decomposition of A, (4) the numerical column rank of A, (5) the numerical fit A %*% x, (6) NULL if lm.fit() is directly called, (7) a list containing the return of QR decomposition performed by a modified PDGEQPF, (8) degrees of freedom of residuals, i.e. n minus the column rank of A.

Details

Solves the linear least squares problem, which is to find an x (possibly non-uniquely) such that || Ax - b ||^2 is minimized, where A is a given n-by-p model matrix, b is a "right hand side" n-by-1 vector (multiple right hand sides can be solved at once, but the solutions are independent, i.e. not simultaneous), and "||" is the l2 norm. Uses level 3 PBLAS and ScaLAPACK routines (modified PDGELS) to get a linear least squares solution, using the 'limited pivoting strategy' from R's DQRDC2 (unsed in DQRLS) routine as a way of dealing with (possibly) rank deficient model matrices. A model matrix with many dependent columns will likely experience poor performance, especially at scale, due to all the data swapping that must occur to handle rank deficiency.

Examples

Run this code

# Save code in a file "demo.r" and run with 2 processors by
# > mpiexec -np 2 Rscript demo.r

library(pbdDMAT, quiet = TRUE)
init.grid()

# don't do this in production code
x <- matrix(rnorm(9), 3)
y <- matrix(rnorm(3))

dx <- as.ddmatrix(x)
dy <- as.ddmatrix(y)

fit <- lm.fit(x=dx, y=dy)

print(fit)

finalize()