fastQR (version 1.1.4)

qr_fast: Fast full QR decomposition

Description

qr_fast provides the fast QR factorization of the matrix \(X\in\mathbb{R}^{n\times p}\) with \(n>p\). The full QR factorization of the matrix \(X\) returns the matrices \(Q\in\mathbb{R}^{n\times p}\) and the upper triangular matrix \(R\in\mathbb{R}^{p\times p}\) such that \(X=QR\). See Golub and Van Loan (2013) for further details on the method.

Usage

qr_fast(X, tol = NULL, pivot = NULL)

Value

A named list containing

qr

a matrix with the same dimensions as \(X\). The upper triangle contains the \(R\) of the decomposition and the lower triangle contains information on the \(Q\) of the decomposition (stored in compact form).

qraux

a vector of length ncol(x) which contains additional information on \(Q\).

rank

the rank of \(X\) as computed by the decomposition.

pivot

information on the pivoting strategy used during the decomposition.

pivoted

a boolean variable returning one if the pivoting has been performed and zero otherwise.

Arguments

X

a \(n\times p\) matrix with \(n>p\).

tol

the tolerance for detecting linear dependencies in the columns of \(X\).

pivot

a logical value indicating whether to pivot the columns of \(X\). Defaults to FALSE, meaning no pivoting is performed.

Details

The QR decomposition plays an important role in many statistical techniques. In particular it can be used to solve the equation \(Ax=b\) for given matrix \(A\in\mathbb{R}^{n\times p}\) and vectors \(x\in\mathbb{R}^{p}\) and \(b\in\mathbb{R}^{n}\). It is useful for computing regression coefficients and in applying the Newton-Raphson algorithm.

References

golub_van_loan.2013fastQR

bjorck.2015fastQR

bjorck.2024fastQR

bernardi_etal.2024fastQR

Examples

## generate sample data
set.seed(1234)
n <- 12
p <- 5
X <- matrix(rnorm(n * p), n, p)

## get the full QR decomposition with pivot
qr_res <- fastQR::qr_fast(X = X,
                          tol = sqrt(.Machine$double.eps),
                          pivot = TRUE)

## reconstruct the reduced Q and R matrices
## reduced Q matrix
Q1 <- qr_Q(qr = qr_res$qr, tau = qr_res$qraux,
           rank = qr_res$rank, complete = FALSE)
Q1

## check the Q matrix (orthogonality)
max(abs(crossprod(Q1)-diag(1, p)))

## complete Q matrix
Q2 <- qr_Q(qr = qr_res$qr, tau = qr_res$qraux,
           rank = NULL, complete = TRUE)
Q2

## check the Q matrix (orthogonality)
max(abs(crossprod(Q2)-diag(1, n)))

## reduced R matrix
R1 <- qr_R(qr = qr_res$qr,
           rank = NULL,
           complete = FALSE)

## check that X^TX = R^TR
## get the permutation matrix
P <- qr_pivot2perm(pivot = qr_res$pivot)
max(abs(crossprod(R1 %*% P) - crossprod(X)))
max(abs(crossprod(R1) - crossprod(X %*% t(P))))

## complete R matrix
R2 <- qr_R(qr = qr_res$qr,
           rank = NULL,
           complete = TRUE)

## check that X^TX = R^TR
## get the permutation matrix
P <- qr_pivot2perm(pivot = qr_res$pivot)
max(abs(crossprod(R2 %*% P) - crossprod(X)))
max(abs(crossprod(R2) - crossprod(X %*% t(P))))

## check that X = Q %*% R
max(abs(Q2 %*% R2 %*% P - X))
max(abs(Q1 %*% R1 %*% P - X))

## create data: n > p
set.seed(1234)
n <- 120
p <- 75
X <- matrix(rnorm(n * p), n, p)

## get the full QR decomposition with pivot
qr_res <- fastQR::qr_fast(X = X, pivot = FALSE)

## reconstruct the reduced Q and R matrices
## reduced Q matrix
Q1 <- qr_Q(qr = qr_res$qr, tau = qr_res$qraux,
           rank = p,
           complete = FALSE)

## check the Q matrix (orthogonality)
max(abs(crossprod(Q1)-diag(1, p)))

## complete Q matrix
Q2 <- qr_Q(qr = qr_res$qr, tau = qr_res$qraux,
           rank = NULL, complete = TRUE)

## check the Q matrix (orthogonality)
max(abs(crossprod(Q2)-diag(1, n)))

## reduced R matrix
R1 <- qr_R(qr = qr_res$qr,
           rank = NULL,
           complete = FALSE)


## check that X^TX = R^TR
max(abs(crossprod(R1) - crossprod(X)))

## complete R matrix
R2 <- qr_R(qr = qr_res$qr,
           rank = NULL,
           complete = TRUE)

## check that X^TX = R^TR
max(abs(crossprod(R2) - crossprod(X)))

## check that X^TX = R^TR
max(abs(crossprod(R2) - crossprod(X)))
max(abs(crossprod(R2) - crossprod(X)))
max(abs(crossprod(R1) - crossprod(X)))

# check that X = Q %*% R
max(abs(Q2 %*% R2 - X))
max(abs(Q1 %*% R1 - X))