Title: | Sparse and Dense Matrix Classes and Methods |
---|---|
Description: | A rich hierarchy of sparse and dense matrix classes, including general, symmetric, triangular, and diagonal matrices with numeric, logical, or pattern entries. Efficient methods for operating on such matrices, often wrapping the 'BLAS', 'LAPACK', and 'SuiteSparse' libraries. |
Authors: | Douglas Bates [aut] |
Maintainer: | Martin Maechler <[email protected]> |
License: | GPL (>= 2) | file LICENCE |
Version: | 1.7-0 |
Built: | 2024-06-15 17:28:00 UTC |
Source: | CRAN |
The "abIndex"
class
, short for “Abstract
Index Vector”, is used for dealing with large index vectors more
efficiently, than using integer (or numeric
) vectors of
the kind 2:1000000
or c(0:1e5, 1000:1e6)
.
Note that the current implementation details are subject to change,
and if you consider working with these classes, please contact the
package maintainers (packageDescription("Matrix")$Maintainer
).
Objects can be created by calls of the form new("abIndex", ...)
,
but more easily and typically either by as(x, "abIndex")
where
x
is an integer (valued) vector, or directly by
abIseq()
and combination c(...)
of such.
kind
:a character
string,
one of ("int32", "double", "rleDiff")
, denoting the
internal structure of the abIndex object.
x
:Object of class "numLike"
; is
used (i.e., not of length 0
) only iff the object is not
compressed, i.e., currently exactly when kind != "rleDiff"
.
rleD
:signature(x = "abIndex")
: ...
signature(x = "abIndex", i = "index", j = "ANY", drop = "ANY")
: ...
signature(from = "numeric", to = "abIndex")
: ...
signature(from = "abIndex", to = "numeric")
: ...
signature(from = "abIndex", to = "integer")
: ...
signature(x = "abIndex")
: ...
signature(e1 = "numeric", e2 = "abIndex")
: These
and the following arithmetic and logic operations are
not yet implemented; see Ops
for a
list of these (S4) group methods.
signature(e1 = "abIndex", e2 = "abIndex")
: ...
signature(e1 = "abIndex", e2 = "numeric")
: ...
signature(x = "abIndex")
: ...
("abIndex")
: simple show
method,
building on show(<rleDiff>)
.
("abIndex")
: works analogously to regular vectors.
("abIndex")
: ditto.
This is currently experimental and not yet used for our own code.
Please contact us (packageDescription("Matrix")$Maintainer
),
if you plan to make use of this class.
Partly builds on ideas and code from Jens Oehlschlaegel, as implemented (around 2008, in the GPL'ed part of) package ff.
rle
(base) which is used here;
numeric
showClass("abIndex") ii <- c(-3:40, 20:70) str(ai <- as(ii, "abIndex"))# note ai # -> show() method stopifnot(identical(-3:20, as(abIseq1(-3,20), "vector")))
showClass("abIndex") ii <- c(-3:40, 20:70) str(ai <- as(ii, "abIndex"))# note ai # -> show() method stopifnot(identical(-3:20, as(abIseq1(-3,20), "vector")))
Generation of abstract index vectors, i.e., objects of class
"abIndex"
.
abIseq()
is designed to work entirely like seq
,
but producing "abIndex"
vectors.abIseq1()
is its basic building block, where
abIseq1(n,m)
corresponds to n:m
.
c(x, ...)
will return an "abIndex"
vector, when x
is one.
abIseq1(from = 1, to = 1) abIseq (from = 1, to = 1, by = ((to - from)/(length.out - 1)), length.out = NULL, along.with = NULL) ## S3 method for class 'abIndex' c(...)
abIseq1(from = 1, to = 1) abIseq (from = 1, to = 1, by = ((to - from)/(length.out - 1)), length.out = NULL, along.with = NULL) ## S3 method for class 'abIndex' c(...)
from , to
|
the starting and (maximal) end value of the sequence. |
by |
number: increment of the sequence. |
length.out |
desired length of the sequence. A
non-negative number, which for |
along.with |
take the length from the length of this argument. |
... |
in general an arbitrary number of R objects; here,
when the first is an |
An abstract index vector, i.e., object of class
"abIndex"
.
the class abIndex
documentation;
rep2abI()
for another constructor;
rle
(base).
stopifnot(identical(-3:20, as(abIseq1(-3,20), "vector"))) try( ## (arithmetic) not yet implemented abIseq(1, 50, by = 3) )
stopifnot(identical(-3:20, as(abIseq1(-3,20), "vector"))) try( ## (arithmetic) not yet implemented abIseq(1, 50, by = 3) )
Methods for function all.equal()
(from R package
base) are defined for all Matrix
classes.
\
\
these three methods are
simply using all.equal.numeric
directly and work via
as.vector()
.
There are more methods, notably also for
"sparseVector"
's, see showMethods("all.equal")
.
showMethods("all.equal") (A <- spMatrix(3,3, i= c(1:3,2:1), j=c(3:1,1:2), x = 1:5)) ex <- expand(lu. <- lu(A)) stopifnot( all.equal(as(A[lu.@p + 1L, lu.@q + 1L], "CsparseMatrix"), lu.@L %*% lu.@U), with(ex, all.equal(as(P %*% A %*% t(Q), "CsparseMatrix"), L %*% U)), with(ex, all.equal(as(A, "CsparseMatrix"), t(P) %*% L %*% U %*% Q)))
showMethods("all.equal") (A <- spMatrix(3,3, i= c(1:3,2:1), j=c(3:1,1:2), x = 1:5)) ex <- expand(lu. <- lu(A)) stopifnot( all.equal(as(A[lu.@p + 1L, lu.@q + 1L], "CsparseMatrix"), lu.@L %*% lu.@U), with(ex, all.equal(as(P %*% A %*% t(Q), "CsparseMatrix"), L %*% U)), with(ex, all.equal(as(A, "CsparseMatrix"), t(P) %*% L %*% U %*% Q)))
Detect or standardize a TsparseMatrix
with
unsorted or duplicated pairs.
anyDuplicatedT(x, ...) isUniqueT(x, byrow = FALSE, isT = is(x, "TsparseMatrix")) asUniqueT(x, byrow = FALSE, isT = is(x, "TsparseMatrix")) aggregateT(x)
anyDuplicatedT(x, ...) isUniqueT(x, byrow = FALSE, isT = is(x, "TsparseMatrix")) asUniqueT(x, byrow = FALSE, isT = is(x, "TsparseMatrix")) aggregateT(x)
x |
an R object. |
... |
optional arguments passed to the default method for
generic function |
byrow |
a logical indicating if |
isT |
a logical indicating if |
anyDuplicatedT(x)
returns the index of the first duplicated
pair in
x
(0 if there are no duplicated pairs).
isUniqueT(x)
returns TRUE
if x
is a
TsparseMatrix
with sorted, nonduplicated
pairs and
FALSE
otherwise.
asUniqueT(x)
returns the unique
TsparseMatrix
representation of x
with
sorted, nonduplicated pairs. Values corresponding to
identical
pairs are aggregated by addition, where in the
logical case “addition” refers to logical OR.
aggregateT(x)
aggregates without sorting.
Virtual class TsparseMatrix
.
example("dgTMatrix-class", echo=FALSE) ## -> 'T2' with (i,j,x) slots of length 5 each T2u <- asUniqueT(T2) stopifnot(## They "are" the same (and print the same): all.equal(T2, T2u, tol=0), ## but not internally: anyDuplicatedT(T2) == 2, anyDuplicatedT(T2u) == 0, length(T2 @x) == 5, length(T2u@x) == 3) isUniqueT(T2 ) # FALSE isUniqueT(T2u) # TRUE T3 <- T2u T3[1, c(1,3)] <- 10; T3[2, c(1,5)] <- 20 T3u <- asUniqueT(T3) str(T3u) # sorted in 'j', and within j, sorted in i stopifnot(isUniqueT(T3u)) ## Logical l.TMatrix and n.TMatrix : (L2 <- T2 > 0) validObject(L2u <- asUniqueT(L2)) (N2 <- as(L2, "nMatrix")) validObject(N2u <- asUniqueT(N2)) stopifnot(N2u@i == L2u@i, L2u@i == T2u@i, N2@i == L2@i, L2@i == T2@i, N2u@j == L2u@j, L2u@j == T2u@j, N2@j == L2@j, L2@j == T2@j) # now with a nasty NA [partly failed in Matrix 1.1-5]: L.0N <- L.1N <- L2 L.0N@x[1:2] <- c(FALSE, NA) L.1N@x[1:2] <- c(TRUE, NA) validObject(L.0N) validObject(L.1N) (m.0N <- as.matrix(L.0N)) (m.1N <- as.matrix(L.1N)) stopifnot(identical(10L, which(is.na(m.0N))), !anyNA(m.1N)) symnum(m.0N) symnum(m.1N)
example("dgTMatrix-class", echo=FALSE) ## -> 'T2' with (i,j,x) slots of length 5 each T2u <- asUniqueT(T2) stopifnot(## They "are" the same (and print the same): all.equal(T2, T2u, tol=0), ## but not internally: anyDuplicatedT(T2) == 2, anyDuplicatedT(T2u) == 0, length(T2 @x) == 5, length(T2u@x) == 3) isUniqueT(T2 ) # FALSE isUniqueT(T2u) # TRUE T3 <- T2u T3[1, c(1,3)] <- 10; T3[2, c(1,5)] <- 20 T3u <- asUniqueT(T3) str(T3u) # sorted in 'j', and within j, sorted in i stopifnot(isUniqueT(T3u)) ## Logical l.TMatrix and n.TMatrix : (L2 <- T2 > 0) validObject(L2u <- asUniqueT(L2)) (N2 <- as(L2, "nMatrix")) validObject(N2u <- asUniqueT(N2)) stopifnot(N2u@i == L2u@i, L2u@i == T2u@i, N2@i == L2@i, L2@i == T2@i, N2u@j == L2u@j, L2u@j == T2u@j, N2@j == L2@j, L2@j == T2@j) # now with a nasty NA [partly failed in Matrix 1.1-5]: L.0N <- L.1N <- L2 L.0N@x[1:2] <- c(FALSE, NA) L.1N@x[1:2] <- c(TRUE, NA) validObject(L.0N) validObject(L.1N) (m.0N <- as.matrix(L.0N)) (m.1N <- as.matrix(L.1N)) stopifnot(identical(10L, which(is.na(m.0N))), !anyNA(m.1N)) symnum(m.0N) symnum(m.1N)
Return the matrix obtained by setting to zero elements below a diagonal
(triu
), above a diagonal (tril
), or outside of a general
band (band
).
band(x, k1, k2, ...) triu(x, k = 0L, ...) tril(x, k = 0L, ...)
band(x, k1, k2, ...) triu(x, k = 0L, ...) tril(x, k = 0L, ...)
x |
a matrix-like object |
k , k1 , k2
|
integers specifying the diagonals that are not set to
zero, |
... |
optional arguments passed to methods, currently unused by package Matrix. |
triu(x, k)
is equivalent to band(x, k, dim(x)[2])
.
Similarly,
tril(x, k)
is equivalent to band(x, -dim(x)[1], k)
.
An object of a suitable matrix class, inheriting from
triangularMatrix
where appropriate.
It inherits from sparseMatrix
if
and only if x
does.
method for compressed, sparse, column-oriented matrices.
method for compressed, sparse, row-oriented matrices.
method for sparse matrices in triplet format.
method for diagonal matrices.
method for dense matrices in packed or unpacked format.
method for traditional matrices
of implicit class matrix
.
bandSparse
for the construction of a
banded sparse matrix directly from its non-zero diagonals.
## A random sparse matrix : set.seed(7) m <- matrix(0, 5, 5) m[sample(length(m), size = 14)] <- rep(1:9, length=14) (mm <- as(m, "CsparseMatrix")) tril(mm) # lower triangle tril(mm, -1) # strict lower triangle triu(mm, 1) # strict upper triangle band(mm, -1, 2) # general band (m5 <- Matrix(rnorm(25), ncol = 5)) tril(m5) # lower triangle tril(m5, -1) # strict lower triangle triu(m5, 1) # strict upper triangle band(m5, -1, 2) # general band (m65 <- Matrix(rnorm(30), ncol = 5)) # not square triu(m65) # result not "dtrMatrix" unless square (sm5 <- crossprod(m65)) # symmetric band(sm5, -1, 1)# "dsyMatrix": symmetric band preserves symmetry property as(band(sm5, -1, 1), "sparseMatrix")# often preferable (sm <- round(crossprod(triu(mm/2)))) # sparse symmetric ("dsC*") band(sm, -1,1) # remains "dsC", *however* band(sm, -2,1) # -> "dgC"
## A random sparse matrix : set.seed(7) m <- matrix(0, 5, 5) m[sample(length(m), size = 14)] <- rep(1:9, length=14) (mm <- as(m, "CsparseMatrix")) tril(mm) # lower triangle tril(mm, -1) # strict lower triangle triu(mm, 1) # strict upper triangle band(mm, -1, 2) # general band (m5 <- Matrix(rnorm(25), ncol = 5)) tril(m5) # lower triangle tril(m5, -1) # strict lower triangle triu(m5, 1) # strict upper triangle band(m5, -1, 2) # general band (m65 <- Matrix(rnorm(30), ncol = 5)) # not square triu(m65) # result not "dtrMatrix" unless square (sm5 <- crossprod(m65)) # symmetric band(sm5, -1, 1)# "dsyMatrix": symmetric band preserves symmetry property as(band(sm5, -1, 1), "sparseMatrix")# often preferable (sm <- round(crossprod(triu(mm/2)))) # sparse symmetric ("dsC*") band(sm, -1,1) # remains "dsC", *however* band(sm, -2,1) # -> "dgC"
Construct a sparse banded matrix by specifying its non-zero sup- and super-diagonals.
bandSparse(n, m = n, k, diagonals, symmetric = FALSE, repr = "C", giveCsparse = (repr == "C"))
bandSparse(n, m = n, k, diagonals, symmetric = FALSE, repr = "C", giveCsparse = (repr == "C"))
n , m
|
the matrix dimension |
k |
integer vector of “diagonal numbers”, with identical
meaning as in |
diagonals |
optional list of sub-/super- diagonals; if missing,
the result will be a pattern matrix, i.e., inheriting from
class
|
symmetric |
logical; if true the result will be symmetric
(inheriting from class |
repr |
|
giveCsparse |
(deprecated, replaced with |
a sparse matrix (of class
CsparseMatrix
) of dimension
with diagonal “bands” as specified.
band
, for extraction of matrix bands;
bdiag
, diag
,
sparseMatrix
,
Matrix
.
diags <- list(1:30, 10*(1:20), 100*(1:20)) s1 <- bandSparse(13, k = -c(0:2, 6), diag = c(diags, diags[2]), symm=TRUE) s1 s2 <- bandSparse(13, k = c(0:2, 6), diag = c(diags, diags[2]), symm=TRUE) stopifnot(identical(s1, t(s2)), is(s1,"dsCMatrix")) ## a pattern Matrix of *full* (sub-)diagonals: bk <- c(0:4, 7,9) (s3 <- bandSparse(30, k = bk, symm = TRUE)) ## If you want a pattern matrix, but with "sparse"-diagonals, ## you currently need to go via logical sparse: lLis <- lapply(list(rpois(20, 2), rpois(20, 1), rpois(20, 3))[c(1:3, 2:3, 3:2)], as.logical) (s4 <- bandSparse(20, k = bk, symm = TRUE, diag = lLis)) (s4. <- as(drop0(s4), "nsparseMatrix")) n <- 1e4 bk <- c(0:5, 7,11) bMat <- matrix(1:8, n, 8, byrow=TRUE) bLis <- as.data.frame(bMat) B <- bandSparse(n, k = bk, diag = bLis) Bs <- bandSparse(n, k = bk, diag = bLis, symmetric=TRUE) B [1:15, 1:30] Bs[1:15, 1:30] ## can use a list *or* a matrix for specifying the diagonals: stopifnot(identical(B, bandSparse(n, k = bk, diag = bMat)), identical(Bs, bandSparse(n, k = bk, diag = bMat, symmetric=TRUE)) , inherits(B, "dtCMatrix") # triangular! )
diags <- list(1:30, 10*(1:20), 100*(1:20)) s1 <- bandSparse(13, k = -c(0:2, 6), diag = c(diags, diags[2]), symm=TRUE) s1 s2 <- bandSparse(13, k = c(0:2, 6), diag = c(diags, diags[2]), symm=TRUE) stopifnot(identical(s1, t(s2)), is(s1,"dsCMatrix")) ## a pattern Matrix of *full* (sub-)diagonals: bk <- c(0:4, 7,9) (s3 <- bandSparse(30, k = bk, symm = TRUE)) ## If you want a pattern matrix, but with "sparse"-diagonals, ## you currently need to go via logical sparse: lLis <- lapply(list(rpois(20, 2), rpois(20, 1), rpois(20, 3))[c(1:3, 2:3, 3:2)], as.logical) (s4 <- bandSparse(20, k = bk, symm = TRUE, diag = lLis)) (s4. <- as(drop0(s4), "nsparseMatrix")) n <- 1e4 bk <- c(0:5, 7,11) bMat <- matrix(1:8, n, 8, byrow=TRUE) bLis <- as.data.frame(bMat) B <- bandSparse(n, k = bk, diag = bLis) Bs <- bandSparse(n, k = bk, diag = bLis, symmetric=TRUE) B [1:15, 1:30] Bs[1:15, 1:30] ## can use a list *or* a matrix for specifying the diagonals: stopifnot(identical(B, bandSparse(n, k = bk, diag = bMat)), identical(Bs, bandSparse(n, k = bk, diag = bMat, symmetric=TRUE)) , inherits(B, "dtCMatrix") # triangular! )
Build a block diagonal matrix given several building block matrices.
bdiag(...) .bdiag(lst)
bdiag(...) .bdiag(lst)
... |
individual matrices or a |
lst |
non-empty |
For non-trivial argument list, bdiag()
calls .bdiag()
.
The latter maybe useful to programmers.
A sparse matrix obtained by combining the arguments into a block diagonal matrix.
The value of bdiag()
inherits from class
CsparseMatrix
, whereas
.bdiag()
returns a TsparseMatrix
.
This function has been written and is efficient for the case of relatively few block matrices which are typically sparse themselves.
It is currently inefficient for the case of many small dense
block matrices.
For the case of many dense matrices,
the
bdiag_m()
function in the ‘Examples’ is an order of
magnitude faster.
Martin Maechler, built on a version posted by Berton Gunter to
R-help; earlier versions have been posted by other authors, notably
Scott Chasalow to S-news. Doug Bates's faster implementation builds
on TsparseMatrix
objects.
Diagonal
for constructing matrices of
class diagonalMatrix
, or kronecker
which also works for "Matrix"
inheriting matrices.
bandSparse
constructs a banded sparse matrix from
its non-zero sub-/super - diagonals.
Note that other CRAN R packages have own versions of bdiag()
which return traditional matrices.
bdiag(matrix(1:4, 2), diag(3)) ## combine "Matrix" class and traditional matrices: bdiag(Diagonal(2), matrix(1:3, 3,4), diag(3:2)) mlist <- list(1, 2:3, diag(x=5:3), 27, cbind(1,3:6), 100:101) bdiag(mlist) stopifnot(identical(bdiag(mlist), bdiag(lapply(mlist, as.matrix)))) ml <- c(as(matrix((1:24)%% 11 == 0, 6,4),"nMatrix"), rep(list(Diagonal(2, x=TRUE)), 3)) mln <- c(ml, Diagonal(x = 1:3)) stopifnot(is(bdiag(ml), "lsparseMatrix"), is(bdiag(mln),"dsparseMatrix") ) ## random (diagonal-)block-triangular matrices: rblockTri <- function(nb, max.ni, lambda = 3) { .bdiag(replicate(nb, { n <- sample.int(max.ni, 1) tril(Matrix(rpois(n * n, lambda = lambda), n, n)) })) } (T4 <- rblockTri(4, 10, lambda = 1)) image(T1 <- rblockTri(12, 20)) ##' Fast version of Matrix :: .bdiag() -- for the case of *many* (k x k) matrices: ##' @param lmat list(<mat1>, <mat2>, ....., <mat_N>) where each mat_j is a k x k 'matrix' ##' @return a sparse (N*k x N*k) matrix of class \code{"\linkS4class{dgCMatrix}"}. bdiag_m <- function(lmat) { ## Copyright (C) 2016 Martin Maechler, ETH Zurich if(!length(lmat)) return(new("dgCMatrix")) stopifnot(is.list(lmat), is.matrix(lmat[[1]]), (k <- (d <- dim(lmat[[1]]))[1]) == d[2], # k x k all(vapply(lmat, dim, integer(2)) == k)) # all of them N <- length(lmat) if(N * k > .Machine$integer.max) stop("resulting matrix too large; would be M x M, with M=", N*k) M <- as.integer(N * k) ## result: an M x M matrix new("dgCMatrix", Dim = c(M,M), ## 'i :' maybe there's a faster way (w/o matrix indexing), but elegant? i = as.vector(matrix(0L:(M-1L), nrow=k)[, rep(seq_len(N), each=k)]), p = k * 0L:M, x = as.double(unlist(lmat, recursive=FALSE, use.names=FALSE))) } l12 <- replicate(12, matrix(rpois(16, lambda = 6.4), 4, 4), simplify=FALSE) dim(T12 <- bdiag_m(l12))# 48 x 48 T12[1:20, 1:20]
bdiag(matrix(1:4, 2), diag(3)) ## combine "Matrix" class and traditional matrices: bdiag(Diagonal(2), matrix(1:3, 3,4), diag(3:2)) mlist <- list(1, 2:3, diag(x=5:3), 27, cbind(1,3:6), 100:101) bdiag(mlist) stopifnot(identical(bdiag(mlist), bdiag(lapply(mlist, as.matrix)))) ml <- c(as(matrix((1:24)%% 11 == 0, 6,4),"nMatrix"), rep(list(Diagonal(2, x=TRUE)), 3)) mln <- c(ml, Diagonal(x = 1:3)) stopifnot(is(bdiag(ml), "lsparseMatrix"), is(bdiag(mln),"dsparseMatrix") ) ## random (diagonal-)block-triangular matrices: rblockTri <- function(nb, max.ni, lambda = 3) { .bdiag(replicate(nb, { n <- sample.int(max.ni, 1) tril(Matrix(rpois(n * n, lambda = lambda), n, n)) })) } (T4 <- rblockTri(4, 10, lambda = 1)) image(T1 <- rblockTri(12, 20)) ##' Fast version of Matrix :: .bdiag() -- for the case of *many* (k x k) matrices: ##' @param lmat list(<mat1>, <mat2>, ....., <mat_N>) where each mat_j is a k x k 'matrix' ##' @return a sparse (N*k x N*k) matrix of class \code{"\linkS4class{dgCMatrix}"}. bdiag_m <- function(lmat) { ## Copyright (C) 2016 Martin Maechler, ETH Zurich if(!length(lmat)) return(new("dgCMatrix")) stopifnot(is.list(lmat), is.matrix(lmat[[1]]), (k <- (d <- dim(lmat[[1]]))[1]) == d[2], # k x k all(vapply(lmat, dim, integer(2)) == k)) # all of them N <- length(lmat) if(N * k > .Machine$integer.max) stop("resulting matrix too large; would be M x M, with M=", N*k) M <- as.integer(N * k) ## result: an M x M matrix new("dgCMatrix", Dim = c(M,M), ## 'i :' maybe there's a faster way (w/o matrix indexing), but elegant? i = as.vector(matrix(0L:(M-1L), nrow=k)[, rep(seq_len(N), each=k)]), p = k * 0L:M, x = as.double(unlist(lmat, recursive=FALSE, use.names=FALSE))) } l12 <- replicate(12, matrix(rpois(16, lambda = 6.4), 4, 4), simplify=FALSE) dim(T12 <- bdiag_m(l12))# 48 x 48 T12[1:20, 1:20]
%&%
and MethodsFor boolean or “pattern” matrices, i.e., R objects of
class nMatrix
, it is natural to allow matrix
products using boolean instead of numerical arithmetic.
In package Matrix, we use the binary operator %&%
(aka
“infix”) function) for this and provide methods for all our
matrices and the traditional R matrices (see matrix
).
a pattern matrix, i.e., inheriting from "nMatrix"
,
or an "ldiMatrix"
in case of a diagonal matrix.
We provide methods for both the “traditional” (R base) matrices
and numeric vectors and conceptually all matrices and
sparseVector
s in package Matrix.
signature(x = "ANY", y = "ANY")
signature(x = "ANY", y = "Matrix")
signature(x = "Matrix", y = "ANY")
signature(x = "nMatrix", y = "nMatrix")
signature(x = "nMatrix", y = "nsparseMatrix")
signature(x = "nsparseMatrix", y = "nMatrix")
signature(x = "nsparseMatrix", y = "nsparseMatrix")
signature(x = "sparseVector", y = "sparseVector")
These boolean arithmetic matrix products had been newly introduced for Matrix 1.2.0 (March 2015). Its implementation has still not been tested extensively.
Originally, it was left unspecified how non-structural zeros, i.e., 0
's
as part of the M@x
slot should be treated for numeric
("dMatrix"
) and logical ("lMatrix"
)
sparse matrices. We now specify that boolean matrix products should behave as if
applied to drop0(M)
, i.e., as if dropping such zeros from
the matrix before using it.
Equivalently, for all matrices M
, boolean arithmetic should work as if
applied to M != 0
(or M != FALSE
).
The current implementation ends up coercing both x
and y
to
(virtual) class nsparseMatrix
which may be quite inefficient
for dense matrices. A future implementation may well return a matrix
with different class, but the “same” content, i.e., the
same matrix entries .
%*%
, crossprod()
, or tcrossprod()
,
for (regular) matrix product methods.
set.seed(7) L <- Matrix(rnorm(20) > 1, 4,5) (N <- as(L, "nMatrix")) L. <- L; L.[1:2,1] <- TRUE; L.@x[1:2] <- FALSE; L. # has "zeros" to drop0() D <- Matrix(round(rnorm(30)), 5,6) # -> values in -1:1 (for this seed) L %&% D stopifnot(identical(L %&% D, N %&% D), all(L %&% D == as((L %*% abs(D)) > 0, "sparseMatrix"))) ## cross products , possibly with boolArith = TRUE : crossprod(N) # -> sparse patter'n' (TRUE/FALSE : boolean arithmetic) crossprod(N +0) # -> numeric Matrix (with same "pattern") stopifnot(all(crossprod(N) == t(N) %&% N), identical(crossprod(N), crossprod(N +0, boolArith=TRUE)), identical(crossprod(L), crossprod(N , boolArith=FALSE))) crossprod(D, boolArith = TRUE) # pattern: "nsCMatrix" crossprod(L, boolArith = TRUE) # ditto crossprod(L, boolArith = FALSE) # numeric: "dsCMatrix"
set.seed(7) L <- Matrix(rnorm(20) > 1, 4,5) (N <- as(L, "nMatrix")) L. <- L; L.[1:2,1] <- TRUE; L.@x[1:2] <- FALSE; L. # has "zeros" to drop0() D <- Matrix(round(rnorm(30)), 5,6) # -> values in -1:1 (for this seed) L %&% D stopifnot(identical(L %&% D, N %&% D), all(L %&% D == as((L %*% abs(D)) > 0, "sparseMatrix"))) ## cross products , possibly with boolArith = TRUE : crossprod(N) # -> sparse patter'n' (TRUE/FALSE : boolean arithmetic) crossprod(N +0) # -> numeric Matrix (with same "pattern") stopifnot(all(crossprod(N) == t(N) %&% N), identical(crossprod(N), crossprod(N +0, boolArith=TRUE)), identical(crossprod(L), crossprod(N , boolArith=FALSE))) crossprod(D, boolArith = TRUE) # pattern: "nsCMatrix" crossprod(L, boolArith = TRUE) # ditto crossprod(L, boolArith = FALSE) # numeric: "dsCMatrix"
Classes BunchKaufman
and pBunchKaufman
represent
Bunch-Kaufman factorizations of real,
symmetric matrices
, having the general form
where
and
are symmetric, block diagonal
matrices composed of
and
or
diagonal blocks;
is the product of
row-permuted unit upper triangular
matrices, each having nonzero entries above the diagonal in 1 or 2 columns;
and
is the product of
row-permuted unit lower triangular
matrices, each having nonzero entries below the diagonal in 1 or 2 columns.
These classes store the nonzero entries of the
or
factors,
which are individually sparse,
in a dense format as a vector of length
(
BunchKaufman
) or
(
pBunchKaufman
),
the latter giving the “packed” representation.
Dim
, Dimnames
inherited from virtual class
MatrixFactorization
.
uplo
a string, either "U"
or "L"
,
indicating which triangle (upper or lower) of the factorized
symmetric matrix was used to compute the factorization and
in turn how the x
slot is partitioned.
x
a numeric vector of length n*n
(BunchKaufman
) or n*(n+1)/2
(pBunchKaufman
),
where n=Dim[1]
.
The details of the representation are specified by the manual
for LAPACK routines dsytrf
and dsptrf
.
perm
an integer vector of length n=Dim[1]
specifying row and column interchanges as described in the manual
for LAPACK routines dsytrf
and dsptrf
.
Class BunchKaufmanFactorization
, directly.
Class MatrixFactorization
, by class
BunchKaufmanFactorization
, distance 2.
Objects can be generated directly by calls of the form
new("BunchKaufman", ...)
or new("pBunchKaufman", ...)
,
but they are more typically obtained as the value of
BunchKaufman(x)
for x
inheriting from
dsyMatrix
or dspMatrix
.
coerce
signature(from = "BunchKaufman", to = "dtrMatrix")
:
returns a dtrMatrix
, useful for inspecting
the internal representation of the factorization; see ‘Note’.
coerce
signature(from = "pBunchKaufman", to = "dtpMatrix")
:
returns a dtpMatrix
, useful for inspecting
the internal representation of the factorization; see ‘Note’.
determinant
signature(from = "p?BunchKaufman", logarithm = "logical")
:
computes the determinant of the factorized matrix
or its logarithm.
expand1
signature(x = "p?BunchKaufman")
:
see expand1-methods
.
expand2
signature(x = "p?BunchKaufman")
:
see expand2-methods
.
solve
signature(a = "p?BunchKaufman", b = .)
:
see solve-methods
.
In Matrix < 1.6-0
, class BunchKaufman
extended
dtrMatrix
and class pBunchKaufman
extended
dtpMatrix
, reflecting the fact that the internal
representation of the factorization is fundamentally triangular:
there are “parameters”, and these
can be arranged systematically to form an
triangular matrix.
Matrix
1.6-0
removed these extensions so that methods
would no longer be inherited from dtrMatrix
and dtpMatrix
.
The availability of such methods gave the wrong impression that
BunchKaufman
and pBunchKaufman
represent a (singular)
matrix, when in fact they represent an ordered set of matrix factors.
The coercions as(., "dtrMatrix")
and as(., "dtpMatrix")
are provided for users who understand the caveats.
The LAPACK source code, including documentation; see https://netlib.org/lapack/double/dsytrf.f and https://netlib.org/lapack/double/dsptrf.f.
Golub, G. H., & Van Loan, C. F. (2013). Matrix computations (4th ed.). Johns Hopkins University Press. doi:10.56021/9781421407944
Class dsyMatrix
and its packed counterpart.
Generic functions BunchKaufman
,
expand1
, and expand2
.
showClass("BunchKaufman") set.seed(1) n <- 6L (A <- forceSymmetric(Matrix(rnorm(n * n), n, n))) ## With dimnames, to see that they are propagated : dimnames(A) <- rep.int(list(paste0("x", seq_len(n))), 2L) (bk.A <- BunchKaufman(A)) str(e.bk.A <- expand2(bk.A, complete = FALSE), max.level = 2L) str(E.bk.A <- expand2(bk.A, complete = TRUE), max.level = 2L) ## Underlying LAPACK representation (m.bk.A <- as(bk.A, "dtrMatrix")) stopifnot(identical(as(m.bk.A, "matrix"), `dim<-`(bk.A@x, bk.A@Dim))) ## Number of factors is 2*b+1, b <= n, which can be nontrivial ... (b <- (length(E.bk.A) - 1L) %/% 2L) ae1 <- function(a, b, ...) all.equal(as(a, "matrix"), as(b, "matrix"), ...) ae2 <- function(a, b, ...) ae1(unname(a), unname(b), ...) ## A ~ U DU U', U := prod(Pk Uk) in floating point stopifnot(exprs = { identical(names(e.bk.A), c("U", "DU", "U.")) identical(e.bk.A[["U" ]], Reduce(`%*%`, E.bk.A[seq_len(b)])) identical(e.bk.A[["U."]], t(e.bk.A[["U"]])) ae1(A, with(e.bk.A, U %*% DU %*% U.)) }) ## Factorization handled as factorized matrix b <- rnorm(n) stopifnot(identical(det(A), det(bk.A)), identical(solve(A, b), solve(bk.A, b)))
showClass("BunchKaufman") set.seed(1) n <- 6L (A <- forceSymmetric(Matrix(rnorm(n * n), n, n))) ## With dimnames, to see that they are propagated : dimnames(A) <- rep.int(list(paste0("x", seq_len(n))), 2L) (bk.A <- BunchKaufman(A)) str(e.bk.A <- expand2(bk.A, complete = FALSE), max.level = 2L) str(E.bk.A <- expand2(bk.A, complete = TRUE), max.level = 2L) ## Underlying LAPACK representation (m.bk.A <- as(bk.A, "dtrMatrix")) stopifnot(identical(as(m.bk.A, "matrix"), `dim<-`(bk.A@x, bk.A@Dim))) ## Number of factors is 2*b+1, b <= n, which can be nontrivial ... (b <- (length(E.bk.A) - 1L) %/% 2L) ae1 <- function(a, b, ...) all.equal(as(a, "matrix"), as(b, "matrix"), ...) ae2 <- function(a, b, ...) ae1(unname(a), unname(b), ...) ## A ~ U DU U', U := prod(Pk Uk) in floating point stopifnot(exprs = { identical(names(e.bk.A), c("U", "DU", "U.")) identical(e.bk.A[["U" ]], Reduce(`%*%`, E.bk.A[seq_len(b)])) identical(e.bk.A[["U."]], t(e.bk.A[["U"]])) ae1(A, with(e.bk.A, U %*% DU %*% U.)) }) ## Factorization handled as factorized matrix b <- rnorm(n) stopifnot(identical(det(A), det(bk.A)), identical(solve(A, b), solve(bk.A, b)))
Computes the Bunch-Kaufman factorization of an
real, symmetric matrix
, which has the general form
where
and
are symmetric, block diagonal
matrices composed of
and
or
diagonal blocks;
is the product of
row-permuted unit upper triangular
matrices, each having nonzero entries above the diagonal in 1 or 2 columns;
and
is the product of
row-permuted unit lower triangular
matrices, each having nonzero entries below the diagonal in 1 or 2 columns.
Methods are built on LAPACK routines dsytrf
and dsptrf
.
BunchKaufman(x, ...) ## S4 method for signature 'dsyMatrix' BunchKaufman(x, warnSing = TRUE, ...) ## S4 method for signature 'dspMatrix' BunchKaufman(x, warnSing = TRUE, ...) ## S4 method for signature 'matrix' BunchKaufman(x, uplo = "U", ...)
BunchKaufman(x, ...) ## S4 method for signature 'dsyMatrix' BunchKaufman(x, warnSing = TRUE, ...) ## S4 method for signature 'dspMatrix' BunchKaufman(x, warnSing = TRUE, ...) ## S4 method for signature 'matrix' BunchKaufman(x, uplo = "U", ...)
x |
a finite symmetric matrix or
|
warnSing |
a logical indicating if a warning should
be signaled for singular |
uplo |
a string, either |
... |
further arguments passed to or from methods. |
An object representing the factorization, inheriting from
virtual class BunchKaufmanFactorization
.
The specific class is BunchKaufman
unless
x
inherits from virtual class packedMatrix
,
in which case it is pBunchKaufman
.
The LAPACK source code, including documentation; see https://netlib.org/lapack/double/dsytrf.f and https://netlib.org/lapack/double/dsptrf.f.
Golub, G. H., & Van Loan, C. F. (2013). Matrix computations (4th ed.). Johns Hopkins University Press. doi:10.56021/9781421407944
Classes BunchKaufman
and
pBunchKaufman
and their methods.
Classes dsyMatrix
and
dspMatrix
.
Generic functions expand1
and expand2
,
for constructing matrix factors from the result.
Generic functions Cholesky
, Schur
,
lu
, and qr
,
for computing other factorizations.
showMethods("BunchKaufman", inherited = FALSE) set.seed(0) data(CAex, package = "Matrix") class(CAex) # dgCMatrix isSymmetric(CAex) # symmetric, but not formally A <- as(CAex, "symmetricMatrix") class(A) # dsCMatrix ## Have methods for denseMatrix (unpacked and packed), ## but not yet sparseMatrix ... ## Not run: (bk.A <- BunchKaufman(A)) ## End(Not run) (bk.A <- BunchKaufman(as(A, "unpackedMatrix"))) ## A ~ U DU U' in floating point str(e.bk.A <- expand2(bk.A), max.level = 2L) stopifnot(all.equal(as(A, "matrix"), as(Reduce(`%*%`, e.bk.A), "matrix")))
showMethods("BunchKaufman", inherited = FALSE) set.seed(0) data(CAex, package = "Matrix") class(CAex) # dgCMatrix isSymmetric(CAex) # symmetric, but not formally A <- as(CAex, "symmetricMatrix") class(A) # dsCMatrix ## Have methods for denseMatrix (unpacked and packed), ## but not yet sparseMatrix ... ## Not run: (bk.A <- BunchKaufman(A)) ## End(Not run) (bk.A <- BunchKaufman(as(A, "unpackedMatrix"))) ## A ~ U DU U' in floating point str(e.bk.A <- expand2(bk.A), max.level = 2L) stopifnot(all.equal(as(A, "matrix"), as(Reduce(`%*%`, e.bk.A), "matrix")))
An example of a sparse matrix for which eigen()
seemed
to be difficult, an unscaled version of this has been posted to the
web, accompanying an E-mail to R-help
(https://stat.ethz.ch/mailman/listinfo/r-help), by
Casper J Albers, Open University, UK.
data(CAex)
data(CAex)
This is a symmetric matrix with 216
non-zero entries in five bands, stored as sparse matrix of class
dgCMatrix
.
Historical note (2006-03-30):
In earlier versions of R, eigen(CAex)
fell into an
infinite loop whereas eigen(CAex, EISPACK=TRUE)
had been okay.
data(CAex, package = "Matrix") str(CAex) # of class "dgCMatrix" image(CAex)# -> it's a simple band matrix with 5 bands ## and the eigen values are basically 1 (42 times) and 0 (30 x): zapsmall(ev <- eigen(CAex, only.values=TRUE)$values) ## i.e., the matrix is symmetric, hence sCA <- as(CAex, "symmetricMatrix") ## and stopifnot(class(sCA) == "dsCMatrix", as(sCA, "matrix") == as(CAex, "matrix"))
data(CAex, package = "Matrix") str(CAex) # of class "dgCMatrix" image(CAex)# -> it's a simple band matrix with 5 bands ## and the eigen values are basically 1 (42 times) and 0 (30 x): zapsmall(ev <- eigen(CAex, only.values=TRUE)$values) ## i.e., the matrix is symmetric, hence sCA <- as(CAex, "symmetricMatrix") ## and stopifnot(class(sCA) == "dsCMatrix", as(sCA, "matrix") == as(CAex, "matrix"))
The base functions cbind
and rbind
are
defined for an arbitrary number of arguments and hence have the first
formal argument ...
. Now, when S4 objects are found among the arguments,
base cbind()
and rbind()
internally “dispatch”
recursively, calling cbind2
or rbind2
respectively, where these have methods defined and so should dispatch
appropriately.
cbind2()
and rbind2()
are from the
methods package, i.e., standard R, and have been provided for
binding together two matrices, where in Matrix, we have
defined methods for these and the 'Matrix'
matrices.
## cbind(..., deparse.level = 1) ## rbind(..., deparse.level = 1) ## S4 method for signature 'Matrix,Matrix' cbind2(x, y, ...) ## S4 method for signature 'Matrix,Matrix' rbind2(x, y, ...)
## cbind(..., deparse.level = 1) ## rbind(..., deparse.level = 1) ## S4 method for signature 'Matrix,Matrix' cbind2(x, y, ...) ## S4 method for signature 'Matrix,Matrix' rbind2(x, y, ...)
... |
for |
deparse.level |
integer controlling the construction of labels
in the case of non-matrix-like arguments; see |
x , y
|
vector- or matrix-like R objects to be bound together. |
typically a ‘matrix-like’ object of a similar
class
as the first argument in ...
.
Note that sometimes by default, the result is a
sparseMatrix
if one of the arguments is (even in
the case where this is not efficient). In other cases,
the result is chosen to be sparse when there are more zero entries is
than non-zero ones (as the default sparse
in
Matrix()
).
Martin Maechler
Our class definition help pages mentioning cbind2()
and
rbind2()
methods:
"denseMatrix"
,
"diagonalMatrix"
,
"indMatrix"
.
(a <- matrix(c(2:1,1:2), 2,2)) (M1 <- cbind(0, rbind(a, 7))) # a traditional matrix D <- Diagonal(2) (M2 <- cbind(4, a, D, -1, D, 0)) # a sparse Matrix stopifnot(validObject(M2), inherits(M2, "sparseMatrix"), dim(M2) == c(2,9))
(a <- matrix(c(2:1,1:2), 2,2)) (M1 <- cbind(0, rbind(a, 7))) # a traditional matrix D <- Diagonal(2) (M2 <- cbind(4, a, D, -1, D, 0)) # a sparse Matrix stopifnot(validObject(M2), inherits(M2, "sparseMatrix"), dim(M2) == c(2,9))
CHMfactor
is the virtual class of sparse Cholesky
factorizations of real, symmetric
matrices
, having the general form
or (equivalently)
where
is a permutation matrix,
is a unit lower triangular matrix,
is a diagonal matrix, and
.
The second equalities hold only for positive semidefinite
,
for which the diagonal entries of
are non-negative
and
is well-defined.
The implementation of class CHMfactor
is based on
CHOLMOD's C-level cholmod_factor_struct
. Virtual
subclasses CHMsimpl
and CHMsuper
separate
the simplicial and supernodal variants. These have nonvirtual
subclasses [dn]CHMsimpl
and [dn]CHMsuper
,
where prefix ‘d’ and prefix ‘n’ are reserved
for numeric and symbolic factorizations, respectively.
isLDL(x)
isLDL(x)
x |
an object inheriting from virtual class |
isLDL(x)
returns TRUE
or FALSE
:
TRUE
if x
stores the lower triangular entries
of ,
FALSE
if x
stores the lower triangular entries
of .
Of CHMfactor
:
Dim
, Dimnames
inherited from virtual class
MatrixFactorization
.
colcount
an integer vector of length Dim[1]
giving an estimate of the number of nonzero entries in
each column of the lower triangular Cholesky factor.
If symbolic analysis was performed prior to factorization,
then the estimate is exact.
perm
a 0-based integer vector of length Dim[1]
specifying the permutation applied to the rows and columns
of the factorized matrix. perm
of length 0 is valid and
equivalent to the identity permutation, implying no pivoting.
type
an integer vector of length 6 specifying
details of the factorization. The elements correspond to
members ordering
, is_ll
, is_super
,
is_monotonic
, maxcsize
, and maxesize
of the original cholmod_factor_struct
.
Simplicial and supernodal factorizations are distinguished
by is_super
. Simplicial factorizations do not use
maxcsize
or maxesize
. Supernodal factorizations
do not use is_ll
or is_monotonic
.
Of CHMsimpl
(all unused by nCHMsimpl
):
nz
an integer vector of length Dim[1]
giving the number of nonzero entries in each column of the
lower triangular Cholesky factor. There is at least one
nonzero entry in each column, because the diagonal elements
of the factor are stored explicitly.
p
an integer vector of length Dim[1]+1
.
Row indices of nonzero entries in column j
of the
lower triangular Cholesky factor are obtained as
i[p[j]+seq_len(nz[j])]+1
.
i
an integer vector of length greater than or equal
to sum(nz)
containing the row indices of nonzero entries
in the lower triangular Cholesky factor. These are grouped by
column and sorted within columns, but the columns themselves
need not be ordered monotonically. Columns may be overallocated,
i.e., the number of elements of i
reserved for column
j
may exceed nz[j]
.
prv
, nxt
integer vectors of length
Dim[1]+2
indicating the order in which the columns of
the lower triangular Cholesky factor are stored in i
and x
.
Starting from j <- Dim[1]+2
,
the recursion j <- nxt[j+1]+1
traverses the columns
in forward order and terminates when nxt[j+1] = -1
.
Starting from j <- Dim[1]+1
,
the recursion j <- prv[j+1]+1
traverses the columns
in backward order and terminates when prv[j+1] = -1
.
Of dCHMsimpl
:
x
a numeric vector parallel to i
containing
the corresponding nonzero entries of the lower triangular
Cholesky factor or (if and only if
type[2]
is 0) of the lower triangular matrix .
Of CHMsuper
:
super
, pi
, px
integer vectors of
length nsuper+1
, where nsuper
is the number of
supernodes. super[j]+1
is the index of the leftmost
column of supernode j
. The row indices of supernode
j
are obtained as s[pi[j]+seq_len(pi[j+1]-pi[j])]+1
.
The numeric entries of supernode j
are obtained as
x[px[j]+seq_len(px[j+1]-px[j])]+1
(if slot x
is available).
s
an integer vector of length greater than or equal
to Dim[1]
containing the row indices of the supernodes.
s
may contain duplicates, but not within a supernode,
where the row indices must be increasing.
Of dCHMsuper
:
x
a numeric vector of length less than or equal to
prod(Dim)
containing the numeric entries of the supernodes.
Class MatrixFactorization
, directly.
Objects can be generated directly by calls of the form
new("dCHMsimpl", ...)
, etc., but dCHMsimpl
and
dCHMsuper
are more typically obtained as the value of
Cholesky(x, ...)
for x
inheriting from
sparseMatrix
(often dsCMatrix
).
There is currently no API outside of calls to new
for generating nCHMsimpl
and nCHMsuper
. These
classes are vestigial and may be formally deprecated in a future
version of Matrix.
coerce
signature(from = "CHMsimpl", to = "dtCMatrix")
:
returns a dtCMatrix
representing
the lower triangular Cholesky factor or
the lower triangular matrix
,
the latter if and only if
from@type[2]
is 0.
coerce
signature(from = "CHMsuper", to = "dgCMatrix")
:
returns a dgCMatrix
representing
the lower triangular Cholesky factor . Note that,
for supernodes spanning two or more columns, the supernodal
algorithm by design stores non-structural zeros above
the main diagonal, hence
dgCMatrix
is
indeed more appropriate than dtCMatrix
as a coercion target.
determinant
signature(from = "CHMfactor", logarithm = "logical")
:
behaves according to an optional argument sqrt
.
If sqrt = FALSE
, then this method computes the determinant
of the factorized matrix or its logarithm.
If
sqrt = TRUE
, then this method computes the determinant
of the factor or
its logarithm, giving
NaN
for the modulus when
has negative diagonal elements. For backwards compatibility,
the default value of
sqrt
is TRUE
, but that can
be expected change in a future version of Matrix, hence
defensive code will always set sqrt
(to TRUE
,
if the code must remain backwards compatible with Matrix
< 1.6-0
). Calls to this method not setting sqrt
may warn about the pending change. The warnings can be disabled
with options(Matrix.warnSqrtDefault = 0)
.
diag
signature(x = "CHMfactor")
:
returns a numeric vector of length containing the diagonal
elements of
, which (if they are all non-negative)
are the squared diagonal elements of
.
expand
signature(x = "CHMfactor")
:
see expand-methods
.
expand1
signature(x = "CHMsimpl")
:
see expand1-methods
.
expand1
signature(x = "CHMsuper")
:
see expand1-methods
.
expand2
signature(x = "CHMsimpl")
:
see expand2-methods
.
expand2
signature(x = "CHMsuper")
:
see expand2-methods
.
image
signature(x = "CHMfactor")
:
see image-methods
.
nnzero
signature(x = "CHMfactor")
:
see nnzero-methods
.
solve
signature(a = "CHMfactor", b = .)
:
see solve-methods
.
update
signature(object = "CHMfactor")
:
returns a copy of object
with the same nonzero pattern
but with numeric entries updated according to additional
arguments parent
and mult
, where parent
is (coercible to) a dsCMatrix
or a
dgCMatrix
and mult
is a numeric
vector of positive length.
The numeric entries are updated with those of the Cholesky
factor of F(parent) + mult[1] * I
, i.e.,
F(parent)
plus mult[1]
times the identity matrix,
where F = identity
for symmetric parent
and F = tcrossprod
for other parent
.
The nonzero pattern of F(parent)
must match
that of S
if object = Cholesky(S, ...)
.
updown
signature(update = ., C = ., object = "CHMfactor")
:
see updown-methods
.
The CHOLMOD source code; see
https://github.com/DrTimothyAldenDavis/SuiteSparse,
notably the header file ‘CHOLMOD/Include/cholmod.h’
defining cholmod_factor_struct
.
Chen, Y., Davis, T. A., Hager, W. W., & Rajamanickam, S. (2008). Algorithm 887: CHOLMOD, supernodal sparse Cholesky factorization and update/downdate. ACM Transactions on Mathematical Software, 35(3), Article 22, 1-14. doi:10.1145/1391989.1391995
Amestoy, P. R., Davis, T. A., & Duff, I. S. (2004). Algorithm 837: AMD, an approximate minimum degree ordering algorithm. ACM Transactions on Mathematical Software, 17(4), 886-905. doi:10.1145/1024074.1024081
Golub, G. H., & Van Loan, C. F. (2013). Matrix computations (4th ed.). Johns Hopkins University Press. doi:10.56021/9781421407944
Class dsCMatrix
.
Generic functions Cholesky
, updown
,
expand1
and expand2
.
showClass("dCHMsimpl") showClass("dCHMsuper") set.seed(2) m <- 1000L n <- 200L M <- rsparsematrix(m, n, 0.01) A <- crossprod(M) ## With dimnames, to see that they are propagated : dimnames(A) <- dn <- rep.int(list(paste0("x", seq_len(n))), 2L) (ch.A <- Cholesky(A)) # pivoted, by default str(e.ch.A <- expand2(ch.A, LDL = TRUE), max.level = 2L) str(E.ch.A <- expand2(ch.A, LDL = FALSE), max.level = 2L) ae1 <- function(a, b, ...) all.equal(as(a, "matrix"), as(b, "matrix"), ...) ae2 <- function(a, b, ...) ae1(unname(a), unname(b), ...) ## A ~ P1' L1 D L1' P1 ~ P1' L L' P1 in floating point stopifnot(exprs = { identical(names(e.ch.A), c("P1.", "L1", "D", "L1.", "P1")) identical(names(E.ch.A), c("P1.", "L" , "L." , "P1")) identical(e.ch.A[["P1"]], new("pMatrix", Dim = c(n, n), Dimnames = c(list(NULL), dn[2L]), margin = 2L, perm = invertPerm(ch.A@perm, 0L, 1L))) identical(e.ch.A[["P1."]], t(e.ch.A[["P1"]])) identical(e.ch.A[["L1."]], t(e.ch.A[["L1"]])) identical(E.ch.A[["L." ]], t(E.ch.A[["L" ]])) identical(e.ch.A[["D"]], Diagonal(x = diag(ch.A))) all.equal(E.ch.A[["L"]], with(e.ch.A, L1 %*% sqrt(D))) ae1(A, with(e.ch.A, P1. %*% L1 %*% D %*% L1. %*% P1)) ae1(A, with(E.ch.A, P1. %*% L %*% L. %*% P1)) ae2(A[ch.A@perm + 1L, ch.A@perm + 1L], with(e.ch.A, L1 %*% D %*% L1.)) ae2(A[ch.A@perm + 1L, ch.A@perm + 1L], with(E.ch.A, L %*% L. )) }) ## Factorization handled as factorized matrix ## (in some cases only optionally, depending on arguments) b <- rnorm(n) stopifnot(identical(det(A), det(ch.A, sqrt = FALSE)), identical(solve(A, b), solve(ch.A, b, system = "A"))) u1 <- update(ch.A, A , mult = sqrt(2)) u2 <- update(ch.A, t(M), mult = sqrt(2)) # updating with crossprod(M), not M stopifnot(all.equal(u1, u2, tolerance = 1e-14))
showClass("dCHMsimpl") showClass("dCHMsuper") set.seed(2) m <- 1000L n <- 200L M <- rsparsematrix(m, n, 0.01) A <- crossprod(M) ## With dimnames, to see that they are propagated : dimnames(A) <- dn <- rep.int(list(paste0("x", seq_len(n))), 2L) (ch.A <- Cholesky(A)) # pivoted, by default str(e.ch.A <- expand2(ch.A, LDL = TRUE), max.level = 2L) str(E.ch.A <- expand2(ch.A, LDL = FALSE), max.level = 2L) ae1 <- function(a, b, ...) all.equal(as(a, "matrix"), as(b, "matrix"), ...) ae2 <- function(a, b, ...) ae1(unname(a), unname(b), ...) ## A ~ P1' L1 D L1' P1 ~ P1' L L' P1 in floating point stopifnot(exprs = { identical(names(e.ch.A), c("P1.", "L1", "D", "L1.", "P1")) identical(names(E.ch.A), c("P1.", "L" , "L." , "P1")) identical(e.ch.A[["P1"]], new("pMatrix", Dim = c(n, n), Dimnames = c(list(NULL), dn[2L]), margin = 2L, perm = invertPerm(ch.A@perm, 0L, 1L))) identical(e.ch.A[["P1."]], t(e.ch.A[["P1"]])) identical(e.ch.A[["L1."]], t(e.ch.A[["L1"]])) identical(E.ch.A[["L." ]], t(E.ch.A[["L" ]])) identical(e.ch.A[["D"]], Diagonal(x = diag(ch.A))) all.equal(E.ch.A[["L"]], with(e.ch.A, L1 %*% sqrt(D))) ae1(A, with(e.ch.A, P1. %*% L1 %*% D %*% L1. %*% P1)) ae1(A, with(E.ch.A, P1. %*% L %*% L. %*% P1)) ae2(A[ch.A@perm + 1L, ch.A@perm + 1L], with(e.ch.A, L1 %*% D %*% L1.)) ae2(A[ch.A@perm + 1L, ch.A@perm + 1L], with(E.ch.A, L %*% L. )) }) ## Factorization handled as factorized matrix ## (in some cases only optionally, depending on arguments) b <- rnorm(n) stopifnot(identical(det(A), det(ch.A, sqrt = FALSE)), identical(solve(A, b), solve(ch.A, b, system = "A"))) u1 <- update(ch.A, A , mult = sqrt(2)) u2 <- update(ch.A, t(M), mult = sqrt(2)) # updating with crossprod(M), not M stopifnot(all.equal(u1, u2, tolerance = 1e-14))
Computes the upper triangular Cholesky factor of an
real, symmetric, positive semidefinite
matrix
, optionally after pivoting.
That is the factor
in
or (equivalently)
where
is a permutation matrix.
Methods for denseMatrix
are built on
LAPACK routines dpstrf
, dpotrf
, and dpptrf
,
The latter two do not permute rows or columns,
so that is an identity matrix.
Methods for sparseMatrix
are built on
CHOLMOD routines cholmod_analyze
and cholmod_factorize_p
.
chol(x, ...) ## S4 method for signature 'dsyMatrix' chol(x, pivot = FALSE, tol = -1, ...) ## S4 method for signature 'dspMatrix' chol(x, ...) ## S4 method for signature 'dsCMatrix' chol(x, pivot = FALSE, ...) ## S4 method for signature 'ddiMatrix' chol(x, ...) ## S4 method for signature 'generalMatrix' chol(x, uplo = "U", ...) ## S4 method for signature 'triangularMatrix' chol(x, uplo = "U", ...)
chol(x, ...) ## S4 method for signature 'dsyMatrix' chol(x, pivot = FALSE, tol = -1, ...) ## S4 method for signature 'dspMatrix' chol(x, ...) ## S4 method for signature 'dsCMatrix' chol(x, pivot = FALSE, ...) ## S4 method for signature 'ddiMatrix' chol(x, ...) ## S4 method for signature 'generalMatrix' chol(x, uplo = "U", ...) ## S4 method for signature 'triangularMatrix' chol(x, uplo = "U", ...)
x |
a finite, symmetric, positive
semidefinite matrix or |
pivot |
a logical indicating if the rows and columns
of |
tol |
a finite numeric tolerance,
used only if |
uplo |
a string, either |
... |
further arguments passed to or from methods. |
For x
inheriting from diagonalMatrix
,
the diagonal result is computed directly and without pivoting,
i.e., bypassing CHOLMOD.
For all other x
, chol(x, pivot = value)
calls
Cholesky(x, perm = value, ...)
under the hood.
If you must know the permutation in addition
to the Cholesky factor
, then call
Cholesky
directly, as the result of chol(x, pivot = TRUE)
specifies
but not
.
A matrix, triangularMatrix
,
or diagonalMatrix
representing
the upper triangular Cholesky factor .
The result is a traditional matrix if
x
is a
traditional matrix, dense if x
is dense, and
sparse if x
is sparse.
The LAPACK source code, including documentation; see https://netlib.org/lapack/double/dpstrf.f, https://netlib.org/lapack/double/dpotrf.f, and https://netlib.org/lapack/double/dpptrf.f.
The CHOLMOD source code; see
https://github.com/DrTimothyAldenDavis/SuiteSparse,
notably the header file ‘CHOLMOD/Include/cholmod.h’
defining cholmod_factor_struct
.
Chen, Y., Davis, T. A., Hager, W. W., & Rajamanickam, S. (2008). Algorithm 887: CHOLMOD, supernodal sparse Cholesky factorization and update/downdate. ACM Transactions on Mathematical Software, 35(3), Article 22, 1-14. doi:10.1145/1391989.1391995
Amestoy, P. R., Davis, T. A., & Duff, I. S. (2004). Algorithm 837: AMD, an approximate minimum degree ordering algorithm. ACM Transactions on Mathematical Software, 17(4), 886-905. doi:10.1145/1024074.1024081
Golub, G. H., & Van Loan, C. F. (2013). Matrix computations (4th ed.). Johns Hopkins University Press. doi:10.56021/9781421407944
The default method from base, chol
,
called for traditional matrices x
.
Generic function Cholesky
, for more flexibility
notably when computing the Cholesky factorization and
not only the factor .
showMethods("chol", inherited = FALSE) set.seed(0) ## ---- Dense ---------------------------------------------------------- ## chol(x, pivot = value) wrapping Cholesky(x, perm = value) selectMethod("chol", "dsyMatrix") ## Except in packed cases where pivoting is not yet available selectMethod("chol", "dspMatrix") ## .... Positive definite .............................................. (A1 <- new("dsyMatrix", Dim = c(2L, 2L), x = c(1, 2, 2, 5))) (R1.nopivot <- chol(A1)) (R1 <- chol(A1, pivot = TRUE)) ## In 2-by-2 cases, we know that the permutation is 1:2 or 2:1, ## even if in general 'chol' does not say ... stopifnot(exprs = { all.equal( A1 , as(crossprod(R1.nopivot), "dsyMatrix")) all.equal(t(A1[2:1, 2:1]), as(crossprod(R1 ), "dsyMatrix")) identical(Cholesky(A1)@perm, 2:1) # because 5 > 1 }) ## .... Positive semidefinite but not positive definite ................ (A2 <- new("dpoMatrix", Dim = c(2L, 2L), x = c(1, 2, 2, 4))) try(R2.nopivot <- chol(A2)) # fails as not positive definite (R2 <- chol(A2, pivot = TRUE)) # returns, with a warning and ... stopifnot(exprs = { all.equal(t(A2[2:1, 2:1]), as(crossprod(R2), "dsyMatrix")) identical(Cholesky(A2)@perm, 2:1) # because 4 > 1 }) ## .... Not positive semidefinite ...................................... (A3 <- new("dsyMatrix", Dim = c(2L, 2L), x = c(1, 2, 2, 3))) try(R3.nopivot <- chol(A3)) # fails as not positive definite (R3 <- chol(A3, pivot = TRUE)) # returns, with a warning and ... ## _Not_ equal: see details and examples in help("Cholesky") all.equal(t(A3[2:1, 2:1]), as(crossprod(R3), "dsyMatrix")) ## ---- Sparse --------------------------------------------------------- ## chol(x, pivot = value) wrapping ## Cholesky(x, perm = value, LDL = FALSE, super = FALSE) selectMethod("chol", "dsCMatrix") ## Except in diagonal cases which are handled "directly" selectMethod("chol", "ddiMatrix") (A4 <- toeplitz(as(c(10, 0, 1, 0, 3), "sparseVector"))) (ch.A4.nopivot <- Cholesky(A4, perm = FALSE, LDL = FALSE, super = FALSE)) (ch.A4 <- Cholesky(A4, perm = TRUE, LDL = FALSE, super = FALSE)) (R4.nopivot <- chol(A4)) (R4 <- chol(A4, pivot = TRUE)) det4 <- det(A4) b4 <- rnorm(5L) x4 <- solve(A4, b4) stopifnot(exprs = { identical(R4.nopivot, expand1(ch.A4.nopivot, "L.")) identical(R4, expand1(ch.A4, "L.")) all.equal(A4, crossprod(R4.nopivot)) all.equal(A4[ch.A4@perm + 1L, ch.A4@perm + 1L], crossprod(R4)) all.equal(diag(R4.nopivot), sqrt(diag(ch.A4.nopivot))) all.equal(diag(R4), sqrt(diag(ch.A4))) all.equal(sqrt(det4), det(R4.nopivot)) all.equal(sqrt(det4), det(R4)) all.equal(det4, det(ch.A4.nopivot, sqrt = FALSE)) all.equal(det4, det(ch.A4, sqrt = FALSE)) all.equal(x4, solve(R4.nopivot, solve(t(R4.nopivot), b4))) all.equal(x4, solve(ch.A4.nopivot, b4)) all.equal(x4, solve(ch.A4, b4)) })
showMethods("chol", inherited = FALSE) set.seed(0) ## ---- Dense ---------------------------------------------------------- ## chol(x, pivot = value) wrapping Cholesky(x, perm = value) selectMethod("chol", "dsyMatrix") ## Except in packed cases where pivoting is not yet available selectMethod("chol", "dspMatrix") ## .... Positive definite .............................................. (A1 <- new("dsyMatrix", Dim = c(2L, 2L), x = c(1, 2, 2, 5))) (R1.nopivot <- chol(A1)) (R1 <- chol(A1, pivot = TRUE)) ## In 2-by-2 cases, we know that the permutation is 1:2 or 2:1, ## even if in general 'chol' does not say ... stopifnot(exprs = { all.equal( A1 , as(crossprod(R1.nopivot), "dsyMatrix")) all.equal(t(A1[2:1, 2:1]), as(crossprod(R1 ), "dsyMatrix")) identical(Cholesky(A1)@perm, 2:1) # because 5 > 1 }) ## .... Positive semidefinite but not positive definite ................ (A2 <- new("dpoMatrix", Dim = c(2L, 2L), x = c(1, 2, 2, 4))) try(R2.nopivot <- chol(A2)) # fails as not positive definite (R2 <- chol(A2, pivot = TRUE)) # returns, with a warning and ... stopifnot(exprs = { all.equal(t(A2[2:1, 2:1]), as(crossprod(R2), "dsyMatrix")) identical(Cholesky(A2)@perm, 2:1) # because 4 > 1 }) ## .... Not positive semidefinite ...................................... (A3 <- new("dsyMatrix", Dim = c(2L, 2L), x = c(1, 2, 2, 3))) try(R3.nopivot <- chol(A3)) # fails as not positive definite (R3 <- chol(A3, pivot = TRUE)) # returns, with a warning and ... ## _Not_ equal: see details and examples in help("Cholesky") all.equal(t(A3[2:1, 2:1]), as(crossprod(R3), "dsyMatrix")) ## ---- Sparse --------------------------------------------------------- ## chol(x, pivot = value) wrapping ## Cholesky(x, perm = value, LDL = FALSE, super = FALSE) selectMethod("chol", "dsCMatrix") ## Except in diagonal cases which are handled "directly" selectMethod("chol", "ddiMatrix") (A4 <- toeplitz(as(c(10, 0, 1, 0, 3), "sparseVector"))) (ch.A4.nopivot <- Cholesky(A4, perm = FALSE, LDL = FALSE, super = FALSE)) (ch.A4 <- Cholesky(A4, perm = TRUE, LDL = FALSE, super = FALSE)) (R4.nopivot <- chol(A4)) (R4 <- chol(A4, pivot = TRUE)) det4 <- det(A4) b4 <- rnorm(5L) x4 <- solve(A4, b4) stopifnot(exprs = { identical(R4.nopivot, expand1(ch.A4.nopivot, "L.")) identical(R4, expand1(ch.A4, "L.")) all.equal(A4, crossprod(R4.nopivot)) all.equal(A4[ch.A4@perm + 1L, ch.A4@perm + 1L], crossprod(R4)) all.equal(diag(R4.nopivot), sqrt(diag(ch.A4.nopivot))) all.equal(diag(R4), sqrt(diag(ch.A4))) all.equal(sqrt(det4), det(R4.nopivot)) all.equal(sqrt(det4), det(R4)) all.equal(det4, det(ch.A4.nopivot, sqrt = FALSE)) all.equal(det4, det(ch.A4, sqrt = FALSE)) all.equal(x4, solve(R4.nopivot, solve(t(R4.nopivot), b4))) all.equal(x4, solve(ch.A4.nopivot, b4)) all.equal(x4, solve(ch.A4, b4)) })
Given formally upper and lower triangular matrices
and
, compute
and
, respectively.
This function can be seen as way to compute the inverse of a
symmetric positive definite matrix given its Cholesky factor.
Equivalently, it can be seen as a way to compute
given the
part of the
QR factorization of
, if
is constrained to have
positive diagonal entries.
chol2inv(x, ...) ## S4 method for signature 'dtrMatrix' chol2inv(x, ...) ## S4 method for signature 'dtCMatrix' chol2inv(x, ...) ## S4 method for signature 'generalMatrix' chol2inv(x, uplo = "U", ...)
chol2inv(x, ...) ## S4 method for signature 'dtrMatrix' chol2inv(x, ...) ## S4 method for signature 'dtCMatrix' chol2inv(x, ...) ## S4 method for signature 'generalMatrix' chol2inv(x, uplo = "U", ...)
x |
a square matrix or |
uplo |
a string, either |
... |
further arguments passed to or from methods. |
A matrix, symmetricMatrix
,
or diagonalMatrix
representing
the inverse of the positive definite matrix whose
Cholesky factor is x
.
The result is a traditional matrix if x
is a
traditional matrix, dense if x
is dense, and
sparse if x
is sparse.
The default method from base, chol2inv
,
called for traditional matrices x
.
Generic function chol
, for computing the upper
triangular Cholesky factor of a symmetric positive
semidefinite matrix.
Generic function solve
, for solving linear systems
and (as a corollary) for computing inverses more generally.
(A <- Matrix(cbind(c(1, 1, 1), c(1, 2, 4), c(1, 4, 16)))) (R <- chol(A)) (L <- t(R)) (R2i <- chol2inv(R)) (L2i <- chol2inv(R)) stopifnot(exprs = { all.equal(R2i, tcrossprod(solve(R))) all.equal(L2i, crossprod(solve(L))) all.equal(as(R2i %*% A, "matrix"), diag(3L)) # the identity all.equal(as(L2i %*% A, "matrix"), diag(3L)) # ditto })
(A <- Matrix(cbind(c(1, 1, 1), c(1, 2, 4), c(1, 4, 16)))) (R <- chol(A)) (L <- t(R)) (R2i <- chol2inv(R)) (L2i <- chol2inv(R)) stopifnot(exprs = { all.equal(R2i, tcrossprod(solve(R))) all.equal(L2i, crossprod(solve(L))) all.equal(as(R2i %*% A, "matrix"), diag(3L)) # the identity all.equal(as(L2i %*% A, "matrix"), diag(3L)) # ditto })
Classes Cholesky
and pCholesky
represent
dense, pivoted Cholesky factorizations of
real, symmetric, positive semidefinite matrices
,
having the general form
or (equivalently)
where
is a permutation matrix,
is a unit lower triangular matrix,
is a non-negative diagonal matrix, and
.
These classes store the entries of the Cholesky factor
or its transpose
in a dense format as
a vector of length
(
Cholesky
) or
(
pCholesky
), the latter
giving the “packed” representation.
Dim
, Dimnames
inherited from virtual class
MatrixFactorization
.
uplo
a string, either "U"
or "L"
,
indicating which triangle (upper or lower) of the factorized
symmetric matrix was used to compute the factorization and
in turn whether x
stores or
.
x
a numeric vector of length n*n
(Cholesky
) or n*(n+1)/2
(pCholesky
),
where n=Dim[1]
, listing the entries of the Cholesky
factor or its transpose
in column-major
order.
perm
a 1-based integer vector of length Dim[1]
specifying the permutation applied to the rows and columns
of the factorized matrix. perm
of length 0 is valid and
equivalent to the identity permutation, implying no pivoting.
Class CholeskyFactorization
, directly.
Class MatrixFactorization
, by class
CholeskyFactorization
, distance 2.
Objects can be generated directly by calls of the form
new("Cholesky", ...)
or new("pCholesky", ...)
,
but they are more typically obtained as the value of
Cholesky(x)
for x
inheriting from
dsyMatrix
or dspMatrix
(often the subclasses of those reserved for positive
semidefinite matrices, namely dpoMatrix
and dppMatrix
).
coerce
signature(from = "Cholesky", to = "dtrMatrix")
:
returns a dtrMatrix
representing
the Cholesky factor or its transpose
;
see ‘Note’.
coerce
signature(from = "pCholesky", to = "dtpMatrix")
:
returns a dtpMatrix
representing
the Cholesky factor or its transpose
;
see ‘Note’.
determinant
signature(from = "p?Cholesky", logarithm = "logical")
:
computes the determinant of the factorized matrix
or its logarithm.
diag
signature(x = "p?Cholesky")
:
returns a numeric vector of length containing the diagonal
elements of
, which are the squared diagonal elements of
.
expand1
signature(x = "p?Cholesky")
:
see expand1-methods
.
expand2
signature(x = "p?Cholesky")
:
see expand2-methods
.
solve
signature(a = "p?Cholesky", b = .)
:
see solve-methods
.
In Matrix < 1.6-0
, class Cholesky
extended
dtrMatrix
and class pCholesky
extended
dtpMatrix
, reflecting the fact that the factor
is indeed a triangular matrix.
Matrix
1.6-0
removed these extensions so that methods
would no longer be inherited from dtrMatrix
and dtpMatrix
.
The availability of such methods gave the wrong impression that
Cholesky
and pCholesky
represent a (singular)
matrix, when in fact they represent an ordered set of matrix factors.
The coercions as(., "dtrMatrix")
and as(., "dtpMatrix")
are provided for users who understand the caveats.
The LAPACK source code, including documentation; see https://netlib.org/lapack/double/dpstrf.f, https://netlib.org/lapack/double/dpotrf.f, and https://netlib.org/lapack/double/dpptrf.f.
Lucas, C. (2004). LAPACK-style codes for level 2 and 3 pivoted Cholesky factorizations. LAPACK Working Note, Number 161. https://www.netlib.org/lapack/lawnspdf/lawn161.pdf
Golub, G. H., & Van Loan, C. F. (2013). Matrix computations (4th ed.). Johns Hopkins University Press. doi:10.56021/9781421407944
Class CHMfactor
for sparse Cholesky factorizations.
Classes dpoMatrix
and dppMatrix
.
Generic functions Cholesky
,
expand1
and expand2
.
showClass("Cholesky") set.seed(1) m <- 30L n <- 6L (A <- crossprod(Matrix(rnorm(m * n), m, n))) ## With dimnames, to see that they are propagated : dimnames(A) <- dn <- rep.int(list(paste0("x", seq_len(n))), 2L) (ch.A <- Cholesky(A)) # pivoted, by default str(e.ch.A <- expand2(ch.A, LDL = TRUE), max.level = 2L) str(E.ch.A <- expand2(ch.A, LDL = FALSE), max.level = 2L) ## Underlying LAPACK representation (m.ch.A <- as(ch.A, "dtrMatrix")) # which is L', not L, because A@uplo == "U" stopifnot(identical(as(m.ch.A, "matrix"), `dim<-`(ch.A@x, ch.A@Dim))) ae1 <- function(a, b, ...) all.equal(as(a, "matrix"), as(b, "matrix"), ...) ae2 <- function(a, b, ...) ae1(unname(a), unname(b), ...) ## A ~ P1' L1 D L1' P1 ~ P1' L L' P1 in floating point stopifnot(exprs = { identical(names(e.ch.A), c("P1.", "L1", "D", "L1.", "P1")) identical(names(E.ch.A), c("P1.", "L" , "L." , "P1")) identical(e.ch.A[["P1"]], new("pMatrix", Dim = c(n, n), Dimnames = c(list(NULL), dn[2L]), margin = 2L, perm = invertPerm(ch.A@perm))) identical(e.ch.A[["P1."]], t(e.ch.A[["P1"]])) identical(e.ch.A[["L1."]], t(e.ch.A[["L1"]])) identical(E.ch.A[["L." ]], t(E.ch.A[["L" ]])) identical(e.ch.A[["D"]], Diagonal(x = diag(ch.A))) all.equal(E.ch.A[["L"]], with(e.ch.A, L1 %*% sqrt(D))) ae1(A, with(e.ch.A, P1. %*% L1 %*% D %*% L1. %*% P1)) ae1(A, with(E.ch.A, P1. %*% L %*% L. %*% P1)) ae2(A[ch.A@perm, ch.A@perm], with(e.ch.A, L1 %*% D %*% L1.)) ae2(A[ch.A@perm, ch.A@perm], with(E.ch.A, L %*% L. )) }) ## Factorization handled as factorized matrix b <- rnorm(n) all.equal(det(A), det(ch.A), tolerance = 0) all.equal(solve(A, b), solve(ch.A, b), tolerance = 0) ## For identical results, we need the _unpivoted_ factorization ## computed by det(A) and solve(A, b) (ch.A.nopivot <- Cholesky(A, perm = FALSE)) stopifnot(identical(det(A), det(ch.A.nopivot)), identical(solve(A, b), solve(ch.A.nopivot, b)))
showClass("Cholesky") set.seed(1) m <- 30L n <- 6L (A <- crossprod(Matrix(rnorm(m * n), m, n))) ## With dimnames, to see that they are propagated : dimnames(A) <- dn <- rep.int(list(paste0("x", seq_len(n))), 2L) (ch.A <- Cholesky(A)) # pivoted, by default str(e.ch.A <- expand2(ch.A, LDL = TRUE), max.level = 2L) str(E.ch.A <- expand2(ch.A, LDL = FALSE), max.level = 2L) ## Underlying LAPACK representation (m.ch.A <- as(ch.A, "dtrMatrix")) # which is L', not L, because A@uplo == "U" stopifnot(identical(as(m.ch.A, "matrix"), `dim<-`(ch.A@x, ch.A@Dim))) ae1 <- function(a, b, ...) all.equal(as(a, "matrix"), as(b, "matrix"), ...) ae2 <- function(a, b, ...) ae1(unname(a), unname(b), ...) ## A ~ P1' L1 D L1' P1 ~ P1' L L' P1 in floating point stopifnot(exprs = { identical(names(e.ch.A), c("P1.", "L1", "D", "L1.", "P1")) identical(names(E.ch.A), c("P1.", "L" , "L." , "P1")) identical(e.ch.A[["P1"]], new("pMatrix", Dim = c(n, n), Dimnames = c(list(NULL), dn[2L]), margin = 2L, perm = invertPerm(ch.A@perm))) identical(e.ch.A[["P1."]], t(e.ch.A[["P1"]])) identical(e.ch.A[["L1."]], t(e.ch.A[["L1"]])) identical(E.ch.A[["L." ]], t(E.ch.A[["L" ]])) identical(e.ch.A[["D"]], Diagonal(x = diag(ch.A))) all.equal(E.ch.A[["L"]], with(e.ch.A, L1 %*% sqrt(D))) ae1(A, with(e.ch.A, P1. %*% L1 %*% D %*% L1. %*% P1)) ae1(A, with(E.ch.A, P1. %*% L %*% L. %*% P1)) ae2(A[ch.A@perm, ch.A@perm], with(e.ch.A, L1 %*% D %*% L1.)) ae2(A[ch.A@perm, ch.A@perm], with(E.ch.A, L %*% L. )) }) ## Factorization handled as factorized matrix b <- rnorm(n) all.equal(det(A), det(ch.A), tolerance = 0) all.equal(solve(A, b), solve(ch.A, b), tolerance = 0) ## For identical results, we need the _unpivoted_ factorization ## computed by det(A) and solve(A, b) (ch.A.nopivot <- Cholesky(A, perm = FALSE)) stopifnot(identical(det(A), det(ch.A.nopivot)), identical(solve(A, b), solve(ch.A.nopivot, b)))
Computes the pivoted Cholesky factorization of an
real, symmetric matrix
,
which has the general form
or (equivalently)
where
is a permutation matrix,
is a unit lower triangular matrix,
is a diagonal matrix, and
.
The second equalities hold only for positive semidefinite
,
for which the diagonal entries of
are non-negative
and
is well-defined.
Methods for denseMatrix
are built on
LAPACK routines dpstrf
, dpotrf
, and dpptrf
.
The latter two do not permute rows or columns,
so that is an identity matrix.
Methods for sparseMatrix
are built on
CHOLMOD routines cholmod_analyze
and cholmod_factorize_p
.
Cholesky(A, ...) ## S4 method for signature 'dsyMatrix' Cholesky(A, perm = TRUE, tol = -1, ...) ## S4 method for signature 'dspMatrix' Cholesky(A, ...) ## S4 method for signature 'dsCMatrix' Cholesky(A, perm = TRUE, LDL = !super, super = FALSE, Imult = 0, ...) ## S4 method for signature 'ddiMatrix' Cholesky(A, ...) ## S4 method for signature 'generalMatrix' Cholesky(A, uplo = "U", ...) ## S4 method for signature 'triangularMatrix' Cholesky(A, uplo = "U", ...) ## S4 method for signature 'matrix' Cholesky(A, uplo = "U", ...)
Cholesky(A, ...) ## S4 method for signature 'dsyMatrix' Cholesky(A, perm = TRUE, tol = -1, ...) ## S4 method for signature 'dspMatrix' Cholesky(A, ...) ## S4 method for signature 'dsCMatrix' Cholesky(A, perm = TRUE, LDL = !super, super = FALSE, Imult = 0, ...) ## S4 method for signature 'ddiMatrix' Cholesky(A, ...) ## S4 method for signature 'generalMatrix' Cholesky(A, uplo = "U", ...) ## S4 method for signature 'triangularMatrix' Cholesky(A, uplo = "U", ...) ## S4 method for signature 'matrix' Cholesky(A, uplo = "U", ...)
A |
a finite, symmetric matrix or
|
perm |
a logical indicating if the rows and columns
of |
tol |
a finite numeric tolerance,
used only if |
LDL |
a logical indicating if the simplicial factorization
should be computed as
|
super |
a logical indicating if the factorization should
use the supernodal algorithm. The alternative is the simplicial
algorithm. Setting |
Imult |
a finite number. The matrix
that is factorized is |
uplo |
a string, either |
... |
further arguments passed to or from methods. |
Note that the result of a call to Cholesky
inherits
from CholeskyFactorization
but not
Matrix
. Users who just want a matrix
should consider using chol
, whose methods are
simple wrappers around Cholesky
returning just the
upper triangular Cholesky factor ,
typically as a
triangularMatrix
.
However, a more principled approach would be to construct
factors as needed from the CholeskyFactorization
object,
e.g., with expand1(x, "L")
, if x
is the
object.
The behaviour of Cholesky(A, perm = TRUE)
for dense A
is somewhat exceptional, in that it expects without checking
that A
is positive semidefinite. By construction, if
is positive semidefinite and the exact algorithm encounters a zero
pivot, then the unfactorized trailing submatrix is the zero matrix,
and there is nothing left to do. Hence when the finite precision
algorithm encounters a pivot less than
tol
, it signals a
warning instead of an error and zeros the trailing submatrix in
order to guarantee that is positive
semidefinite even if
is not. It follows that one way to
test for positive semidefiniteness of
in the event of a
warning is to analyze the error
See the examples and LAPACK Working Note (“LAWN”) 161 for details.
An object representing the factorization, inheriting from
virtual class CholeskyFactorization
.
For a traditional matrix A
, the specific class is
Cholesky
.
For A
inheriting from
unpackedMatrix
,
packedMatrix
, and
sparseMatrix
,
the specific class is
Cholesky
,
pCholesky
, and
dCHMsimpl
or dCHMsuper
,
respectively.
The LAPACK source code, including documentation; see https://netlib.org/lapack/double/dpstrf.f, https://netlib.org/lapack/double/dpotrf.f, and https://netlib.org/lapack/double/dpptrf.f.
The CHOLMOD source code; see
https://github.com/DrTimothyAldenDavis/SuiteSparse,
notably the header file ‘CHOLMOD/Include/cholmod.h’
defining cholmod_factor_struct
.
Lucas, C. (2004). LAPACK-style codes for level 2 and 3 pivoted Cholesky factorizations. LAPACK Working Note, Number 161. https://www.netlib.org/lapack/lawnspdf/lawn161.pdf
Chen, Y., Davis, T. A., Hager, W. W., & Rajamanickam, S. (2008). Algorithm 887: CHOLMOD, supernodal sparse Cholesky factorization and update/downdate. ACM Transactions on Mathematical Software, 35(3), Article 22, 1-14. doi:10.1145/1391989.1391995
Amestoy, P. R., Davis, T. A., & Duff, I. S. (2004). Algorithm 837: AMD, an approximate minimum degree ordering algorithm. ACM Transactions on Mathematical Software, 17(4), 886-905. doi:10.1145/1024074.1024081
Golub, G. H., & Van Loan, C. F. (2013). Matrix computations (4th ed.). Johns Hopkins University Press. doi:10.56021/9781421407944
Classes Cholesky
, pCholesky
,
dCHMsimpl
and dCHMsuper
and their methods.
Classes dpoMatrix
, dppMatrix
,
and dsCMatrix
.
Generic function chol
,
for obtaining the upper triangular Cholesky factor as a
matrix or
Matrix
.
Generic functions expand1
and expand2
,
for constructing matrix factors from the result.
Generic functions BunchKaufman
, Schur
,
lu
, and qr
,
for computing other factorizations.
showMethods("Cholesky", inherited = FALSE) set.seed(0) ## ---- Dense ---------------------------------------------------------- ## .... Positive definite .............................................. n <- 6L (A1 <- crossprod(Matrix(rnorm(n * n), n, n))) (ch.A1.nopivot <- Cholesky(A1, perm = FALSE)) (ch.A1 <- Cholesky(A1)) stopifnot(exprs = { length(ch.A1@perm) == ncol(A1) isPerm(ch.A1@perm) is.unsorted(ch.A1@perm) # typically not the identity permutation length(ch.A1.nopivot@perm) == 0L }) ## A ~ P1' L D L' P1 ~ P1' L L' P1 in floating point str(e.ch.A1 <- expand2(ch.A1, LDL = TRUE), max.level = 2L) str(E.ch.A1 <- expand2(ch.A1, LDL = FALSE), max.level = 2L) stopifnot(exprs = { all.equal(as(A1, "matrix"), as(Reduce(`%*%`, e.ch.A1), "matrix")) all.equal(as(A1, "matrix"), as(Reduce(`%*%`, E.ch.A1), "matrix")) }) ## .... Positive semidefinite but not positive definite ................ A2 <- A1 A2[1L, ] <- A2[, 1L] <- 0 A2 try(Cholesky(A2, perm = FALSE)) # fails as not positive definite ch.A2 <- Cholesky(A2) # returns, with a warning and ... A2.hat <- Reduce(`%*%`, expand2(ch.A2, LDL = FALSE)) norm(A2 - A2.hat, "2") / norm(A2, "2") # 7.670858e-17 ## .... Not positive semidefinite ...................................... A3 <- A1 A3[1L, ] <- A3[, 1L] <- -1 A3 try(Cholesky(A3, perm = FALSE)) # fails as not positive definite ch.A3 <- Cholesky(A3) # returns, with a warning and ... A3.hat <- Reduce(`%*%`, expand2(ch.A3, LDL = FALSE)) norm(A3 - A3.hat, "2") / norm(A3, "2") # 1.781568 ## Indeed, 'A3' is not positive semidefinite, but 'A3.hat' _is_ ch.A3.hat <- Cholesky(A3.hat) A3.hat.hat <- Reduce(`%*%`, expand2(ch.A3.hat, LDL = FALSE)) norm(A3.hat - A3.hat.hat, "2") / norm(A3.hat, "2") # 1.777944e-16 ## ---- Sparse --------------------------------------------------------- ## Really just three cases modulo permutation : ## ## type factorization minors of P1 A P1' ## 1 simplicial P1 A P1' = L1 D L1' nonzero ## 2 simplicial P1 A P1' = L L ' positive ## 3 supernodal P1 A P2' = L L ' positive data(KNex, package = "Matrix") A4 <- crossprod(KNex[["mm"]]) ch.A4 <- list(pivoted = list(simpl1 = Cholesky(A4, perm = TRUE, super = FALSE, LDL = TRUE), simpl0 = Cholesky(A4, perm = TRUE, super = FALSE, LDL = FALSE), super0 = Cholesky(A4, perm = TRUE, super = TRUE )), unpivoted = list(simpl1 = Cholesky(A4, perm = FALSE, super = FALSE, LDL = TRUE), simpl0 = Cholesky(A4, perm = FALSE, super = FALSE, LDL = FALSE), super0 = Cholesky(A4, perm = FALSE, super = TRUE ))) ch.A4 s <- simplify2array rapply2 <- function(object, f, ...) rapply(object, f, , , how = "list", ...) s(rapply2(ch.A4, isLDL)) s(m.ch.A4 <- rapply2(ch.A4, expand1, "L")) # giving L = L1 sqrt(D) ## By design, the pivoted and simplicial factorizations ## are more sparse than the unpivoted and supernodal ones ... s(rapply2(m.ch.A4, object.size)) ## Which is nicely visualized by lattice-based methods for 'image' inm <- c("pivoted", "unpivoted") jnm <- c("simpl1", "simpl0", "super0") for(i in 1:2) for(j in 1:3) print(image(m.ch.A4[[c(i, j)]], main = paste(inm[i], jnm[j])), split = c(j, i, 3L, 2L), more = i * j < 6L) simpl1 <- ch.A4[[c("pivoted", "simpl1")]] stopifnot(exprs = { length(simpl1@perm) == ncol(A4) isPerm(simpl1@perm, 0L) is.unsorted(simpl1@perm) # typically not the identity permutation }) ## One can expand with and without D regardless of isLDL(.), ## but "without" requires L = L1 sqrt(D), which is conditional ## on min(diag(D)) >= 0, hence "with" is the default isLDL(simpl1) stopifnot(min(diag(simpl1)) >= 0) str(e.ch.A4 <- expand2(simpl1, LDL = TRUE), max.level = 2L) # default str(E.ch.A4 <- expand2(simpl1, LDL = FALSE), max.level = 2L) stopifnot(exprs = { all.equal(E.ch.A4[["L" ]], e.ch.A4[["L1" ]] %*% sqrt(e.ch.A4[["D"]])) all.equal(E.ch.A4[["L."]], sqrt(e.ch.A4[["D"]]) %*% e.ch.A4[["L1."]]) all.equal(A4, as(Reduce(`%*%`, e.ch.A4), "symmetricMatrix")) all.equal(A4, as(Reduce(`%*%`, E.ch.A4), "symmetricMatrix")) }) ## The "same" permutation matrix with "alternate" representation ## [i, perm[i]] {margin=1} <-> [invertPerm(perm)[j], j] {margin=2} alt <- function(P) { P@margin <- 1L + !(P@margin - 1L) # 1 <-> 2 P@perm <- invertPerm(P@perm) P } ## Expansions are elegant but inefficient (transposes are redundant) ## hence programmers should consider methods for 'expand1' and 'diag' stopifnot(exprs = { identical(expand1(simpl1, "P1"), alt(e.ch.A4[["P1"]])) identical(expand1(simpl1, "L"), E.ch.A4[["L"]]) identical(Diagonal(x = diag(simpl1)), e.ch.A4[["D"]]) }) ## chol(A, pivot = value) is a simple wrapper around ## Cholesky(A, perm = value, LDL = FALSE, super = FALSE), ## returning L' = sqrt(D) L1' _but_ giving no information ## about the permutation P1 selectMethod("chol", "dsCMatrix") stopifnot(all.equal(chol(A4, pivot = TRUE), E.ch.A4[["L."]])) ## Now a symmetric matrix with positive _and_ negative eigenvalues, ## hence _not_ positive semidefinite A5 <- new("dsCMatrix", Dim = c(7L, 7L), p = c(0:1, 3L, 6:7, 10:11, 15L), i = c(0L, 0:1, 0:3, 2:5, 3:6), x = c(1, 6, 38, 10, 60, 103, -4, 6, -32, -247, -2, -16, -128, -2, -67)) (ev <- eigen(A5, only.values = TRUE)$values) (t.ev <- table(factor(sign(ev), -1:1))) # the matrix "inertia" ch.A5 <- Cholesky(A5) isLDL(ch.A5) (d.A5 <- diag(ch.A5)) # diag(D) is partly negative ## Sylvester's law of inertia holds here, but not in general ## in finite precision arithmetic stopifnot(identical(table(factor(sign(d.A5), -1:1)), t.ev)) try(expand1(ch.A5, "L")) # unable to compute L = L1 sqrt(D) try(expand2(ch.A5, LDL = FALSE)) # ditto try(chol(A5, pivot = TRUE)) # ditto ## The default expansion is "square root free" and still works here str(e.ch.A5 <- expand2(ch.A5, LDL = TRUE), max.level = 2L) stopifnot(all.equal(A5, as(Reduce(`%*%`, e.ch.A5), "symmetricMatrix"))) ## Version of the SuiteSparse library, which includes CHOLMOD Mv <- Matrix.Version() Mv[["suitesparse"]]
showMethods("Cholesky", inherited = FALSE) set.seed(0) ## ---- Dense ---------------------------------------------------------- ## .... Positive definite .............................................. n <- 6L (A1 <- crossprod(Matrix(rnorm(n * n), n, n))) (ch.A1.nopivot <- Cholesky(A1, perm = FALSE)) (ch.A1 <- Cholesky(A1)) stopifnot(exprs = { length(ch.A1@perm) == ncol(A1) isPerm(ch.A1@perm) is.unsorted(ch.A1@perm) # typically not the identity permutation length(ch.A1.nopivot@perm) == 0L }) ## A ~ P1' L D L' P1 ~ P1' L L' P1 in floating point str(e.ch.A1 <- expand2(ch.A1, LDL = TRUE), max.level = 2L) str(E.ch.A1 <- expand2(ch.A1, LDL = FALSE), max.level = 2L) stopifnot(exprs = { all.equal(as(A1, "matrix"), as(Reduce(`%*%`, e.ch.A1), "matrix")) all.equal(as(A1, "matrix"), as(Reduce(`%*%`, E.ch.A1), "matrix")) }) ## .... Positive semidefinite but not positive definite ................ A2 <- A1 A2[1L, ] <- A2[, 1L] <- 0 A2 try(Cholesky(A2, perm = FALSE)) # fails as not positive definite ch.A2 <- Cholesky(A2) # returns, with a warning and ... A2.hat <- Reduce(`%*%`, expand2(ch.A2, LDL = FALSE)) norm(A2 - A2.hat, "2") / norm(A2, "2") # 7.670858e-17 ## .... Not positive semidefinite ...................................... A3 <- A1 A3[1L, ] <- A3[, 1L] <- -1 A3 try(Cholesky(A3, perm = FALSE)) # fails as not positive definite ch.A3 <- Cholesky(A3) # returns, with a warning and ... A3.hat <- Reduce(`%*%`, expand2(ch.A3, LDL = FALSE)) norm(A3 - A3.hat, "2") / norm(A3, "2") # 1.781568 ## Indeed, 'A3' is not positive semidefinite, but 'A3.hat' _is_ ch.A3.hat <- Cholesky(A3.hat) A3.hat.hat <- Reduce(`%*%`, expand2(ch.A3.hat, LDL = FALSE)) norm(A3.hat - A3.hat.hat, "2") / norm(A3.hat, "2") # 1.777944e-16 ## ---- Sparse --------------------------------------------------------- ## Really just three cases modulo permutation : ## ## type factorization minors of P1 A P1' ## 1 simplicial P1 A P1' = L1 D L1' nonzero ## 2 simplicial P1 A P1' = L L ' positive ## 3 supernodal P1 A P2' = L L ' positive data(KNex, package = "Matrix") A4 <- crossprod(KNex[["mm"]]) ch.A4 <- list(pivoted = list(simpl1 = Cholesky(A4, perm = TRUE, super = FALSE, LDL = TRUE), simpl0 = Cholesky(A4, perm = TRUE, super = FALSE, LDL = FALSE), super0 = Cholesky(A4, perm = TRUE, super = TRUE )), unpivoted = list(simpl1 = Cholesky(A4, perm = FALSE, super = FALSE, LDL = TRUE), simpl0 = Cholesky(A4, perm = FALSE, super = FALSE, LDL = FALSE), super0 = Cholesky(A4, perm = FALSE, super = TRUE ))) ch.A4 s <- simplify2array rapply2 <- function(object, f, ...) rapply(object, f, , , how = "list", ...) s(rapply2(ch.A4, isLDL)) s(m.ch.A4 <- rapply2(ch.A4, expand1, "L")) # giving L = L1 sqrt(D) ## By design, the pivoted and simplicial factorizations ## are more sparse than the unpivoted and supernodal ones ... s(rapply2(m.ch.A4, object.size)) ## Which is nicely visualized by lattice-based methods for 'image' inm <- c("pivoted", "unpivoted") jnm <- c("simpl1", "simpl0", "super0") for(i in 1:2) for(j in 1:3) print(image(m.ch.A4[[c(i, j)]], main = paste(inm[i], jnm[j])), split = c(j, i, 3L, 2L), more = i * j < 6L) simpl1 <- ch.A4[[c("pivoted", "simpl1")]] stopifnot(exprs = { length(simpl1@perm) == ncol(A4) isPerm(simpl1@perm, 0L) is.unsorted(simpl1@perm) # typically not the identity permutation }) ## One can expand with and without D regardless of isLDL(.), ## but "without" requires L = L1 sqrt(D), which is conditional ## on min(diag(D)) >= 0, hence "with" is the default isLDL(simpl1) stopifnot(min(diag(simpl1)) >= 0) str(e.ch.A4 <- expand2(simpl1, LDL = TRUE), max.level = 2L) # default str(E.ch.A4 <- expand2(simpl1, LDL = FALSE), max.level = 2L) stopifnot(exprs = { all.equal(E.ch.A4[["L" ]], e.ch.A4[["L1" ]] %*% sqrt(e.ch.A4[["D"]])) all.equal(E.ch.A4[["L."]], sqrt(e.ch.A4[["D"]]) %*% e.ch.A4[["L1."]]) all.equal(A4, as(Reduce(`%*%`, e.ch.A4), "symmetricMatrix")) all.equal(A4, as(Reduce(`%*%`, E.ch.A4), "symmetricMatrix")) }) ## The "same" permutation matrix with "alternate" representation ## [i, perm[i]] {margin=1} <-> [invertPerm(perm)[j], j] {margin=2} alt <- function(P) { P@margin <- 1L + !(P@margin - 1L) # 1 <-> 2 P@perm <- invertPerm(P@perm) P } ## Expansions are elegant but inefficient (transposes are redundant) ## hence programmers should consider methods for 'expand1' and 'diag' stopifnot(exprs = { identical(expand1(simpl1, "P1"), alt(e.ch.A4[["P1"]])) identical(expand1(simpl1, "L"), E.ch.A4[["L"]]) identical(Diagonal(x = diag(simpl1)), e.ch.A4[["D"]]) }) ## chol(A, pivot = value) is a simple wrapper around ## Cholesky(A, perm = value, LDL = FALSE, super = FALSE), ## returning L' = sqrt(D) L1' _but_ giving no information ## about the permutation P1 selectMethod("chol", "dsCMatrix") stopifnot(all.equal(chol(A4, pivot = TRUE), E.ch.A4[["L."]])) ## Now a symmetric matrix with positive _and_ negative eigenvalues, ## hence _not_ positive semidefinite A5 <- new("dsCMatrix", Dim = c(7L, 7L), p = c(0:1, 3L, 6:7, 10:11, 15L), i = c(0L, 0:1, 0:3, 2:5, 3:6), x = c(1, 6, 38, 10, 60, 103, -4, 6, -32, -247, -2, -16, -128, -2, -67)) (ev <- eigen(A5, only.values = TRUE)$values) (t.ev <- table(factor(sign(ev), -1:1))) # the matrix "inertia" ch.A5 <- Cholesky(A5) isLDL(ch.A5) (d.A5 <- diag(ch.A5)) # diag(D) is partly negative ## Sylvester's law of inertia holds here, but not in general ## in finite precision arithmetic stopifnot(identical(table(factor(sign(d.A5), -1:1)), t.ev)) try(expand1(ch.A5, "L")) # unable to compute L = L1 sqrt(D) try(expand2(ch.A5, LDL = FALSE)) # ditto try(chol(A5, pivot = TRUE)) # ditto ## The default expansion is "square root free" and still works here str(e.ch.A5 <- expand2(ch.A5, LDL = TRUE), max.level = 2L) stopifnot(all.equal(A5, as(Reduce(`%*%`, e.ch.A5), "symmetricMatrix"))) ## Version of the SuiteSparse library, which includes CHOLMOD Mv <- Matrix.Version() Mv[["suitesparse"]]
Since 2005, package Matrix has supported coercions to and
from class graph
from package
graph.
Since 2013, this functionality has been exposed via functions
T2graph
and graph2T
, which, unlike methods for
as(from, "<Class>")
, support optional arguments.
graph2T(from, use.weights = ) T2graph(from, need.uniq = !isUniqueT(from), edgemode = NULL)
graph2T(from, use.weights = ) T2graph(from, need.uniq = !isUniqueT(from), edgemode = NULL)
from |
for |
use.weights |
logical indicating if weights should be used, i.e.,
equivalently the result will be numeric, i.e. of class
|
need.uniq |
a logical indicating if |
edgemode |
one of |
For graph2T()
, a sparse matrix inheriting from
"TsparseMatrix"
.
For T2graph()
an R object of class "graph"
.
Package igraph, which provides similar coercions
to and from its class igraph
via functions
graph_from_adjacency_matrix
and as_adjacency_matrix
.
if(requireNamespace("graph")) { n4 <- LETTERS[1:4]; dns <- list(n4,n4) show(a1 <- sparseMatrix(i= c(1:4), j=c(2:4,1), x = 2, dimnames=dns)) show(g1 <- as(a1, "graph")) # directed unlist(graph::edgeWeights(g1)) # all '2' show(a2 <- sparseMatrix(i= c(1:4,4), j=c(2:4,1:2), x = TRUE, dimnames=dns)) show(g2 <- as(a2, "graph")) # directed # now if you want it undirected: show(g3 <- T2graph(as(a2,"TsparseMatrix"), edgemode="undirected")) show(m3 <- as(g3,"Matrix")) show( graph2T(g3) ) # a "pattern Matrix" (nsTMatrix) a. <- sparseMatrix(i=4:1, j=1:4, dimnames=list(n4, n4), repr="T") # no 'x' show(a.) # "ngTMatrix" show(g. <- as(a., "graph")) }
if(requireNamespace("graph")) { n4 <- LETTERS[1:4]; dns <- list(n4,n4) show(a1 <- sparseMatrix(i= c(1:4), j=c(2:4,1), x = 2, dimnames=dns)) show(g1 <- as(a1, "graph")) # directed unlist(graph::edgeWeights(g1)) # all '2' show(a2 <- sparseMatrix(i= c(1:4,4), j=c(2:4,1:2), x = TRUE, dimnames=dns)) show(g2 <- as(a2, "graph")) # directed # now if you want it undirected: show(g3 <- T2graph(as(a2,"TsparseMatrix"), edgemode="undirected")) show(m3 <- as(g3,"Matrix")) show( graph2T(g3) ) # a "pattern Matrix" (nsTMatrix) a. <- sparseMatrix(i=4:1, j=1:4, dimnames=list(n4, n4), repr="T") # no 'x' show(a.) # "ngTMatrix" show(g. <- as(a., "graph")) }
Methods for coercion from and to sparse matrices from package SparseM
are provided here, for ease of porting functionality to the
Matrix package, and comparing functionality of the two
packages. All these work via the usual as(., "<class>")
coercion,
as(from, Class)
...
...
...
...
...
...
...
...
...
...
...
...
The documentation in CRAN package SparseM, such as
SparseM.ontology
, and one important class,
matrix.csr
.
Form row and column sums and means for
objects, for sparseMatrix
the result may
optionally be sparse (sparseVector
), too.
Row or column names are kept respectively as for base matrices
and colSums
methods, when the result is
numeric
vector.
colSums(x, na.rm = FALSE, dims = 1L, ...) rowSums(x, na.rm = FALSE, dims = 1L, ...) colMeans(x, na.rm = FALSE, dims = 1L, ...) rowMeans(x, na.rm = FALSE, dims = 1L, ...) ## S4 method for signature 'CsparseMatrix' colSums(x, na.rm = FALSE, dims = 1L, sparseResult = FALSE, ...) ## S4 method for signature 'CsparseMatrix' rowSums(x, na.rm = FALSE, dims = 1L, sparseResult = FALSE, ...) ## S4 method for signature 'CsparseMatrix' colMeans(x, na.rm = FALSE, dims = 1L, sparseResult = FALSE, ...) ## S4 method for signature 'CsparseMatrix' rowMeans(x, na.rm = FALSE, dims = 1L, sparseResult = FALSE, ...)
colSums(x, na.rm = FALSE, dims = 1L, ...) rowSums(x, na.rm = FALSE, dims = 1L, ...) colMeans(x, na.rm = FALSE, dims = 1L, ...) rowMeans(x, na.rm = FALSE, dims = 1L, ...) ## S4 method for signature 'CsparseMatrix' colSums(x, na.rm = FALSE, dims = 1L, sparseResult = FALSE, ...) ## S4 method for signature 'CsparseMatrix' rowSums(x, na.rm = FALSE, dims = 1L, sparseResult = FALSE, ...) ## S4 method for signature 'CsparseMatrix' colMeans(x, na.rm = FALSE, dims = 1L, sparseResult = FALSE, ...) ## S4 method for signature 'CsparseMatrix' rowMeans(x, na.rm = FALSE, dims = 1L, sparseResult = FALSE, ...)
x |
a Matrix, i.e., inheriting from |
na.rm |
logical. Should missing values (including |
dims |
completely ignored by the |
... |
potentially further arguments, for method |
sparseResult |
logical indicating if the result should be sparse,
i.e., inheriting from class |
returns a numeric vector if sparseResult
is FALSE
as per
default. Otherwise, returns a sparseVector
.
dimnames(x)
are only kept (as names(v)
)
when the resulting v
is numeric
, since
sparseVector
s do not have names.
colSums
and the
sparseVector
classes.
(M <- bdiag(Diagonal(2), matrix(1:3, 3,4), diag(3:2))) # 7 x 8 colSums(M) d <- Diagonal(10, c(0,0,10,0,2,rep(0,5))) MM <- kronecker(d, M) dim(MM) # 70 80 length(MM@x) # 160, but many are '0' ; drop those: MM <- drop0(MM) length(MM@x) # 32 cm <- colSums(MM) (scm <- colSums(MM, sparseResult = TRUE)) stopifnot(is(scm, "sparseVector"), identical(cm, as.numeric(scm))) rowSums (MM, sparseResult = TRUE) # 14 of 70 are not zero colMeans(MM, sparseResult = TRUE) # 16 of 80 are not zero ## Since we have no 'NA's, these two are equivalent : stopifnot(identical(rowMeans(MM, sparseResult = TRUE), rowMeans(MM, sparseResult = TRUE, na.rm = TRUE)), rowMeans(Diagonal(16)) == 1/16, colSums(Diagonal(7)) == 1) ## dimnames(x) --> names( <value> ) : dimnames(M) <- list(paste0("r", 1:7), paste0("V",1:8)) M colSums(M) rowMeans(M) ## Assertions : stopifnot(exprs = { all.equal(colSums(M), structure(c(1,1,6,6,6,6,3,2), names = colnames(M))) all.equal(rowMeans(M), structure(c(1,1,4,8,12,3,2)/8, names = paste0("r", 1:7))) })
(M <- bdiag(Diagonal(2), matrix(1:3, 3,4), diag(3:2))) # 7 x 8 colSums(M) d <- Diagonal(10, c(0,0,10,0,2,rep(0,5))) MM <- kronecker(d, M) dim(MM) # 70 80 length(MM@x) # 160, but many are '0' ; drop those: MM <- drop0(MM) length(MM@x) # 32 cm <- colSums(MM) (scm <- colSums(MM, sparseResult = TRUE)) stopifnot(is(scm, "sparseVector"), identical(cm, as.numeric(scm))) rowSums (MM, sparseResult = TRUE) # 14 of 70 are not zero colMeans(MM, sparseResult = TRUE) # 16 of 80 are not zero ## Since we have no 'NA's, these two are equivalent : stopifnot(identical(rowMeans(MM, sparseResult = TRUE), rowMeans(MM, sparseResult = TRUE, na.rm = TRUE)), rowMeans(Diagonal(16)) == 1/16, colSums(Diagonal(7)) == 1) ## dimnames(x) --> names( <value> ) : dimnames(M) <- list(paste0("r", 1:7), paste0("V",1:8)) M colSums(M) rowMeans(M) ## Assertions : stopifnot(exprs = { all.equal(colSums(M), structure(c(1,1,6,6,6,6,3,2), names = colnames(M))) all.equal(rowMeans(M), structure(c(1,1,4,8,12,3,2)/8, names = paste0("r", 1:7))) })
“Estimate”, i.e. compute approximately the CONDition number of
a (potentially large, often sparse) matrix A
.
It works by apply a fast randomized approximation of the 1-norm,
norm(A,"1")
, through onenormest(.)
.
condest(A, t = min(n, 5), normA = norm(A, "1"), silent = FALSE, quiet = TRUE) onenormest(A, t = min(n, 5), A.x, At.x, n, silent = FALSE, quiet = silent, iter.max = 10, eps = 4 * .Machine$double.eps)
condest(A, t = min(n, 5), normA = norm(A, "1"), silent = FALSE, quiet = TRUE) onenormest(A, t = min(n, 5), A.x, At.x, n, silent = FALSE, quiet = silent, iter.max = 10, eps = 4 * .Machine$double.eps)
A |
a square matrix, optional for |
t |
number of columns to use in the iterations. |
normA |
number; (an estimate of) the 1-norm of |
silent |
logical indicating if warning and (by default) convergence messages should be displayed. |
quiet |
logical indicating if convergence messages should be displayed. |
A.x , At.x
|
when |
n |
|
iter.max |
maximal number of iterations for the 1-norm estimator. |
eps |
the relative change that is deemed irrelevant. |
condest()
calls lu(A)
, and subsequently
onenormest(A.x = , At.x = )
to compute an approximate norm of
the inverse of A
, , in a way which
keeps using sparse matrices efficiently when
A
is sparse.
Note that onenormest()
uses random vectors and hence
both functions' results are random, i.e., depend on the random
seed, see, e.g., set.seed()
.
Both functions return a list
;
condest()
with components,
est |
a number |
v |
the maximal |
The function onenormest()
returns a list with components,
est |
a number |
v |
0-1 integer vector length |
w |
numeric vector, the largest |
iter |
the number of iterations used. |
This is based on octave's condest()
and
onenormest()
implementations with original author
Jason Riedy, U Berkeley; translation to R and
adaption by Martin Maechler.
Nicholas J. Higham and Françoise Tisseur (2000). A Block Algorithm for Matrix 1-Norm Estimation, with an Application to 1-Norm Pseudospectra. SIAM J. Matrix Anal. Appl. 21, 4, 1185–1201.
William W. Hager (1984). Condition Estimates. SIAM J. Sci. Stat. Comput. 5, 311–316.
data(KNex, package = "Matrix") mtm <- with(KNex, crossprod(mm)) system.time(ce <- condest(mtm)) sum(abs(ce$v)) ## || v ||_1 == 1 ## Prove that || A v || = || A || / est (as ||v|| = 1): stopifnot(all.equal(norm(mtm %*% ce$v), norm(mtm) / ce$est)) ## reciprocal 1 / ce$est system.time(rc <- rcond(mtm)) # takes ca 3 x longer rc all.equal(rc, 1/ce$est) # TRUE -- the approximation was good one <- onenormest(mtm) str(one) ## est = 12.3 ## the maximal column: which(one$v == 1) # mostly 4, rarely 1, depending on random seed
data(KNex, package = "Matrix") mtm <- with(KNex, crossprod(mm)) system.time(ce <- condest(mtm)) sum(abs(ce$v)) ## || v ||_1 == 1 ## Prove that || A v || = || A || / est (as ||v|| = 1): stopifnot(all.equal(norm(mtm %*% ce$v), norm(mtm) / ce$est)) ## reciprocal 1 / ce$est system.time(rc <- rcond(mtm)) # takes ca 3 x longer rc all.equal(rc, 1/ce$est) # TRUE -- the approximation was good one <- onenormest(mtm) str(one) ## est = 12.3 ## the maximal column: which(one$v == 1) # mostly 4, rarely 1, depending on random seed
The "CsparseMatrix"
class is the virtual class of
all sparse matrices coded in sorted compressed column-oriented form.
Since it is a virtual class, no objects may be created from it. See
showClass("CsparseMatrix")
for its subclasses.
i
:Object of class "integer"
of length nnzero
(number of non-zero elements). These are the 0-based row numbers for
each non-zero element in the matrix, i.e., i
must be in
0:(nrow(.)-1)
.
p
:integer
vector for providing pointers, one
for each column, to the initial (zero-based) index of elements in
the column. .@p
is of length ncol(.) + 1
, with
p[1] == 0
and p[length(p)] == nnzero
, such that in
fact, diff(.@p)
are the number of non-zero elements for
each column.
In other words, m@p[1:ncol(m)]
contains the indices of
those elements in m@x
that are the first elements in the
respective column of m
.
Dim
, Dimnames
:inherited from
the superclass, see the sparseMatrix
class.
Class "sparseMatrix"
, directly.
Class "Matrix"
, by class "sparseMatrix"
.
matrix products %*%
,
crossprod()
and tcrossprod()
,
several solve
methods,
and other matrix methods available:
signature(e1 = "CsparseMatrix", e2 = "numeric")
: ...
signature(e1 = "numeric", e2 = "CsparseMatrix")
: ...
signature(x = "CsparseMatrix")
: ...
signature(x = "CsparseMatrix")
: ...
signature(e1 = "CsparseMatrix", e2 = "numeric")
: ...
signature(e1 = "numeric", e2 = "CsparseMatrix")
: ...
signature(e1 = "CsparseMatrix", e2 = "numeric")
: ...
signature(e1 = "numeric", e2 = "CsparseMatrix")
: ...
signature(from = "CsparseMatrix", to = "TsparseMatrix")
: ...
signature(from = "CsparseMatrix", to = "denseMatrix")
: ...
signature(from = "CsparseMatrix", to = "matrix")
: ...
signature(from = "TsparseMatrix", to = "CsparseMatrix")
: ...
signature(from = "denseMatrix", to = "CsparseMatrix")
: ...
signature(x = "CsparseMatrix")
: ...
signature(x = "CsparseMatrix")
: ...
signature(x = "CsparseMatrix")
: ...
signature(x = "CsparseMatrix")
: ...
signature(x = "CsparseMatrix")
: ...
signature(x = "CsparseMatrix")
: ...
signature(x = "CsparseMatrix")
: ...
All classes extending CsparseMatrix
have a common validity
(see validObject
) check function. That function
additionally checks the i
slot for each column to contain
increasing row numbers.
In earlier versions of Matrix (<= 0.999375-16
),
validObject
automatically re-sorted the entries when
necessary, and hence new()
calls with somewhat permuted
i
and x
slots worked, as new(...)
(with slot arguments) automatically checks the validity.
Now, you have to use sparseMatrix
to achieve the same
functionality or know how to use .validateCsparse()
to do so.
colSums
, kronecker
, and other such methods
with own help pages.
Further, the super class of CsparseMatrix
,
sparseMatrix
, and, e.g.,
class dgCMatrix
for the links to other classes.
getClass("CsparseMatrix") ## The common validity check function (based on C code): getValidity(getClass("CsparseMatrix"))
getClass("CsparseMatrix") ## The common validity check function (based on C code): getValidity(getClass("CsparseMatrix"))
This is the virtual class of all dense numeric (i.e., double, hence “ddense”) S4 matrices.
Its most important subclass is the dgeMatrix
class.
Class "dMatrix"
directly;
class "Matrix"
, by the above.
the same slots at its subclass dgeMatrix
, see
there.
Most methods are implemented via as(*, "generalMatrix")
and are
mainly used as “fallbacks” when the subclass doesn't need its
own specialized method.
Use showMethods(class = "ddenseMatrix", where =
"package:Matrix")
for an overview.
The virtual classes Matrix
,
dMatrix
, and dsparseMatrix
.
showClass("ddenseMatrix") showMethods(class = "ddenseMatrix", where = "package:Matrix")
showClass("ddenseMatrix") showMethods(class = "ddenseMatrix", where = "package:Matrix")
The class "ddiMatrix"
of numerical diagonal matrices.
Note that diagonal matrices now extend sparseMatrix
, whereas
they did extend dense matrices earlier.
Objects can be created by calls of the form new("ddiMatrix", ...)
but typically rather via Diagonal
.
x
:numeric vector. For an
matrix, the
x
slot is of length or
0
,
depending on the diag
slot:
diag
:"character"
string, either "U"
or
"N"
where "U"
denotes unit-diagonal, i.e., identity
matrices.
Dim
,Dimnames
:matrix dimension and
dimnames
, see the Matrix
class
description.
Class "diagonalMatrix"
, directly.
Class "dMatrix"
, directly.
Class "sparseMatrix"
, indirectly, see
showClass("ddiMatrix")
.
signature(x = "ddiMatrix", y = "ddiMatrix")
: ...
Class diagonalMatrix
and function Diagonal
.
(d2 <- Diagonal(x = c(10,1))) str(d2) ## slightly larger in internal size: str(as(d2, "sparseMatrix")) M <- Matrix(cbind(1,2:4)) M %*% d2 #> `fast' multiplication chol(d2) # trivial stopifnot(is(cd2 <- chol(d2), "ddiMatrix"), all.equal(cd2@x, c(sqrt(10),1)))
(d2 <- Diagonal(x = c(10,1))) str(d2) ## slightly larger in internal size: str(as(d2, "sparseMatrix")) M <- Matrix(cbind(1,2:4)) M %*% d2 #> `fast' multiplication chol(d2) # trivial stopifnot(is(cd2 <- chol(d2), "ddiMatrix"), all.equal(cd2@x, c(sqrt(10),1)))
denseLU
is the class of dense, row-pivoted LU factorizations
of real matrices
,
having the general form
or (equivalently)
where
is an
permutation matrix,
is an
unit lower trapezoidal matrix, and
is a
upper trapezoidal matrix. If
, then the factors
and
are triangular.
Dim
, Dimnames
inherited from virtual class
MatrixFactorization
.
x
a numeric vector of length prod(Dim)
storing
the triangular and
factors together in a packed
format. The details of the representation are specified by the
manual for LAPACK routine
dgetrf
.
perm
an integer vector of length min(Dim)
specifying the permutation as a product of
transpositions. The corresponding permutation vector can
be obtained as
asPerm(perm)
.
Class LU
, directly.
Class MatrixFactorization
, by class
LU
, distance 2.
Objects can be generated directly by calls of the form
new("denseLU", ...)
, but they are more typically obtained
as the value of lu(x)
for x
inheriting from
denseMatrix
(often dgeMatrix
).
coerce
signature(from = "denseLU", to = "dgeMatrix")
:
returns a dgeMatrix
with the dimensions
of the factorized matrix , equal to
below the
diagonal and equal to
on and above the diagonal.
determinant
signature(from = "denseLU", logarithm = "logical")
:
computes the determinant of the factorized matrix
or its logarithm.
expand
signature(x = "denseLU")
:
see expand-methods
.
expand1
signature(x = "denseLU")
:
see expand1-methods
.
expand2
signature(x = "denseLU")
:
see expand2-methods
.
solve
signature(a = "denseLU", b = "missing")
:
see solve-methods
.
The LAPACK source code, including documentation; see https://netlib.org/lapack/double/dgetrf.f.
Golub, G. H., & Van Loan, C. F. (2013). Matrix computations (4th ed.). Johns Hopkins University Press. doi:10.56021/9781421407944
Class sparseLU
for sparse LU factorizations.
Class dgeMatrix
.
Generic functions lu
,
expand1
and expand2
.
showClass("denseLU") set.seed(1) n <- 3L (A <- Matrix(round(rnorm(n * n), 2L), n, n)) ## With dimnames, to see that they are propagated : dimnames(A) <- dn <- list(paste0("r", seq_len(n)), paste0("c", seq_len(n))) (lu.A <- lu(A)) str(e.lu.A <- expand2(lu.A), max.level = 2L) ## Underlying LAPACK representation (m.lu.A <- as(lu.A, "dgeMatrix")) # which is L and U interlaced stopifnot(identical(as(m.lu.A, "matrix"), `dim<-`(lu.A@x, lu.A@Dim))) ae1 <- function(a, b, ...) all.equal(as(a, "matrix"), as(b, "matrix"), ...) ae2 <- function(a, b, ...) ae1(unname(a), unname(b), ...) ## A ~ P1' L U in floating point stopifnot(exprs = { identical(names(e.lu.A), c("P1.", "L", "U")) identical(e.lu.A[["P1."]], new( "pMatrix", Dim = c(n, n), Dimnames = c(dn[1L], list(NULL)), margin = 1L, perm = invertPerm(asPerm(lu.A@perm)))) identical(e.lu.A[["L"]], new("dtrMatrix", Dim = c(n, n), Dimnames = list(NULL, NULL), uplo = "L", diag = "U", x = lu.A@x)) identical(e.lu.A[["U"]], new("dtrMatrix", Dim = c(n, n), Dimnames = c(list(NULL), dn[2L]), uplo = "U", diag = "N", x = lu.A@x)) ae1(A, with(e.lu.A, P1. %*% L %*% U)) ae2(A[asPerm(lu.A@perm), ], with(e.lu.A, L %*% U)) }) ## Factorization handled as factorized matrix b <- rnorm(n) stopifnot(identical(det(A), det(lu.A)), identical(solve(A, b), solve(lu.A, b)))
showClass("denseLU") set.seed(1) n <- 3L (A <- Matrix(round(rnorm(n * n), 2L), n, n)) ## With dimnames, to see that they are propagated : dimnames(A) <- dn <- list(paste0("r", seq_len(n)), paste0("c", seq_len(n))) (lu.A <- lu(A)) str(e.lu.A <- expand2(lu.A), max.level = 2L) ## Underlying LAPACK representation (m.lu.A <- as(lu.A, "dgeMatrix")) # which is L and U interlaced stopifnot(identical(as(m.lu.A, "matrix"), `dim<-`(lu.A@x, lu.A@Dim))) ae1 <- function(a, b, ...) all.equal(as(a, "matrix"), as(b, "matrix"), ...) ae2 <- function(a, b, ...) ae1(unname(a), unname(b), ...) ## A ~ P1' L U in floating point stopifnot(exprs = { identical(names(e.lu.A), c("P1.", "L", "U")) identical(e.lu.A[["P1."]], new( "pMatrix", Dim = c(n, n), Dimnames = c(dn[1L], list(NULL)), margin = 1L, perm = invertPerm(asPerm(lu.A@perm)))) identical(e.lu.A[["L"]], new("dtrMatrix", Dim = c(n, n), Dimnames = list(NULL, NULL), uplo = "L", diag = "U", x = lu.A@x)) identical(e.lu.A[["U"]], new("dtrMatrix", Dim = c(n, n), Dimnames = c(list(NULL), dn[2L]), uplo = "U", diag = "N", x = lu.A@x)) ae1(A, with(e.lu.A, P1. %*% L %*% U)) ae2(A[asPerm(lu.A@perm), ], with(e.lu.A, L %*% U)) }) ## Factorization handled as factorized matrix b <- rnorm(n) stopifnot(identical(det(A), det(lu.A)), identical(solve(A, b), solve(lu.A, b)))
This is the virtual class of all dense (S4) matrices.
It partitions into two subclasses
packedMatrix
and
unpackedMatrix
.
Alternatively into the (currently) three subclasses
ddenseMatrix
,
ldenseMatrix
, and
ndenseMatrix
.
denseMatrix
is (hence) the direct superclass of these () classes.
class "Matrix"
directly.
exactly those of its superclass "Matrix"
, i.e.,
"Dim"
and "Dimnames"
.
Use showMethods(class = "denseMatrix", where =
"package:Matrix")
for an overview of methods.
Extraction ("["
) methods,
see [-methods
.
colSums
, kronecker
, and other such methods
with own help pages.
Its superclass Matrix
, and main subclasses,
ddenseMatrix
and sparseMatrix
.
showClass("denseMatrix")
showClass("denseMatrix")
The dgCMatrix
class is a class of sparse numeric
matrices in the compressed, sparse, column-oriented format. In this
implementation the non-zero elements in the columns are sorted into
increasing row order. dgCMatrix
is the
“standard” class for sparse numeric matrices in the
Matrix package.
Objects can be created by calls of the form new("dgCMatrix",
...)
, more typically via as(*, "CsparseMatrix")
or similar.
Often however, more easily via Matrix(*, sparse = TRUE)
,
or most efficiently via sparseMatrix()
.
x
:Object of class "numeric"
- the non-zero
elements of the matrix.
all other slots are inherited from the superclass
"CsparseMatrix"
.
Matrix products (e.g., crossprod-methods), and (among other)
signature(from = "matrix", to = "dgCMatrix")
signature(x = "dgCMatrix")
: returns the diagonal
of x
signature(x = "dgCMatrix")
: returns the dimensions
of x
signature(x = "dgCMatrix")
: plots an image of
x
using the levelplot
function
signature(a = "dgCMatrix", b = "...")
:
see solve-methods
, notably the extra argument
sparse
.
signature(x = "dgCMatrix")
: computes the LU
decomposition of a square dgCMatrix
object
Classes dsCMatrix
,
dtCMatrix
, lu
(m <- Matrix(c(0,0,2:0), 3,5)) str(m) m[,1]
(m <- Matrix(c(0,0,2:0), 3,5)) str(m) m[,1]
A general numeric dense matrix in the S4 Matrix
representation. dgeMatrix
is the “standard”
class for dense numeric matrices in the Matrix package.
Objects can be created by calls of the form new("dgeMatrix", ...)
or, more commonly, by coercion from the Matrix
class (see
Matrix) or by Matrix(..)
.
x
:Object of class "numeric"
- the numeric
values contained in the matrix, in column-major order.
Dim
:Object of class "integer"
- the dimensions
of the matrix - must be an integer vector with exactly two
non-negative values.
Dimnames
:a list of length two - inherited from class
Matrix
.
factors
:Object of class "list"
- a list
of factorizations of the matrix.
The are group methods (see, e.g., Arith
)
signature(e1 = "dgeMatrix", e2 = "dgeMatrix")
: ...
signature(e1 = "dgeMatrix", e2 = "numeric")
: ...
signature(e1 = "numeric", e2 = "dgeMatrix")
: ...
signature(x = "dgeMatrix")
: ...
signature(x = "dgeMatrix", digits = "numeric")
: ...
matrix products %*%
,
crossprod()
and tcrossprod()
,
several solve
methods,
and other matrix methods available:
signature(x = "dgeMatrix", vectors = "logical")
: ...
signature(x = "dgeMatrix", vectors = "missing")
: ...
signature(x = "dgeMatrix")
: see chol
.
signature(x = "dgeMatrix")
: columnwise means (averages)
signature(x = "dgeMatrix")
: columnwise sums
signature(x = "dgeMatrix")
: ...
signature(x = "dgeMatrix")
: ...
signature(x = "dgeMatrix")
: ...
signature(x = "dgeMatrix", only.values= "logical")
: ...
signature(x = "dgeMatrix", only.values= "missing")
: ...
signature(x = "dgeMatrix", type = "character")
: ...
signature(x = "dgeMatrix", type = "missing")
: ...
signature(x = "dgeMatrix", norm = "character")
or norm = "missing"
:
the reciprocal condition number, rcond()
.
signature(x = "dgeMatrix")
: rowwise means (averages)
signature(x = "dgeMatrix")
: rowwise sums
signature(x = "dgeMatrix")
: matrix transpose
Classes Matrix
,
dtrMatrix
, and dsyMatrix
.
The dgRMatrix
class is a class of sparse numeric
matrices in the compressed, sparse, row-oriented format. In this
implementation the non-zero elements in the rows are sorted into
increasing column order.
Note: The column-oriented sparse classes, e.g.,
dgCMatrix
, are preferred and better supported in
the Matrix package.
Objects can be created by calls of the form new("dgRMatrix", ...)
.
j
:Object of class "integer"
of length nnzero
(number of non-zero elements). These are the column numbers for
each non-zero element in the matrix.
p
:Object of class "integer"
of pointers, one
for each row, to the initial (zero-based) index of elements in
the row.
x
:Object of class "numeric"
- the non-zero
elements of the matrix.
Dim
:Object of class "integer"
- the dimensions
of the matrix.
signature(x = "dgRMatrix")
: returns the diagonal
of x
signature(x = "dgRMatrix")
: returns the dimensions
of x
signature(x = "dgRMatrix")
: plots an image of
x
using the levelplot
function
the RsparseMatrix
class, the virtual class of all
sparse compressed row-oriented matrices, with its methods.
The dgCMatrix
class (column compressed
sparse) is really preferred.
The "dgTMatrix"
class is the class of sparse
matrices stored as (possibly redundant) triplets. The internal
representation is not at all unique, contrary to the one for class
dgCMatrix
.
Objects can be created by calls of the form
new("dgTMatrix", ...)
, but more typically via
spMatrix()
or sparseMatrix(*, repr = "T")
.
i
:integer
row indices of non-zero
entries in 0-base, i.e., must be in 0:(nrow(.)-1)
.
j
:integer
column indices of non-zero
entries. Must be the same length as slot i
and
0-based as well, i.e., in 0:(ncol(.)-1)
.
x
:numeric
vector - the (non-zero)
entry at position (i,j)
. Must be the same length as slot
i
. If an index pair occurs more than once, the corresponding
values of slot x
are added to form the element of the matrix.
Dim
:Object of class "integer"
of length 2 -
the dimensions of the matrix.
signature(e1 = "dgTMatrix", e2 = "dgTMatrix")
signature(x = "dgTMatrix")
: plots an image of
x
using the levelplot
function
signature(x = "dgTMatrix")
: returns the transpose of
x
Triplet matrices are a convenient form in which to construct sparse
matrices after which they can be coerced to
dgCMatrix
objects.
Note that both new(.)
and spMatrix
constructors
for "dgTMatrix"
(and other "TsparseMatrix"
classes) implicitly add 's that belong to identical
pairs.
However this means that a matrix typically can be stored in more than
one possible "TsparseMatrix"
representations.
Use asUniqueT()
in order to ensure uniqueness of the
internal representation of such a matrix.
Class dgCMatrix
or the superclasses
dsparseMatrix
and
TsparseMatrix
; asUniqueT
.
m <- Matrix(0+1:28, nrow = 4) m[-3,c(2,4:5,7)] <- m[ 3, 1:4] <- m[1:3, 6] <- 0 (mT <- as(m, "TsparseMatrix")) str(mT) mT[1,] mT[4, drop = FALSE] stopifnot(identical(mT[lower.tri(mT)], m [lower.tri(m) ])) mT[lower.tri(mT,diag=TRUE)] <- 0 mT ## Triplet representation with repeated (i,j) entries ## *adds* the corresponding x's: T2 <- new("dgTMatrix", i = as.integer(c(1,1,0,3,3)), j = as.integer(c(2,2,4,0,0)), x=10*1:5, Dim=4:5) str(T2) # contains (i,j,x) slots exactly as above, but T2 ## has only three non-zero entries, as for repeated (i,j)'s, ## the corresponding x's are "implicitly" added stopifnot(nnzero(T2) == 3)
m <- Matrix(0+1:28, nrow = 4) m[-3,c(2,4:5,7)] <- m[ 3, 1:4] <- m[1:3, 6] <- 0 (mT <- as(m, "TsparseMatrix")) str(mT) mT[1,] mT[4, drop = FALSE] stopifnot(identical(mT[lower.tri(mT)], m [lower.tri(m) ])) mT[lower.tri(mT,diag=TRUE)] <- 0 mT ## Triplet representation with repeated (i,j) entries ## *adds* the corresponding x's: T2 <- new("dgTMatrix", i = as.integer(c(1,1,0,3,3)), j = as.integer(c(2,2,4,0,0)), x=10*1:5, Dim=4:5) str(T2) # contains (i,j,x) slots exactly as above, but T2 ## has only three non-zero entries, as for repeated (i,j)'s, ## the corresponding x's are "implicitly" added stopifnot(nnzero(T2) == 3)
Construct a formally diagonal Matrix
,
i.e., an object inheriting from virtual class
diagonalMatrix
(or, if desired, a mathematically diagonal
CsparseMatrix
).
Diagonal(n, x = NULL, names = FALSE) .sparseDiagonal(n, x = NULL, uplo = "U", shape = "t", unitri = TRUE, kind, cols) .trDiagonal(n, x = NULL, uplo = "U", unitri = TRUE, kind) .symDiagonal(n, x = NULL, uplo = "U", kind)
Diagonal(n, x = NULL, names = FALSE) .sparseDiagonal(n, x = NULL, uplo = "U", shape = "t", unitri = TRUE, kind, cols) .trDiagonal(n, x = NULL, uplo = "U", unitri = TRUE, kind) .symDiagonal(n, x = NULL, uplo = "U", kind)
n |
integer indicating the dimension of the (square) matrix.
If missing, then |
x |
numeric or logical vector listing values for the diagonal
entries, to be recycled as necessary. If |
names |
either |
uplo |
one of |
shape |
one of |
unitri |
logical indicating if a formally triangular result with
ones on the diagonal should be formally unit triangular, i.e.,
with |
kind |
one of |
cols |
optional integer vector with values in |
Diagonal()
returns an object inheriting from virtual class
diagonalMatrix
.
.sparseDiagonal()
returns a CsparseMatrix
representation of Diagonal(n, x)
or, if cols
is given,
of Diagonal(n, x)[, cols+1]
. The precise class of the result
depends on shape
and kind
.
.trDiagonal()
and .symDiagonal()
are simple wrappers,
for .sparseDiagonal(shape = "t")
and
.sparseDiagonal(shape = "s")
, respectively.
.sparseDiagonal()
exists primarily to leverage efficient
C-level methods available for CsparseMatrix
.
Martin Maechler
the generic function diag
for extraction
of the diagonal from a matrix works for all “Matrices”.
bandSparse
constructs a banded sparse matrix from
its non-zero sub-/super - diagonals. band(A)
returns a
band matrix containing some sub-/super - diagonals of A
.
Matrix
for general matrix construction;
further, class diagonalMatrix
.
Diagonal(3) Diagonal(x = 10^(3:1)) Diagonal(x = (1:4) >= 2)#-> "ldiMatrix" ## Use Diagonal() + kronecker() for "repeated-block" matrices: M1 <- Matrix(0+0:5, 2,3) (M <- kronecker(Diagonal(3), M1)) (S <- crossprod(Matrix(rbinom(60, size=1, prob=0.1), 10,6))) (SI <- S + 10*.symDiagonal(6)) # sparse symmetric still stopifnot(is(SI, "dsCMatrix")) (I4 <- .sparseDiagonal(4, shape="t"))# now (2012-10) unitriangular stopifnot(I4@diag == "U", all(I4 == diag(4)))
Diagonal(3) Diagonal(x = 10^(3:1)) Diagonal(x = (1:4) >= 2)#-> "ldiMatrix" ## Use Diagonal() + kronecker() for "repeated-block" matrices: M1 <- Matrix(0+0:5, 2,3) (M <- kronecker(Diagonal(3), M1)) (S <- crossprod(Matrix(rbinom(60, size=1, prob=0.1), 10,6))) (SI <- S + 10*.symDiagonal(6)) # sparse symmetric still stopifnot(is(SI, "dsCMatrix")) (I4 <- .sparseDiagonal(4, shape="t"))# now (2012-10) unitriangular stopifnot(I4@diag == "U", all(I4 == diag(4)))
Class "diagonalMatrix" is the virtual class of all diagonal matrices.
A virtual Class: No objects may be created from it.
diag
:character
string, either "U"
or
"N"
, where "U"
means ‘unit-diagonal’.
Dim
:matrix dimension, and
Dimnames
:the dimnames
, a
list
, see the Matrix
class
description. Typically list(NULL,NULL)
for diagonal matrices.
Class "sparseMatrix"
, directly.
These are just a subset of the signature for which defined methods. Currently, there are (too) many explicit methods defined in order to ensure efficient methods for diagonal matrices.
signature(from = "matrix", to = "diagonalMatrix")
: ...
signature(from = "Matrix", to = "diagonalMatrix")
: ...
signature(from = "diagonalMatrix", to = "generalMatrix")
: ...
signature(from = "diagonalMatrix", to = "triangularMatrix")
: ...
signature(from = "diagonalMatrix", to = "nMatrix")
: ...
signature(from = "diagonalMatrix", to = "matrix")
: ...
signature(from = "diagonalMatrix", to = "sparseVector")
: ...
signature(x = "diagonalMatrix")
: ...
and many more methods
signature(a = "diagonalMatrix", b, ...)
: is
trivially implemented, of course; see also solve-methods
.
signature(x = "nMatrix")
, semantically
equivalent to base function which(x, arr.ind)
.
signature(x = "diagonalMatrix")
: all these
group methods return a "diagonalMatrix"
, apart from
cumsum()
etc which return a vector also for
base matrix
.
signature(e1 = "ddiMatrix", e2="denseMatrix")
:
arithmetic and other operators from the Ops
group have a few dozen explicit method definitions, in order to
keep the results diagonal in many cases, including the following:
signature(e1 = "ddiMatrix", e2="denseMatrix")
:
the result is from class ddiMatrix
which is
typically very desirable. Note that when e2
contains
off-diagonal zeros or NA
s, we implicitly use , hence
differing from traditional R arithmetic (where
), in order to preserve sparsity.
(object = "diagonalMatrix")
: Returns
an object of S3 class "diagSummary"
which is the summary of
the vector object@x
plus a simple heading, and an
appropriate print
method.
Diagonal()
as constructor of these matrices, and
isDiagonal
.
ddiMatrix
and ldiMatrix
are
“actual” classes extending "diagonalMatrix"
.
I5 <- Diagonal(5) D5 <- Diagonal(x = 10*(1:5)) ## trivial (but explicitly defined) methods: stopifnot(identical(crossprod(I5), I5), identical(tcrossprod(I5), I5), identical(crossprod(I5, D5), D5), identical(tcrossprod(D5, I5), D5), identical(solve(D5), solve(D5, I5)), all.equal(D5, solve(solve(D5)), tolerance = 1e-12) ) solve(D5)# efficient as is diagonal # an unusual way to construct a band matrix: rbind2(cbind2(I5, D5), cbind2(D5, I5))
I5 <- Diagonal(5) D5 <- Diagonal(x = 10*(1:5)) ## trivial (but explicitly defined) methods: stopifnot(identical(crossprod(I5), I5), identical(tcrossprod(I5), I5), identical(crossprod(I5, D5), D5), identical(tcrossprod(D5, I5), D5), identical(solve(D5), solve(D5, I5)), all.equal(D5, solve(solve(D5)), tolerance = 1e-12) ) solve(D5)# efficient as is diagonal # an unusual way to construct a band matrix: rbind2(cbind2(I5, D5), cbind2(D5, I5))
Transform a triangular matrix x
, i.e., of class
triangularMatrix
,
from (internally!) unit triangular (“unitriangular”) to
“general” triangular (diagU2N(x)
) or back (diagN2U(x)
).
Note that the latter, diagN2U(x)
, also sets the diagonal to one
in cases where diag(x)
was not all one.
.diagU2N(x)
and .diagN2U(x)
assume without
checking that x
is a triangularMatrix
with
suitable diag
slot ("U"
and "N"
, respectively),
hence they should be used with care.
diagU2N(x, cl = getClassDef(class(x)), checkDense = FALSE) diagN2U(x, cl = getClassDef(class(x)), checkDense = FALSE) .diagU2N(x, cl = getClassDef(class(x)), checkDense = FALSE) .diagN2U(x, cl = getClassDef(class(x)), checkDense = FALSE)
diagU2N(x, cl = getClassDef(class(x)), checkDense = FALSE) diagN2U(x, cl = getClassDef(class(x)), checkDense = FALSE) .diagU2N(x, cl = getClassDef(class(x)), checkDense = FALSE) .diagN2U(x, cl = getClassDef(class(x)), checkDense = FALSE)
x |
a |
cl |
(optional, for speedup only:) class (definition) of |
checkDense |
logical indicating if dense (see
|
The concept of unit triangular matrices with a diag
slot of
"U"
stems from LAPACK.
a triangular matrix of the same class
but with a
different diag
slot. For diagU2N
(semantically) with
identical entries as x
, whereas in diagN2U(x)
, the
off-diagonal entries are unchanged and the diagonal is set to all
1
even if it was not previously.
Such internal storage details should rarely be of relevance to the user. Hence, these functions really are rather internal utilities.
"triangularMatrix"
,
"dtCMatrix"
.
(T <- Diagonal(7) + triu(Matrix(rpois(49, 1/4), 7, 7), k = 1)) (uT <- diagN2U(T)) # "unitriangular" (t.u <- diagN2U(10*T))# changes the diagonal! stopifnot(all(T == uT), diag(t.u) == 1, identical(T, diagU2N(uT))) T[upper.tri(T)] <- 5 # still "dtC" T <- diagN2U(as(T,"triangularMatrix")) dT <- as(T, "denseMatrix") # (unitriangular) dT.n <- diagU2N(dT, checkDense = TRUE) sT.n <- diagU2N(dT) stopifnot(is(dT.n, "denseMatrix"), is(sT.n, "sparseMatrix"), dT@diag == "U", dT.n@diag == "N", sT.n@diag == "N", all(dT == dT.n), all(dT == sT.n))
(T <- Diagonal(7) + triu(Matrix(rpois(49, 1/4), 7, 7), k = 1)) (uT <- diagN2U(T)) # "unitriangular" (t.u <- diagN2U(10*T))# changes the diagonal! stopifnot(all(T == uT), diag(t.u) == 1, identical(T, diagU2N(uT))) T[upper.tri(T)] <- 5 # still "dtC" T <- diagN2U(as(T,"triangularMatrix")) dT <- as(T, "denseMatrix") # (unitriangular) dT.n <- diagU2N(dT, checkDense = TRUE) sT.n <- diagU2N(dT) stopifnot(is(dT.n, "denseMatrix"), is(sT.n, "sparseMatrix"), dT@diag == "U", dT.n@diag == "N", sT.n@diag == "N", all(dT == dT.n), all(dT == sT.n))
dimScale
, rowScale
, and colScale
implement
D1 %*% x %*% D2
, D %*% x
, and x %*% D
for diagonal matrices D1
, D2
, and D
with
diagonal entries d1
, d2
, and d
, respectively.
Unlike the explicit products, these functions preserve dimnames(x)
and symmetry where appropriate.
dimScale(x, d1 = sqrt(1/diag(x, names = FALSE)), d2 = d1) rowScale(x, d) colScale(x, d)
dimScale(x, d1 = sqrt(1/diag(x, names = FALSE)), d2 = d1) rowScale(x, d) colScale(x, d)
x |
a matrix, possibly inheriting from virtual class
|
d1 , d2 , d
|
numeric vectors giving factors by which to scale
the rows or columns of |
dimScale(x)
(with d1
and d2
unset) is only
roughly equivalent to cov2cor(x)
. cov2cor
sets the diagonal entries of the result to 1 (exactly);
dimScale
does not.
The result of scaling x
, currently always inheriting from
virtual class dMatrix
.
It inherits from triangularMatrix
if and only
if x
does. In the special case of dimScale(x, d1, d2)
with identical d1
and d2
, it inherits from
symmetricMatrix
if and only if x
does.
Mikael Jagan
n <- 6L (x <- forceSymmetric(matrix(1, n, n))) dimnames(x) <- rep.int(list(letters[seq_len(n)]), 2L) d <- seq_len(n) (D <- Diagonal(x = d)) (scx <- dimScale(x, d)) # symmetry and 'dimnames' kept (mmx <- D %*% x %*% D) # symmetry and 'dimnames' lost stopifnot(identical(unname(as(scx, "generalMatrix")), mmx)) rowScale(x, d) colScale(x, d)
n <- 6L (x <- forceSymmetric(matrix(1, n, n))) dimnames(x) <- rep.int(list(letters[seq_len(n)]), 2L) d <- seq_len(n) (D <- Diagonal(x = d)) (scx <- dimScale(x, d)) # symmetry and 'dimnames' kept (mmx <- D %*% x %*% D) # symmetry and 'dimnames' lost stopifnot(identical(unname(as(scx, "generalMatrix")), mmx)) rowScale(x, d) colScale(x, d)
The dMatrix
class is a virtual class contained by all actual
classes of numeric matrices in the Matrix package. Similarly,
all the actual classes of logical matrices inherit from the
lMatrix
class.
Common to all matrix object in the package:
Dim
:Object of class "integer"
- the dimensions
of the matrix - must be an integer vector with exactly two
non-negative values.
Dimnames
:list of length two; each component
containing NULL or a character
vector length
equal the corresponding Dim
element.
There are (relatively simple) group methods (see, e.g., Arith
)
signature(e1 = "dMatrix", e2 = "dMatrix")
: ...
signature(e1 = "dMatrix", e2 = "numeric")
: ...
signature(e1 = "numeric", e2 = "dMatrix")
: ...
signature(x = "dMatrix")
: ...
signature(x = "dMatrix", digits = "numeric")
:
this group contains round()
and signif()
.
signature(e1 = "numeric", e2 = "dMatrix")
: ...
signature(e1 = "dMatrix", e2 = "numeric")
: ...
signature(e1 = "dMatrix", e2 = "dMatrix")
: ...
signature(x = "dMatrix")
: The "Summary"
group contains the seven functions
max()
, min()
, range()
,
prod()
, sum()
,
any()
, and all()
.
The following methods are also defined for all double matrices:
signature(x = "dMatrix")
: computes the
“Matrix Exponential”, see expm
.
The following methods are defined for all logical matrices:
signature(x = "lsparseMatrix")
and many other
subclasses of "lMatrix"
: as the base function
which(x, arr.ind)
returns the indices of the
TRUE
entries in x
; if arr.ind
is true,
as a 2-column matrix of row and column indices. Since Matrix
version 1.2-9, if useNames
is true, as by default, with
dimnames
, the same as base::which
.
The nonzero-pattern matrix class nMatrix
, which
can be used to store non-NA
logical
matrices even more compactly.
The numeric matrix classes dgeMatrix
,
dgCMatrix
, and Matrix
.
drop0(x, tol=1e-10)
is sometimes preferable to (and
more efficient than) zapsmall(x, digits=10)
.
showClass("dMatrix") set.seed(101) round(Matrix(rnorm(28), 4,7), 2) M <- Matrix(rlnorm(56, sd=10), 4,14) (M. <- zapsmall(M)) table(as.logical(M. == 0))
showClass("dMatrix") set.seed(101) round(Matrix(rnorm(28), 4,7), 2) M <- Matrix(rlnorm(56, sd=10), 4,14) (M. <- zapsmall(M)) table(as.logical(M. == 0))
For any (typically) sparse matrix
x
compute the Dulmage-Mendelsohn row and columns permutations which at
first splits the rows and
m
columns into coarse partitions
each; and then a finer one, reordering rows and columns such that the
permutated matrix is “as upper triangular” as possible.
dmperm(x, nAns = 6L, seed = 0L)
dmperm(x, nAns = 6L, seed = 0L)
x |
a typically sparse matrix; internally coerced to either
|
nAns |
an integer specifying the |
seed |
an integer code in -1,0,1; determining the (initial)
permutation; by default, |
See the book section by Tim Davis; page 122–127, in the References.
a named list
with (by default) 6 components,
p |
integer vector with the permutation |
q |
integer vector with the permutation |
r |
integer vector of length |
s |
integer vector of length |
rr5 |
integer vector of length 5, defining the coarse row decomposition. |
cc5 |
integer vector of length 5, defining the coarse column decomposition. |
Martin Maechler, with a lot of “encouragement” by Mauricio Vargas.
Section 7.4 Dulmage-Mendelsohn decomposition, pp. 122 ff of
Timothy A. Davis (2006)
Direct Methods for Sparse Linear Systems, SIAM Series
“Fundamentals of Algorithms”.
Schur
, the class of permutation matrices; "pMatrix"
.
set.seed(17) (S9 <- rsparsematrix(9, 9, nnz = 10, symmetric=TRUE)) # dsCMatrix str( dm9 <- dmperm(S9) ) (S9p <- with(dm9, S9[p, q])) ## looks good, but *not* quite upper triangular; these, too: str( dm9.0 <- dmperm(S9, seed=-1)) # non-random too. str( dm9_1 <- dmperm(S9, seed= 1)) # a random one ## The last two permutations differ, but have the same effect! (S9p0 <- with(dm9.0, S9[p, q])) # .. hmm .. stopifnot(all.equal(S9p0, S9p))# same as as default, but different from the random one set.seed(11) (M <- triu(rsparsematrix(9,11, 1/4))) dM <- dmperm(M); with(dM, M[p, q]) (Mp <- M[sample.int(nrow(M)), sample.int(ncol(M))]) dMp <- dmperm(Mp); with(dMp, Mp[p, q]) set.seed(7) (n7 <- rsparsematrix(5, 12, nnz = 10, rand.x = NULL)) str( dm.7 <- dmperm(n7) ) stopifnot(exprs = { lengths(dm.7[1:2]) == dim(n7) identical(dm.7, dmperm(as(n7, "dMatrix"))) identical(dm.7[1:4], dmperm(n7, nAns=4)) identical(dm.7[1:2], dmperm(n7, nAns=2)) })
set.seed(17) (S9 <- rsparsematrix(9, 9, nnz = 10, symmetric=TRUE)) # dsCMatrix str( dm9 <- dmperm(S9) ) (S9p <- with(dm9, S9[p, q])) ## looks good, but *not* quite upper triangular; these, too: str( dm9.0 <- dmperm(S9, seed=-1)) # non-random too. str( dm9_1 <- dmperm(S9, seed= 1)) # a random one ## The last two permutations differ, but have the same effect! (S9p0 <- with(dm9.0, S9[p, q])) # .. hmm .. stopifnot(all.equal(S9p0, S9p))# same as as default, but different from the random one set.seed(11) (M <- triu(rsparsematrix(9,11, 1/4))) dM <- dmperm(M); with(dM, M[p, q]) (Mp <- M[sample.int(nrow(M)), sample.int(ncol(M))]) dMp <- dmperm(Mp); with(dMp, Mp[p, q]) set.seed(7) (n7 <- rsparsematrix(5, 12, nnz = 10, rand.x = NULL)) str( dm.7 <- dmperm(n7) ) stopifnot(exprs = { lengths(dm.7[1:2]) == dim(n7) identical(dm.7, dmperm(as(n7, "dMatrix"))) identical(dm.7[1:4], dmperm(n7, nAns=4)) identical(dm.7[1:2], dmperm(n7, nAns=2)) })
The "dpoMatrix"
class is the class of
positive-semidefinite symmetric matrices in nonpacked storage.
The "dppMatrix"
class is the same except in packed
storage. Only the upper triangle or the lower triangle is
required to be available.
The "corMatrix"
and "copMatrix"
classes
represent correlation matrices. They extend "dpoMatrix"
and "dppMatrix"
, respectively, with an additional slot
sd
allowing restoration of the original covariance matrix.
Objects can be created by calls of the
form new("dpoMatrix", ...)
or from crossprod
applied to
an "dgeMatrix"
object.
uplo
:Object of class "character"
. Must be
either "U", for upper triangular, and "L", for lower triangular.
x
:Object of class "numeric"
. The numeric
values that constitute the matrix, stored in column-major order.
Dim
:Object of class "integer"
. The dimensions
of the matrix which must be a two-element vector of non-negative
integers.
Dimnames
:inherited from class "Matrix"
factors
:Object of class "list"
. A named
list of factorizations that have been computed for the matrix.
sd
:(for "corMatrix"
and "copMatrix"
)
a numeric
vector of length n
containing the
(original) entries which allow
reconstruction of a covariance matrix from the correlation matrix.
Class "dsyMatrix"
, directly.
Classes "dgeMatrix"
, "symmetricMatrix"
, and many more
by class "dsyMatrix"
.
signature(x = "dpoMatrix")
:
Returns (and stores) the Cholesky decomposition of x
, see
chol
.
signature(x = "dpoMatrix")
:
Returns the determinant
of x
, via
chol(x)
, see above.
signature(x = "dpoMatrix", norm = "character")
:
Returns (and stores) the reciprocal of the condition number of
x
. The norm
can be "O"
for the
one-norm (the default) or "I"
for the infinity-norm. For
symmetric matrices the result does not depend on the norm.
signature(a = "dpoMatrix", b = "....")
, and
signature(a = "dppMatrix", b = "....")
work
via the Cholesky composition, see also the Matrix solve-methods
.
signature(e1 = "dpoMatrix", e2 = "numeric")
(and
quite a few other signatures): The result of (“elementwise”
defined) arithmetic operations is typically not
positive-definite anymore. The only exceptions, currently, are
multiplications, divisions or additions with positive
length(.) == 1
numbers (or logical
s).
Currently the validity methods for these classes such as
getValidity(getClass("dpoMatrix"))
for efficiency reasons
only check the diagonal entries of the matrix – they may not be negative.
This is only necessary but not sufficient for a symmetric matrix to be
positive semi-definite.
A more reliable (but often more expensive) check for positive
semi-definiteness would look at the signs of diag(BunchKaufman(.))
(with some tolerance for very small negative values), and for (strict)
positive definiteness at something like
!inherits(tryCatch(chol(.), error=identity), "error")
.
Indeed, when coercing to these classes, a version
of Cholesky()
or chol()
is
typically used, e.g., see selectMethod("coerce",
c(from="dsyMatrix", to="dpoMatrix"))
.
Classes dsyMatrix
and dgeMatrix
;
further, Matrix
, rcond
,
chol
, solve
, crossprod
.
h6 <- Hilbert(6) rcond(h6) str(h6) h6 * 27720 # is ``integer'' solve(h6) str(hp6 <- pack(h6)) ### Note that as(*, "corMatrix") *scales* the matrix (ch6 <- as(h6, "corMatrix")) stopifnot(all.equal(as(h6 * 27720, "dsyMatrix"), round(27720 * h6), tolerance = 1e-14), all.equal(ch6@sd^(-2), 2*(1:6)-1, tolerance = 1e-12)) chch <- Cholesky(ch6, perm = FALSE) stopifnot(identical(chch, ch6@factors$Cholesky), all(abs(crossprod(as(chch, "dtrMatrix")) - ch6) < 1e-10))
h6 <- Hilbert(6) rcond(h6) str(h6) h6 * 27720 # is ``integer'' solve(h6) str(hp6 <- pack(h6)) ### Note that as(*, "corMatrix") *scales* the matrix (ch6 <- as(h6, "corMatrix")) stopifnot(all.equal(as(h6 * 27720, "dsyMatrix"), round(27720 * h6), tolerance = 1e-14), all.equal(ch6@sd^(-2), 2*(1:6)-1, tolerance = 1e-12)) chch <- Cholesky(ch6, perm = FALSE) stopifnot(identical(chch, ch6@factors$Cholesky), all(abs(crossprod(as(chch, "dtrMatrix")) - ch6) < 1e-10))
Deletes “non-structural” zeros (i.e., zeros stored explicitly, in memory) from a sparse matrix and returns the result.
drop0(x, tol = 0, is.Csparse = NA, give.Csparse = TRUE)
drop0(x, tol = 0, is.Csparse = NA, give.Csparse = TRUE)
x |
a |
tol |
a non-negative number. If |
is.Csparse |
a logical used only if |
give.Csparse |
a logical indicating if the result must
inherit from virtual class |
A sparseMatrix
, the result of deleting
non-structural zeros from x
, possibly after coercion.
drop0
is sometimes called in conjunction with
zapsmall
, e.g., when dealing with sparse
matrix products; see the example.
Function sparseMatrix
, for constructing objects
inheriting from virtual class sparseMatrix
;
nnzero
.
(m <- sparseMatrix(i = 1:8, j = 2:9, x = c(0:2, 3:-1), dims = c(10L, 20L))) drop0(m) ## A larger example: t5 <- new("dtCMatrix", Dim = c(5L, 5L), uplo = "L", x = c(10, 1, 3, 10, 1, 10, 1, 10, 10), i = c(0L,2L,4L, 1L, 3L,2L,4L, 3L, 4L), p = c(0L, 3L, 5L, 7:9)) TT <- kronecker(t5, kronecker(kronecker(t5, t5), t5)) IT <- solve(TT) I. <- TT %*% IT ; nnzero(I.) # 697 ( == 625 + 72 ) I.0 <- drop0(zapsmall(I.)) ## which actually can be more efficiently achieved by I.. <- drop0(I., tol = 1e-15) stopifnot(all(I.0 == Diagonal(625)), nnzero(I..) == 625)
(m <- sparseMatrix(i = 1:8, j = 2:9, x = c(0:2, 3:-1), dims = c(10L, 20L))) drop0(m) ## A larger example: t5 <- new("dtCMatrix", Dim = c(5L, 5L), uplo = "L", x = c(10, 1, 3, 10, 1, 10, 1, 10, 10), i = c(0L,2L,4L, 1L, 3L,2L,4L, 3L, 4L), p = c(0L, 3L, 5L, 7:9)) TT <- kronecker(t5, kronecker(kronecker(t5, t5), t5)) IT <- solve(TT) I. <- TT %*% IT ; nnzero(I.) # 697 ( == 625 + 72 ) I.0 <- drop0(zapsmall(I.)) ## which actually can be more efficiently achieved by I.. <- drop0(I., tol = 1e-15) stopifnot(all(I.0 == Diagonal(625)), nnzero(I..) == 625)
The dsCMatrix
class is a class of symmetric, sparse
numeric matrices in the compressed, column-oriented format. In
this implementation the non-zero elements in the columns are sorted
into increasing row order.
The dsTMatrix
class is the class of symmetric, sparse numeric
matrices in triplet format.
Objects can be created by calls of the form new("dsCMatrix",
...)
or new("dsTMatrix", ...)
, or automatically via e.g.,
as(*, "symmetricMatrix")
, or (for dsCMatrix
) also
from Matrix(.)
.
Creation “from scratch” most efficiently happens via
sparseMatrix(*, symmetric=TRUE)
.
uplo
:A character object indicating if the upper
triangle ("U"
) or the lower triangle ("L"
) is stored.
i
:Object of class "integer"
of length nnZ
(half number of non-zero elements). These are the row
numbers for each non-zero element in the lower triangle of the matrix.
p
:(only in class "dsCMatrix"
:) an
integer
vector for providing pointers, one for each
column, see the detailed description in CsparseMatrix
.
j
:(only in class "dsTMatrix"
:) Object of
class "integer"
of length nnZ (as i
). These are the
column numbers for each non-zero element in the lower triangle of
the matrix.
x
:Object of class "numeric"
of length nnZ –
the non-zero elements of the matrix (to be duplicated for full matrix).
factors
:Object of class "list"
- a list
of factorizations of the matrix.
Dim
:Object of class "integer"
- the dimensions
of the matrix - must be an integer vector with exactly two
non-negative values.
Both classes extend classes and symmetricMatrix
dsparseMatrix
directly;
dsCMatrix
further directly extends
CsparseMatrix
, where
dsTMatrix
does TsparseMatrix
.
signature(a = "dsCMatrix", b = "....")
: x
<- solve(a,b)
solves for
; see
solve-methods
.
signature(x = "dsCMatrix", pivot = "logical")
:
Returns (and stores) the Cholesky decomposition of x
, see
chol
.
signature(A = "dsCMatrix",...)
:
Computes more flexibly Cholesky decompositions,
see Cholesky
.
signature(x = "dsCMatrix", logarithm =
"missing")
: Evaluate the determinant of x
on the
logarithm scale. This creates and stores the Cholesky factorization.
signature(x = "dsCMatrix", logarithm =
"logical")
: Evaluate the determinant of x
on the
logarithm scale or not, according to the logarithm
argument. This creates and stores the Cholesky factorization.
signature(x = "dsCMatrix")
: Transpose. As for all
symmetric matrices, a matrix for which the upper triangle is
stored produces a matrix for which the lower triangle is stored
and vice versa, i.e., the uplo
slot is swapped, and the row
and column indices are interchanged.
signature(x = "dsTMatrix")
: Transpose. The
uplo
slot is swapped from "U"
to "L"
or vice
versa, as for a "dsCMatrix"
, see above.
Classes dgCMatrix
, dgTMatrix
,
dgeMatrix
and those mentioned above.
mm <- Matrix(toeplitz(c(10, 0, 1, 0, 3)), sparse = TRUE) mm # automatically dsCMatrix str(mm) mT <- as(as(mm, "generalMatrix"), "TsparseMatrix") ## Either (symM <- as(mT, "symmetricMatrix")) # dsT (symC <- as(symM, "CsparseMatrix")) # dsC ## or sT <- Matrix(mT, sparse=TRUE, forceCheck=TRUE) # dsT sym2 <- as(symC, "TsparseMatrix") ## --> the same as 'symM', a "dsTMatrix"
mm <- Matrix(toeplitz(c(10, 0, 1, 0, 3)), sparse = TRUE) mm # automatically dsCMatrix str(mm) mT <- as(as(mm, "generalMatrix"), "TsparseMatrix") ## Either (symM <- as(mT, "symmetricMatrix")) # dsT (symC <- as(symM, "CsparseMatrix")) # dsC ## or sT <- Matrix(mT, sparse=TRUE, forceCheck=TRUE) # dsT sym2 <- as(symC, "TsparseMatrix") ## --> the same as 'symM', a "dsTMatrix"
The Class "dsparseMatrix"
is the virtual (super) class of
all numeric sparse matrices.
Dim
:the matrix dimension, see class "Matrix"
.
Dimnames
:see the "Matrix"
class.
x
:a numeric
vector containing the
(non-zero) matrix entries.
Class "dMatrix"
and "sparseMatrix"
, directly.
Class "Matrix"
, by the above classes.
the documentation of the (non virtual) sub classes, see
showClass("dsparseMatrix")
; in particular,
dgTMatrix, dgCMatrix, and
dgRMatrix.
showClass("dsparseMatrix")
showClass("dsparseMatrix")
The dsRMatrix
class is a class of symmetric, sparse
matrices in the compressed, row-oriented format. In this
implementation the non-zero elements in the rows are sorted into
increasing column order.
These "..RMatrix"
classes are currently still mostly unimplemented!
Objects can be created by calls of the form new("dsRMatrix", ...)
.
uplo
:A character object indicating if the upper
triangle ("U"
) or the lower triangle ("L"
) is
stored. At present only the lower triangle form is allowed.
j
:Object of class "integer"
of length
nnzero
(number of non-zero elements). These are the row
numbers for each non-zero element in the matrix.
p
:Object of class "integer"
of pointers, one
for each row, to the initial (zero-based) index of elements in
the row.
factors
:Object of class "list"
- a list
of factorizations of the matrix.
x
:Object of class "numeric"
- the non-zero
elements of the matrix.
Dim
:Object of class "integer"
- the dimensions
of the matrix - must be an integer vector with exactly two
non-negative values.
Dimnames
:List of length two, see Matrix
.
Classes dsparseMatrix
,
symmetricMatrix
, and
RsparseMatrix
, directly.
Class "dMatrix"
, by class "dsparseMatrix"
;
class "sparseMatrix"
, by classes "dsparseMatrix"
and
"RsparseMatrix"
.
signature(x = "dsRMatrix", uplo = "missing")
:
a trivial method just returning x
signature(x = "dsRMatrix", uplo = "character")
:
if uplo == x@uplo
, this trivially returns x
;
otherwise t(x)
.
the classes dgCMatrix
,
dgTMatrix
, and dgeMatrix
.
(m0 <- new("dsRMatrix")) m2 <- new("dsRMatrix", Dim = c(2L,2L), x = c(3,1), j = c(1L,1L), p = 0:2) m2 stopifnot(colSums(as(m2, "TsparseMatrix")) == 3:4) str(m2) (ds2 <- forceSymmetric(diag(2))) # dsy* dR <- as(ds2, "RsparseMatrix") dR # dsRMatrix
(m0 <- new("dsRMatrix")) m2 <- new("dsRMatrix", Dim = c(2L,2L), x = c(3,1), j = c(1L,1L), p = 0:2) m2 stopifnot(colSums(as(m2, "TsparseMatrix")) == 3:4) str(m2) (ds2 <- forceSymmetric(diag(2))) # dsy* dR <- as(ds2, "RsparseMatrix") dR # dsRMatrix
The "dsyMatrix"
class is the class of symmetric, dense matrices
in non-packed storage and
"dspMatrix"
is the class of symmetric dense matrices in
packed storage, see pack()
. Only the upper
triangle or the lower triangle is stored.
Objects can be created by calls of the form new("dsyMatrix",
...)
or new("dspMatrix", ...)
, respectively.
uplo
:Object of class "character"
. Must be
either "U", for upper triangular, and "L", for lower triangular.
x
:Object of class "numeric"
. The numeric
values that constitute the matrix, stored in column-major order.
Dim
,Dimnames
:The dimension (a length-2
"integer"
) and corresponding names (or NULL
), see the
Matrix
.
factors
:Object of class "list"
. A named
list of factorizations that have been computed for the matrix.
"dsyMatrix"
extends class "dgeMatrix"
, directly, whereas"dspMatrix"
extends class "ddenseMatrix"
, directly.
Both extend class "symmetricMatrix"
, directly,
and class "Matrix"
and others, indirectly, use
showClass("dsyMatrix")
, e.g., for details.
signature(x = "dspMatrix", type = "character")
, or
x = "dsyMatrix"
or type = "missing"
: Computes the
matrix norm of the desired type, see, norm
.
signature(x = "dspMatrix", type = "character")
, or
x = "dsyMatrix"
or type = "missing"
: Computes the
reciprocal condition number, rcond()
.
signature(a = "dspMatrix", b = "....")
, and
signature(a = "dsyMatrix", b = "....")
: x
<- solve(a,b)
solves for
; see
solve-methods
.
signature(x = "dsyMatrix")
: Transpose; swaps from
upper triangular to lower triangular storage, i.e., the uplo slot
from "U"
to "L"
or vice versa, the same as for all
symmetric matrices.
The positive (Semi-)definite dense (packed or non-packed
numeric matrix classes dpoMatrix
,
dppMatrix
and corMatrix
,
Classes dgeMatrix
and Matrix
;
solve
, norm
, rcond
,
t
## Only upper triangular part matters (when uplo == "U" as per default) (sy2 <- new("dsyMatrix", Dim = as.integer(c(2,2)), x = c(14, NA,32,77))) str(t(sy2)) # uplo = "L", and the lower tri. (i.e. NA is replaced). chol(sy2) #-> "Cholesky" matrix (sp2 <- pack(sy2)) # a "dspMatrix" ## Coercing to dpoMatrix gives invalid object: sy3 <- new("dsyMatrix", Dim = as.integer(c(2,2)), x = c(14, -1, 2, -7)) try(as(sy3, "dpoMatrix")) # -> error: not positive definite ## 4x4 example m <- matrix(0,4,4); m[upper.tri(m)] <- 1:6 (sym <- m+t(m)+diag(11:14, 4)) (S1 <- pack(sym)) (S2 <- t(S1)) stopifnot(all(S1 == S2)) # equal "seen as matrix", but differ internally : str(S1) S2@x
## Only upper triangular part matters (when uplo == "U" as per default) (sy2 <- new("dsyMatrix", Dim = as.integer(c(2,2)), x = c(14, NA,32,77))) str(t(sy2)) # uplo = "L", and the lower tri. (i.e. NA is replaced). chol(sy2) #-> "Cholesky" matrix (sp2 <- pack(sy2)) # a "dspMatrix" ## Coercing to dpoMatrix gives invalid object: sy3 <- new("dsyMatrix", Dim = as.integer(c(2,2)), x = c(14, -1, 2, -7)) try(as(sy3, "dpoMatrix")) # -> error: not positive definite ## 4x4 example m <- matrix(0,4,4); m[upper.tri(m)] <- 1:6 (sym <- m+t(m)+diag(11:14, 4)) (S1 <- pack(sym)) (S2 <- t(S1)) stopifnot(all(S1 == S2)) # equal "seen as matrix", but differ internally : str(S1) S2@x
The "dtCMatrix"
class is a class of triangular, sparse
matrices in the compressed, column-oriented format. In this
implementation the non-zero elements in the columns are sorted into
increasing row order.
The "dtTMatrix"
class is a class of triangular, sparse matrices
in triplet format.
Objects can be created by calls of the form new("dtCMatrix",
...)
or calls of the form new("dtTMatrix", ...)
,
but more typically automatically via Matrix()
or coercions such as as(x, "triangularMatrix")
.
uplo
:Object of class "character"
. Must be
either "U", for upper triangular, and "L", for lower triangular.
diag
:Object of class "character"
. Must be
either "U"
, for unit triangular (diagonal is all ones), or
"N"
; see triangularMatrix
.
p
:(only present in "dtCMatrix"
:) an
integer
vector for providing pointers, one for each
column, see the detailed description in CsparseMatrix
.
i
:Object of class "integer"
of length nnzero
(number of non-zero elements). These are the row numbers for
each non-zero element in the matrix.
j
:Object of class "integer"
of length nnzero
(number of non-zero elements). These are the column numbers for
each non-zero element in the matrix. (Only present in the
dtTMatrix
class.)
x
:Object of class "numeric"
- the non-zero
elements of the matrix.
Dim
,Dimnames
:The dimension (a length-2
"integer"
) and corresponding names (or NULL
),
inherited from the Matrix
, see there.
Class "dgCMatrix"
, directly.
Class "triangularMatrix"
, directly.
Class "dMatrix"
, "sparseMatrix"
, and more by class
"dgCMatrix"
etc, see the examples.
signature(a = "dtCMatrix", b = "....")
:
sparse triangular solve (aka “backsolve” or
“forwardsolve”), see solve-methods
.
signature(x = "dtCMatrix")
: returns the transpose of
x
signature(x = "dtTMatrix")
: returns the transpose of
x
Classes dgCMatrix
, dgTMatrix
,
dgeMatrix
, and dtrMatrix
.
showClass("dtCMatrix") showClass("dtTMatrix") t1 <- new("dtTMatrix", x= c(3,7), i= 0:1, j=3:2, Dim= as.integer(c(4,4))) t1 ## from 0-diagonal to unit-diagonal {low-level step}: tu <- t1 ; tu@diag <- "U" tu (cu <- as(tu, "CsparseMatrix")) str(cu)# only two entries in @i and @x stopifnot(cu@i == 1:0, all(2 * symmpart(cu) == Diagonal(4) + forceSymmetric(cu))) t1[1,2:3] <- -1:-2 diag(t1) <- 10*c(1:2,3:2) t1 # still triangular (it1 <- solve(t1)) t1. <- solve(it1) all(abs(t1 - t1.) < 10 * .Machine$double.eps) ## 2nd example U5 <- new("dtCMatrix", i= c(1L, 0:3), p=c(0L,0L,0:2, 5L), Dim = c(5L, 5L), x = rep(1, 5), diag = "U") U5 (iu <- solve(U5)) # contains one '0' validObject(iu2 <- solve(U5, Diagonal(5)))# failed in earlier versions I5 <- iu %*% U5 # should equal the identity matrix i5 <- iu2 %*% U5 m53 <- matrix(1:15, 5,3, dimnames=list(NULL,letters[1:3])) asDiag <- function(M) as(drop0(M), "diagonalMatrix") stopifnot( all.equal(Diagonal(5), asDiag(I5), tolerance=1e-14) , all.equal(Diagonal(5), asDiag(i5), tolerance=1e-14) , identical(list(NULL, dimnames(m53)[[2]]), dimnames(solve(U5, m53))) )
showClass("dtCMatrix") showClass("dtTMatrix") t1 <- new("dtTMatrix", x= c(3,7), i= 0:1, j=3:2, Dim= as.integer(c(4,4))) t1 ## from 0-diagonal to unit-diagonal {low-level step}: tu <- t1 ; tu@diag <- "U" tu (cu <- as(tu, "CsparseMatrix")) str(cu)# only two entries in @i and @x stopifnot(cu@i == 1:0, all(2 * symmpart(cu) == Diagonal(4) + forceSymmetric(cu))) t1[1,2:3] <- -1:-2 diag(t1) <- 10*c(1:2,3:2) t1 # still triangular (it1 <- solve(t1)) t1. <- solve(it1) all(abs(t1 - t1.) < 10 * .Machine$double.eps) ## 2nd example U5 <- new("dtCMatrix", i= c(1L, 0:3), p=c(0L,0L,0:2, 5L), Dim = c(5L, 5L), x = rep(1, 5), diag = "U") U5 (iu <- solve(U5)) # contains one '0' validObject(iu2 <- solve(U5, Diagonal(5)))# failed in earlier versions I5 <- iu %*% U5 # should equal the identity matrix i5 <- iu2 %*% U5 m53 <- matrix(1:15, 5,3, dimnames=list(NULL,letters[1:3])) asDiag <- function(M) as(drop0(M), "diagonalMatrix") stopifnot( all.equal(Diagonal(5), asDiag(I5), tolerance=1e-14) , all.equal(Diagonal(5), asDiag(i5), tolerance=1e-14) , identical(list(NULL, dimnames(m53)[[2]]), dimnames(solve(U5, m53))) )
The "dtpMatrix"
class is the class of triangular,
dense, numeric matrices in packed storage. The "dtrMatrix"
class is the same except in nonpacked storage.
Objects can be created by calls of the form new("dtpMatrix",
...)
or by coercion from other classes of matrices.
uplo
:Object of class "character"
. Must be
either "U", for upper triangular, and "L", for lower triangular.
diag
:Object of class "character"
. Must be
either "U"
, for unit triangular (diagonal is all ones), or
"N"
; see triangularMatrix
.
x
:Object of class "numeric"
. The numeric
values that constitute the matrix, stored in column-major order.
For a packed square matrix of dimension ,
length(x)
is of length (also when
diag == "U"
!).
Dim
,Dimnames
:The dimension (a length-2
"integer"
) and corresponding names (or NULL
),
inherited from the Matrix
, see there.
Class "ddenseMatrix"
, directly.
Class "triangularMatrix"
, directly.
Class "dMatrix"
and more by class "ddenseMatrix"
etc, see
the examples.
signature(x = "dtpMatrix", y = "dgeMatrix")
:
Matrix multiplication; ditto for several other signature
combinations, see showMethods("%*%", class = "dtpMatrix")
.
signature(x = "dtpMatrix", logarithm = "logical")
:
the determinant(x)
trivially is
prod(diag(x))
, but computed on log scale to prevent over-
and underflow.
signature(x = "dtpMatrix")
: ...
signature(x = "dtpMatrix", type = "character")
: ...
signature(x = "dtpMatrix", norm = "character")
: ...
signature(a = "dtpMatrix", b = "...")
:
efficiently using internal backsolve or forwardsolve, see
solve-methods
.
signature(x = "dtpMatrix")
: t(x)
remains
a "dtpMatrix"
, lower triangular if x
is upper
triangular, and vice versa.
Class dtrMatrix
showClass("dtrMatrix") example("dtrMatrix-class", echo=FALSE) (p1 <- pack(T2)) str(p1) (pp <- pack(T)) ip1 <- solve(p1) stopifnot(length(p1@x) == 3, length(pp@x) == 3, p1 @ uplo == T2 @ uplo, pp @ uplo == T @ uplo, identical(t(pp), p1), identical(t(p1), pp), all((l.d <- p1 - T2) == 0), is(l.d, "dtpMatrix"), all((u.d <- pp - T ) == 0), is(u.d, "dtpMatrix"), l.d@uplo == T2@uplo, u.d@uplo == T@uplo, identical(t(ip1), solve(pp)), is(ip1, "dtpMatrix"), all.equal(as(solve(p1,p1), "diagonalMatrix"), Diagonal(2)))
showClass("dtrMatrix") example("dtrMatrix-class", echo=FALSE) (p1 <- pack(T2)) str(p1) (pp <- pack(T)) ip1 <- solve(p1) stopifnot(length(p1@x) == 3, length(pp@x) == 3, p1 @ uplo == T2 @ uplo, pp @ uplo == T @ uplo, identical(t(pp), p1), identical(t(p1), pp), all((l.d <- p1 - T2) == 0), is(l.d, "dtpMatrix"), all((u.d <- pp - T ) == 0), is(u.d, "dtpMatrix"), l.d@uplo == T2@uplo, u.d@uplo == T@uplo, identical(t(ip1), solve(pp)), is(ip1, "dtpMatrix"), all.equal(as(solve(p1,p1), "diagonalMatrix"), Diagonal(2)))
The "dtrMatrix"
class is the class of triangular, dense,
numeric matrices in nonpacked storage. The "dtpMatrix"
class
is the same except in packed storage, see pack()
.
Objects can be created by calls of the form new("dtrMatrix", ...)
.
uplo
:Object of class "character"
. Must be
either "U", for upper triangular, and "L", for lower triangular.
diag
:Object of class "character"
. Must be
either "U"
, for unit triangular (diagonal is all ones), or
"N"
; see triangularMatrix
.
x
:Object of class "numeric"
. The numeric
values that constitute the matrix, stored in column-major order.
Dim
:Object of class "integer"
. The dimensions
of the matrix which must be a two-element vector of non-negative
integers.
Class "ddenseMatrix"
, directly.
Class "triangularMatrix"
, directly.
Class "Matrix"
and others, by class "ddenseMatrix"
.
Among others (such as matrix products, e.g. ?crossprod-methods
),
signature(x = "dtrMatrix", type = "character")
: ..
signature(x = "dtrMatrix", norm = "character")
: ..
signature(a = "dtrMatrix", b = "....")
: efficiently
use a “forwardsolve” or backsolve
for a lower or
upper triangular matrix, respectively, see also
solve-methods
.
all the Ops
group
methods are available. When applied to two triangular matrices,
these return a triangular matrix when easily possible.
Classes ddenseMatrix
, dtpMatrix
,
triangularMatrix
(m <- rbind(2:3, 0:-1)) (M <- as(m, "generalMatrix")) (T <- as(M, "triangularMatrix")) # formally upper triangular (T2 <- as(t(M), "triangularMatrix")) stopifnot(T@uplo == "U", T2@uplo == "L", identical(T2, t(T))) m <- matrix(0,4,4); m[upper.tri(m)] <- 1:6 (t1 <- Matrix(m+diag(,4))) str(t1p <- pack(t1)) (t1pu <- diagN2U(t1p)) stopifnot(exprs = { inherits(t1 , "dtrMatrix"); validObject(t1) inherits(t1p, "dtpMatrix"); validObject(t1p) inherits(t1pu,"dtCMatrix"); validObject(t1pu) t1pu@x == 1:6 all(t1pu == t1p) identical((t1pu - t1)@x, numeric())# sparse all-0 })
(m <- rbind(2:3, 0:-1)) (M <- as(m, "generalMatrix")) (T <- as(M, "triangularMatrix")) # formally upper triangular (T2 <- as(t(M), "triangularMatrix")) stopifnot(T@uplo == "U", T2@uplo == "L", identical(T2, t(T))) m <- matrix(0,4,4); m[upper.tri(m)] <- 1:6 (t1 <- Matrix(m+diag(,4))) str(t1p <- pack(t1)) (t1pu <- diagN2U(t1p)) stopifnot(exprs = { inherits(t1 , "dtrMatrix"); validObject(t1) inherits(t1p, "dtpMatrix"); validObject(t1p) inherits(t1pu,"dtCMatrix"); validObject(t1pu) t1pu@x == 1:6 all(t1pu == t1p) identical((t1pu - t1)@x, numeric())# sparse all-0 })
The dtRMatrix
class is a class of triangular, sparse
matrices in the compressed, row-oriented format. In this
implementation the non-zero elements in the rows are sorted into
increasing columnd order.
This class is currently still mostly unimplemented!
Objects can be created by calls of the form new("dtRMatrix", ...)
.
uplo
:Object of class "character"
. Must be
either "U", for upper triangular, and "L", for lower triangular.
At present only the lower triangle form is allowed.
diag
:Object of class "character"
. Must be
either "U"
, for unit triangular (diagonal is all ones), or
"N"
; see triangularMatrix
.
j
:Object of class "integer"
of length
nnzero(.)
(number of non-zero elements). These are
the row numbers for each non-zero element in the matrix.
p
:Object of class "integer"
of pointers, one
for each row, to the initial (zero-based) index of elements in
the row. (Only present in the dsRMatrix
class.)
x
:Object of class "numeric"
- the non-zero
elements of the matrix.
Dim
:The dimension (a length-2 "integer"
)
Dimnames
:corresponding names (or NULL
),
inherited from the Matrix
, see there.
Class "dgRMatrix"
, directly.
Class "dsparseMatrix"
, by class "dgRMatrix"
.
Class "dMatrix"
, by class "dgRMatrix"
.
Class "sparseMatrix"
, by class "dgRMatrix"
.
Class "Matrix"
, by class "dgRMatrix"
.
No methods currently with class "dsRMatrix" in the signature.
Classes dgCMatrix
, dgTMatrix
,
dgeMatrix
(m0 <- new("dtRMatrix")) (m2 <- new("dtRMatrix", Dim = c(2L,2L), x = c(5, 1:2), p = c(0L,2:3), j= c(0:1,1L))) str(m2) (m3 <- as(Diagonal(2), "RsparseMatrix"))# --> dtRMatrix
(m0 <- new("dtRMatrix")) (m2 <- new("dtRMatrix", Dim = c(2L,2L), x = c(5, 1:2), p = c(0L,2:3), j= c(0:1,1L))) str(m2) (m3 <- as(Diagonal(2), "RsparseMatrix"))# --> dtRMatrix
expand1
and expand2
construct matrix factors from
objects specifying matrix factorizations. Such objects typically
do not store the factors explicitly, employing instead a compact
representation to save memory.
expand1(x, which, ...) expand2(x, ...) expand (x, ...)
expand1(x, which, ...) expand2(x, ...) expand (x, ...)
x |
a matrix factorization, typically inheriting from
virtual class |
which |
a character string indicating a matrix factor. |
... |
further arguments passed to or from methods. |
Methods for expand
are retained only for backwards
compatibility with Matrix < 1.6-0
. New code
should use expand1
and expand2
, whose methods
provide more control and behave more consistently. Notably,
expand2
obeys the rule that the product of the matrix
factors in the returned list should reproduce
(within some tolerance) the factorized matrix,
including its dimnames
.
Hence if x
is a matrix and y
is its factorization,
then
all.equal(as(x, "matrix"), as(Reduce(`%*%`, expand2(y)), "matrix"))
should in most cases return TRUE
.
expand1
returns an object inheriting from virtual class
Matrix
, representing the factor indicated
by which
, always without row and column names.
expand2
returns a list of factors, typically with names
using conventional notation, as in list(L=, U=)
.
The first and last factors get the row and column names of the
factorized matrix, which are preserved in the Dimnames
slot of x
.
The following table lists methods for expand1
together with
allowed values of argument which
.
class(x) |
which
|
Schur |
c("Q", "T", "Q.")
|
denseLU |
c("P1", "P1.", "L", "U")
|
sparseLU |
c("P1", "P1.", "P2", "P2.", "L", "U")
|
sparseQR |
c("P1", "P1.", "P2", "P2.", "Q", "Q1", "R", "R1")
|
BunchKaufman , pBunchKaufman |
c("U", "DU", "U.", "L", "DL", "L.")
|
Cholesky , pCholesky |
c("P1", "P1.", "L1", "D", "L1.", "L", "L.")
|
CHMsimpl , CHMsimpl |
c("P1", "P1.", "L1", "D", "L1.", "L", "L.")
|
Methods for expand2
and expand
are described
below. Factor names and classes apply also to expand1
.
expand2
signature(x = "CHMsimpl")
:
expands the factorization
as
list(P1., L1, D, L1., P1)
(the default)
or as list(P1., L, L., P1)
,
depending on optional logical argument LDL
.
P1
and P1.
are pMatrix
,
L1
, L1.
, L
, and L.
are
dtCMatrix
,
and D
is a ddiMatrix
.
expand2
signature(x = "CHMsuper")
:
as CHMsimpl
, but the triangular factors are
stored as dgCMatrix
.
expand2
signature(x = "p?Cholesky")
:
expands the factorization
as
list(L1, D, L1.)
(the default) or as list(L, L.)
,
depending on optional logical argument LDL
.
L1
, L1.
, L
, and L.
are
dtrMatrix
or dtpMatrix
,
and D
is a ddiMatrix
.
expand2
signature(x = "p?BunchKaufman")
:
expands the factorization
where
and
as
list(U, DU, U.)
or list(L, DL, L.)
,
depending on x@uplo
. If optional argument complete
is TRUE
, then an unnamed list giving the full expansion
with or
matrix
factors is returned instead.
are represented as
pMatrix
,
and
are represented as
dtCMatrix
, and
and
are represented as
dsCMatrix
.
expand2
signature(x = "Schur")
:
expands the factorization
as
list(Q, T, Q.)
.
Q
and Q.
are x@Q
and t(x@Q)
modulo Dimnames
, and T
is x@T
.
expand2
signature(x = "sparseLU")
:
expands the factorization
as
list(P1., L, U, P2.)
.
P1.
and P2.
are pMatrix
,
and L
and U
are dtCMatrix
.
expand2
signature(x = "denseLU")
:
expands the factorization
as
list(P1., L, U)
.
P1.
is a pMatrix
,
and L
and U
are dtrMatrix
if square and dgeMatrix
otherwise.
expand2
signature(x = "sparseQR")
:
expands the factorization
as
list(P1., Q, R, P2.)
or list(P1., Q1, R1, P2.)
,
depending on optional logical argument complete
.
P1.
and P2.
are pMatrix
,
Q
and Q1
are dgeMatrix
,
R
is a dgCMatrix
,
and R1
is a dtCMatrix
.
expand
signature(x = "CHMfactor")
:
as expand2
, but returning list(P, L)
.
expand(x)[["P"]]
and expand2(x)[["P1"]]
represent the same permutation matrix
but have opposite
margin
slots and inverted
perm
slots. The components of expand(x)
do not preserve x@Dimnames
.
expand
signature(x = "sparseLU")
:
as expand2
, but returning list(P, L, U, Q)
.
expand(x)[["Q"]]
and expand2(x)[["P2."]]
represent the same permutation matrix
but have opposite
margin
slots and inverted
perm
slots. expand(x)[["P"]]
represents
the permutation matrix rather than its
transpose
; it is
expand2(x)[["P1."]]
with an inverted perm
slot. expand(x)[["L"]]
and expand2(x)[["L"]]
represent the same unit lower
triangular matrix , but with
diag
slot equal
to "N"
and "U"
, respectively.
expand(x)[["L"]]
and expand(x)[["U"]]
store the permuted first and second components of
x@Dimnames
in their Dimnames
slots.
expand
signature(x = "denseLU")
:
as expand2
, but returning list(L, U, P)
.
expand(x)[["P"]]
and expand2(x)[["P1."]]
are identical modulo Dimnames
. The components
of expand(x)
do not preserve x@Dimnames
.
The virtual class MatrixFactorization
of matrix factorizations.
Generic functions Cholesky
, BunchKaufman
,
Schur
, lu
, and qr
for
computing factorizations.
showMethods("expand1", inherited = FALSE) showMethods("expand2", inherited = FALSE) set.seed(0) (A <- Matrix(rnorm(9L, 0, 10), 3L, 3L)) (lu.A <- lu(A)) (e.lu.A <- expand2(lu.A)) stopifnot(exprs = { is.list(e.lu.A) identical(names(e.lu.A), c("P1.", "L", "U")) all(sapply(e.lu.A, is, "Matrix")) all.equal(as(A, "matrix"), as(Reduce(`%*%`, e.lu.A), "matrix")) }) ## 'expand1' and 'expand2' give equivalent results modulo ## dimnames and representation of permutation matrices; ## see also function 'alt' in example("Cholesky-methods") (a1 <- sapply(names(e.lu.A), expand1, x = lu.A, simplify = FALSE)) all.equal(a1, e.lu.A) ## see help("denseLU-class") and others for more examples
showMethods("expand1", inherited = FALSE) showMethods("expand2", inherited = FALSE) set.seed(0) (A <- Matrix(rnorm(9L, 0, 10), 3L, 3L)) (lu.A <- lu(A)) (e.lu.A <- expand2(lu.A)) stopifnot(exprs = { is.list(e.lu.A) identical(names(e.lu.A), c("P1.", "L", "U")) all(sapply(e.lu.A, is, "Matrix")) all.equal(as(A, "matrix"), as(Reduce(`%*%`, e.lu.A), "matrix")) }) ## 'expand1' and 'expand2' give equivalent results modulo ## dimnames and representation of permutation matrices; ## see also function 'alt' in example("Cholesky-methods") (a1 <- sapply(names(e.lu.A), expand1, x = lu.A, simplify = FALSE)) all.equal(a1, e.lu.A) ## see help("denseLU-class") and others for more examples
Compute the exponential of a matrix.
expm(x)
expm(x)
x |
a matrix, typically inheriting from the
|
The exponential of a matrix is defined as the infinite Taylor
series expm(A) = I + A + A^2/2! + A^3/3! + ...
(although this is
definitely not the way to compute it). The method for the
dgeMatrix
class uses Ward's diagonal Pade' approximation with
three step preconditioning, a recommendation from
Moler & Van Loan (1978) “Nineteen dubious ways...”.
The matrix exponential of x
.
This is a translation of the implementation of the corresponding Octave function contributed to the Octave project by A. Scottedward Hodel [email protected]. A bug in there has been fixed by Martin Maechler.
https://en.wikipedia.org/wiki/Matrix_exponential
Cleve Moler and Charles Van Loan (2003) Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later. SIAM Review 45, 1, 3–49. doi:10.1137/S00361445024180
for historical reference mostly:
Moler, C. and Van Loan, C. (1978)
Nineteen dubious ways to compute the exponential of a matrix.
SIAM Review 20, 4, 801–836.
doi:10.1137/1020098
Eric W. Weisstein et al. (1999) Matrix Exponential. From MathWorld, https://mathworld.wolfram.com/MatrixExponential.html
Package expm, which provides newer (in some cases
faster, more accurate) algorithms for computing the matrix
exponential via its own (non-generic) function expm()
.
expm also implements logm()
, sqrtm()
, etc.
Generic function Schur
.
(m1 <- Matrix(c(1,0,1,1), ncol = 2)) (e1 <- expm(m1)) ; e <- exp(1) stopifnot(all.equal(e1@x, c(e,0,e,e), tolerance = 1e-15)) (m2 <- Matrix(c(-49, -64, 24, 31), ncol = 2)) (e2 <- expm(m2)) (m3 <- Matrix(cbind(0,rbind(6*diag(3),0))))# sparse! (e3 <- expm(m3)) # upper triangular
(m1 <- Matrix(c(1,0,1,1), ncol = 2)) (e1 <- expm(m1)) ; e <- exp(1) stopifnot(all.equal(e1@x, c(e,0,e,e), tolerance = 1e-15)) (m2 <- Matrix(c(-49, -64, 24, 31), ncol = 2)) (e2 <- expm(m2)) (m3 <- Matrix(cbind(0,rbind(6*diag(3),0))))# sparse! (e3 <- expm(m3)) # upper triangular
Read matrices stored in the Harwell-Boeing or MatrixMarket formats
or write sparseMatrix
objects to one of these
formats.
readHB(file) readMM(file) writeMM(obj, file, ...)
readHB(file) readMM(file) writeMM(obj, file, ...)
obj |
a real sparse matrix |
file |
for Alternatively, |
... |
optional additional arguments. Currently none are used in any methods. |
The readHB
and readMM
functions return an object that
inherits from the "Matrix"
class. Methods for the
writeMM
generic functions usually return
NULL
and, as a side effect, the matrix obj
is
written to file
in the MatrixMarket format (writeMM).
The Harwell-Boeing format is older and less flexible than the
MatrixMarket format. The function writeHB
was deprecated and
has now been removed. Please use writeMM
instead.
Note that these formats do not know anything about
dimnames
, hence these are dropped by writeMM()
.
A very simple way to export small sparse matrices S
, is to use
summary(S)
which returns a data.frame
with
columns i
, j
, and possibly x
, see summary
in
sparseMatrix-class
, and an example below.
https://math.nist.gov/MatrixMarket/
str(pores <- readMM(system.file("external/pores_1.mtx", package = "Matrix"))) str(utm <- readHB(system.file("external/utm300.rua" , package = "Matrix"))) str(lundA <- readMM(system.file("external/lund_a.mtx" , package = "Matrix"))) str(lundA <- readHB(system.file("external/lund_a.rsa" , package = "Matrix"))) ## https://math.nist.gov/MatrixMarket/data/Harwell-Boeing/counterx/counterx.htm str(jgl <- readMM(system.file("external/jgl009.mtx" , package = "Matrix"))) ## NOTE: The following examples take quite some time ## ---- even on a fast internet connection: if(FALSE) { ## The URL has been corrected, but we need an untar step: u. <- url("https://www.cise.ufl.edu/research/sparse/RB/Boeing/msc00726.tar.gz") str(sm <- readHB(gzcon(u.))) } data(KNex, package = "Matrix") ## Store as MatrixMarket (".mtx") file, here inside temporary dir./folder: (MMfile <- file.path(tempdir(), "mmMM.mtx")) writeMM(KNex$mm, file=MMfile) file.info(MMfile)[,c("size", "ctime")] # (some confirmation of the file's) ## very simple export - in triplet format - to text file: data(CAex, package = "Matrix") s.CA <- summary(CAex) s.CA # shows (i, j, x) [columns of a data frame] message("writing to ", outf <- tempfile()) write.table(s.CA, file = outf, row.names=FALSE) ## and read it back -- showing off sparseMatrix(): str(dd <- read.table(outf, header=TRUE)) ## has columns (i, j, x) -> we can use via do.call() as arguments to sparseMatrix(): mm <- do.call(sparseMatrix, dd) stopifnot(all.equal(mm, CAex, tolerance=1e-15))
str(pores <- readMM(system.file("external/pores_1.mtx", package = "Matrix"))) str(utm <- readHB(system.file("external/utm300.rua" , package = "Matrix"))) str(lundA <- readMM(system.file("external/lund_a.mtx" , package = "Matrix"))) str(lundA <- readHB(system.file("external/lund_a.rsa" , package = "Matrix"))) ## https://math.nist.gov/MatrixMarket/data/Harwell-Boeing/counterx/counterx.htm str(jgl <- readMM(system.file("external/jgl009.mtx" , package = "Matrix"))) ## NOTE: The following examples take quite some time ## ---- even on a fast internet connection: if(FALSE) { ## The URL has been corrected, but we need an untar step: u. <- url("https://www.cise.ufl.edu/research/sparse/RB/Boeing/msc00726.tar.gz") str(sm <- readHB(gzcon(u.))) } data(KNex, package = "Matrix") ## Store as MatrixMarket (".mtx") file, here inside temporary dir./folder: (MMfile <- file.path(tempdir(), "mmMM.mtx")) writeMM(KNex$mm, file=MMfile) file.info(MMfile)[,c("size", "ctime")] # (some confirmation of the file's) ## very simple export - in triplet format - to text file: data(CAex, package = "Matrix") s.CA <- summary(CAex) s.CA # shows (i, j, x) [columns of a data frame] message("writing to ", outf <- tempfile()) write.table(s.CA, file = outf, row.names=FALSE) ## and read it back -- showing off sparseMatrix(): str(dd <- read.table(outf, header=TRUE)) ## has columns (i, j, x) -> we can use via do.call() as arguments to sparseMatrix(): mm <- do.call(sparseMatrix, dd) stopifnot(all.equal(mm, CAex, tolerance=1e-15))
Multiplies a matrix or vector on the left or right by a factor from a matrix factorization or its transpose.
facmul(x, factor, y, trans = FALSE, left = TRUE, ...)
facmul(x, factor, y, trans = FALSE, left = TRUE, ...)
x |
a |
factor |
a character string indicating a factor in the
factorization represented by |
y |
a matrix or vector to be multiplied on the left or right by the factor or its transpose. |
trans |
a logical indicating if the transpose of the factor should be used, rather than the factor itself. |
left |
a logical indicating if the |
... |
further arguments passed to or from methods. |
facmul
is experimental and currently no methods are
exported from Matrix.
The value of op(M) %*% y
or y %*% op(M)
,
depending on left
, where M
is the factor
(always without dimnames
) and op(M)
is M
or t(M)
, depending on trans
.
## Conceptually, methods for 'facmul' _would_ behave as follows ... ## Not run: n <- 3L x <- lu(Matrix(rnorm(n * n), n, n)) y <- rnorm(n) L <- unname(expand2(x)[[nm <- "L"]]) stopifnot(exprs = { all.equal(facmul(x, nm, y, trans = FALSE, left = TRUE), L %*% y) all.equal(facmul(x, nm, y, trans = FALSE, left = FALSE), y %*% L) all.equal(facmul(x, nm, y, trans = TRUE, left = TRUE), crossprod(L, y)) all.equal(facmul(x, nm, y, trans = TRUE, left = FALSE), tcrossprod(y, L)) }) ## End(Not run)
## Conceptually, methods for 'facmul' _would_ behave as follows ... ## Not run: n <- 3L x <- lu(Matrix(rnorm(n * n), n, n)) y <- rnorm(n) L <- unname(expand2(x)[[nm <- "L"]]) stopifnot(exprs = { all.equal(facmul(x, nm, y, trans = FALSE, left = TRUE), L %*% y) all.equal(facmul(x, nm, y, trans = FALSE, left = FALSE), y %*% L) all.equal(facmul(x, nm, y, trans = TRUE, left = TRUE), crossprod(L, y)) all.equal(facmul(x, nm, y, trans = TRUE, left = FALSE), tcrossprod(y, L)) }) ## End(Not run)
“Semi-API” functions used internally by Matrix,
often to bypass S4 dispatch and avoid the associated overhead.
These are exported to provide this capability to expert users.
Typical users should continue to rely on S4 generic functions
to dispatch suitable methods, by calling,
e.g., as(., <class>)
for coercions.
.M2kind(from, kind = ".", sparse = NA) .M2gen(from, kind = ".") .M2sym(from, ...) .M2tri(from, ...) .M2diag(from) .M2v(from) .M2m(from) .M2unpacked(from) .M2packed(from) .M2C(from) .M2R(from) .M2T(from) .M2V(from) .m2V(from, kind = ".") .sparse2dense(from, packed = FALSE) .diag2dense(from, kind = ".", shape = "t", packed = FALSE, uplo = "U") .ind2dense(from, kind = "n") .m2dense(from, class = ".ge", uplo = "U", diag = "N", trans = FALSE) .dense2sparse(from, repr = "C") .diag2sparse(from, kind = ".", shape = "t", repr = "C", uplo = "U") .ind2sparse(from, kind = "n", repr = ".") .m2sparse(from, class = ".gC", uplo = "U", diag = "N", trans = FALSE) .tCRT(x, lazy = TRUE) .diag.dsC(x, Chx = Cholesky(x, LDL = TRUE), res.kind = "diag") .solve.dgC.lu (a, b, tol = .Machine$double.eps, check = TRUE) .solve.dgC.qr (a, b, order = 3L, check = TRUE) .solve.dgC.chol(a, b, check = TRUE) .updateCHMfactor(object, parent, mult = 0)
.M2kind(from, kind = ".", sparse = NA) .M2gen(from, kind = ".") .M2sym(from, ...) .M2tri(from, ...) .M2diag(from) .M2v(from) .M2m(from) .M2unpacked(from) .M2packed(from) .M2C(from) .M2R(from) .M2T(from) .M2V(from) .m2V(from, kind = ".") .sparse2dense(from, packed = FALSE) .diag2dense(from, kind = ".", shape = "t", packed = FALSE, uplo = "U") .ind2dense(from, kind = "n") .m2dense(from, class = ".ge", uplo = "U", diag = "N", trans = FALSE) .dense2sparse(from, repr = "C") .diag2sparse(from, kind = ".", shape = "t", repr = "C", uplo = "U") .ind2sparse(from, kind = "n", repr = ".") .m2sparse(from, class = ".gC", uplo = "U", diag = "N", trans = FALSE) .tCRT(x, lazy = TRUE) .diag.dsC(x, Chx = Cholesky(x, LDL = TRUE), res.kind = "diag") .solve.dgC.lu (a, b, tol = .Machine$double.eps, check = TRUE) .solve.dgC.qr (a, b, order = 3L, check = TRUE) .solve.dgC.chol(a, b, check = TRUE) .updateCHMfactor(object, parent, mult = 0)
from , x , a , b
|
a |
kind |
a string ( |
shape |
a string ( |
repr |
a string ( |
packed |
a logical indicating if the result should
inherit from |
sparse |
a logical indicating if the result should inherit
from |
uplo |
a string ( |
diag |
a string ( |
trans |
a logical indicating if the result should be a 1-row
matrix rather than a 1-column matrix where |
class |
a string whose first three characters specify the class
of the result. It should match the pattern
|
... |
optional arguments passed to |
lazy |
a logical indicating if the transpose should be constructed with minimal allocation, but possibly without preserving representation. |
Chx |
optionally, the |
res.kind |
a string in |
tol |
see |
order |
see |
check |
a logical indicating if the first argument should be
tested for inheritance from |
object |
a Cholesky factorization inheriting from virtual class
|
parent |
|
mult |
a numeric vector of postive length. Only the first element is used, and that must be finite. |
Functions with names of the form .<A>2<B>
implement coercions
from virtual class A to the “nearest” non-virtual subclass of
virtual class B, where the virtual classes are abbreviated as follows:
M
V
m
matrix
v
vector
dense
unpacked
packed
sparse
C
R
T
gen
sym
tri
diag
ind
Abbreviations should be seen as a guide, rather than as an
exact description of behaviour. Notably, .m2dense
,
.m2sparse
, and .m2V
accept vectors that are
not matrices.
.tCRT(x)
If lazy = TRUE
, then .tCRT
constructs the transpose
of x
using the most efficient representation,
which for ‘CRT’ is ‘RCT’. If lazy = FALSE
,
then .tCRT
preserves the representation of x
,
behaving as the corresponding methods for generic function t
.
.diag.dsC(x)
.diag.dsC
computes (or uses if Chx
is supplied)
the Cholesky factorization of x
as in order
to calculate one of several possible statistics from the diagonal
entries of
. See
res.kind
under ‘Arguments’.
.solve.dgC.*(a, b)
.solve.dgC.lu(a, b)
needs a square matrix a
.
.solve.dgC.qr(a, b)
needs a “long” matrix a
,
with nrow(a) >= ncol(a)
.
.solve.dgC.chol(a, b)
needs a “wide” matrix a
,
with nrow(a) <= ncol(a)
.
All three may be used to solve sparse linear systems directly.
Only .solve.dgC.qr
and .solve.dgC.chol
be used
to solve sparse least squares problems.
.updateCHMfactor(object, parent, mult)
.updateCHMfactor
updates object
with the result
of Cholesky factorizing
F(parent) + mult[1] * diag(nrow(parent))
,
i.e., F(parent)
plus mult[1]
times the identity matrix,
where F = identity
if parent
is a dsCMatrix
and F = tcrossprod
if parent
is a dgCMatrix
.
The nonzero pattern of F(parent)
must match
that of S
if object = Cholesky(S, ...)
.
D. <- diag(x = c(1, 1, 2, 3, 5, 8)) D.0 <- Diagonal(x = c(0, 0, 0, 3, 5, 8)) S. <- toeplitz(as.double(1:6)) C. <- new("dgCMatrix", Dim = c(3L, 4L), p = c(0L, 1L, 1L, 1L, 3L), i = c(1L, 0L, 2L), x = c(-8, 2, 3)) stopifnot(exprs = { identical(.M2tri (D.), as(D., "triangularMatrix")) identical(.M2sym (D.), as(D., "symmetricMatrix")) identical(.M2diag(D.), as(D., "diagonalMatrix")) identical(.M2kind(C., "l"), as(C., "lMatrix")) identical(.M2kind(.sparse2dense(C.), "l"), as(as(C., "denseMatrix"), "lMatrix")) identical(.diag2sparse(D.0, ".", "t", "C"), .dense2sparse(.diag2dense(D.0, ".", "t", TRUE), "C")) identical(.M2gen(.diag2dense(D.0, ".", "s", FALSE)), .sparse2dense(.M2gen(.diag2sparse(D.0, ".", "s", "T")))) identical(S., .M2m(.m2sparse(S., ".sR"))) identical(S. * lower.tri(S.) + diag(1, 6L), .M2m(.m2dense (S., ".tr", "L", "U"))) identical(.M2R(C.), .M2R(.M2T(C.))) identical(.tCRT(C.), .M2R(t(C.))) }) A <- tcrossprod(C.)/6 + Diagonal(3, 1/3); A[1,2] <- 3; A stopifnot(exprs = { is.numeric( x. <- c(2.2, 0, -1.2) ) all.equal(x., .solve.dgC.lu(A, c(1,0,0), check=FALSE)) all.equal(x., .solve.dgC.qr(A, c(1,0,0), check=FALSE)) }) ## Solving sparse least squares: X <- rbind(A, Diagonal(3)) # design matrix X (for L.S.) Xt <- t(X) # *transposed* X (for L.S.) (y <- drop(crossprod(Xt, 1:3)) + c(-1,1)/1000) # small rand.err. str(solveCh <- .solve.dgC.chol(Xt, y, check=FALSE)) # Xt *is* dgC.. stopifnot(exprs = { all.equal(solveCh$coef, 1:3, tol = 1e-3)# rel.err ~ 1e-4 all.equal(solveCh$coef, drop(solve(tcrossprod(Xt), Xt %*% y))) all.equal(solveCh$coef, .solve.dgC.qr(X, y, check=FALSE)) })
D. <- diag(x = c(1, 1, 2, 3, 5, 8)) D.0 <- Diagonal(x = c(0, 0, 0, 3, 5, 8)) S. <- toeplitz(as.double(1:6)) C. <- new("dgCMatrix", Dim = c(3L, 4L), p = c(0L, 1L, 1L, 1L, 3L), i = c(1L, 0L, 2L), x = c(-8, 2, 3)) stopifnot(exprs = { identical(.M2tri (D.), as(D., "triangularMatrix")) identical(.M2sym (D.), as(D., "symmetricMatrix")) identical(.M2diag(D.), as(D., "diagonalMatrix")) identical(.M2kind(C., "l"), as(C., "lMatrix")) identical(.M2kind(.sparse2dense(C.), "l"), as(as(C., "denseMatrix"), "lMatrix")) identical(.diag2sparse(D.0, ".", "t", "C"), .dense2sparse(.diag2dense(D.0, ".", "t", TRUE), "C")) identical(.M2gen(.diag2dense(D.0, ".", "s", FALSE)), .sparse2dense(.M2gen(.diag2sparse(D.0, ".", "s", "T")))) identical(S., .M2m(.m2sparse(S., ".sR"))) identical(S. * lower.tri(S.) + diag(1, 6L), .M2m(.m2dense (S., ".tr", "L", "U"))) identical(.M2R(C.), .M2R(.M2T(C.))) identical(.tCRT(C.), .M2R(t(C.))) }) A <- tcrossprod(C.)/6 + Diagonal(3, 1/3); A[1,2] <- 3; A stopifnot(exprs = { is.numeric( x. <- c(2.2, 0, -1.2) ) all.equal(x., .solve.dgC.lu(A, c(1,0,0), check=FALSE)) all.equal(x., .solve.dgC.qr(A, c(1,0,0), check=FALSE)) }) ## Solving sparse least squares: X <- rbind(A, Diagonal(3)) # design matrix X (for L.S.) Xt <- t(X) # *transposed* X (for L.S.) (y <- drop(crossprod(Xt, 1:3)) + c(-1,1)/1000) # small rand.err. str(solveCh <- .solve.dgC.chol(Xt, y, check=FALSE)) # Xt *is* dgC.. stopifnot(exprs = { all.equal(solveCh$coef, 1:3, tol = 1e-3)# rel.err ~ 1e-4 all.equal(solveCh$coef, drop(solve(tcrossprod(Xt), Xt %*% y))) all.equal(solveCh$coef, .solve.dgC.qr(X, y, check=FALSE)) })
Force a square matrix x
to a symmetricMatrix
,
without a symmetry check as it would be applied for as(x,
"symmetricMatrix")
.
forceSymmetric(x, uplo)
forceSymmetric(x, uplo)
x |
any square matrix (of numbers), either “"traditional"”
( |
uplo |
optional string, |
a square matrix inheriting from class
symmetricMatrix
.
symmpart
for the symmetric part of a matrix, or
the coercions as(x, <symmetricMatrix class>)
.
## Hilbert matrix i <- 1:6 h6 <- 1/outer(i - 1L, i, "+") sd <- sqrt(diag(h6)) hh <- t(h6/sd)/sd # theoretically symmetric isSymmetric(hh, tol=0) # FALSE; hence try( as(hh, "symmetricMatrix") ) # fails, but this works fine: H6 <- forceSymmetric(hh) ## result can be pretty surprising: (M <- Matrix(1:36, 6)) forceSymmetric(M) # symmetric, hence very different in lower triangle (tm <- tril(M)) forceSymmetric(tm)
## Hilbert matrix i <- 1:6 h6 <- 1/outer(i - 1L, i, "+") sd <- sqrt(diag(h6)) hh <- t(h6/sd)/sd # theoretically symmetric isSymmetric(hh, tol=0) # FALSE; hence try( as(hh, "symmetricMatrix") ) # fails, but this works fine: H6 <- forceSymmetric(hh) ## result can be pretty surprising: (M <- Matrix(1:36, 6)) forceSymmetric(M) # symmetric, hence very different in lower triangle (tm <- tril(M)) forceSymmetric(tm)
Utilities for formatting sparse numeric matrices in a flexible way.
These functions are used by the format
and print
methods for sparse matrices and can be applied as well to standard R
matrices. Note that all arguments but the first are optional.
formatSparseM()
is the main “workhorse” of
formatSpMatrix
, the format
method for sparse
matrices.
.formatSparseSimple()
is a simple helper function, also dealing
with (short/empty) column names construction.
formatSparseM(x, zero.print = ".", align = c("fancy", "right"), m = as(x,"matrix"), asLogical=NULL, uniDiag=NULL, digits=NULL, cx, iN0, dn = dimnames(m)) .formatSparseSimple(m, asLogical=FALSE, digits=NULL, col.names, note.dropping.colnames = TRUE, dn=dimnames(m))
formatSparseM(x, zero.print = ".", align = c("fancy", "right"), m = as(x,"matrix"), asLogical=NULL, uniDiag=NULL, digits=NULL, cx, iN0, dn = dimnames(m)) .formatSparseSimple(m, asLogical=FALSE, digits=NULL, col.names, note.dropping.colnames = TRUE, dn=dimnames(m))
x |
an R object inheriting from class |
zero.print |
character which should be used for
structural zeroes. The default |
align |
a string specifying how the |
m |
(optional) a (standard R) |
asLogical |
should the matrix be formatted as a logical matrix
(or rather as a numeric one); mostly for |
uniDiag |
logical indicating if the diagonal entries of a sparse
unit triangular or unit-diagonal matrix should be formatted as
|
digits |
significant digits to use for printing, see
|
cx |
(optional) character matrix; a formatted version of |
iN0 |
(optional) integer vector, specifying the location of the
non-zeroes of |
col.names , note.dropping.colnames
|
see |
dn |
|
a character matrix like cx
, where the zeros have been replaced
with (padded versions of) zero.print
.
As this is a dense matrix, do not use these functions for
really large (really) sparse matrices!
Martin Maechler
formatSpMatrix
which calls formatSparseM()
and is
the format
method for sparse matrices.printSpMatrix
which is used by the (typically
implicitly called) show
and print
methods
for sparse matrices.
m <- suppressWarnings(matrix(c(0, 3.2, 0,0, 11,0,0,0,0,-7,0), 4,9)) fm <- formatSparseM(m) noquote(fm) ## nice, but this is nicer {with "units" vertically aligned}: print(fm, quote=FALSE, right=TRUE) ## and "the same" as : Matrix(m) ## align = "right" is cheaper --> the "." are not aligned: noquote(f2 <- formatSparseM(m,align="r")) stopifnot(f2 == fm | m == 0, dim(f2) == dim(m), (f2 == ".") == (m == 0))
m <- suppressWarnings(matrix(c(0, 3.2, 0,0, 11,0,0,0,0,-7,0), 4,9)) fm <- formatSparseM(m) noquote(fm) ## nice, but this is nicer {with "units" vertically aligned}: print(fm, quote=FALSE, right=TRUE) ## and "the same" as : Matrix(m) ## align = "right" is cheaper --> the "." are not aligned: noquote(f2 <- formatSparseM(m,align="r")) stopifnot(f2 == fm | m == 0, dim(f2) == dim(m), (f2 == ".") == (m == 0))
Virtual class of “general” matrices; i.e., matrices that do not have a known property such as symmetric, triangular, or diagonal.
A virtual Class: No objects may be created from it.
Dim, Dimnames
inherited from virtual class
Matrix
.
factors
a list of
MatrixFactorization
objects caching
factorizations of the matrix. Typically, it is initialized
as an empty list and updated “automagically” whenever
a factorization is computed.
Class "Matrix"
, directly.
Virtual classes
symmetricMatrix
,
triangularMatrix
, and
diagonalMatrix
.
Generate the n
by n
symmetric Hilbert matrix. Because
these matrices are ill-conditioned for moderate to large n
,
they are often used for testing numerical linear algebra code.
Hilbert(n)
Hilbert(n)
n |
a non-negative integer. |
the n
by n
symmetric Hilbert matrix as a
"dpoMatrix"
object.
the class dpoMatrix
Hilbert(6)
Hilbert(6)
Methods for function image
in package
Matrix. An image of a matrix simply color codes all matrix
entries and draws the matrix using an
grid of (colored) rectangles.
The Matrix package image
methods are based on
levelplot()
from package lattice; hence
these methods return an “object” of class "trellis"
,
producing a graphic when (auto-) print()
ed.
## S4 method for signature 'dgTMatrix' image(x, xlim = c(1, di[2]), ylim = c(di[1], 1), aspect = "iso", sub = sprintf("Dimensions: %d x %d", di[1], di[2]), xlab = "Column", ylab = "Row", cuts = 15, useRaster = FALSE, useAbs = NULL, colorkey = !useAbs, col.regions = NULL, lwd = NULL, border.col = NULL, ...)
## S4 method for signature 'dgTMatrix' image(x, xlim = c(1, di[2]), ylim = c(di[1], 1), aspect = "iso", sub = sprintf("Dimensions: %d x %d", di[1], di[2]), xlab = "Column", ylab = "Row", cuts = 15, useRaster = FALSE, useAbs = NULL, colorkey = !useAbs, col.regions = NULL, lwd = NULL, border.col = NULL, ...)
x |
a Matrix object, i.e., fulfilling |
xlim , ylim
|
x- and y-axis limits; may be used to “zoom
into” matrix. Note that |
aspect |
aspect ratio specified as number (y/x) or string;
see |
sub , xlab , ylab
|
axis annotation with sensible defaults;
see |
cuts |
number of levels the range of matrix values would be divided into. |
useRaster |
logical indicating if raster graphics should be used
(instead of the tradition rectangle vector drawing). If true,
Note that using raster graphics may often be faster, but can be slower, depending on the matrix dimensions and the graphics device (dimensions). |
useAbs |
logical indicating if |
colorkey |
logical indicating if a color key aka ‘legend’
should be produced. Default is to draw one, unless |
col.regions |
vector of gradually varying colors; see
|
lwd |
(only used when |
border.col |
color for the border of each rectangle. |
... |
further arguments passed to methods and
|
as all lattice graphics functions, image(<Matrix>)
returns a "trellis"
object, effectively the result of
levelplot()
.
All methods currently end up calling the method for the
dgTMatrix
class.
Use showMethods(image)
to list them all.
levelplot
, and
print.trellis
from package lattice.
showMethods(image) ## And if you want to see the method definitions: showMethods(image, includeDefs = TRUE, inherited = FALSE) data(CAex, package = "Matrix") image(CAex, main = "image(CAex)") -> imgC; imgC stopifnot(!is.null(leg <- imgC$legend), is.list(leg$right)) # failed for 2 days .. image(CAex, useAbs=TRUE, main = "image(CAex, useAbs=TRUE)") cCA <- Cholesky(crossprod(CAex), Imult = .01) ## See ?print.trellis --- place two image() plots side by side: print(image(cCA, main="Cholesky(crossprod(CAex), Imult = .01)"), split=c(x=1,y=1,nx=2, ny=1), more=TRUE) print(image(cCA, useAbs=TRUE), split=c(x=2,y=1,nx=2,ny=1)) data(USCounties, package = "Matrix") image(USCounties)# huge image(sign(USCounties))## just the pattern # how the result looks, may depend heavily on # the device, screen resolution, antialiasing etc # e.g. x11(type="Xlib") may show very differently than cairo-based ## Drawing borders around each rectangle; # again, viewing depends very much on the device: image(USCounties[1:400,1:200], lwd=.1) ## Using (xlim,ylim) has advantage : matrix dimension and (col/row) indices: image(USCounties, c(1,200), c(1,400), lwd=.1) image(USCounties, c(1,300), c(1,200), lwd=.5 ) image(USCounties, c(1,300), c(1,200), lwd=.01) ## These 3 are all equivalent : (I1 <- image(USCounties, c(1,100), c(1,100), useAbs=FALSE)) I2 <- image(USCounties, c(1,100), c(1,100), useAbs=FALSE, border.col=NA) I3 <- image(USCounties, c(1,100), c(1,100), useAbs=FALSE, lwd=2, border.col=NA) stopifnot(all.equal(I1, I2, check.environment=FALSE), all.equal(I2, I3, check.environment=FALSE)) ## using an opaque border color image(USCounties, c(1,100), c(1,100), useAbs=FALSE, lwd=3, border.col = adjustcolor("skyblue", 1/2)) if(interactive() || nzchar(Sys.getenv("R_MATRIX_CHECK_EXTRA"))) { ## Using raster graphics: For PDF this would give a 77 MB file, ## however, for such a large matrix, this is typically considerably ## *slower* (than vector graphics rectangles) in most cases : if(doPNG <- !dev.interactive()) png("image-USCounties-raster.png", width=3200, height=3200) image(USCounties, useRaster = TRUE) # should not suffer from anti-aliasing if(doPNG) dev.off() ## and now look at the *.png image in a viewer you can easily zoom in and out }#only if(doExtras)
showMethods(image) ## And if you want to see the method definitions: showMethods(image, includeDefs = TRUE, inherited = FALSE) data(CAex, package = "Matrix") image(CAex, main = "image(CAex)") -> imgC; imgC stopifnot(!is.null(leg <- imgC$legend), is.list(leg$right)) # failed for 2 days .. image(CAex, useAbs=TRUE, main = "image(CAex, useAbs=TRUE)") cCA <- Cholesky(crossprod(CAex), Imult = .01) ## See ?print.trellis --- place two image() plots side by side: print(image(cCA, main="Cholesky(crossprod(CAex), Imult = .01)"), split=c(x=1,y=1,nx=2, ny=1), more=TRUE) print(image(cCA, useAbs=TRUE), split=c(x=2,y=1,nx=2,ny=1)) data(USCounties, package = "Matrix") image(USCounties)# huge image(sign(USCounties))## just the pattern # how the result looks, may depend heavily on # the device, screen resolution, antialiasing etc # e.g. x11(type="Xlib") may show very differently than cairo-based ## Drawing borders around each rectangle; # again, viewing depends very much on the device: image(USCounties[1:400,1:200], lwd=.1) ## Using (xlim,ylim) has advantage : matrix dimension and (col/row) indices: image(USCounties, c(1,200), c(1,400), lwd=.1) image(USCounties, c(1,300), c(1,200), lwd=.5 ) image(USCounties, c(1,300), c(1,200), lwd=.01) ## These 3 are all equivalent : (I1 <- image(USCounties, c(1,100), c(1,100), useAbs=FALSE)) I2 <- image(USCounties, c(1,100), c(1,100), useAbs=FALSE, border.col=NA) I3 <- image(USCounties, c(1,100), c(1,100), useAbs=FALSE, lwd=2, border.col=NA) stopifnot(all.equal(I1, I2, check.environment=FALSE), all.equal(I2, I3, check.environment=FALSE)) ## using an opaque border color image(USCounties, c(1,100), c(1,100), useAbs=FALSE, lwd=3, border.col = adjustcolor("skyblue", 1/2)) if(interactive() || nzchar(Sys.getenv("R_MATRIX_CHECK_EXTRA"))) { ## Using raster graphics: For PDF this would give a 77 MB file, ## however, for such a large matrix, this is typically considerably ## *slower* (than vector graphics rectangles) in most cases : if(doPNG <- !dev.interactive()) png("image-USCounties-raster.png", width=3200, height=3200) image(USCounties, useRaster = TRUE) # should not suffer from anti-aliasing if(doPNG) dev.off() ## and now look at the *.png image in a viewer you can easily zoom in and out }#only if(doExtras)
Class index
is a virtual class designating index vectors,
or “subscripts”, for (possibly named) vectors and arrays.
It is typically used in signatures of methods for the subscript
and subassignment operators, namely [
and [<-
.
It is implemented as a union of the atomic vector classes
numeric
, logical
,
and character
.
[
, [-methods
, and
[<–methods
.
showClass("index")
showClass("index")
The indMatrix
class is the class of row and column
index matrices, stored as 1-based integer index vectors.
A row (column) index matrix is a matrix whose rows (columns)
are standard unit vectors. Such matrices are useful
when mapping observations to discrete sets of covariate values.
Multiplying a matrix on the left by a row index matrix is equivalent to indexing its rows, i.e., sampling the rows “with replacement”. Analogously, multiplying a matrix on the right by a column index matrix is equivalent to indexing its columns. Indeed, such products are implemented in Matrix as indexing operations; see ‘Details’ below.
A matrix whose rows and columns are standard unit vectors
is called a permutation matrix. This special case is
designated by the pMatrix
class, a direct
subclass of indMatrix
.
The transpose of an index matrix is an index matrix with identical
perm
but opposite margin
. Hence the transpose of a
row index matrix is a column index matrix, and vice versa.
The cross product of a row index matrix R
and itself is a
diagonal matrix whose diagonal entries are the the number of entries
in each column of R
.
Given a row index matrix R
with perm
slot p
,
a column index matrix C
with perm
slot q
,
and a matrix M
with conformable dimensions, we have
|
= | R %*% M |
= | M[p, ]
|
|
= | M %*% C |
= | M[, q]
|
|
= | crossprod(C, M) |
= | M[q, ]
|
|
= | tcrossprod(M, R) |
= | M[, p]
|
|
= | crossprod(R) |
= | Diagonal(x=tabulate(p, ncol(R)))
|
|
= | tcrossprod(C) |
= | Diagonal(x=tabulate(q, nrow(C)))
|
Operations on index matrices that result in index matrices will
accordingly return an indMatrix
. These include products
of two column index matrices and (equivalently) column-indexing
of a column index matrix (when dimensions are not dropped).
Most other operations on indMatrix
treat them as sparse
nonzero pattern matrices (i.e., inheriting from virtual class
nsparseMatrix
). Hence vector-valued subsets
of indMatrix
, such as those given by diag
,
are always of type "logical"
.
Objects can be created explicitly with calls of the form
new("indMatrix", ...)
, but they are more commonly created
by coercing 1-based integer index vectors, with calls of the
form as(., "indMatrix")
; see ‘Methods’ below.
margin
an integer, either 1 or 2, specifying whether the matrix is a row (1) or column (2) index.
perm
a 1-based integer index vector, i.e.,
a vector of length Dim[margin]
with elements
taken from 1:Dim[1+margin%%2]
.
Dim
,Dimnames
inherited from virtual
superclass Matrix
.
Classes "sparseMatrix"
and
"generalMatrix"
, directly.
%*%
signature(x = "indMatrix", y = "Matrix")
and others listed by showMethods("%*%", classes = "indMatrix")
:
matrix products implemented where appropriate as indexing operations.
coerce
signature(from = "numeric", to = "indMatrix")
:
supporting typical indMatrix
construction from
a vector of positive integers. Row indexing is assumed.
coerce
signature(from = "list", to = "indMatrix")
:
supporting indMatrix
construction for row and
column indexing, including index vectors of length 0 and
index vectors whose maximum is less than the number of rows
or columns being indexed.
coerce
signature(from = "indMatrix", to = "matrix")
:
coercion to a traditional matrix
of logical type,
with FALSE
and TRUE
in place of 0 and 1.
t
signature(x = "indMatrix")
:
the transpose, which is an indMatrix
with identical
perm
but opposite margin
.
rowSums
,rowMeans
,colSums
,colMeans
signature(x = "indMatrix")
:
row and column sums and means.
rbind2
,cbind2
signature(x = "indMatrix", y = "indMatrix")
:
row-wise catenation of two row index matrices with equal numbers
of columns and column-wise catenation of two column index matrices
with equal numbers of rows.
signature(X = "indMatrix", Y = "indMatrix")
:
Kronecker product of two row index matrices or two column index
matrices, giving the row or column index matrix corresponding to
their “interaction”.
Fabian Scheipl at ‘uni-muenchen.de’, building on the existing class
pMatrix
after a nice hike's conversation with
Martin Maechler. Methods for crossprod(x, y)
and
kronecker(x, y)
with both arguments inheriting from
indMatrix
were made considerably faster thanks to a suggestion
by Boris Vaillant. Diverse tweaks by Martin Maechler and
Mikael Jagan, notably the latter's implementation of margin
,
prior to which the indMatrix
class was designated only for
row index matrices.
Subclass pMatrix
of permutation matrices,
a special case of index matrices; virtual class
nMatrix
of nonzero pattern matrices,
and its subclasses.
p1 <- as(c(2,3,1), "pMatrix") (sm1 <- as(rep(c(2,3,1), e=3), "indMatrix")) stopifnot(all(sm1 == p1[rep(1:3, each=3),])) ## row-indexing of a <pMatrix> turns it into an <indMatrix>: class(p1[rep(1:3, each=3),]) set.seed(12) # so we know '10' is in sample ## random index matrix for 30 observations and 10 unique values: (s10 <- as(sample(10, 30, replace=TRUE),"indMatrix")) ## Sample rows of a numeric matrix : (mm <- matrix(1:10, nrow=10, ncol=3)) s10 %*% mm set.seed(27) IM1 <- as(sample(1:20, 100, replace=TRUE), "indMatrix") IM2 <- as(sample(1:18, 100, replace=TRUE), "indMatrix") (c12 <- crossprod(IM1,IM2)) ## same as cross-tabulation of the two index vectors: stopifnot(all(c12 - unclass(table(IM1@perm, IM2@perm)) == 0)) # 3 observations, 4 implied values, first does not occur in sample: as(2:4, "indMatrix") # 3 observations, 5 values, first and last do not occur in sample: as(list(2:4, 5), "indMatrix") as(sm1, "nMatrix") s10[1:7, 1:4] # gives an "ngTMatrix" (most economic!) s10[1:4, ] # preserves "indMatrix"-class I1 <- as(c(5:1,6:4,7:3), "indMatrix") I2 <- as(7:1, "pMatrix") (I12 <- rbind(I1, I2)) stopifnot(is(I12, "indMatrix"), identical(I12, rbind(I1, I2)), colSums(I12) == c(2L,2:4,4:2))
p1 <- as(c(2,3,1), "pMatrix") (sm1 <- as(rep(c(2,3,1), e=3), "indMatrix")) stopifnot(all(sm1 == p1[rep(1:3, each=3),])) ## row-indexing of a <pMatrix> turns it into an <indMatrix>: class(p1[rep(1:3, each=3),]) set.seed(12) # so we know '10' is in sample ## random index matrix for 30 observations and 10 unique values: (s10 <- as(sample(10, 30, replace=TRUE),"indMatrix")) ## Sample rows of a numeric matrix : (mm <- matrix(1:10, nrow=10, ncol=3)) s10 %*% mm set.seed(27) IM1 <- as(sample(1:20, 100, replace=TRUE), "indMatrix") IM2 <- as(sample(1:18, 100, replace=TRUE), "indMatrix") (c12 <- crossprod(IM1,IM2)) ## same as cross-tabulation of the two index vectors: stopifnot(all(c12 - unclass(table(IM1@perm, IM2@perm)) == 0)) # 3 observations, 4 implied values, first does not occur in sample: as(2:4, "indMatrix") # 3 observations, 5 values, first and last do not occur in sample: as(list(2:4, 5), "indMatrix") as(sm1, "nMatrix") s10[1:7, 1:4] # gives an "ngTMatrix" (most economic!) s10[1:4, ] # preserves "indMatrix"-class I1 <- as(c(5:1,6:4,7:3), "indMatrix") I2 <- as(7:1, "pMatrix") (I12 <- rbind(I1, I2)) stopifnot(is(I12, "indMatrix"), identical(I12, rbind(I1, I2)), colSums(I12) == c(2L,2:4,4:2))
invertPerm
and signPerm
compute the inverse and sign
of a length-n
permutation vector. isPerm
tests
if a length-n
integer vector is a valid permutation vector.
asPerm
coerces a length-m
transposition vector to a
length-n
permutation vector, where m <= n
.
invertPerm(p, off = 1L, ioff = 1L) signPerm(p, off = 1L) isPerm(p, off = 1L) asPerm(pivot, off = 1L, ioff = 1L, n = length(pivot)) invPerm(p, zero.p = FALSE, zero.res = FALSE)
invertPerm(p, off = 1L, ioff = 1L) signPerm(p, off = 1L) isPerm(p, off = 1L) asPerm(pivot, off = 1L, ioff = 1L, n = length(pivot)) invPerm(p, zero.p = FALSE, zero.res = FALSE)
p |
an integer vector of length |
pivot |
an integer vector of length |
off |
an integer offset, indicating that |
ioff |
an integer offset, indicating that the result
should be a permutation of |
n |
a integer greater than or equal to |
zero.p |
a logical. Equivalent to |
zero.res |
a logical. Equivalent to |
invertPerm(p, off, ioff=1)
is equivalent to
order(p)
or sort.list(p)
for all values of off
. For the default value
off=1
, it returns the value of p
after
p[p] <- seq_along(p)
.
invPerm
is a simple wrapper around invertPerm
,
retained for backwards compatibility.
By default, i.e., with off=1
and ioff=1
:
invertPerm(p)
returns an integer vector of length
length(p)
such that p[invertPerm(p)]
and invertPerm(p)[p]
are both seq_along(p)
,
i.e., the identity permutation.
signPerm(p)
returns 1 if p
is an even permutation
and -1
otherwise (i.e., if p
is odd).
isPerm(p)
returns TRUE
if p
is a
permutation of seq_along(p)
and FALSE
otherwise.
asPerm(pivot)
returns the result of transposing elements
i
and pivot[i]
of a permutation vector initialized
as seq_len(n)
, for i
in seq_along(pivot)
.
Class pMatrix
of permutation matrices.
p <- sample(10L) # a random permutation vector ip <- invertPerm(p) s <- signPerm(p) ## 'p' and 'ip' are indeed inverses: stopifnot(exprs = { isPerm(p) isPerm(ip) identical(s, 1L) || identical(s, -1L) identical(s, signPerm(ip)) identical(p[ip], 1:10) identical(ip[p], 1:10) identical(invertPerm(ip), p) }) ## Product of transpositions (1 2)(2 1)(4 3)(6 8)(10 1) = (3 4)(6 8)(1 10) pivot <- c(2L, 1L, 3L, 3L, 5L, 8L, 7L, 8L, 9L, 1L) q <- asPerm(pivot) stopifnot(exprs = { identical(q, c(10L, 2L, 4L, 3L, 5L, 8L, 7L, 6L, 9L, 1L)) identical(q[q], seq_len(10L)) # because the permutation is odd: signPerm(q) == -1L }) invPerm # a less general version of 'invertPerm'
p <- sample(10L) # a random permutation vector ip <- invertPerm(p) s <- signPerm(p) ## 'p' and 'ip' are indeed inverses: stopifnot(exprs = { isPerm(p) isPerm(ip) identical(s, 1L) || identical(s, -1L) identical(s, signPerm(ip)) identical(p[ip], 1:10) identical(ip[p], 1:10) identical(invertPerm(ip), p) }) ## Product of transpositions (1 2)(2 1)(4 3)(6 8)(10 1) = (3 4)(6 8)(1 10) pivot <- c(2L, 1L, 3L, 3L, 5L, 8L, 7L, 8L, 9L, 1L) q <- asPerm(pivot) stopifnot(exprs = { identical(q, c(10L, 2L, 4L, 3L, 5L, 8L, 7L, 6L, 9L, 1L)) identical(q[q], seq_len(10L)) # because the permutation is odd: signPerm(q) == -1L }) invPerm # a less general version of 'invertPerm'
Methods for generic functions anyNA()
,
is.na()
, is.nan()
,
is.infinite()
, and is.finite()
,
for objects inheriting from virtual class
Matrix
or sparseVector
.
## S4 method for signature 'denseMatrix' is.na(x) ## S4 method for signature 'sparseMatrix' is.na(x) ## S4 method for signature 'diagonalMatrix' is.na(x) ## S4 method for signature 'indMatrix' is.na(x) ## S4 method for signature 'sparseVector' is.na(x) ## ... ## and likewise for anyNA, is.nan, is.infinite, is.finite
## S4 method for signature 'denseMatrix' is.na(x) ## S4 method for signature 'sparseMatrix' is.na(x) ## S4 method for signature 'diagonalMatrix' is.na(x) ## S4 method for signature 'indMatrix' is.na(x) ## S4 method for signature 'sparseVector' is.na(x) ## ... ## and likewise for anyNA, is.nan, is.infinite, is.finite
x |
an R object, here a sparse or dense matrix or vector. |
For is.*()
, an nMatrix
or
nsparseVector
matching the dimensions
of x
and specifying the positions in x
of
(some subset of) NA
, NaN
,
Inf
, and -Inf
.
For anyNA
(), TRUE
if x
contains NA
or NaN
and FALSE
otherwise.
(M <- Matrix(1:6, nrow = 4, ncol = 3, dimnames = list(letters[1:4], LETTERS[1:3]))) stopifnot(!anyNA(M), !any(is.na(M))) M[2:3, 2] <- NA (inM <- is.na(M)) stopifnot(anyNA(M), sum(inM) == 2) (A <- spMatrix(nrow = 10, ncol = 20, i = c(1, 3:8), j = c(2, 9, 6:10), x = 7 * (1:7))) stopifnot(!anyNA(A), !any(is.na(A))) A[2, 3] <- A[1, 2] <- A[5, 5:9] <- NA (inA <- is.na(A)) stopifnot(anyNA(A), sum(inA) == 1 + 1 + 5)
(M <- Matrix(1:6, nrow = 4, ncol = 3, dimnames = list(letters[1:4], LETTERS[1:3]))) stopifnot(!anyNA(M), !any(is.na(M))) M[2:3, 2] <- NA (inM <- is.na(M)) stopifnot(anyNA(M), sum(inM) == 2) (A <- spMatrix(nrow = 10, ncol = 20, i = c(1, 3:8), j = c(2, 9, 6:10), x = 7 * (1:7))) stopifnot(!anyNA(A), !any(is.na(A))) A[2, 3] <- A[1, 2] <- A[5, 5:9] <- NA (inA <- is.na(A)) stopifnot(anyNA(A), sum(inA) == 1 + 1 + 5)
dn
NULL-like ?Are the dimnames
dn
NULL
-like?
is.null.DN(dn)
is less strict than is.null(dn)
,
because it is also true (TRUE
) when the dimnames
dn
are “like” NULL
, or list(NULL,NULL)
, as
they can easily be for the traditional R matrices
(matrix
) which have no formal class
definition, and hence much freedom in how their dimnames
look like.
is.null.DN(dn)
is.null.DN(dn)
dn |
This function is really to be used on “traditional” matrices
rather than those inheriting from Matrix
, as
the latter will always have dimnames list(NULL,NULL)
exactly,
in such a case.
Martin Maechler
m1 <- m2 <- m3 <- m4 <- m <- matrix(round(100 * rnorm(6)), 2, 3) dimnames(m1) <- list(NULL, NULL) dimnames(m2) <- list(NULL, character()) dimnames(m3) <- rev(dimnames(m2)) dimnames(m4) <- rep(list(character()),2) m4 # prints absolutely identically to m c.o <- capture.output cm <- c.o(m) stopifnot(exprs = { m == m1; m == m2; m == m3; m == m4 identical(cm, c.o(m1)); identical(cm, c.o(m2)) identical(cm, c.o(m3)); identical(cm, c.o(m4)) }) hasNoDimnames <- function(.) is.null.DN(dimnames(.)) stopifnot(exprs = { hasNoDimnames(m) hasNoDimnames(m1); hasNoDimnames(m2) hasNoDimnames(m3); hasNoDimnames(m4) hasNoDimnames(Matrix(m) -> M) hasNoDimnames(as(M, "sparseMatrix")) })
m1 <- m2 <- m3 <- m4 <- m <- matrix(round(100 * rnorm(6)), 2, 3) dimnames(m1) <- list(NULL, NULL) dimnames(m2) <- list(NULL, character()) dimnames(m3) <- rev(dimnames(m2)) dimnames(m4) <- rep(list(character()),2) m4 # prints absolutely identically to m c.o <- capture.output cm <- c.o(m) stopifnot(exprs = { m == m1; m == m2; m == m3; m == m4 identical(cm, c.o(m1)); identical(cm, c.o(m2)) identical(cm, c.o(m3)); identical(cm, c.o(m4)) }) hasNoDimnames <- function(.) is.null.DN(dimnames(.)) stopifnot(exprs = { hasNoDimnames(m) hasNoDimnames(m1); hasNoDimnames(m2) hasNoDimnames(m3); hasNoDimnames(m4) hasNoDimnames(Matrix(m) -> M) hasNoDimnames(as(M, "sparseMatrix")) })
isSymmetric
tests whether its argument is a symmetric square
matrix, by default tolerating some numerical fuzz and requiring
symmetric [dD]imnames
in addition to symmetry in the
mathematical sense. isSymmetric
is a generic function in
base, which has a method for traditional
matrices of implicit class
"matrix"
.
Methods are defined here for various proper and virtual classes
in Matrix, so that isSymmetric
works for all objects
inheriting from virtual class "Matrix"
.
## S4 method for signature 'denseMatrix' isSymmetric(object, checkDN = TRUE, ...) ## S4 method for signature 'CsparseMatrix' isSymmetric(object, checkDN = TRUE, ...) ## S4 method for signature 'RsparseMatrix' isSymmetric(object, checkDN = TRUE, ...) ## S4 method for signature 'TsparseMatrix' isSymmetric(object, checkDN = TRUE, ...) ## S4 method for signature 'diagonalMatrix' isSymmetric(object, checkDN = TRUE, ...) ## S4 method for signature 'indMatrix' isSymmetric(object, checkDN = TRUE, ...) ## S4 method for signature 'dgeMatrix' isSymmetric(object, checkDN = TRUE, tol = 100 * .Machine$double.eps, tol1 = 8 * tol, ...) ## S4 method for signature 'dgCMatrix' isSymmetric(object, checkDN = TRUE, tol = 100 * .Machine$double.eps, ...)
## S4 method for signature 'denseMatrix' isSymmetric(object, checkDN = TRUE, ...) ## S4 method for signature 'CsparseMatrix' isSymmetric(object, checkDN = TRUE, ...) ## S4 method for signature 'RsparseMatrix' isSymmetric(object, checkDN = TRUE, ...) ## S4 method for signature 'TsparseMatrix' isSymmetric(object, checkDN = TRUE, ...) ## S4 method for signature 'diagonalMatrix' isSymmetric(object, checkDN = TRUE, ...) ## S4 method for signature 'indMatrix' isSymmetric(object, checkDN = TRUE, ...) ## S4 method for signature 'dgeMatrix' isSymmetric(object, checkDN = TRUE, tol = 100 * .Machine$double.eps, tol1 = 8 * tol, ...) ## S4 method for signature 'dgCMatrix' isSymmetric(object, checkDN = TRUE, tol = 100 * .Machine$double.eps, ...)
object |
a |
checkDN |
a logical indicating whether symmetry of the
|
tol , tol1
|
numerical tolerances allowing approximate
symmetry of numeric (rather than logical) matrices. See also
|
... |
further arguments passed to methods
(typically methods for |
The Dimnames
slot of object
, say dn
,
is considered to be symmetric if and only if
dn[[1]]
and dn[[2]]
are identical or
one is NULL
; and
ndn <- names(dn)
is NULL
or
ndn[1]
and ndn[2]
are identical or
one is the empty string ""
.
Hence list(a=nms, a=nms)
is considered to be symmetric,
and so too are list(a=nms, NULL)
and list(NULL, a=nms)
.
Note that this definition is looser than that employed by
isSymmetric.matrix
, which requires dn[1]
and
dn[2]
to be identical, where dn
is the dimnames
attribute of a traditional matrix.
A logical, either TRUE
or FALSE
(never NA
).
forceSymmetric
;
symmpart
and skewpart
;
virtual class "symmetricMatrix"
and its subclasses.
isSymmetric(Diagonal(4)) # TRUE of course M <- Matrix(c(1,2,2,1), 2,2) isSymmetric(M) # TRUE (*and* of formal class "dsyMatrix") isSymmetric(as(M, "generalMatrix")) # still symmetric, even if not "formally" isSymmetric(triu(M)) # FALSE ## Look at implementations: showMethods("isSymmetric", includeDefs = TRUE) # includes S3 generic from base
isSymmetric(Diagonal(4)) # TRUE of course M <- Matrix(c(1,2,2,1), 2,2) isSymmetric(M) # TRUE (*and* of formal class "dsyMatrix") isSymmetric(as(M, "generalMatrix")) # still symmetric, even if not "formally" isSymmetric(triu(M)) # FALSE ## Look at implementations: showMethods("isSymmetric", includeDefs = TRUE) # includes S3 generic from base
isTriangular
and isDiagonal
test whether their argument
is a triangular or diagonal matrix, respectively. Unlike the analogous
isSymmetric
, these two functions are generically
from Matrix rather than base
. Hence Matrix
defines methods for traditional matrices of implicit class
"matrix"
in addition to matrices inheriting from
virtual class "Matrix"
.
By our definition, triangular and diagonal matrices are square, i.e., they have the same number of rows and columns.
isTriangular(object, upper = NA, ...) isDiagonal(object)
isTriangular(object, upper = NA, ...) isDiagonal(object)
object |
an R object, typically a matrix. |
upper |
a logical, either |
... |
further arguments passed to methods (currently unused by Matrix). |
A logical, either TRUE
or FALSE
(never NA
).
If object
is triangular and upper
is NA
, then
isTriangular
returns TRUE
with an attribute
kind
, either "U"
or "L"
, indicating that
object
is upper or lower triangular, respectively.
Users should not rely on how kind
is determined for diagonal
matrices, which are both upper and lower triangular.
isSymmetric
;
virtual classes "triangularMatrix"
and
"diagonalMatrix"
and their subclasses.
isTriangular(Diagonal(4)) ## is TRUE: a diagonal matrix is also (both upper and lower) triangular (M <- Matrix(c(1,2,0,1), 2,2)) isTriangular(M) # TRUE (*and* of formal class "dtrMatrix") isTriangular(as(M, "generalMatrix")) # still triangular, even if not "formally" isTriangular(crossprod(M)) # FALSE isDiagonal(matrix(c(2,0,0,1), 2,2)) # TRUE ## Look at implementations: showMethods("isTriangular", includeDefs = TRUE) showMethods("isDiagonal", includeDefs = TRUE)
isTriangular(Diagonal(4)) ## is TRUE: a diagonal matrix is also (both upper and lower) triangular (M <- Matrix(c(1,2,0,1), 2,2)) isTriangular(M) # TRUE (*and* of formal class "dtrMatrix") isTriangular(as(M, "generalMatrix")) # still triangular, even if not "formally" isTriangular(crossprod(M)) # FALSE isDiagonal(matrix(c(2,0,0,1), 2,2)) # TRUE ## Look at implementations: showMethods("isTriangular", includeDefs = TRUE) showMethods("isDiagonal", includeDefs = TRUE)
Computes Khatri-Rao products for any kind of matrices.
The Khatri-Rao product is a column-wise Kronecker product. Originally introduced by Khatri and Rao (1968), it has many different applications, see Liu and Trenkler (2008) for a survey. Notably, it is used in higher-dimensional tensor decompositions, see Bader and Kolda (2008).
KhatriRao(X, Y = X, FUN = "*", sparseY = TRUE, make.dimnames = FALSE)
KhatriRao(X, Y = X, FUN = "*", sparseY = TRUE, make.dimnames = FALSE)
X , Y
|
matrices of with the same number of columns. |
FUN |
the (name of the) |
sparseY |
logical specifying if |
make.dimnames |
logical indicating if the result should inherit
|
a "CsparseMatrix"
, say R
, the Khatri-Rao
product of X
() and
Y
(), is of dimension
,
where the j-th column,
R[,j]
is the kronecker product
kronecker(X[,j], Y[,j])
.
The current implementation is efficient for large sparse matrices.
Original by Michael Cysouw, Univ. Marburg; minor tweaks, bug fixes etc, by Martin Maechler.
Khatri, C. G., and Rao, C. Radhakrishna (1968) Solutions to Some Functional Equations and Their Applications to Characterization of Probability Distributions. Sankhya: Indian J. Statistics, Series A 30, 167–180.
Bader, Brett W, and Tamara G Kolda (2008) Efficient MATLAB Computations with Sparse and Factored Tensors. SIAM J. Scientific Computing 30, 205–231.
## Example with very small matrices: m <- matrix(1:12,3,4) d <- diag(1:4) KhatriRao(m,d) KhatriRao(d,m) dimnames(m) <- list(LETTERS[1:3], letters[1:4]) KhatriRao(m,d, make.dimnames=TRUE) KhatriRao(d,m, make.dimnames=TRUE) dimnames(d) <- list(NULL, paste0("D", 1:4)) KhatriRao(m,d, make.dimnames=TRUE) KhatriRao(d,m, make.dimnames=TRUE) dimnames(d) <- list(paste0("d", 10*1:4), paste0("D", 1:4)) (Kmd <- KhatriRao(m,d, make.dimnames=TRUE)) (Kdm <- KhatriRao(d,m, make.dimnames=TRUE)) nm <- as(m, "nsparseMatrix") nd <- as(d, "nsparseMatrix") KhatriRao(nm,nd, make.dimnames=TRUE) KhatriRao(nd,nm, make.dimnames=TRUE) stopifnot(dim(KhatriRao(m,d)) == c(nrow(m)*nrow(d), ncol(d))) ## border cases / checks: zm <- nm; zm[] <- FALSE # all FALSE matrix stopifnot(all(K1 <- KhatriRao(nd, zm) == 0), identical(dim(K1), c(12L, 4L)), all(K2 <- KhatriRao(zm, nd) == 0), identical(dim(K2), c(12L, 4L))) d0 <- d; d0[] <- 0; m0 <- Matrix(d0[-1,]) stopifnot(all(K3 <- KhatriRao(d0, m) == 0), identical(dim(K3), dim(Kdm)), all(K4 <- KhatriRao(m, d0) == 0), identical(dim(K4), dim(Kmd)), all(KhatriRao(d0, d0) == 0), all(KhatriRao(m0, d0) == 0), all(KhatriRao(d0, m0) == 0), all(KhatriRao(m0, m0) == 0), identical(dimnames(KhatriRao(m, d0, make.dimnames=TRUE)), dimnames(Kmd))) ## a matrix with "structural" and non-structural zeros: m01 <- new("dgCMatrix", i = c(0L, 2L, 0L, 1L), p = c(0L, 0L, 0L, 2L, 4L), Dim = 3:4, x = c(1, 0, 1, 0)) D4 <- Diagonal(4, x=1:4) # "as" d DU <- Diagonal(4)# unit-diagonal: uplo="U" (K5 <- KhatriRao( d, m01)) K5d <- KhatriRao( d, m01, sparseY=FALSE) K5Dd <- KhatriRao(D4, m01, sparseY=FALSE) K5Ud <- KhatriRao(DU, m01, sparseY=FALSE) (K6 <- KhatriRao(diag(3), t(m01))) K6D <- KhatriRao(Diagonal(3), t(m01)) K6d <- KhatriRao(diag(3), t(m01), sparseY=FALSE) K6Dd <- KhatriRao(Diagonal(3), t(m01), sparseY=FALSE) stopifnot(exprs = { all(K5 == K5d) identical(cbind(c(7L, 10L), c(3L, 4L)), which(K5 != 0, arr.ind = TRUE, useNames=FALSE)) identical(K5d, K5Dd) identical(K6, K6D) all(K6 == K6d) identical(cbind(3:4, 1L), which(K6 != 0, arr.ind = TRUE, useNames=FALSE)) identical(K6d, K6Dd) })
## Example with very small matrices: m <- matrix(1:12,3,4) d <- diag(1:4) KhatriRao(m,d) KhatriRao(d,m) dimnames(m) <- list(LETTERS[1:3], letters[1:4]) KhatriRao(m,d, make.dimnames=TRUE) KhatriRao(d,m, make.dimnames=TRUE) dimnames(d) <- list(NULL, paste0("D", 1:4)) KhatriRao(m,d, make.dimnames=TRUE) KhatriRao(d,m, make.dimnames=TRUE) dimnames(d) <- list(paste0("d", 10*1:4), paste0("D", 1:4)) (Kmd <- KhatriRao(m,d, make.dimnames=TRUE)) (Kdm <- KhatriRao(d,m, make.dimnames=TRUE)) nm <- as(m, "nsparseMatrix") nd <- as(d, "nsparseMatrix") KhatriRao(nm,nd, make.dimnames=TRUE) KhatriRao(nd,nm, make.dimnames=TRUE) stopifnot(dim(KhatriRao(m,d)) == c(nrow(m)*nrow(d), ncol(d))) ## border cases / checks: zm <- nm; zm[] <- FALSE # all FALSE matrix stopifnot(all(K1 <- KhatriRao(nd, zm) == 0), identical(dim(K1), c(12L, 4L)), all(K2 <- KhatriRao(zm, nd) == 0), identical(dim(K2), c(12L, 4L))) d0 <- d; d0[] <- 0; m0 <- Matrix(d0[-1,]) stopifnot(all(K3 <- KhatriRao(d0, m) == 0), identical(dim(K3), dim(Kdm)), all(K4 <- KhatriRao(m, d0) == 0), identical(dim(K4), dim(Kmd)), all(KhatriRao(d0, d0) == 0), all(KhatriRao(m0, d0) == 0), all(KhatriRao(d0, m0) == 0), all(KhatriRao(m0, m0) == 0), identical(dimnames(KhatriRao(m, d0, make.dimnames=TRUE)), dimnames(Kmd))) ## a matrix with "structural" and non-structural zeros: m01 <- new("dgCMatrix", i = c(0L, 2L, 0L, 1L), p = c(0L, 0L, 0L, 2L, 4L), Dim = 3:4, x = c(1, 0, 1, 0)) D4 <- Diagonal(4, x=1:4) # "as" d DU <- Diagonal(4)# unit-diagonal: uplo="U" (K5 <- KhatriRao( d, m01)) K5d <- KhatriRao( d, m01, sparseY=FALSE) K5Dd <- KhatriRao(D4, m01, sparseY=FALSE) K5Ud <- KhatriRao(DU, m01, sparseY=FALSE) (K6 <- KhatriRao(diag(3), t(m01))) K6D <- KhatriRao(Diagonal(3), t(m01)) K6d <- KhatriRao(diag(3), t(m01), sparseY=FALSE) K6Dd <- KhatriRao(Diagonal(3), t(m01), sparseY=FALSE) stopifnot(exprs = { all(K5 == K5d) identical(cbind(c(7L, 10L), c(3L, 4L)), which(K5 != 0, arr.ind = TRUE, useNames=FALSE)) identical(K5d, K5Dd) identical(K6, K6D) all(K6 == K6d) identical(cbind(3:4, 1L), which(K6 != 0, arr.ind = TRUE, useNames=FALSE)) identical(K6d, K6Dd) })
A model matrix mm
and corresponding response vector y
used in an example by Koenker and Ng. The matrix mm
is a sparse
matrix with 1850 rows and 712 columns but only 8758 non-zero entries.
It is a "dgCMatrix"
object. The vector y
is just
numeric
of length 1850.
data(KNex)
data(KNex)
Roger Koenker and Pin Ng (2003). SparseM: A sparse matrix package for R; J. of Statistical Software, 8 (6), doi:10.18637/jss.v008.i06
data(KNex, package = "Matrix") class(KNex$mm) dim(KNex$mm) image(KNex$mm) str(KNex) system.time( # a fraction of a second sparse.sol <- with(KNex, solve(crossprod(mm), crossprod(mm, y)))) head(round(sparse.sol,3)) ## Compare with QR-based solution ("more accurate, but slightly slower"): system.time( sp.sol2 <- with(KNex, qr.coef(qr(mm), y) )) all.equal(sparse.sol, sp.sol2, tolerance = 1e-13) # TRUE
data(KNex, package = "Matrix") class(KNex$mm) dim(KNex$mm) image(KNex$mm) str(KNex) system.time( # a fraction of a second sparse.sol <- with(KNex, solve(crossprod(mm), crossprod(mm, y)))) head(round(sparse.sol,3)) ## Compare with QR-based solution ("more accurate, but slightly slower"): system.time( sp.sol2 <- with(KNex, qr.coef(qr(mm), y) )) all.equal(sparse.sol, sp.sol2, tolerance = 1e-13) # TRUE
Computes Kronecker products for objects inheriting from
"Matrix"
.
In order to preserver sparseness, we treat 0 * NA
as 0
,
not as NA
as usually in R (and as used for the
base function kronecker
).
signature(X = "Matrix", Y = "ANY")
.......
signature(X = "ANY", Y = "Matrix")
.......
signature(X = "diagonalMatrix", Y = "ANY")
.......
signature(X = "sparseMatrix", Y = "ANY")
.......
signature(X = "TsparseMatrix", Y = "TsparseMatrix")
.......
signature(X = "dgTMatrix", Y = "dgTMatrix")
.......
signature(X = "dtTMatrix", Y = "dtTMatrix")
.......
signature(X = "indMatrix", Y = "indMatrix")
.......
(t1 <- spMatrix(5,4, x= c(3,2,-7,11), i= 1:4, j=4:1)) # 5 x 4 (t2 <- kronecker(Diagonal(3, 2:4), t1)) # 15 x 12 ## should also work with special-cased logical matrices l3 <- upper.tri(matrix(,3,3)) M <- Matrix(l3) (N <- as(M, "nsparseMatrix")) # "ntCMatrix" (upper triangular) N2 <- as(N, "generalMatrix") # (lost "t"riangularity) MM <- kronecker(M,M) NN <- kronecker(N,N) # "dtTMatrix" i.e. did keep NN2 <- kronecker(N2,N2) stopifnot(identical(NN,MM), is(NN2, "sparseMatrix"), all(NN2 == NN), is(NN, "triangularMatrix"))
(t1 <- spMatrix(5,4, x= c(3,2,-7,11), i= 1:4, j=4:1)) # 5 x 4 (t2 <- kronecker(Diagonal(3, 2:4), t1)) # 15 x 12 ## should also work with special-cased logical matrices l3 <- upper.tri(matrix(,3,3)) M <- Matrix(l3) (N <- as(M, "nsparseMatrix")) # "ntCMatrix" (upper triangular) N2 <- as(N, "generalMatrix") # (lost "t"riangularity) MM <- kronecker(M,M) NN <- kronecker(N,N) # "dtTMatrix" i.e. did keep NN2 <- kronecker(N2,N2) stopifnot(identical(NN,MM), is(NN2, "sparseMatrix"), all(NN2 == NN), is(NN, "triangularMatrix"))
ldenseMatrix
is the virtual class of all dense logical
(S4) matrices. It extends both denseMatrix
and lMatrix
directly.
x
:logical vector containing the entries of the matrix.
Dim
, Dimnames
:see Matrix
.
Class "lMatrix"
, directly.
Class "denseMatrix"
, directly.
Class "Matrix"
, by class "lMatrix"
.
Class "Matrix"
, by class "denseMatrix"
.
signature(x = "ldenseMatrix", mode = "missing")
: ...
signature(x = "ndenseMatrix")
, semantically
equivalent to base function which(x, arr.ind)
;
for details, see the lMatrix
class documentation.
Class lgeMatrix
and the other subclasses.
showClass("ldenseMatrix") as(diag(3) > 0, "ldenseMatrix")
showClass("ldenseMatrix") as(diag(3) > 0, "ldenseMatrix")
The class "ldiMatrix"
of logical diagonal matrices.
Objects can be created by calls of the form new("ldiMatrix", ...)
but typically rather via Diagonal
.
x
:"logical"
vector.
diag
:"character"
string, either "U" or "N",
see ddiMatrix
.
Dim
,Dimnames
:matrix dimension and
dimnames
, see the Matrix
class
description.
Class "diagonalMatrix"
and
class "lMatrix"
, directly.
Class "sparseMatrix"
, by class "diagonalMatrix"
.
Classes ddiMatrix