sparseQR-class: Sparse QR Factorizations
In Matrix: Sparse and Dense Matrix Classes and Methods
sparseQR-class R Documentation
Sparse QR Factorizations
Description
sparseQR is the class of sparse, row- and column-pivoted
QR factorizations of m \times n (m \ge n)
real matrices, having the general form
P_1 A P_2 = Q R = \begin{bmatrix} Q_1 & Q_2 \end{bmatrix} \begin{bmatrix} R_1 \\ 0 \end{bmatrix} = Q_1 R_1
or (equivalently)
A = P_1' Q R P_2' = P_1' \begin{bmatrix} Q_1 & Q_2 \end{bmatrix} \begin{bmatrix} R_1 \\ 0 \end{bmatrix} P_2' = P_1' Q_1 R_1 P_2'
where
P_1 and P_2 are permutation matrices,
Q = \prod_{j = 1}^{n} H_j
is an m \times m orthogonal matrix
(Q_1 contains the first n column vectors)
equal to the product of n Householder matrices H_j, and
R is an m \times n upper trapezoidal matrix
(R_1 contains the first n row vectors and is
upper triangular).
Details
The method for qr.Q does not return Q but rather the
(also orthogonal) product P_1' Q. This behaviour
is algebraically consistent with the base implementation
(see qr), which can be seen by noting that
qr.default in base does not pivot rows, constraining
P_1 to be an identity matrix. It follows that
qr.Q(qr.default(x)) also returns P_1' Q.
Similarly, the methods for qr.qy and qr.qty multiply
on the left by P_1' Q and Q' P_1
rather than Q and Q'.
It is wrong to expect the values of qr.Q (or qr.R,
qr.qy, qr.qty) computed from “equivalent”
sparse and dense factorizations
(say, qr(x) and qr(as(x, "matrix")) for x
of class dgCMatrix) to compare equal.
The underlying factorization algorithms are quite different,
notably as they employ different pivoting strategies,
and in general the factorization is not unique even for fixed
P_1 and P_2.
On the other hand, the values of qr.X, qr.coef,
qr.fitted, and qr.resid are well-defined, and
in those cases the sparse and dense computations should
compare equal (within some tolerance).
The method for qr.R is a simple wrapper around qrR,
but not back-permuting by default and never giving row names.
It did not support backPermute = TRUE until Matrix
1.6-0, hence code needing the back-permuted result should
call qrR if Matrix >= 1.6-0 is not known.
Slots
Dim, Dimnamesinherited from virtual class
MatrixFactorization.
betaa numeric vector of length Dim[2],
used to construct Householder matrices; see V below.
Van object of class dgCMatrix
with Dim[2] columns. The number of rows nrow(V)
is at least Dim[1] and at most Dim[1]+Dim[2].
V is lower trapezoidal, and its column vectors generate the
Householder matrices H_j that compose the orthogonal
Q factor. Specifically, H_j is constructed as
diag(Dim[1]) - beta[j] * tcrossprod(V[, j]).
Ran object of class dgCMatrix
with nrow(V) rows and Dim[2] columns.
R is the upper trapezoidal R factor.
p, q0-based integer vectors of length
nrow(V) and Dim[2], respectively,
specifying the permutations applied to the rows and columns of
the factorized matrix. q of length 0 is valid and
equivalent to the identity permutation, implying no column pivoting.
Using R syntax, the matrix P_1 A P_2
is precisely A[p+1, q+1]
(A[p+1, ] when q has length 0).
Extends
Class QR, directly.
Class MatrixFactorization, by class
QR, distance 2.
Instantiation
Objects can be generated directly by calls of the form
new("sparseQR", ...), but they are more typically obtained
as the value of qr(x) for x inheriting from
sparseMatrix (often dgCMatrix).
Methods
determinantsignature(x = "sparseQR", logarithm = "logical"):
computes the determinant of the factorized matrix A
or its logarithm.
expand1signature(x = "sparseQR"):
see expand1-methods.
expand2signature(x = "sparseQR"):
see expand2-methods.
qr.Qsignature(qr = "sparseQR"):
returns as a dgeMatrix either
P_1' Q or P_1' Q_1,
depending on optional argument complete. The default
is FALSE, indicating P_1' Q_1.
qr.Rsignature(qr = "sparseQR"):
qrR returns R, R_1,
R P2', or R_1 P2',
depending on optional arguments complete and
backPermute. The default in both cases is FALSE,
indicating R_1, for compatibility with base.
The class of the result in that case is
dtCMatrix. In the other three cases,
it is dgCMatrix.
qr.Xsignature(qr = "sparseQR"):
returns A as a dgeMatrix,
by default. If m > n and optional argument
ncol is greater than n, then the result
is augmented with P_1' Q J, where
J is composed of columns (n+1) through
ncol of the m \times m identity matrix.
qr.coefsignature(qr = "sparseQR", y = .):
returns as a dgeMatrix or vector
the result of multiplying y on the left by
P_2 R_1^{-1} Q_1' P_1.
qr.fittedsignature(qr = "sparseQR", y = .):
returns as a dgeMatrix or vector
the result of multiplying y on the left by
P_1' Q_1 Q_1' P_1.
qr.residsignature(qr = "sparseQR", y = .):
returns as a dgeMatrix or vector
the result of multiplying y on the left by
P_1' Q_2 Q_2' P_1.
qr.qtysignature(qr = "sparseQR", y = .):
returns as a dgeMatrix or vector
the result of multiplying y on the left by
Q' P_1.
qr.qysignature(qr = "sparseQR", y = .):
returns as a dgeMatrix or vector
the result of multiplying y on the left by
P_1' Q.
solvesignature(a = "sparseQR", b = .):
see solve-methods.
References
Davis, T. A. (2006).
Direct methods for sparse linear systems.
Society for Industrial and Applied Mathematics.
\Sexpr[results=rd]{tools:::Rd_expr_doi("10.1137/1.9780898718881")}
Golub, G. H., & Van Loan, C. F. (2013).
Matrix computations (4th ed.).
Johns Hopkins University Press.
\Sexpr[results=rd]{tools:::Rd_expr_doi("10.56021/9781421407944")}
See Also
Class dgCMatrix.
Generic function qr from base,
whose default method qr.default “defines”
the S3 class qr of dense QR factorizations.
qr-methods for methods defined in Matrix.
Generic functions expand1 and expand2.
The many auxiliary functions for QR factorizations:
qr.Q, qr.R, qr.X,
qr.coef, qr.fitted, qr.resid,
qr.qty, qr.qy, and qr.solve.
Examples
showClass("sparseQR")
set.seed(2)
m <- 300L
n <- 60L
A <- rsparsematrix(m, n, 0.05)
## With dimnames, to see that they are propagated :
dimnames(A) <- dn <- list(paste0("r", seq_len(m)),
paste0("c", seq_len(n)))
(qr.A <- qr(A))
str(e.qr.A <- expand2(qr.A, complete = FALSE), max.level = 2L)
str(E.qr.A <- expand2(qr.A, complete = TRUE), max.level = 2L)
t(sapply(e.qr.A, dim))
t(sapply(E.qr.A, dim))
## Horribly inefficient, but instructive :
slowQ <- function(V, beta) {
d <- dim(V)
Q <- diag(d[1L])
if(d[2L] > 0L) {
for(j in d[2L]:1L) {
cat(j, "\n", sep = "")
Q <- Q - (beta[j] * tcrossprod(V[, j])) %*% Q
}
}
Q
}
ae1 <- function(a, b, ...) all.equal(as(a, "matrix"), as(b, "matrix"), ...)
ae2 <- function(a, b, ...) ae1(unname(a), unname(b), ...)
## A ~ P1' Q R P2' ~ P1' Q1 R1 P2' in floating point
stopifnot(exprs = {
identical(names(e.qr.A), c("P1.", "Q1", "R1", "P2."))
identical(names(E.qr.A), c("P1.", "Q" , "R" , "P2."))
identical(e.qr.A[["P1."]],
new("pMatrix", Dim = c(m, m), Dimnames = c(dn[1L], list(NULL)),
margin = 1L, perm = invertPerm(qr.A@p, 0L, 1L)))
identical(e.qr.A[["P2."]],
new("pMatrix", Dim = c(n, n), Dimnames = c(list(NULL), dn[2L]),
margin = 2L, perm = invertPerm(qr.A@q, 0L, 1L)))
identical(e.qr.A[["R1"]], triu(E.qr.A[["R"]][seq_len(n), ]))
identical(e.qr.A[["Q1"]], E.qr.A[["Q"]][, seq_len(n)] )
identical(E.qr.A[["R"]], qr.A@R)
## ae1(E.qr.A[["Q"]], slowQ(qr.A@V, qr.A@beta))
ae1(crossprod(E.qr.A[["Q"]]), diag(m))
ae1(A, with(e.qr.A, P1. %*% Q1 %*% R1 %*% P2.))
ae1(A, with(E.qr.A, P1. %*% Q %*% R %*% P2.))
ae2(A.perm <- A[qr.A@p + 1L, qr.A@q + 1L], with(e.qr.A, Q1 %*% R1))
ae2(A.perm , with(E.qr.A, Q %*% R ))
})
## More identities
b <- rnorm(m)
stopifnot(exprs = {
ae1(qrX <- qr.X (qr.A ), A)
ae2(qrQ <- qr.Q (qr.A ), with(e.qr.A, P1. %*% Q1))
ae2( qr.R (qr.A ), with(e.qr.A, R1))
ae2(qrc <- qr.coef (qr.A, b), with(e.qr.A, solve(R1 %*% P2., t(qrQ)) %*% b))
ae2(qrf <- qr.fitted(qr.A, b), with(e.qr.A, tcrossprod(qrQ) %*% b))
ae2(qrr <- qr.resid (qr.A, b), b - qrf)
ae2(qrq <- qr.qy (qr.A, b), with(E.qr.A, P1. %*% Q %*% b))
ae2(qr.qty(qr.A, qrq), b)
})
## Sparse and dense computations should agree here
qr.Am <- qr(as(A, "matrix")) # <=> qr.default(A)
stopifnot(exprs = {
ae2(qrX, qr.X (qr.Am ))
ae2(qrc, qr.coef (qr.Am, b))
ae2(qrf, qr.fitted(qr.Am, b))
ae2(qrr, qr.resid (qr.Am, b))
})
Matrix documentation built on Aug. 13, 2024, 3:01 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
| sparseQR-class | R Documentation |
Sparse QR Factorizations
Description
sparseQR is the class of sparse, row- and column-pivoted
QR factorizations of m \times n (m \ge n)
real matrices, having the general form
P_1 A P_2 = Q R = \begin{bmatrix} Q_1 & Q_2 \end{bmatrix} \begin{bmatrix} R_1 \\ 0 \end{bmatrix} = Q_1 R_1
or (equivalently)
A = P_1' Q R P_2' = P_1' \begin{bmatrix} Q_1 & Q_2 \end{bmatrix} \begin{bmatrix} R_1 \\ 0 \end{bmatrix} P_2' = P_1' Q_1 R_1 P_2'
where
P_1 and P_2 are permutation matrices,
Q = \prod_{j = 1}^{n} H_j
is an m \times m orthogonal matrix
(Q_1 contains the first n column vectors)
equal to the product of n Householder matrices H_j, and
R is an m \times n upper trapezoidal matrix
(R_1 contains the first n row vectors and is
upper triangular).
Details
The method for qr.Q does not return Q but rather the
(also orthogonal) product P_1' Q. This behaviour
is algebraically consistent with the base implementation
(see qr), which can be seen by noting that
qr.default in base does not pivot rows, constraining
P_1 to be an identity matrix. It follows that
qr.Q(qr.default(x)) also returns P_1' Q.
Similarly, the methods for qr.qy and qr.qty multiply
on the left by P_1' Q and Q' P_1
rather than Q and Q'.
It is wrong to expect the values of qr.Q (or qr.R,
qr.qy, qr.qty) computed from “equivalent”
sparse and dense factorizations
(say, qr(x) and qr(as(x, "matrix")) for x
of class dgCMatrix) to compare equal.
The underlying factorization algorithms are quite different,
notably as they employ different pivoting strategies,
and in general the factorization is not unique even for fixed
P_1 and P_2.
On the other hand, the values of qr.X, qr.coef,
qr.fitted, and qr.resid are well-defined, and
in those cases the sparse and dense computations should
compare equal (within some tolerance).
The method for qr.R is a simple wrapper around qrR,
but not back-permuting by default and never giving row names.
It did not support backPermute = TRUE until Matrix
1.6-0, hence code needing the back-permuted result should
call qrR if Matrix >= 1.6-0 is not known.
Slots
Dim,Dimnamesinherited from virtual class
MatrixFactorization.betaa numeric vector of length
Dim[2], used to construct Householder matrices; seeVbelow.Van object of class
dgCMatrixwithDim[2]columns. The number of rowsnrow(V)is at leastDim[1]and at mostDim[1]+Dim[2].Vis lower trapezoidal, and its column vectors generate the Householder matricesH_jthat compose the orthogonalQfactor. Specifically,H_jis constructed asdiag(Dim[1]) - beta[j] * tcrossprod(V[, j]).Ran object of class
dgCMatrixwithnrow(V)rows andDim[2]columns.Ris the upper trapezoidalRfactor.p,q0-based integer vectors of length
nrow(V)andDim[2], respectively, specifying the permutations applied to the rows and columns of the factorized matrix.qof length 0 is valid and equivalent to the identity permutation, implying no column pivoting. Using R syntax, the matrixP_1 A P_2is preciselyA[p+1, q+1](A[p+1, ]whenqhas length 0).
Extends
Class QR, directly.
Class MatrixFactorization, by class
QR, distance 2.
Instantiation
Objects can be generated directly by calls of the form
new("sparseQR", ...), but they are more typically obtained
as the value of qr(x) for x inheriting from
sparseMatrix (often dgCMatrix).
Methods
determinantsignature(x = "sparseQR", logarithm = "logical"): computes the determinant of the factorized matrixAor its logarithm.expand1signature(x = "sparseQR"): seeexpand1-methods.expand2signature(x = "sparseQR"): seeexpand2-methods.qr.Qsignature(qr = "sparseQR"): returns as adgeMatrixeitherP_1' QorP_1' Q_1, depending on optional argumentcomplete. The default isFALSE, indicatingP_1' Q_1.qr.Rsignature(qr = "sparseQR"):qrRreturnsR,R_1,R P2', orR_1 P2', depending on optional argumentscompleteandbackPermute. The default in both cases isFALSE, indicatingR_1, for compatibility with base. The class of the result in that case isdtCMatrix. In the other three cases, it isdgCMatrix.qr.Xsignature(qr = "sparseQR"): returnsAas adgeMatrix, by default. Ifm > nand optional argumentncolis greater thann, then the result is augmented withP_1' Q J, whereJis composed of columns(n+1)throughncolof them \times midentity matrix.qr.coefsignature(qr = "sparseQR", y = .): returns as adgeMatrixor vector the result of multiplyingyon the left byP_2 R_1^{-1} Q_1' P_1.qr.fittedsignature(qr = "sparseQR", y = .): returns as adgeMatrixor vector the result of multiplyingyon the left byP_1' Q_1 Q_1' P_1.qr.residsignature(qr = "sparseQR", y = .): returns as adgeMatrixor vector the result of multiplyingyon the left byP_1' Q_2 Q_2' P_1.qr.qtysignature(qr = "sparseQR", y = .): returns as adgeMatrixor vector the result of multiplyingyon the left byQ' P_1.qr.qysignature(qr = "sparseQR", y = .): returns as adgeMatrixor vector the result of multiplyingyon the left byP_1' Q.solvesignature(a = "sparseQR", b = .): seesolve-methods.
References
Davis, T. A. (2006). Direct methods for sparse linear systems. Society for Industrial and Applied Mathematics. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1137/1.9780898718881")}
Golub, G. H., & Van Loan, C. F. (2013). Matrix computations (4th ed.). Johns Hopkins University Press. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.56021/9781421407944")}
See Also
Class dgCMatrix.
Generic function qr from base,
whose default method qr.default “defines”
the S3 class qr of dense QR factorizations.
qr-methods for methods defined in Matrix.
Generic functions expand1 and expand2.
The many auxiliary functions for QR factorizations:
qr.Q, qr.R, qr.X,
qr.coef, qr.fitted, qr.resid,
qr.qty, qr.qy, and qr.solve.
Examples
showClass("sparseQR")
set.seed(2)
m <- 300L
n <- 60L
A <- rsparsematrix(m, n, 0.05)
## With dimnames, to see that they are propagated :
dimnames(A) <- dn <- list(paste0("r", seq_len(m)),
paste0("c", seq_len(n)))
(qr.A <- qr(A))
str(e.qr.A <- expand2(qr.A, complete = FALSE), max.level = 2L)
str(E.qr.A <- expand2(qr.A, complete = TRUE), max.level = 2L)
t(sapply(e.qr.A, dim))
t(sapply(E.qr.A, dim))
## Horribly inefficient, but instructive :
slowQ <- function(V, beta) {
d <- dim(V)
Q <- diag(d[1L])
if(d[2L] > 0L) {
for(j in d[2L]:1L) {
cat(j, "\n", sep = "")
Q <- Q - (beta[j] * tcrossprod(V[, j])) %*% Q
}
}
Q
}
ae1 <- function(a, b, ...) all.equal(as(a, "matrix"), as(b, "matrix"), ...)
ae2 <- function(a, b, ...) ae1(unname(a), unname(b), ...)
## A ~ P1' Q R P2' ~ P1' Q1 R1 P2' in floating point
stopifnot(exprs = {
identical(names(e.qr.A), c("P1.", "Q1", "R1", "P2."))
identical(names(E.qr.A), c("P1.", "Q" , "R" , "P2."))
identical(e.qr.A[["P1."]],
new("pMatrix", Dim = c(m, m), Dimnames = c(dn[1L], list(NULL)),
margin = 1L, perm = invertPerm(qr.A@p, 0L, 1L)))
identical(e.qr.A[["P2."]],
new("pMatrix", Dim = c(n, n), Dimnames = c(list(NULL), dn[2L]),
margin = 2L, perm = invertPerm(qr.A@q, 0L, 1L)))
identical(e.qr.A[["R1"]], triu(E.qr.A[["R"]][seq_len(n), ]))
identical(e.qr.A[["Q1"]], E.qr.A[["Q"]][, seq_len(n)] )
identical(E.qr.A[["R"]], qr.A@R)
## ae1(E.qr.A[["Q"]], slowQ(qr.A@V, qr.A@beta))
ae1(crossprod(E.qr.A[["Q"]]), diag(m))
ae1(A, with(e.qr.A, P1. %*% Q1 %*% R1 %*% P2.))
ae1(A, with(E.qr.A, P1. %*% Q %*% R %*% P2.))
ae2(A.perm <- A[qr.A@p + 1L, qr.A@q + 1L], with(e.qr.A, Q1 %*% R1))
ae2(A.perm , with(E.qr.A, Q %*% R ))
})
## More identities
b <- rnorm(m)
stopifnot(exprs = {
ae1(qrX <- qr.X (qr.A ), A)
ae2(qrQ <- qr.Q (qr.A ), with(e.qr.A, P1. %*% Q1))
ae2( qr.R (qr.A ), with(e.qr.A, R1))
ae2(qrc <- qr.coef (qr.A, b), with(e.qr.A, solve(R1 %*% P2., t(qrQ)) %*% b))
ae2(qrf <- qr.fitted(qr.A, b), with(e.qr.A, tcrossprod(qrQ) %*% b))
ae2(qrr <- qr.resid (qr.A, b), b - qrf)
ae2(qrq <- qr.qy (qr.A, b), with(E.qr.A, P1. %*% Q %*% b))
ae2(qr.qty(qr.A, qrq), b)
})
## Sparse and dense computations should agree here
qr.Am <- qr(as(A, "matrix")) # <=> qr.default(A)
stopifnot(exprs = {
ae2(qrX, qr.X (qr.Am ))
ae2(qrc, qr.coef (qr.Am, b))
ae2(qrf, qr.fitted(qr.Am, b))
ae2(qrr, qr.resid (qr.Am, b))
})
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.