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] , Martin Maechler [aut, cre] , Mikael Jagan [aut] , Timothy A. Davis [ctb] (<https://orcid.org/0000000176146899>, SuiteSparse libraries, notably CHOLMOD and AMD, collaborators listed in dir(pattern="^[AZ]+[.]txt$", full.names=TRUE, system.file("doc", "SuiteSparse", package="Matrix"))), Jens Oehlschlägel [ctb] (initial nearPD()), Jason Riedy [ctb] (<https://orcid.org/0000000243454200>, GNU Octave's condest() and onenormest(), Copyright: Regents of the University of California), R Core Team [ctb] (base R's matrix implementation) 
Maintainer:  Martin Maechler <[email protected]> 
License:  GPL (>= 2)  file LICENCE 
Version:  1.61.1 
Built:  20231101 20:53:46 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")))
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(...)
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
nonnegative 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)
)
The basic R functions all
and any
now
have methods for Matrix
objects and should
behave as for matrix
ones.
signature(x = "Matrix", ..., na.rm = FALSE)
: ...
signature(x = "Matrix", ..., na.rm = FALSE)
: ...
signature(x = "ldenseMatrix", ..., na.rm = FALSE)
: ...
signature(x = "lsparseMatrix", ..., na.rm = FALSE)
: ...
M < Matrix(1:12 +0, 3,4)
all(M >= 1) # TRUE
any(M < 0 ) # FALSE
MN < M; MN[2,3] < NA; MN
all(MN >= 0) # NA
any(MN < 0) # NA
any(MN < 0, na.rm = TRUE) # > FALSE
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)))
The class
"atomicVector"
is a
virtual class containing all atomic vector classes of base R,
as also implicitly defined via is.atomic
.
A virtual Class: No objects may be created from it.
In the Matrix package, the "atomicVector" is used in signatures where typically “oldstyle” "matrix" objects can be used and can be substituted by simple vectors.
The atomic classes
"logical"
, "integer"
, "double"
, "numeric"
,
"complex"
, "raw"
and "character"
are extended
directly. Note that "numeric"
already contains "integer"
and "double"
, but we want all of them to be direct subclasses of
"atomicVector"
.
Martin Maechler
is.atomic
, integer
, numeric
,
complex
, etc.
showClass("atomicVector")
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, ...)
x 
a matrixlike object 
k , k1 , k2

integers specifying the diagonals that are not set to
zero. These are interpreted relative to the main diagonal, which
is 
... 
optional arguments passed 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, columnoriented matrices.
method for compressed, sparse, roworiented 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 nonzero 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"
Construct a sparse banded matrix by specifying its nonzero sup and superdiagonals.
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 $n \times m$
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!
)
Build a block diagonal matrix given several building block matrices.
bdiag(...)
.bdiag(lst)
... 
individual matrices or a 
lst 
nonempty 
For nontrivial 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 $k \times k$
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
Rhelp; earlier versions have been posted by other authors, notably
Scott Chasalow to Snews. 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 nonzero 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)blocktriangular 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:(M1L), 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 = "mMatrix", y = "mMatrix")
signature(x = "nMatrix", y = "nMatrix")
signature(x = "nMatrix", y = "nsparseMatrix")
signature(x = "nsparseMatrix", y = "nMatrix")
signature(x = "nsparseMatrix", y = "nsparseMatrix")
signature(x = "sparseVector", y = "mMatrix")
signature(x = "mMatrix", y = "sparseVector")
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 nonstructural 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 $m_ij$
.
%*%
, 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"
Classes BunchKaufman
and pBunchKaufman
represent
BunchKaufman factorizations of $n \times n$
real,
symmetric matrices $A$
, having the general form
$A = U D_{U} U' = L D_{L} L'$
where
$D_{U}$
and $D_{L}$
are symmetric, block diagonal
matrices composed of $b_{U}$
and $b_{L}$
$1 \times 1$
or $2 \times 2$
diagonal blocks;
$U = \prod_{k = 1}^{b_{U}} P_{k} U_{k}$
is the product of $b_{U}$
rowpermuted unit upper triangular
matrices, each having nonzero entries above the diagonal in 1 or 2 columns;
and
$L = \prod_{k = 1}^{b_{L}} P_{k} L_{k}$
is the product of $b_{L}$
rowpermuted unit lower triangular
matrices, each having nonzero entries below the diagonal in 1 or 2 columns.
These classes store the nonzero entries of the
$2 b_{U} + 1$
or $2 b_{L} + 1$
factors,
which are individually sparse,
in a dense format as a vector of length
$nn$
(BunchKaufman
) or
$n(n+1)/2$
(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 $A$
or its logarithm.
expand1
signature(x = "p?BunchKaufman")
:
see expand1methods
.
expand2
signature(x = "p?BunchKaufman")
:
see expand2methods
.
solve
signature(a = "p?BunchKaufman", b = .)
:
see solvemethods
.
In Matrix < 1.60
, class BunchKaufman
extended
dtrMatrix
and class pBunchKaufman
extended
dtpMatrix
, reflecting the fact that the internal
representation of the factorization is fundamentally triangular:
there are $n(n+1)/2$
“parameters”, and these
can be arranged systematically to form an $n \times n$
triangular matrix.
Matrix 1.60
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)))
Computes the BunchKaufman factorization of an $n \times n$
real, symmetric matrix $A$
, which has the general form
$A = U D_{U} U' = L D_{L} L'$
where
$D_{U}$
and $D_{L}$
are symmetric, block diagonal
matrices composed of $b_{U}$
and $b_{L}$
$1 \times 1$
or $2 \times 2$
diagonal blocks;
$U = \prod_{k = 1}^{b_{U}} P_{k} U_{k}$
is the product of $b_{U}$
rowpermuted unit upper triangular
matrices, each having nonzero entries above the diagonal in 1 or 2 columns;
and
$L = \prod_{k = 1}^{b_{L}} P_{k} L_{k}$
is the product of $b_{L}$
rowpermuted 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", ...)
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")))
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 Email to Rhelp
(https://stat.ethz.ch/mailman/listinfo/rhelp), by
Casper J Albers, Open University, UK.
data(CAex)
This is a $72 \times 72$
symmetric matrix with 216
nonzero entries in five bands, stored as sparse matrix of class
dgCMatrix
.
Historical note (20060330):
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"))
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 'denseMatrix,sparseMatrix'
cbind2(x, y, sparse = NA, ...)
## S4 method for signature 'sparseMatrix,denseMatrix'
cbind2(x, y, sparse = NA, ...)
## S4 method for signature 'denseMatrix,sparseMatrix'
rbind2(x, y, sparse = NA, ...)
## S4 method for signature 'sparseMatrix,denseMatrix'
rbind2(x, y, sparse = NA, ...)
... 
for 
deparse.level 
integer controlling the construction of labels
in the case of nonmatrixlike arguments; see 
x , y

vector or matrixlike R objects to be bound together. 
sparse 

typically a ‘matrixlike’ 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 nonzero 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))
CHMfactor
is the virtual class of sparse Cholesky
factorizations of $n \times n$
real, symmetric
matrices $A$
, having the general form
$P_1 A P_1' = L_1 D L_1' \overset{D_{jj} \ge 0}{=} L L'$
or (equivalently)
$A = P_1' L_1 D L_1' P_1 \overset{D_{jj} \ge 0}{=} P_1' L L' P_1$
where
$P_1$
is a permutation matrix,
$L_1$
is a unit lower triangular matrix,
$D$
is a diagonal matrix, and
$L = L_1 \sqrt{D}$
.
The second equalities hold only for positive semidefinite $A$
,
for which the diagonal entries of $D$
are nonnegative
and $\sqrt{D}$
is welldefined.
The implementation of class CHMfactor
is based on
CHOLMOD's Clevel 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)
x 
an object inheriting from virtual class 
isLDL(x)
returns TRUE
or FALSE
:
TRUE
if x
stores the lower triangular entries
of $L_1I+D$
,
FALSE
if x
stores the lower triangular entries
of $L$
.
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 0based 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 $L$
or (if and only if type[2]
is 0) of the lower triangular matrix $L_1I+D$
.
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 $L$
or
the lower triangular matrix $L_1I+D$
,
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 $L$
. Note that,
for supernodes spanning two or more columns, the supernodal
algorithm by design stores nonstructural 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 $A$
or its logarithm.
If sqrt = TRUE
, then this method computes the determinant
of the factor $L = L_1 sqrt(D)$
or
its logarithm, giving NaN
for the modulus when $D$
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.60
). 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 $n$
containing the diagonal
elements of $D$
, which (if they are all nonnegative)
are the squared diagonal elements of $L$
.
expand
signature(x = "CHMfactor")
:
see expandmethods
.
expand1
signature(x = "CHMsimpl")
:
see expand1methods
.
expand1
signature(x = "CHMsuper")
:
see expand1methods
.
expand2
signature(x = "CHMsimpl")
:
see expand2methods
.
expand2
signature(x = "CHMsuper")
:
see expand2methods
.
image
signature(x = "CHMfactor")
:
see imagemethods
.
nnzero
signature(x = "CHMfactor")
:
see nnzeromethods
.
solve
signature(a = "CHMfactor", b = .)
:
see solvemethods
.
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 updownmethods
.
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, 114. 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), 886905. 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 = 1e14))
Computes the upper triangular Cholesky factor of an
$n \times n$
real, symmetric, positive semidefinite
matrix $A$
, optionally after pivoting.
That is the factor $L'$
in
$P_{1} A P_{1}' = L L'$
or (equivalently)
$A = P_{1}' L L' P_{1}$
where
$P_{1}$
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 $P_{1}$
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", ...)
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 $P_{1}$
in addition
to the Cholesky factor $L'$
, then call Cholesky
directly, as the result of chol(x, pivot = TRUE)
specifies
$L'$
but not $P_{1}$
.
A matrix, triangularMatrix
,
or diagonalMatrix
representing
the upper triangular Cholesky factor $L'$
.
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, 114. 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), 886905. 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 $L'$
.
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 2by2 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
$U$
and $L$
, compute $(U' U)^{1}$
and $(L L')^{1}$
, 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
$(X' X)^{1}$
given the $R$
part
of the QR factorization of $X$
.
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 $L'$
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
})
Classes Cholesky
and pCholesky
represent
dense, pivoted Cholesky factorizations of $n \times n$
real, symmetric, positive semidefinite matrices $A$
,
having the general form
$P_{1} A P_{1}' = L_{1} D L_{1}' = L L'$
or (equivalently)
$A = P_{1}' L_{1} D L_{1}' P_{1} = P_{1}' L L' P_{1}$
where
$P_{1}$
is a permutation matrix,
$L_{1}$
is a unit lower triangular matrix,
$D$
is a nonnegative diagonal matrix, and
$L = L_{1} \sqrt{D}$
.
These classes store the entries of the Cholesky factor
$L$
or its transpose $L'$
in a dense format as
a vector of length $nn$
(Cholesky
) or
$n(n+1)/2$
(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 $L'$
or $L$
.
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 $L$
or its transpose $L'$
in columnmajor
order.
perm
a 1based 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 $L$
or its transpose $L'$
;
see ‘Note’.
coerce
signature(from = "pCholesky", to = "dtpMatrix")
:
returns a dtpMatrix
representing
the Cholesky factor $L$
or its transpose $L'$
;
see ‘Note’.
determinant
signature(from = "p?Cholesky", logarithm = "logical")
:
computes the determinant of the factorized matrix $A$
or its logarithm.
diag
signature(x = "p?Cholesky")
:
returns a numeric vector of length $n$
containing the diagonal
elements of $D$
, which are the squared diagonal elements of
$L$
.
expand1
signature(x = "p?Cholesky")
:
see expand1methods
.
expand2
signature(x = "p?Cholesky")
:
see expand2methods
.
solve
signature(a = "p?Cholesky", b = .)
:
see solvemethods
.
In Matrix < 1.60
, class Cholesky
extended
dtrMatrix
and class pCholesky
extended
dtpMatrix
, reflecting the fact that the factor
$L$
is indeed a triangular matrix.
Matrix 1.60
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). LAPACKstyle 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)))
Computes the pivoted Cholesky factorization of an
$n \times n$
real, symmetric matrix $A$
,
which has the general form
$P_1 A P_1' = L_1 D L_1' \overset{D_{jj} \ge 0}{=} L L'$
or (equivalently)
$A = P_1' L_1 D L_1' P_1 \overset{D_{jj} \ge 0}{=} P_1' L L' P_1$
where
$P_1$
is a permutation matrix,
$L_1$
is a unit lower triangular matrix,
$D$
is a diagonal matrix, and
$L = L_1 \sqrt{D}$
.
The second equalities hold only for positive semidefinite $A$
,
for which the diagonal entries of $D$
are nonnegative
and $\sqrt{D}$
is welldefined.
Methods for denseMatrix
are built on
LAPACK routines dpstrf
, dpotrf
, and dpptrf
.
The latter two do not permute rows or columns,
so that $P_1$
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", ...)
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 $L'$
,
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 $A$
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 $P' L L' P$
is positive
semidefinite even if $A$
is not. It follows that one way to
test for positive semidefiniteness of $A$
in the event of a
warning is to analyze the error
$\frac{\lVert A  P' L L' P \rVert}{\lVert A \rVert}\,.$
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). LAPACKstyle 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, 114. 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), 886905. 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 $L'$
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.670858e17
## .... 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.777944e16
##  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 latticebased 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
.SuiteSparse_version()
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 = is_not_uniqT(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"))
}
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, ...)
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)))
})
Virtual class of composite matrices; i.e., matrices that can be factorized, typically as a product of simpler matrices.
A virtual Class: No objects may be created from it.
factors
:Object of class "list"
 a list of
factorizations of the matrix. Note that this is typically empty,
i.e., list()
, initially and is updated
automagically whenever a matrix factorization is computed.
Dim
, Dimnames
:inherited from the
Matrix
class, see there.
Class "Matrix"
, directly.
signature(x = "compMatrix", value = "list")
: set
the dimnames
to a list
of length 2, see
dimnames<
. The factors
slot is currently reset to
empty, as the factorization dimnames
would have to be adapted, too.
The matrix factorization classes
"MatrixFactorization"
and their generators,
lu()
,
qr()
,
chol()
and Cholesky()
,
BunchKaufman()
, Schur()
.
“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 1norm,
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)
A 
a square matrix, optional for 
t 
number of columns to use in the iterations. 
normA 
number; (an estimate of) the 1norm 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 1norm 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
, $A^{1}$
, 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 
01 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 1Norm Estimation, with an Application to 1Norm 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 approxmation 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 columnoriented 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 nonzero elements). These are the 0based row numbers for
each nonzero 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 (zerobased) 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 nonzero 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.99937516
),
validObject
automatically resorted 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"))
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")
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 $n \times n$
matrix, the x
slot is of length $n$
or 0
,
depending on the diag
slot:
diag
:"character"
string, either "U"
or
"N"
where "U"
denotes unitdiagonal, 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)))
denseLU
is the class of dense, rowpivoted LU factorizations
of $m \times n$
real matrices $A$
,
having the general form
$P_{1} A = L U$
or (equivalently)
$A = P_{1}' L U$
where
$P_{1}$
is an $m \times m$
permutation matrix,
$L$
is an $m \times \text{min}(m,n)$
unit lower trapezoidal matrix, and
$U$
is a $\text{min}(m,n) \times n$
upper trapezoidal matrix. If $m = n$
, then the factors
$L$
and $U$
are triangular.
Dim
, Dimnames
inherited from virtual class
MatrixFactorization
.
x
a numeric vector of length prod(Dim)
storing
the triangular $L$
and $U$
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 $P_{1}$
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 $A$
, equal to $L$
below the
diagonal and equal to $U$
on and above the diagonal.
determinant
signature(from = "denseLU", logarithm = "logical")
:
computes the determinant of the factorized matrix $A$
or its logarithm.
expand
signature(x = "denseLU")
:
see expandmethods
.
expand1
signature(x = "denseLU")
:
see expand1methods
.
expand2
signature(x = "denseLU")
:
see expand2methods
.
solve
signature(a = "denseLU", b = "missing")
:
see solvemethods
.
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)))
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 ($2+3 = 5$
) 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")
The dgCMatrix
class is a class of sparse numeric
matrices in the compressed, sparse, columnoriented format. In this
implementation the nonzero 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 nonzero
elements of the matrix.
all other slots are inherited from the superclass
"CsparseMatrix"
.
Matrix products (e.g., crossprodmethods), 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 solvemethods
, 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]
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 columnmajor order.
Dim
:Object of class "integer"
 the dimensions
of the matrix  must be an integer vector with exactly two
nonnegative 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, roworiented format. In this
implementation the nonzero elements in the rows are sorted into
increasing column order.
Note: The columnoriented 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 nonzero elements). These are the column numbers for
each nonzero element in the matrix.
p
:Object of class "integer"
of pointers, one
for each row, to the initial (zerobased) index of elements in
the row.
x
:Object of class "numeric"
 the nonzero
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 roworiented 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 nonzero
entries in 0base, i.e., must be in 0:(nrow(.)1)
.
j
:integer
column indices of nonzero
entries. Must be the same length as slot i
and
0based as well, i.e., in 0:(ncol(.)1)
.
x
:numeric
vector  the (nonzero)
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 $x_k$
's that belong to identical
$(i_k, j_k)$
pairs.
However this means that a matrix typically can be stored in more than
one possible "TsparseMatrix"
representations.
Use uniqTsparse()
in order to ensure uniqueness of the
internal representation of such a matrix.
Class dgCMatrix
or the superclasses
dsparseMatrix
and
TsparseMatrix
; uniqTsparse
.
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 nonzero 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)
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
Clevel 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 nonzero 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 "repeatedblock" 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 (201210) 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 ‘unitdiagonal’.
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 solvemethods
.
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
offdiagonal zeros or NA
s, we implicitly use $0 / x = 0$
, hence
differing from traditional R arithmetic (where $0 / 0
\mapsto \mbox{NaN}$
), 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 = 1e12)
)
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)
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
offdiagonal 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))
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)
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)
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
nonnegative 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
.
signature(x = "dMatrix")
: ...
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 2column matrix of row and column indices. Since Matrix
version 1.29, if useNames
is true, as by default, with
dimnames
, the same as base::which
.
The nonzeropattern matrix class nMatrix
, which
can be used to store nonNA
logical
matrices even more compactly.
The numeric matrix classes dgeMatrix
,
dgCMatrix
, and Matrix
.
drop0(x, tol=1e10)
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))
For any $n \times m$
(typically) sparse matrix x
compute the DulmageMendelsohn row and columns permutations which at
first splits the $n$
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)
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 DulmageMendelsohn 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)) # nonrandom 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
positivesemidefinite 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 "pcorMatrix"
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 columnmajor order.
Dim
:Object of class "integer"
. The dimensions
of the matrix which must be a twoelement vector of nonnegative
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 "pcorMatrix"
)
a numeric
vector of length n
containing the
(original) $\sqrt{var(.)}$
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
onenorm (the default) or "I"
for the infinitynorm. 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 solvemethods
.
signature(e1 = "dpoMatrix", e2 = "numeric")
(and
quite a few other signatures): The result of (“elementwise”
defined) arithmetic operations is typically not
positivedefinite 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 semidefinite.
A more reliable (but often more expensive) check for positive
semidefiniteness 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 < as(h6, "dppMatrix"))
### Note that as(*, "corMatrix") *scales* the matrix
(ch6 < as(h6, "corMatrix"))
stopifnot(all.equal(h6 * 27720, round(27720 * h6), tolerance = 1e14),
all.equal(ch6@sd^(2), 2*(1:6)1, tolerance= 1e12))
chch < Cholesky(ch6, perm = FALSE)
stopifnot(identical(chch, ch6@factors$Cholesky),
all(abs(crossprod(as(chch, "dtrMatrix"))  ch6) < 1e10))
Deletes “nonstructural” 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)
x 
a 
tol 
a nonnegative 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
nonstructural 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 = 1e15)
stopifnot(all(I.0 == Diagonal(625)), nnzero(I..) == 625)
The dsCMatrix
class is a class of symmetric, sparse
numeric matrices in the compressed, columnoriented format. In
this implementation the nonzero 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 nonzero elements). These are the row
numbers for each nonzero 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 nonzero element in the lower triangle of
the matrix.
x
:Object of class "numeric"
of length nnZ –
the nonzero 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
nonnegative 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 $A x = b$
for $x$
; see
solvemethods
.
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"
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
(nonzero) 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")
The dsRMatrix
class is a class of symmetric, sparse
matrices in the compressed, roworiented format. In this
implementation the nonzero 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 nonzero elements). These are the row
numbers for each nonzero element in the matrix.
p
:Object of class "integer"
of pointers, one
for each row, to the initial (zerobased) index of elements in
the row.
factors
:Object of class "list"
 a list
of factorizations of the matrix.
x
:Object of class "numeric"
 the nonzero
elements of the matrix.
Dim
:Object of class "integer"
 the dimensions
of the matrix  must be an integer vector with exactly two
nonnegative values.
Dimnames
:List of length two, see Matrix
.
Classes RsparseMatrix
, dsparseMatrix
and
symmetricMatrix
, directly.
Class "dMatrix"
, by class "dsparseMatrix"
,
class "sparseMatrix"
, by class "dsparseMatrix"
or
"RsparseMatrix"
;
class "compMatrix"
by class "symmetricMatrix"
and of course,
class "Matrix"
.
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
The "dsyMatrix"
class is the class of symmetric, dense matrices
in nonpacked 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 columnmajor order.
Dim
,Dimnames
:The dimension (a length2
"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 $A x = b$
for $x$
; see
solvemethods
.
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 nonpacked
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
The "dtCMatrix"
class is a class of triangular, sparse
matrices in the compressed, columnoriented format. In this
implementation the nonzero 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 nonzero elements). These are the row numbers for
each nonzero element in the matrix.
j
:Object of class "integer"
of length nnzero
(number of nonzero elements). These are the column numbers for
each nonzero element in the matrix. (Only present in the
dtTMatrix
class.)
x
:Object of class "numeric"
 the nonzero
elements of the matrix.
Dim
,Dimnames
:The dimension (a length2
"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 solvemethods
.
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 0diagonal to unitdiagonal {lowlevel 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=1e14) ,
all.equal(Diagonal(5), asDiag(i5), tolerance=1e14) ,
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 columnmajor order.
For a packed square matrix of dimension $d \times d$
,
length(x)
is of length $d(d+1)/2$
(also when
diag == "U"
!).
Dim
,Dimnames
:The dimension (a length2
"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
solvemethods
.
signature(x = "dtpMatrix")
: t(x)
remains
a "dtpMatrix"
, lower triangular if x
is upper
triangular, and vice versa.
Class dtrMatrix
showClass("dtrMatrix")
example("dtrMatrixclass", 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 columnmajor order.
Dim
:Object of class "integer"
. The dimensions
of the matrix which must be a twoelement vector of nonnegative
integers.
Class "ddenseMatrix"
, directly.
Class "triangularMatrix"
, directly.
Class "Matrix"
and others, by class "ddenseMatrix"
.
Among others (such as matrix products, e.g. ?crossprodmethods
),
signature(x = "dtrMatrix", type = "character")
signature(x = "dtrMatrix", norm = "character")
signature(a = "dtrMatrix", b = "....")
efficientely
use a “forwardsolve” or backsolve
for a lower or
upper triangular matrix, respectively, see also
solvemethods
.
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 all0
})
The dtRMatrix
class is a class of triangular, sparse
matrices in the compressed, roworiented format. In this
implementation the nonzero 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 nonzero elements). These are
the row numbers for each nonzero element in the matrix.
p
:Object of class "integer"
of pointers, one
for each row, to the initial (zerobased) index of elements in
the row. (Only present in the dsRMatrix
class.)
x
:Object of class "numeric"
 the nonzero
elements of the matrix.
Dim
:The dimension (a length2 "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
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, ...)
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.60
. 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
$A = P_{1}' L_{1} D L_{1}' P_{1} = P_{1}' L L' P_{1}$
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
$A = L_{1} D L_{1}' = L L'$
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
$A = U D_{U} U' = L D_{L} L'$
where
$U = \prod_{k = 1}^{b_{U}} P_{k} U_{k}$
and
$L = \prod_{k = 1}^{b_{L}} P_{k} L_{k}$
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 $2 b_{U} + 1$
or $2 b_{L} + 1$
matrix
factors is returned instead.
$P_{k}$
are represented as pMatrix
,
$U_{k}$
and $L_{k}$
are represented as
dtCMatrix
, and
$D_{U}$
and $D_{L}$
are represented as
dsCMatrix
.
expand2
signature(x = "Schur")
:
expands the factorization
$A = Q T Q'$
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
$A = P_{1}' L U P_{2}'$
as list(P1., L, U, P2.)
.
P1.
and P2.
are pMatrix
,
and L
and U
are dtCMatrix
.
expand2
signature(x = "denseLU")
:
expands the factorization
$A = P_{1}' L U$
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
$A = P_{1}' Q R P_{2}' = P_{1}' Q_{1} R_{1} P_{2}'$
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 $P_{1}$
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 $P_{2}'$
but have opposite margin
slots and inverted
perm
slots. expand(x)[["P"]]
represents
the permutation matrix $P_{1}$
rather than its
transpose $P_{1}'$
; it is expand2(x)[["P1."]]
with an inverted perm
slot. expand(x)[["L"]]
and expand2(x)[["L"]]
represent the same unit lower
triangular matrix $L$
, 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 compMatrix
of factorizable matrices.
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("Choleskymethods")
(a1 < sapply(names(e.lu.A), expand1, x = lu.A, simplify = FALSE))
all.equal(a1, e.lu.A)
## see help("denseLUclass") and others for more examples
Compute the exponential of a matrix.
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, twentyfive 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 (nongeneric) 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 = 1e15))
(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 HarwellBoeing or MatrixMarket formats
or write sparseMatrix
objects to one of these
formats.
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 HarwellBoeing 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
sparseMatrixclass
, 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/HarwellBoeing/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=1e15))
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, ...)
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)
“SemiAPI” 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)
.sparse2dense(from, packed = FALSE)
.diag2dense(from, shape = "t", packed = FALSE, uplo = "U")
.ind2dense(from, kind = "n")
.m2dense(from, class, uplo = "U", diag = "N")
.dense2sparse(from, repr = "C")
.diag2sparse(from, shape = "t", repr = "C", uplo = "U")
.ind2sparse(from, kind = "n", repr = ".")
.m2sparse(from, class, uplo = "U", diag = "N")
.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 ( 
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” nonvirtual subclass of
virtual class B, where the virtual classes are abbreviated as follows:
M
m
matrix or vector
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. For example, .m2dense
and .m2sparse
accept vectors in addition to 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 $L D L'$
in order
to calculate one of several possible statistics from the diagonal
entries of $D$
. 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 = 1e3)# rel.err ~ 1e4
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)
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)
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))
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 unitdiagonal 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
nonzeroes 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))
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.
factors
,
Dim
,
Dimnames
:all slots inherited from
compMatrix
; see its description.
Class "compMatrix"
, directly.
Class "Matrix"
, by class "compMatrix"
.
Classes compMatrix
, and the nongeneral virtual
classes:
symmetricMatrix
,
triangularMatrix
,
diagonalMatrix
.
Generate the n
by n
symmetric Hilbert matrix. Because
these matrices are illconditioned for moderate to large n
,
they are often used for testing numerical linear algebra code.
Hilbert(n)
n 
a nonnegative integer. 
the n
by n
symmetric Hilbert matrix as a
"dpoMatrix"
object.
the class dpoMatrix
Hilbert(6)
Methods for function image
in package
Matrix. An image of a matrix simply color codes all matrix
entries and draws the $n\times m$
matrix using an
$n\times m$
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, ...)
x 
a Matrix object, i.e., fulfilling 
xlim , ylim

x and yaxis 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 cairobased
## 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("imageUSCountiesraster.png", width=3200, height=3200)
image(USCounties, useRaster = TRUE) # should not suffer from antialiasing
if(doPNG)
dev.off()
## and now look at the *.png image in a viewer you can easily zoom in and out
}#only if(doExtras)
The class "index"
is a virtual class used for
indices (in signatures) for matrix indexing and subassignment of
Matrix matrices.
In fact, it is currently implemented as a simple class union
(setClassUnion
) of
"numeric"
, "logical"
and "character"
.
Since it is a virtual Class, no objects may be created from it.
[methods
, and
Subassignmethods
, also for examples.
showClass("index")
The indMatrix
class is the class of row and column
index matrices, stored as 1based 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$ 
=  R %*% M 
=  M[p, ]

$M C$ 
=  M %*% C 
=  M[, q]

$C'M$ 
=  crossprod(C, M) 
=  M[q, ]

$MR'$ 
=  tcrossprod(M, R) 
=  M[, p]

$R'R$ 
=  crossprod(R) 
=  Diagonal(x=tabulate(p, ncol(R)))

$CC'$ 
=  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) columnindexing
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 vectorvalued 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 1based 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 1based 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")
:
rowwise catenation of two row index matrices with equal numbers
of columns and columnwise 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 and Uni Muenchen, 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),]))
## rowindexing 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 crosstabulation 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 lengthn
permutation vector. isPerm
tests
if a lengthn
integer vector is a valid permutation vector.
asPerm
coerces a lengthm
transposition vector to a
lengthn
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)
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'
Methods for generic functions is.na()
,
is.nan()
, is.finite()
,
is.infinite()
, and anyNA()
,
for objects inheriting from virtual class
Matrix
or sparseVector
.
## S4 method for signature 'dsparseMatrix'
is.na(x)
## S4 method for signature 'dsparseMatrix'
is.nan(x)
## S4 method for signature 'dsparseMatrix'
is.finite(x)
## S4 method for signature 'dsparseMatrix'
is.infinite(x)
## S4 method for signature 'dsparseMatrix'
anyNA(x)
## ...
## and for other classes
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)
dn
NULLlike ?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)
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"))
})
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 'symmetricMatrix'
isSymmetric(object, ...)
## S4 method for signature 'triangularMatrix'
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, tol = 100 * .Machine$double.eps, tol1 = 8 * tol, checkDN = TRUE, ...)
## S4 method for signature 'lgeMatrix'
isSymmetric(object, checkDN = TRUE, ...)
## S4 method for signature 'ngeMatrix'
isSymmetric(object, checkDN = TRUE, ...)
## S4 method for signature 'dgCMatrix'
isSymmetric(object, tol = 100 * .Machine$double.eps, checkDN = TRUE, ...)
## S4 method for signature 'lgCMatrix'
isSymmetric(object, checkDN = TRUE, ...)
## S4 method for signature 'ngCMatrix'
isSymmetric(object, checkDN = TRUE, ...)
object 
a 
tol , tol1

numerical tolerances allowing approximate
symmetry of numeric (rather than logical) matrices. See also

checkDN 
a logical indicating whether symmetry of the

... 
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
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)
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)
Computes KhatriRao products for any kind of matrices.
The KhatriRao product is a columnwise 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 higherdimensional tensor decompositions, see Bader and Kolda (2008).
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 KhatriRao
product of X
($n \times k$
) and Y
($m
\times k$
), is of dimension $(n\cdot m) \times k$
,
where the jth 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 nonstructural 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)# unitdiagonal: 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 nonzero entries.
It is a "dgCMatrix"
object. The vector y
is just
numeric
of length 1850.
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 QRbased solution ("more accurate, but slightly slower"):
system.time(
sp.sol2 < with(KNex, qr.coef(qr(mm), y) ))
all.equal(sparse.sol, sp.sol2, tolerance = 1e13) # 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 specialcased 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")
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
and
diagonalMatrix
; function Diagonal
.
(lM < Diagonal(x = c(TRUE,FALSE,FALSE)))
str(lM)#> gory details (slots)
crossprod(lM) # numeric
(nM < as(lM, "nMatrix"))# > sparse (not formally ``diagonal'')
crossprod(nM) # logical sparse
This is the class of general dense logical
matrices.
x
:Object of class "logical"
. The logical
values that constitute the matrix, stored in columnmajor order.
Dim
,Dimnames
:The dimension (a length2
"integer"
) and corresponding names (or NULL
), see the
Matrix
class.
factors
:Object of class "list"
. A named
list of factorizations that have been computed for the matrix.
Class "ldenseMatrix"
, directly.
Class "lMatrix"
, by class "ldenseMatrix"
.
Class "denseMatrix"
, by class "ldenseMatrix"
.
Class "Matrix"
, by class "ldenseMatrix"
.
Class "Matrix"
, by class "ldenseMatrix"
.
Currently, mainly t()
and coercion methods (for
as(.)
); use, e.g.,
showMethods(class="lgeMatrix")
for details.
Nongeneral logical dense matrix classes such as
ltrMatrix
, or lsyMatrix
;
sparse logical classes such as lgCMatrix
.
showClass("lgeMatrix")
str(new("lgeMatrix"))
set.seed(1)
(lM < Matrix(matrix(rnorm(28), 4,7) > 0))# a simple random lgeMatrix
set.seed(11)
(lC < Matrix(matrix(rnorm(28), 4,7) > 0))# a simple random lgCMatrix
as(lM, "CsparseMatrix")
The lsparseMatrix
class is a virtual class of sparse
matrices with TRUE
/FALSE
or NA
entries. Only the
positions of the elements that are TRUE
are stored.
These can be stored in the “triplet” form (class
TsparseMatrix
, subclasses lgTMatrix
,
lsTMatrix
, and ltTMatrix
) or in compressed
columnoriented form (class CsparseMatrix
,
subclasses lgCMatrix
, lsCMatrix
, and ltCMatrix
)
or–rarely–in compressed roworiented form (class
RsparseMatrix
, subclasses lgRMatrix
,
lsRMatrix
, and ltRMatrix
). The second letter in the
name of these nonvirtual classes indicates g
eneral,
s
ymmetric, or t
riangular.
Note that triplet stored (TsparseMatrix
) matrices
such as lgTMatrix
may contain duplicated pairs of indices
$(i,j)$
as for the corresponding numeric class
dgTMatrix
where for such pairs, the corresponding
x
slot entries are added. For logical matrices, the x
entries corresponding to duplicated index pairs $(i,j)$
are
“added” as well if the addition is defined as logical $or$
,
i.e., “TRUE + TRUE > TRUE
” and
“TRUE + FALSE > TRUE
”.
Note the use of uniqTsparse()
for getting an internally
unique representation without duplicated $(i,j)$
entries.
Objects can be created by calls of the form new("lgCMatrix",
...)
and so on. More frequently objects are created by coercion of
a numeric sparse matrix to the logical form, e.g. in an expression
x != 0
.
The logical form is also used in the symbolic analysis phase of an algorithm involving sparse matrices. Such algorithms often involve two phases: a symbolic phase wherein the positions of the nonzeros in the result are determined and a numeric phase wherein the actual results are calculated. During the symbolic phase only the positions of the nonzero elements in any operands are of interest, hence any numeric sparse matrices can be treated as logical sparse matrices.
x
:Object of class "logical"
, i.e., either
TRUE
, NA
, or FALSE
.
uplo
:Object of class "character"
. Must be
either "U", for upper triangular, and "L", for lower
triangular. Present in the triangular and symmetric classes but not
in the general class.
diag
:Object of class "character"
. Must be
either "U"
, for unit triangular (diagonal is all ones), or
"N"
for nonunit. The implicit diagonal elements are not
explicitly stored when diag
is "U"
. Present in the
triangular classes only.
p
:Object of class "integer"
of pointers, one
for each column (row), to the initial (zerobased) index of elements in
the column. Present in compressed columnoriented and compressed
roworiented forms only.
i
:Object of class "integer"
of length nnzero
(number of nonzero elements). These are the row numbers for
each TRUE element in the matrix. All other elements are FALSE.
Present in triplet and compressed columnoriented forms only.
j
:Object of class "integer"
of length nnzero
(number of nonzero elements). These are the column numbers for
each TRUE element in the matrix. All other elements are FALSE.
Present in triplet and compressed roworiented forms only.
Dim
:Object of class "integer"
 the dimensions
of the matrix.
signature(from = "dgCMatrix", to = "lgCMatrix")
signature(x = "lgCMatrix")
: returns the transpose
of x
signature(x = "lsparseMatrix")
, semantically
equivalent to base function which(x, arr.ind)
;
for details, see the lMatrix
class documentation.
the class dgCMatrix
and dgTMatrix
(m < Matrix(c(0,0,2:0), 3,5, dimnames=list(LETTERS[1:3],NULL)))
(lm < (m > 1)) # lgC
!lm # no longer sparse
stopifnot(is(lm,"lsparseMatrix"),
identical(!lm, m <= 1))
data(KNex, package = "Matrix")
str(mmG.1 < (KNex $ mm) > 0.1)# "lgC..."
table(mmG.1@x)# however with many ``nonstructural zeros''
## from logical to nz_pattern  okay when there are no NA's :
nmG.1 < as(mmG.1, "nMatrix") # <<< has "TRUE" also where mmG.1 had FALSE
## from logical to "double"
dmG.1 < as(mmG.1, "dMatrix") # has '0' and back:
lmG.1 < as(dmG.1, "lMatrix")
stopifnot(identical(nmG.1, as((KNex $ mm) != 0,"nMatrix")),
validObject(lmG.1),
identical(lmG.1, mmG.1))
class(xnx < crossprod(nmG.1))# "nsC.."
class(xlx < crossprod(mmG.1))# "dsC.." : numeric
is0 < (xlx == 0)
mean(as.vector(is0))# 99.3% zeros: quite sparse, but
table(xlx@x == 0)# more than half of the entries are (nonstructural!) 0
stopifnot(isSymmetric(xlx), isSymmetric(xnx),
## compare xnx and xlx : have the *same* nonstructural 0s :
sapply(slotNames(xnx),
function(n) identical(slot(xnx, n), slot(xlx, n))))
The "lsyMatrix"
class is the class of symmetric, dense logical
matrices in nonpacked storage and "lspMatrix"
is the class of
of these in packed storage. In the packed form, only the upper
triangle or the lower triangle is stored.
Objects can be created by calls of the form new("lsyMatrix", ...)
.
uplo
:Object of class "character"
. Must be
either "U", for upper triangular, and "L", for lower triangular.
x
:Object of class "logical"
. The logical
values that constitute the matrix, stored in columnmajor order.
Dim
,Dimnames
:The dimension (a length2
"integer"
) and corresponding names (or NULL
), see the
Matrix
class.
factors
:Object of class "list"
. A named
list of factorizations that have been computed for the matrix.
Both extend classes "ldenseMatrix"
and
"symmetricMatrix"
, directly; further, class
"Matrix"
and others, indirectly. Use
showClass("lsyMatrix")
, e.g., for details.
Currently, mainly t()
and coercion methods (for
as(.)
; use, e.g.,
showMethods(class="lsyMatrix")
for details.
(M2 < Matrix(c(TRUE, NA, FALSE, FALSE), 2, 2)) # logical dense (ltr)
str(M2)
# can
(sM < M2  t(M2)) # "lge"
as(sM, "symmetricMatrix")
str(sM < as(sM, "packedMatrix")) # packed symmetric
The "ltrMatrix"
class is the class of triangular, dense,
logical matrices in nonpacked storage. The "ltpMatrix"
class
is the same except in packed storage.
x
:Object of class "logical"
. The logical
values that constitute the matrix, stored in columnmajor order.
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
.
Dim
,Dimnames
:The dimension (a length2
"integer"
) and corresponding names (or NULL
), see the
Matrix
class.
factors
:Object of class "list"
. A named
list of factorizations that have been computed for the matrix.
Both extend classes "ldenseMatrix"
and
"triangularMatrix"
, directly; further, class
"Matrix"
, "lMatrix"
and others,
indirectly. Use showClass("ltrMatrix")
, e.g.,
for details.
Currently, mainly t()
and coercion methods (for
as(.)
; use, e.g.,
showMethods(class="ltrMatrix")
for details.
Classes lgeMatrix
, Matrix
;
function t
showClass("ltrMatrix")
str(new("ltpMatrix"))
(lutr < as(upper.tri(matrix(, 4, 4)), "ldenseMatrix"))
str(lutp < pack(lutr)) # packed matrix: only 10 = 4*(4+1)/2 entries
!lutp # the logical negation (is *not* logical triangular !)
## but this one is:
stopifnot(all.equal(lutp, pack(!!lutp)))
Computes the pivoted LU factorization of an $m \times n$
real matrix $A$
, which has the general form
$P_{1} A P_{2} = L U$
or (equivalently)
$A = P_{1}' L U P_{2}'$
where
$P_{1}$
is an $m \times m$
permutation matrix,
$P_{2}$
is an $n \times n$
permutation matrix,
$L$
is an $m \times \text{min}(m,n)$
unit lower trapezoidal matrix, and
$U$
is a $\text{min}(m,n) \times n$
upper trapezoidal matrix.
Methods for denseMatrix
are built on
LAPACK routine dgetrf
, which does not permute columns,
so that $P_{2}$
is an identity matrix.
Methods for sparseMatrix
are built on
CSparse routine cs_lu
, which requires $m = n$
,
so that $L$
and $U$
are triangular matrices.
lu(x, ...)
## S4 method for signature 'dgeMatrix'
lu(x, warnSing = TRUE, ...)
## S4 method for signature 'dgCMatrix'
lu(x, errSing = TRUE, order = NA_integer_,
tol = 1, ...)
## S4 method for signature 'dsyMatrix'
lu(x, cache = TRUE, ...)
## S4 method for signature 'dsCMatrix'
lu(x, cache = TRUE, ...)
## S4 method for signature 'matrix'
lu(x, ...)
x 
a finite matrix or

warnSing 
a logical indicating if a warning should
be signaled for singular 
errSing 
a logical indicating if an error should
be signaled for singular 
order 
an integer in 
tol 
a number. The original pivot element is used
if its absolute value exceeds 
cache 
a logical indicating if the result should be
cached in 
... 
further arguments passed to or from methods. 
What happens when x
is determined to be nearsingular
differs by method. The method for class dgeMatrix
completes the factorization, warning if warnSing = TRUE
and in any case returning a valid denseLU
object. Users of this method can detect singular x
with
a suitable warning handler; see tryCatch
.
In contrast, the method for class dgCMatrix
abandons further computation, throwing an error if errSing = TRUE
and otherwise returning NA
. Users of this method can
detect singular x
with an error handler or by setting
errSing = FALSE
and testing for a formal result with
is(., "sparseLU")
.
An object representing the factorization, inheriting from
virtual class LU
. The specific class
is denseLU
unless x
inherits
from virtual class sparseMatrix
,
in which case it is sparseLU
.
The LAPACK source code, including documentation; see https://netlib.org/lapack/double/dgetrf.f.
Davis, T. A. (2006). Direct methods for sparse linear systems. Society for Industrial and Applied Mathematics. doi:10.1137/1.9780898718881
Golub, G. H., & Van Loan, C. F. (2013). Matrix computations (4th ed.). Johns Hopkins University Press. doi:10.56021/9781421407944
Classes denseLU
and
sparseLU
and their methods.
Classes dgeMatrix
and
dgCMatrix
.
Generic functions expand1
and expand2
,
for constructing matrix factors from the result.
Generic functions Cholesky
, BunchKaufman
,
Schur
, and qr
,
for computing other factorizations.
showMethods("lu", inherited = FALSE)
set.seed(0)
##  Dense 
(A1 < Matrix(rnorm(9L), 3L, 3L))
(lu.A1 < lu(A1))
(A2 < round(10 * A1[, 3L]))
(lu.A2 < lu(A2))
## A ~ P1' L U in floating point
str(e.lu.A2 < expand2(lu.A2), max.level = 2L)
stopifnot(all.equal(A2, Reduce(`%*%`, e.lu.A2)))
##  Sparse 
A3 < as(readMM(system.file("external/pores_1.mtx", package = "Matrix")),
"CsparseMatrix")
(lu.A3 < lu(A3))
## A ~ P1' L U P2' in floating point
str(e.lu.A3 < expand2(lu.A3), max.level = 2L)
stopifnot(all.equal(A3, Reduce(`%*%`, e.lu.A3)))
From an R object coercible to "TsparseMatrix"
,
typically a (sparse) matrix, produce its triplet representation which may
collapse to a “Duplet” in the case of binary aka pattern, such as
"nMatrix"
objects.
mat2triplet(x, uniqT = FALSE)
x 
any R object for which 
uniqT 

A list
, typically with three components,
i 
vector of row indices for all nonzero entries of 
i 
vector of columns indices for all nonzero entries of 
x 
vector of all nonzero entries of 
Note that the order
of the entries is determined by the
coercion to "TsparseMatrix"
and hence typically
with increasing j
(and increasing i
within ties of j
).
The mat2triplet()
utility was created to be a more efficient and
more predictable substitute for summary(<sparseMatrix>)
.
UseRs have wrongly expected the latter to return a data frame with
columns i
and j
which however is wrong for a
"diagonalMatrix"
.
The summary()
method for "sparseMatrix"
,
summary,sparseMatrixmethod
.
mat2triplet()
is conceptually the inverse function of
spMatrix
and (one case of) sparseMatrix
.
mat2triplet # simple definition
i < c(1,3:8); j < c(2,9,6:10); x < 7 * (1:7)
(Ax < sparseMatrix(i, j, x = x)) ## 8 x 10 "dgCMatrix"
str(trA < mat2triplet(Ax))
stopifnot(i == sort(trA$i), sort(j) == trA$j, x == sort(trA$x))
D < Diagonal(x=4:2)
summary(D)
str(mat2triplet(D))
The basic matrix product, %*%
is implemented for all our
Matrix
and also for
sparseVector
classes, fully analogously to R's
base matrix
and vector objects.
The functions crossprod
and tcrossprod
are
matrix products or “cross products”, ideally implemented
efficiently without computing t(.)
's unnecessarily.
They also return symmetricMatrix
classed
matrices when easily detectable, e.g., in crossprod(m)
, the one
argument case.
tcrossprod()
takes the crossproduct of the transpose of a matrix.
tcrossprod(x)
is formally equivalent to, but faster than, the
call x %*% t(x)
, and so is tcrossprod(x, y)
instead of
x %*% t(y)
.
Boolean matrix products are computed via either
%&%
or boolArith = TRUE
.
## S4 method for signature 'CsparseMatrix,diagonalMatrix'
x %*% y
## S4 method for signature 'dgeMatrix,missing'
crossprod(x, y = NULL, boolArith = NA, ...)
## S4 method for signature 'CsparseMatrix,diagonalMatrix'
crossprod(x, y = NULL, boolArith = NA, ...)
## .... and for many more signatures
## S4 method for signature 'CsparseMatrix,ddenseMatrix'
tcrossprod(x