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/0000-0001-7614-6899>,
SuiteSparse libraries, collaborators listed in
dir(system.file("doc", "SuiteSparse", package="Matrix"),
pattern="License", full.names=TRUE, recursive=TRUE)),
George Karypis [ctb] (<https://orcid.org/0000-0003-2753-1437>, METIS
library, Copyright: Regents of the University of Minnesota),
Jason Riedy [ctb] (<https://orcid.org/0000-0002-4345-4200>, GNU
Octave's condest() and onenormest(), Copyright: Regents of the
University of California),
Jens Oehlschlägel [ctb] (initial nearPD()),
R Core Team [ctb] (base R's matrix implementation)
Classes BunchKaufman and pBunchKaufman represent
Bunch-Kaufman factorizations of n×n real,
symmetric matrices A, having the general form
A=UDUU′=LDLL′
where
DU and DL are symmetric, block diagonal
matrices composed of bU and bL1×1 or 2×2 diagonal blocks;
U=∏k=1bUPkUk
is the product of bU row-permuted unit upper triangular
matrices, each having nonzero entries above the diagonal in 1 or 2 columns;
and
L=∏k=1bLPkLk
is the product of bL row-permuted unit lower triangular
matrices, each having nonzero entries below the diagonal in 1 or 2 columns.
These classes store the nonzero entries of the
2bU+1 or 2bL+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.
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.
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.
Methods
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.
signature(a = "p?BunchKaufman", b = .):
see solve-methods.
Note
In Matrix< 1.6-0, class BunchKaufman extended
dtrMatrix and class pBunchKaufman extended
dtpMatrix, reflecting the fact that the internal
representation of the factorization is fundamentally triangular:
there are n(n+1)/2 “parameters”, and these
can be arranged systematically to form an n×n
triangular matrix.
Matrix1.6-0 removed these extensions so that methods
would no longer be inherited from dtrMatrix and dtpMatrix.
The availability of such methods gave the wrong impression that
BunchKaufman and pBunchKaufman represent a (singular)
matrix, when in fact they represent an ordered set of matrix factors.
The coercions as(., "dtrMatrix") and as(., "dtpMatrix")
are provided for users who understand the caveats.
showClass("BunchKaufman")
set.seed(1)
n <-6L(A <- forceSymmetric(Matrix(rnorm(n * n), n, n)))## With dimnames, to see that they are propagated :
dimnames(A)<- rep.int(list(paste0("x", seq_len(n))),2L)(bk.A <- BunchKaufman(A))
str(e.bk.A <- expand2(bk.A, complete =FALSE), max.level =2L)
str(E.bk.A <- expand2(bk.A, complete =TRUE), max.level =2L)## Underlying LAPACK representation(m.bk.A <- as(bk.A,"dtrMatrix"))
stopifnot(identical(as(m.bk.A,"matrix"), `dim<-`(bk.A@x, bk.A@Dim)))## Number of factors is 2*b+1, b <= n, which can be nontrivial ...(b <-(length(E.bk.A)-1L)%/%2L)
ae1 <-function(a, b,...) all.equal(as(a,"matrix"), as(b,"matrix"),...)
ae2 <-function(a, b,...) ae1(unname(a), unname(b),...)## A ~ U DU U', U := prod(Pk Uk) in floating point
stopifnot(exprs ={
identical(names(e.bk.A), c("U","DU","U."))
identical(e.bk.A[["U"]], Reduce(`%*%`, E.bk.A[seq_len(b)]))
identical(e.bk.A[["U."]], t(e.bk.A[["U"]]))
ae1(A, with(e.bk.A, U %*% DU %*% U.))})## Factorization handled as factorized matrix
b <- rnorm(n)
stopifnot(identical(det(A), det(bk.A)),
identical(solve(A, b), solve(bk.A, b)))
Computes the Bunch-Kaufman factorization of an n×n
real, symmetric matrix A, which has the general form
A=UDUU′=LDLL′
where
DU and DL are symmetric, block diagonal
matrices composed of bU and bL1×1 or 2×2 diagonal blocks;
U=∏k=1bUPkUk
is the product of bU row-permuted unit upper triangular
matrices, each having nonzero entries above the diagonal in 1 or 2 columns;
and
L=∏k=1bLPkLk
is the product of bL row-permuted unit lower triangular
matrices, each having nonzero entries below the diagonal in 1 or 2 columns.
Methods are built on LAPACK routines dsytrf and dsptrf.
Usage
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",...)
Arguments
x
a finite symmetric matrix or
Matrix to be factorized.
If x is square but not symmetric, then it will be
treated as symmetric; see uplo.
warnSing
a logical indicating if a warning should
be signaled for singular x.
uplo
a string, either "U" or "L",
indicating which triangle of x should be used
to compute the factorization.
Generic functions expand1 and expand2,
for constructing matrix factors from the result.
Generic functions Cholesky, Schur,
lu, and qr,
for computing other factorizations.
Examples
showMethods("BunchKaufman", inherited =FALSE)
set.seed(0)
data(CAex, package ="Matrix")
class(CAex)# dgCMatrix
isSymmetric(CAex)# symmetric, but not formally
A <- as(CAex,"symmetricMatrix")
class(A)# dsCMatrix## Have methods for denseMatrix (unpacked and packed),## but not yet sparseMatrix ...## Not run: (bk.A <- BunchKaufman(A))## End(Not run)(bk.A <- BunchKaufman(as(A,"unpackedMatrix")))## A ~ U DU U' in floating point
str(e.bk.A <- expand2(bk.A), max.level =2L)
stopifnot(all.equal(as(A,"matrix"), as(Reduce(`%*%`, e.bk.A),"matrix")))
An example of a sparse matrix for which eigen() seemed
to be difficult, an unscaled version of this has been posted to the
web, accompanying an E-mail to R-help
(https://stat.ethz.ch/mailman/listinfo/r-help), by
Casper J Albers, Open University, UK.
Usage
data(CAex)
Format
This is a 72×72 symmetric matrix with 216
non-zero entries in five bands, stored as sparse matrix of class
dgCMatrix.
Details
Historical note (2006-03-30):
In earlier versions of R, eigen(CAex) fell into an
infinite loop whereas eigen(CAex, EISPACK=TRUE) had been okay.
Examples
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"))
CHMfactor is the virtual class of sparse Cholesky
factorizations of n×n real, symmetric
matrices A, having the general form
P1AP1′=L1DL1′=Djj≥0LL′
or (equivalently)
A=P1′L1DL1′P1=Djj≥0P1′LL′P1
where
P1 is a permutation matrix,
L1 is a unit lower triangular matrix,
D is a diagonal matrix, and
L=L1D.
The second equalities hold only for positive semidefinite A,
for which the diagonal entries of D are non-negative
and D is well-defined.
The implementation of class CHMfactor is based on
CHOLMOD's C-level cholmod_factor_struct. Virtual
subclasses CHMsimpl and CHMsuper separate
the simplicial and supernodal variants. These have nonvirtual
subclasses [dn]CHMsimpl and [dn]CHMsuper,
where prefix ‘d’ and prefix ‘n’ are reserved
for numeric and symbolic factorizations, respectively.
Usage
isLDL(x)
Arguments
x
an object inheriting from virtual class CHMfactor,
almost always the result of a call to generic function
Cholesky.
Value
isLDL(x) returns TRUE or FALSE:
TRUE if x stores the lower triangular entries
of L1−I+D,
FALSE if x stores the lower triangular entries
of L.
an integer vector of length Dim[1]
giving an estimate of the number of nonzero entries in
each column of the lower triangular Cholesky factor.
If symbolic analysis was performed prior to factorization,
then the estimate is exact.
perm
a 0-based integer vector of length Dim[1]
specifying the permutation applied to the rows and columns
of the factorized matrix. perm of length 0 is valid and
equivalent to the identity permutation, implying no pivoting.
type
an integer vector of length 6 specifying
details of the factorization. The elements correspond to
members ordering, is_ll, is_super,
is_monotonic, maxcsize, and maxesize
of the original cholmod_factor_struct.
Simplicial and supernodal factorizations are distinguished
by is_super. Simplicial factorizations do not use
maxcsize or maxesize. Supernodal factorizations
do not use is_ll or is_monotonic.
Of CHMsimpl (all unused by nCHMsimpl):
nz
an integer vector of length Dim[1]
giving the number of nonzero entries in each column of the
lower triangular Cholesky factor. There is at least one
nonzero entry in each column, because the diagonal elements
of the factor are stored explicitly.
p
an integer vector of length Dim[1]+1.
Row indices of nonzero entries in column j of the
lower triangular Cholesky factor are obtained as
i[p[j]+seq_len(nz[j])]+1.
i
an integer vector of length greater than or equal
to sum(nz) containing the row indices of nonzero entries
in the lower triangular Cholesky factor. These are grouped by
column and sorted within columns, but the columns themselves
need not be ordered monotonically. Columns may be overallocated,
i.e., the number of elements of i reserved for column
j may exceed nz[j].
prv, nxt
integer vectors of length
Dim[1]+2 indicating the order in which the columns of
the lower triangular Cholesky factor are stored in i
and x.
Starting from j <- Dim[1]+2,
the recursion j <- nxt[j+1]+1 traverses the columns
in forward order and terminates when nxt[j+1] = -1.
Starting from j <- Dim[1]+1,
the recursion j <- prv[j+1]+1 traverses the columns
in backward order and terminates when prv[j+1] = -1.
Of dCHMsimpl:
x
a numeric vector parallel to i containing
the corresponding nonzero entries of the lower triangular
Cholesky factor Lor (if and only if type[2]
is 0) of the lower triangular matrix L1−I+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.
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.
Methods
coerce
signature(from = "CHMsimpl", to = "dtCMatrix"):
returns a dtCMatrix representing
the lower triangular Cholesky factor Lor
the lower triangular matrix L1−I+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 non-structural zeros above
the main diagonal, hence dgCMatrix is
indeed more appropriate than dtCMatrix
as a coercion target.
determinant
signature(from = "CHMfactor", logarithm = "logical"):
behaves according to an optional argument sqrt.
If sqrt = FALSE, then this method computes the determinant
of the factorized matrix A or its logarithm.
If sqrt = TRUE, then this method computes the determinant
of the factor L=L1sqrt(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.6-0). Calls to this method not setting sqrt
may warn about the pending change. The warnings can be disabled
with options(Matrix.warnSqrtDefault = 0).
diag
signature(x = "CHMfactor"):
returns a numeric vector of length n containing the diagonal
elements of D, which (if they are all non-negative)
are the squared diagonal elements of L.
signature(a = "CHMfactor", b = .):
see solve-methods.
update
signature(object = "CHMfactor"):
returns a copy of object with the same nonzero pattern
but with numeric entries updated according to additional
arguments parent and mult, where parent
is (coercible to) a dsCMatrix or a
dgCMatrix and mult is a numeric
vector of positive length.
The numeric entries are updated with those of the Cholesky
factor of F(parent) + mult[1] * I, i.e.,
F(parent) plus mult[1] times the identity matrix,
where F = identity for symmetric parent
and F = tcrossprod for other parent.
The nonzero pattern of F(parent) must match
that of S if object = Cholesky(S, ...).
updown
signature(update = ., C = ., object = "CHMfactor"):
see updown-methods.
Chen, Y., Davis, T. A., Hager, W. W., & Rajamanickam, S. (2008).
Algorithm 887: CHOLMOD, supernodal sparse Cholesky factorization
and update/downdate.
ACM Transactions on Mathematical Software,
35(3), Article 22, 1-14.
doi:10.1145/1391989.1391995
Amestoy, P. R., Davis, T. A., & Duff, I. S. (2004).
Algorithm 837: AMD, an approximate minimum degree ordering algorithm.
ACM Transactions on Mathematical Software,
17(4), 886-905.
doi:10.1145/1024074.1024081
Golub, G. H., & Van Loan, C. F. (2013).
Matrix computations (4th ed.).
Johns Hopkins University Press.
doi:10.56021/9781421407944
Classes Cholesky and pCholesky represent
dense, pivoted Cholesky factorizations of n×n
real, symmetric, positive semidefinite matrices A,
having the general form
P1AP1′=L1DL1′=LL′
or (equivalently)
A=P1′L1DL1′P1=P1′LL′P1
where
P1 is a permutation matrix,
L1 is a unit lower triangular matrix,
D is a non-negative diagonal matrix, and
L=L1D.
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.
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 column-major
order.
perm
a 1-based integer vector of length Dim[1]
specifying the permutation applied to the rows and columns
of the factorized matrix. perm of length 0 is valid and
equivalent to the identity permutation, implying no pivoting.
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).
Methods
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.
signature(a = "p?Cholesky", b = .):
see solve-methods.
Note
In Matrix< 1.6-0, class Cholesky extended
dtrMatrix and class pCholesky extended
dtpMatrix, reflecting the fact that the factor
L is indeed a triangular matrix.
Matrix1.6-0 removed these extensions so that methods
would no longer be inherited from dtrMatrix and dtpMatrix.
The availability of such methods gave the wrong impression that
Cholesky and pCholesky represent a (singular)
matrix, when in fact they represent an ordered set of matrix factors.
The coercions as(., "dtrMatrix") and as(., "dtpMatrix")
are provided for users who understand the caveats.
Computes the pivoted Cholesky factorization of an
n×n real, symmetric matrix A,
which has the general form
P1AP1′=L1DL1′=Djj≥0LL′
or (equivalently)
A=P1′L1DL1′P1=Djj≥0P1′LL′P1
where
P1 is a permutation matrix,
L1 is a unit lower triangular matrix,
D is a diagonal matrix, and
L=L1D.
The second equalities hold only for positive semidefinite A,
for which the diagonal entries of D are non-negative
and D is well-defined.
Methods for denseMatrix are built on
LAPACK routines dpstrf, dpotrf, and dpptrf.
The latter two do not permute rows or columns,
so that P1 is an identity matrix.
Methods for sparseMatrix are built on
CHOLMOD routines cholmod_analyze and cholmod_factorize_p.
Usage
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",...)
Arguments
A
a finite, symmetric matrix or
Matrix to be factorized. If A
is square but not symmetric, then it will be treated
as symmetric; see uplo.
Methods for dense A require positive definiteness
when perm = FALSE and positive semidefiniteness
when perm = TRUE.
Methods for sparse A require positive definiteness
when LDL = TRUE and nonzero leading principal minors
(after pivoting) when LDL = FALSE.
Methods for sparse, diagonalA are an exception,
requiring positive semidefiniteness unconditionally.
perm
a logical indicating if the rows and columns
of A should be pivoted. Methods for sparse A
employ the approximate minimum degree (AMD) algorithm
in order to reduce fill-in, i.e., without regard for
numerical stability.
Pivoting for sparsity may introduce nonpositive leading
principal minors, causing the factorization to fail, in
which case it may be necessary to set perm = FALSE.
tol
a finite numeric tolerance,
used only if perm = TRUE.
The factorization algorithm stops if the pivot is less than
or equal to tol. Negative tol is equivalent
to nrow(A) * .Machine$double.eps * max(diag(A)).
LDL
a logical indicating if the simplicial factorization
should be computed as
P1′L1DL1′P1,
such that the result stores the lower triangular entries
of L1−I+D.
The alternative is P1′LL′P1,
such that the result stores the lower triangular entries
of L=L1D. This argument
is ignored if super = TRUE (or if super = NA
and the supernodal algorithm is chosen), as the supernodal
code does not yet support the LDL = TRUE variant.
super
a logical indicating if the factorization should
use the supernodal algorithm. The alternative is the simplicial
algorithm. Setting super = NA leaves the choice to
a CHOLMOD-internal heuristic.
Imult
a finite number. The matrix
that is factorized is A + Imult * diag(nrow(A)),
i.e., A plus Imult times the identity matrix.
This argument is useful for symmetric, indefinite A,
as Imult > max(rowSums(abs(A)) - diag(abs(A))) ensures
that A + Imult * diag(nrow(A)) is diagonally dominant.
(Symmetric, diagonally dominant matrices are positive definite.)
uplo
a string, either "U" or "L",
indicating which triangle of A should be used
to compute the factorization. The default is "U",
even for lower triangular A, to be consistent with
chol from base.
...
further arguments passed to or from methods.
Details
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′LL′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
∥A∥∥A−P′LL′P∥.
See the examples and LAPACK Working Note (“LAWN”) 161
for details.
Chen, Y., Davis, T. A., Hager, W. W., & Rajamanickam, S. (2008).
Algorithm 887: CHOLMOD, supernodal sparse Cholesky factorization
and update/downdate.
ACM Transactions on Mathematical Software,
35(3), Article 22, 1-14.
doi:10.1145/1391989.1391995
Amestoy, P. R., Davis, T. A., & Duff, I. S. (2004).
Algorithm 837: AMD, an approximate minimum degree ordering algorithm.
ACM Transactions on Mathematical Software,
17(4), 886-905.
doi:10.1145/1024074.1024081
Golub, G. H., & Van Loan, C. F. (2013).
Matrix computations (4th ed.).
Johns Hopkins University Press.
doi:10.56021/9781421407944
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.
Examples
showMethods("Cholesky", inherited =FALSE)
set.seed(0)## ---- Dense ----------------------------------------------------------## .... Positive definite ..............................................
n <-6L(A1 <- crossprod(Matrix(rnorm(n * n), n, n)))(ch.A1.nopivot <- Cholesky(A1, perm =FALSE))(ch.A1 <- Cholesky(A1))
stopifnot(exprs ={
length(ch.A1@perm)== ncol(A1)
isPerm(ch.A1@perm)
is.unsorted(ch.A1@perm)# typically not the identity permutation
length(ch.A1.nopivot@perm)==0L})## A ~ P1' L D L' P1 ~ P1' L L' P1 in floating point
str(e.ch.A1 <- expand2(ch.A1, LDL =TRUE), max.level =2L)
str(E.ch.A1 <- expand2(ch.A1, LDL =FALSE), max.level =2L)
stopifnot(exprs ={
all.equal(as(A1,"matrix"), as(Reduce(`%*%`, e.ch.A1),"matrix"))
all.equal(as(A1,"matrix"), as(Reduce(`%*%`, E.ch.A1),"matrix"))})## .... Positive semidefinite but not positive definite ................
A2 <- A1
A2[1L,]<- A2[,1L]<-0
A2
try(Cholesky(A2, perm =FALSE))# fails as not positive definite
ch.A2 <- Cholesky(A2)# returns, with a warning and ...
A2.hat <- Reduce(`%*%`, expand2(ch.A2, LDL =FALSE))
norm(A2 - A2.hat,"2")/ norm(A2,"2")# 7.670858e-17## .... Not positive semidefinite ......................................
A3 <- A1
A3[1L,]<- A3[,1L]<--1
A3
try(Cholesky(A3, perm =FALSE))# fails as not positive definite
ch.A3 <- Cholesky(A3)# returns, with a warning and ...
A3.hat <- Reduce(`%*%`, expand2(ch.A3, LDL =FALSE))
norm(A3 - A3.hat,"2")/ norm(A3,"2")# 1.781568## Indeed, 'A3' is not positive semidefinite, but 'A3.hat' _is_
ch.A3.hat <- Cholesky(A3.hat)
A3.hat.hat <- Reduce(`%*%`, expand2(ch.A3.hat, LDL =FALSE))
norm(A3.hat - A3.hat.hat,"2")/ norm(A3.hat,"2")# 1.777944e-16## ---- Sparse ---------------------------------------------------------## Really just three cases modulo permutation :#### type factorization minors of P1 A P1'## 1 simplicial P1 A P1' = L1 D L1' nonzero## 2 simplicial P1 A P1' = L L ' positive## 3 supernodal P1 A P2' = L L ' positive
data(KNex, package ="Matrix")
A4 <- crossprod(KNex[["mm"]])
ch.A4 <-
list(pivoted =
list(simpl1 = Cholesky(A4, perm =TRUE, super =FALSE, LDL =TRUE),
simpl0 = Cholesky(A4, perm =TRUE, super =FALSE, LDL =FALSE),
super0 = Cholesky(A4, perm =TRUE, super =TRUE)),
unpivoted =
list(simpl1 = Cholesky(A4, perm =FALSE, super =FALSE, LDL =TRUE),
simpl0 = Cholesky(A4, perm =FALSE, super =FALSE, LDL =FALSE),
super0 = Cholesky(A4, perm =FALSE, super =TRUE)))
ch.A4
s <- simplify2array
rapply2 <-function(object, f,...) rapply(object, f,,, how ="list",...)
s(rapply2(ch.A4, isLDL))
s(m.ch.A4 <- rapply2(ch.A4, expand1,"L"))# giving L = L1 sqrt(D)## By design, the pivoted and simplicial factorizations## are more sparse than the unpivoted and supernodal ones ...
s(rapply2(m.ch.A4, object.size))## Which is nicely visualized by lattice-based methods for 'image'
inm <- c("pivoted","unpivoted")
jnm <- c("simpl1","simpl0","super0")for(i in1:2)for(j in1:3)
print(image(m.ch.A4[[c(i, j)]], main = paste(inm[i], jnm[j])),
split = c(j, i,3L,2L), more = i * j <6L)
simpl1 <- ch.A4[[c("pivoted","simpl1")]]
stopifnot(exprs ={
length(simpl1@perm)== ncol(A4)
isPerm(simpl1@perm,0L)
is.unsorted(simpl1@perm)# typically not the identity permutation})## One can expand with and without D regardless of isLDL(.),## but "without" requires L = L1 sqrt(D), which is conditional## on min(diag(D)) >= 0, hence "with" is the default
isLDL(simpl1)
stopifnot(min(diag(simpl1))>=0)
str(e.ch.A4 <- expand2(simpl1, LDL =TRUE), max.level =2L)# default
str(E.ch.A4 <- expand2(simpl1, LDL =FALSE), max.level =2L)
stopifnot(exprs ={
all.equal(E.ch.A4[["L"]], e.ch.A4[["L1"]]%*% sqrt(e.ch.A4[["D"]]))
all.equal(E.ch.A4[["L."]], sqrt(e.ch.A4[["D"]])%*% e.ch.A4[["L1."]])
all.equal(A4, as(Reduce(`%*%`, e.ch.A4),"symmetricMatrix"))
all.equal(A4, as(Reduce(`%*%`, E.ch.A4),"symmetricMatrix"))})## The "same" permutation matrix with "alternate" representation## [i, perm[i]] {margin=1} <-> [invertPerm(perm)[j], j] {margin=2}
alt <-function(P){
P@margin <-1L+!(P@margin -1L)# 1 <-> 2
P@perm <- invertPerm(P@perm)
P
}## Expansions are elegant but inefficient (transposes are redundant)## hence programmers should consider methods for 'expand1' and 'diag'
stopifnot(exprs ={
identical(expand1(simpl1,"P1"), alt(e.ch.A4[["P1"]]))
identical(expand1(simpl1,"L"), E.ch.A4[["L"]])
identical(Diagonal(x = diag(simpl1)), e.ch.A4[["D"]])})## chol(A, pivot = value) is a simple wrapper around## Cholesky(A, perm = value, LDL = FALSE, super = FALSE),## returning L' = sqrt(D) L1' _but_ giving no information## about the permutation P1
selectMethod("chol","dsCMatrix")
stopifnot(all.equal(chol(A4, pivot =TRUE), E.ch.A4[["L."]]))## Now a symmetric matrix with positive _and_ negative eigenvalues,## hence _not_ positive semidefinite
A5 <- new("dsCMatrix",
Dim = c(7L,7L),
p = c(0:1,3L,6:7,10:11,15L),
i = c(0L,0:1,0:3,2:5,3:6),
x = c(1,6,38,10,60,103,-4,6,-32,-247,-2,-16,-128,-2,-67))(ev <- eigen(A5, only.values =TRUE)$values)(t.ev <- table(factor(sign(ev),-1:1)))# the matrix "inertia"
ch.A5 <- Cholesky(A5)
isLDL(ch.A5)(d.A5 <- diag(ch.A5))# diag(D) is partly negative## Sylvester's law of inertia holds here, but not in general## in finite precision arithmetic
stopifnot(identical(table(factor(sign(d.A5),-1:1)), t.ev))
try(expand1(ch.A5,"L"))# unable to compute L = L1 sqrt(D)
try(expand2(ch.A5, LDL =FALSE))# ditto
try(chol(A5, pivot =TRUE))# ditto## The default expansion is "square root free" and still works here
str(e.ch.A5 <- expand2(ch.A5, LDL =TRUE), max.level =2L)
stopifnot(all.equal(A5, as(Reduce(`%*%`, e.ch.A5),"symmetricMatrix")))## Version of the SuiteSparse library, which includes CHOLMOD
Mv <- Matrix.Version()
Mv[["suitesparse"]]
The "CsparseMatrix" class is the virtual class of
all sparse matrices coded in sorted compressed column-oriented form.
Since it is a virtual class, no objects may be created from it. See
showClass("CsparseMatrix") for its subclasses.
Slots
i:
Object of class "integer" of length nnzero
(number of non-zero elements). These are the 0-based row numbers for
each non-zero element in the matrix, i.e., i must be in
0:(nrow(.)-1).
p:
integer vector for providing pointers, one
for each column, to the initial (zero-based) index of elements in
the column. .@p is of length ncol(.) + 1, with
p[1] == 0 and p[length(p)] == nnzero, such that in
fact, diff(.@p) are the number of non-zero elements for
each column.
In other words, m@p[1:ncol(m)] contains the indices of
those elements in m@x that are the first elements in the
respective column of m.
Dim, Dimnames:
inherited from
the superclass, see the sparseMatrix class.
Extends
Class "sparseMatrix", directly.
Class "Matrix", by class "sparseMatrix".
Methods
matrix products %*%,
crossprod() and tcrossprod(),
several solve methods,
and other matrix methods available:
signature(from = "CsparseMatrix", to = "TsparseMatrix"): ...
coerce
signature(from = "CsparseMatrix", to = "denseMatrix"): ...
coerce
signature(from = "CsparseMatrix", to = "matrix"): ...
coerce
signature(from = "TsparseMatrix", to = "CsparseMatrix"): ...
coerce
signature(from = "denseMatrix", to = "CsparseMatrix"): ...
diag
signature(x = "CsparseMatrix"): ...
gamma
signature(x = "CsparseMatrix"): ...
lgamma
signature(x = "CsparseMatrix"): ...
log
signature(x = "CsparseMatrix"): ...
t
signature(x = "CsparseMatrix"): ...
tril
signature(x = "CsparseMatrix"): ...
triu
signature(x = "CsparseMatrix"): ...
Note
All classes extending CsparseMatrix have a common validity
(see validObject) check function. That function
additionally checks the i slot for each column to contain
increasing row numbers.
In earlier versions of Matrix (<= 0.999375-16),
validObject automatically re-sorted the entries when
necessary, and hence new() calls with somewhat permuted
i and x slots worked, as new(...)
(with slot arguments) automatically checks the validity.
Now, you have to use sparseMatrix to achieve the same
functionality or know how to use .validateCsparse() to do so.
Construct a formally diagonal Matrix,
i.e., an object inheriting from virtual class
diagonalMatrix
(or, if desired, a mathematically diagonal
CsparseMatrix).
Usage
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)
Arguments
n
integer indicating the dimension of the (square) matrix.
If missing, then length(x) is used.
x
numeric or logical vector listing values for the diagonal
entries, to be recycled as necessary. If NULL (the default),
then the result is a unit diagonal matrix. .sparseDiagonal()
and friends ignore non-NULLx when kind = "n".
names
either logicalTRUE or FALSE or
then a character vector of lengthn. If true andnames(x) is not
NULL, use that as both row and column names for the resulting
matrix. When a character vector, use it for both dimnames.
uplo
one of c("U","L"), specifying the uplo slot
of the result if the result is formally triangular of symmetric.
shape
one of c("t","s","g"), indicating if the result
should be formally triangular, symmetric, or “general”.
The result will inherit from virtual class
triangularMatrix,
symmetricMatrix, or
generalMatrix, respectively.
unitri
logical indicating if a formally triangular result with
ones on the diagonal should be formally unit triangular, i.e.,
with diag slot equal to "U" rather than "N".
kind
one of c("d","l","n"), indicating the “mode”
of the result: numeric, logical, or pattern.
The result will inherit from virtual class
dsparseMatrix,
lsparseMatrix, or
nsparseMatrix, respectively.
Values other than "n" are ignored when x is
non-NULL; in that case the mode is determined by
typeof(x).
cols
optional integer vector with values in 0:(n-1),
indexing columns of the specified diagonal matrix. If specified,
then the result is (mathematically) D[, cols+1] rather
than D, where D = Diagonal(n, x), and it is always
“general” (i.e., shape is ignored).
Value
Diagonal() returns an object inheriting from virtual class
diagonalMatrix.
.sparseDiagonal() returns a CsparseMatrix
representation of Diagonal(n, x) or, if cols is given,
of Diagonal(n, x)[, cols+1]. The precise class of the result
depends on shape and kind.
.trDiagonal() and .symDiagonal() are simple wrappers,
for .sparseDiagonal(shape = "t") and
.sparseDiagonal(shape = "s"), respectively.
.sparseDiagonal() exists primarily to leverage efficient
C-level methods available for CsparseMatrix.
Author(s)
Martin Maechler
See Also
the generic function diag for extraction
of the diagonal from a matrix works for all “Matrices”.
bandSparse constructs a banded sparse matrix from
its non-zero sub-/super - diagonals. band(A) returns a
band matrix containing some sub-/super - diagonals of A.
Generate the n by n symmetric Hilbert matrix. Because
these matrices are ill-conditioned for moderate to large n,
they are often used for testing numerical linear algebra code.
Usage
Hilbert(n)
Arguments
n
a non-negative integer.
Value
the n by n symmetric Hilbert matrix as a
"dpoMatrix" object.
A model matrix mm and corresponding response vector y
used in an example by Koenker and Ng. The matrix mm is a sparse
matrix with 1850 rows and 712 columns but only 8758 non-zero entries.
It is a "dgCMatrix" object. The vector y is just
numeric of length 1850.
Usage
data(KNex)
References
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
Examples
data(KNex, package ="Matrix")
class(KNex$mm)
dim(KNex$mm)
image(KNex$mm)
str(KNex)
system.time(# a fraction of a second
sparse.sol <- with(KNex, solve(crossprod(mm), crossprod(mm, y))))
head(round(sparse.sol,3))## Compare with QR-based solution ("more accurate, but slightly slower"):
system.time(
sp.sol2 <- with(KNex, qr.coef(qr(mm), y)))
all.equal(sparse.sol, sp.sol2, tolerance =1e-13)# TRUE
Computes Khatri-Rao products for any kind of matrices.
The Khatri-Rao product is a column-wise Kronecker product. Originally
introduced by Khatri and Rao (1968), it has many different applications,
see Liu and Trenkler (2008) for a survey. Notably, it is used in
higher-dimensional tensor decompositions, see Bader and Kolda (2008).
Usage
KhatriRao(X, Y = X, FUN ="*", sparseY =TRUE, make.dimnames =FALSE)
Arguments
X, Y
matrices of with the same number of columns.
FUN
the (name of the) function to be used for the
column-wise Kronecker products, see kronecker,
defaulting to the usual multiplication.
sparseY
logical specifying if Y should be coerced and
treated as sparseMatrix. Set this to
FALSE, e.g., to distinguish structural zeros from zero entries.
make.dimnames
logical indicating if the result should inherit
dimnames from X and Y in a simple way.
Value
a "CsparseMatrix", say R, the Khatri-Rao
product of X (n×k) and Y (m×k), is of dimension (n⋅m)×k,
where the j-th column, R[,j] is the kronecker product
kronecker(X[,j], Y[,j]).
Note
The current implementation is efficient for large sparse matrices.
Author(s)
Original by Michael Cysouw, Univ. Marburg;
minor tweaks, bug fixes etc, by Martin Maechler.
References
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 A30, 167–180.
Bader, Brett W, and Tamara G Kolda (2008)
Efficient MATLAB Computations with Sparse and Factored Tensors.
SIAM J. Scientific Computing30, 205–231.
when data is not a matrix, the desired number of rows
ncol
when data is not a matrix, the desired number of columns
byrow
logical. If FALSE (the default) the matrix is
filled by columns, otherwise the matrix is filled by rows.
dimnames
a dimnames attribute for the matrix: a
list of two character components. They are set if not
NULL (as per default).
sparse
logical or NULL, specifying if the result should
be sparse or not. By default, it is made sparse when more than half
of the entries are 0.
doDiag
logical indicating if a diagonalMatrix
object should be returned when the resulting matrix is diagonal
(mathematically). As class diagonalMatrixextendssparseMatrix, this is a
natural default for all values of sparse.
Otherwise, if doDiag is false, a dense or sparse (depending on
sparse) symmetric matrix will be returned.
forceCheck
logical indicating if the checks for structure
should even happen when data is already a "Matrix"
object.
Details
If either of nrow or ncol is not given, an attempt is
made to infer it from the length of data and the other
parameter.
Further, Matrix() makes efforts to keep logical
matrices logical, i.e., inheriting from class lMatrix,
and to determine specially structured matrices such as symmetric,
triangular or diagonal ones. Note that a symmetric matrix also
needs symmetric dimnames, e.g., by specifying
dimnames = list(NULL,NULL), see the examples.
Most of the time, the function works via a traditional (full)
matrix. However, Matrix(0, nrow,ncol) directly
constructs an “empty” sparseMatrix, as does
Matrix(FALSE, *).
Although it is sometime possible to mix unclassed matrices (created
with matrix) with ones of class "Matrix", it is much
safer to always use carefully constructed ones of class
"Matrix".
Value
Returns matrix of a class that inherits from "Matrix".
Only if data is not a matrix and does not already inherit
from class Matrix are the arguments
nrow, ncol and byrow made use of.
The Matrix class is a class contained by all actual
classes in the Matrix package. It is a “virtual” class.
Slots
Dim
an integer vector of length 2 giving the
dimensions of the matrix.
Dimnames
a list of length 2. Each element must
be NULL or a character vector of length equal to the
corresponding element of Dim.
Methods
determinant
signature(x = "Matrix", logarithm = "missing"): and
determinant
signature(x = "Matrix", logarithm = "logical"):
compute the (log) determinant of x. The method
chosen depends on the actual Matrix class of x. Note that
det also works for all our matrices, calling the
appropriate determinant() method. The Matrix::det
is an exact copy of base::det, but in the correct
namespace, and hence calling the S4-aware version of
determinant().).
diff
signature(x = "Matrix"): As diff()
for traditional matrices, i.e., applying diff() to each
column.
signature(x = "Matrix", value = "ANY"): where
value is integer of length 2. Allows to reshape
Matrix objects, but only when prod(value) == prod(dim(x)).
signature(x = "Matrix", value = "list"): set
the dimnames to a list of length 2, see
dimnames<-.
length
signature(x = "Matrix"): simply defined as
prod(dim(x)) (and hence of mode "double").
show
signature(object = "Matrix"): show
method for printing. For printing sparse
matrices, see printSpMatrix.
zapsmall
signature(x = "Matrix"): typically used for
"dMatrix": round() matrix entries
such that (relatively) very small entries become zero exactly.
image
signature(object = "Matrix"): draws an
image of the matrix entries, using
levelplot() from package lattice.
head
signature(object = "Matrix"): return only the
“head”, i.e., the first few rows.
tail
signature(object = "Matrix"): return only the
“tail”, i.e., the last few rows of the respective matrix.
as.matrix, as.array
signature(x = "Matrix"): the same as
as(x, "matrix"); see also the note below.
as.vector
signature(x = "Matrix", mode = "missing"):
as.vector(m) should be identical to as.vector(as(m,
"matrix")), implemented more efficiently for some subclasses.
as(x, "vector"), as(x, "numeric")
etc, similarly.
coerce
signature(from = "ANY", to = "Matrix"): This
relies on a correct as.matrix() method for from.
There are many more methods that (conceptually should) work for all
"Matrix" objects, e.g., colSums,
rowMeans. Even base functions may work
automagically (if they first call as.matrix() on their
principal argument), e.g., apply, eigen,
svd or kappa all do work via coercion to a
“traditional” (dense) matrix.
Note
Loading the Matrix namespace “overloads”
as.matrix and as.array in the base
namespace by the equivalent of function(x) as(x, "matrix").
Consequently, as.matrix(m) or as.array(m) will properly
work when m inherits from the "Matrix" class —
also for functions in package base and other packages.
E.g., apply or outer can therefore be applied
to "Matrix" matrices.
slotNames("Matrix")
cl <- getClass("Matrix")
names(cl@subclasses)# more than 40 ..
showClass("Matrix")#> output with slots and all subclasses(M <- Matrix(c(0,1,0,0),6,4))
dim(M)
diag(M)
cm <- M[1:4,]+10*Diagonal(4)
diff(M)## can reshape it even :
dim(M)<- c(2,12)
M
stopifnot(identical(M, Matrix(c(0,1,0,0),2,12)),
all.equal(det(cm),
determinant(as(cm,"matrix"), log=FALSE)$modulus,
check.attributes=FALSE))
MatrixFactorization is the virtual class of
factorizations of m×n matrices A,
having the general form
P1AP2=A1⋯Ap
or (equivalently)
A=P1′A1⋯ApP2′
where P1 and P2 are permutation matrices.
Factorizations requiring symmetric A have the constraint
P2=P1′, and factorizations without row
or column pivoting have the constraints
P1=Im and P2=In,
where Im and In are the
m×m and n×n identity matrices.
CholeskyFactorization, BunchKaufmanFactorization,
SchurFactorization, LU, and QR are the virtual
subclasses of MatrixFactorization containing all Cholesky,
Bunch-Kaufman, Schur, LU, and QR factorizations, respectively.
Slots
Dim
an integer vector of length 2 giving the
dimensions of the factorized matrix.
Dimnames
a list of length 2 preserving the
dimnames of the factorized matrix. Each element
must be NULL or a character vector of length equal
to the corresponding element of Dim.
Methods
determinant
signature(x = "MatrixFactorization", logarithm = "missing"):
sets logarithm = TRUE and recalls the generic function.
The "RsparseMatrix" class is the virtual class of
all sparse matrices coded in sorted compressed row-oriented form.
Since it is a virtual class, no objects may be created from it. See
showClass("RsparseMatrix") for its subclasses.
Slots
j:
Object of class "integer" of length nnzero
(number of non-zero elements). These are the row numbers for
each non-zero element in the matrix.
p:
Object of class "integer" of pointers, one
for each row, to the initial (zero-based) index of elements in
the row.
Class "sparseMatrix", directly.
Class "Matrix", by class "sparseMatrix".
Methods
Originally, few methods were defined on purpose, as we
rather use the CsparseMatrix in Matrix.
Then, more methods were added but beware that these
typically do not return "RsparseMatrix" results, but
rather Csparse* or Tsparse* ones; e.g., R[i, j] <- v for an
"RsparseMatrix"R works, but after the assignment, R
is a (triplet) "TsparseMatrix".
t
signature(x = "RsparseMatrix"): ...
coerce
signature(from = "RsparseMatrix", to = "CsparseMatrix"): ...
coerce
signature(from = "RsparseMatrix", to = "TsparseMatrix"): ...
See Also
its superclass, sparseMatrix, and, e.g., class
dgRMatrix for the links to other classes.
Schur is the class of Schur factorizations of
n×n real matrices A,
having the general form
A=QTQ′
where
Q is an orthogonal matrix and
T is a block upper triangular matrix with
1×1 or 2×2 diagonal blocks
specifying the real and complex conjugate eigenvalues of A.
The column vectors of Q are the Schur vectors of A,
and T is the Schur form of A.
The Schur factorization generalizes the spectral decomposition
of normal matrices A, whose Schur form is block diagonal,
to arbitrary square matrices.
Details
The matrix A and its Schur form T are similar
and thus have the same spectrum. The eigenvalues are computed
trivially as the eigenvalues of the diagonal blocks of T.
an orthogonal matrix,
inheriting from virtual class Matrix.
T
a block upper triangular matrix,
inheriting from virtual class Matrix.
The diagonal blocks have dimensions 1-by-1 or 2-by-2.
EValues
a numeric or complex vector containing
the eigenvalues of the diagonal blocks of T, which are
the eigenvalues of T and consequently of the factorized
matrix.
Objects can be generated directly by calls of the form
new("Schur", ...), but they are more typically obtained
as the value of Schur(x) for x inheriting from
Matrix (often dgeMatrix).
Methods
determinant
signature(from = "Schur", logarithm = "logical"):
computes the determinant of the factorized matrix A
or its logarithm.
showClass("Schur")
set.seed(0)
n <-4L(A <- Matrix(rnorm(n * n), n, n))## With dimnames, to see that they are propagated :
dimnames(A)<- list(paste0("r", seq_len(n)),
paste0("c", seq_len(n)))(sch.A <- Schur(A))
str(e.sch.A <- expand2(sch.A), max.level =2L)## A ~ Q T Q' in floating point
stopifnot(exprs ={
identical(names(e.sch.A), c("Q","T","Q."))
all.equal(A, with(e.sch.A, Q %*% T %*% Q.))})## Factorization handled as factorized matrix
b <- rnorm(n)
stopifnot(all.equal(det(A), det(sch.A)),
all.equal(solve(A, b), solve(sch.A, b)))## One of the non-general cases:
Schur(Diagonal(6L))
Computes the Schur factorization of an n×n
real matrix A, which has the general form
A=QTQ′
where
Q is an orthogonal matrix and
T is a block upper triangular matrix with
1×1 and 2×2 diagonal blocks
specifying the real and complex conjugate eigenvalues of A.
The column vectors of Q are the Schur vectors of A,
and T is the Schur form of A.
Methods are built on LAPACK routine dgees.
Usage
Schur(x, vectors =TRUE,...)
Arguments
x
a finite square matrix or
Matrix to be factorized.
vectors
a logical. If TRUE (the default),
then Schur vectors are computed in addition to the Schur form.
...
further arguments passed to or from methods.
Value
An object representing the factorization, inheriting
from virtual class SchurFactorization
if vectors = TRUE. Currently, the specific class
is always Schur in that case.
An exception is if x is a traditional matrix,
in which case the result is a named list containing
Q, T, and EValues slots of the
Schur object.
If vectors = FALSE, then the result is the same
named list but without Q.
The "TsparseMatrix" class is the virtual class of
all sparse matrices coded in triplet form. Since it is a virtual class,
no objects may be created from it. See
showClass("TsparseMatrix") for its subclasses.
Object of class "integer" - the row indices
of non-zero entries in 0-base, i.e., must be in
0:(nrow(.)-1).
j:
Object of class "integer" - the column
indices of non-zero entries. Must be the same length as slot
i and 0-based as well, i.e., in
0:(ncol(.)-1). For numeric Tsparse matrices, (i,j)
pairs can occur more than once, see dgTMatrix.
Extends
Class "sparseMatrix", directly.
Class "Matrix", by class "sparseMatrix".
Most operations with sparse matrices are performed using the
compressed, column-oriented or CsparseMatrix
representation. The triplet representation is convenient for
creating a sparse matrix or for reading and writing such
matrices. Once it is created, however, the matrix is generally
coerced to a CsparseMatrix for further
operations.
Note that all new(.), spMatrix and
sparseMatrix(*, repr="T") constructors
for "TsparseMatrix" classes implicitly add (i.e., “sum up”)
xk's that belong to identical (ik,jk) pairs, see, the
example below, or also "dgTMatrix".
For convenience, methods for some operations such as %*%
and crossprod are defined for
TsparseMatrix objects. These methods simply
coerce the TsparseMatrix object to a
CsparseMatrix object then perform the
operation.
See Also
its superclass, sparseMatrix, and the
dgTMatrix class, for the links to other classes.
Examples
showClass("TsparseMatrix")## or just the subclasses' names
names(getClass("TsparseMatrix")@subclasses)
T3 <- spMatrix(3,4, i=c(1,3:1), j=c(2,4:2), x=1:4)
T3 # only 3 non-zero entries, 5 = 1+4 !
This matrix gives the contiguities of 3111 U.S. counties,
using the queen criterion of at least one shared vertex
or edge.
Usage
data(USCounties)
Format
A 3111×3111 sparse, symmetric
matrix of class dsCMatrix, with 9101
nonzero entries.
Source
GAL lattice file ‘usc_q.GAL’
(retrieved in 2008 from
‘http://sal.uiuc.edu/weights/zips/usc.zip’
with permission from Luc Anselin for use and distribution)
was read into R using function read.gal
from package spdep.
Neighbour lists were augmented with row-standardized
(and then symmetrized) spatial weights, using functions
nb2listw and similar.listw from packages
spdep and spatialreg.
The resulting listw object was coerced to class
dsTMatrix
using as_dsTMatrix_listw from spatialreg,
and subsequently to class dsCMatrix.
References
Ord, J. K. (1975).
Estimation methods for models of spatial interaction.
Journal of the American Statistical Association,
70(349), 120-126.
doi:10.2307/2285387
Examples
data(USCounties, package ="Matrix")(n <- ncol(USCounties))
I <- .symDiagonal(n)
set.seed(1)
r <-50L
rho <-1/ runif(r,0,0.5)
system.time(MJ0 <- sapply(rho,function(mult)
determinant(USCounties + mult * I, logarithm =TRUE)$modulus))## Can be done faster by updating the Cholesky factor:
C1 <- Cholesky(USCounties, Imult =2)
system.time(MJ1 <- sapply(rho,function(mult)
determinant(update(C1, USCounties, mult), sqrt =FALSE)$modulus))
stopifnot(all.equal(MJ0, MJ1))
C2 <- Cholesky(USCounties, super =TRUE, Imult =2)
system.time(MJ2 <- sapply(rho,function(mult)
determinant(update(C2, USCounties, mult), sqrt =FALSE)$modulus))
stopifnot(all.equal(MJ0, MJ2))
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 from the Class
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.
Slots
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:
object of class "rleDiff",
used for compression via rle.
Methods
as.numeric, as.integer, as.vector
signature(x = "abIndex"): ...
[
signature(x = "abIndex", i = "index", j = "ANY", drop = "ANY"): ...
coerce
signature(from = "numeric", to = "abIndex"): ...
coerce
signature(from = "abIndex", to = "numeric"): ...
coerce
signature(from = "abIndex", to = "integer"): ...
length
signature(x = "abIndex"): ...
Ops
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.
Ops
signature(e1 = "abIndex", e2 = "abIndex"): ...
Ops
signature(e1 = "abIndex", e2 = "numeric"): ...
Summary
signature(x = "abIndex"): ...
show
("abIndex"): simple show method,
building on show(<rleDiff>).
is.na
("abIndex"): works analogously to regular vectors.
is.finite, is.infinite
("abIndex"): ditto.
Note
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.
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.
Usage
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(...)
Arguments
from, to
the starting and (maximal) end value of the sequence.
by
number: increment of the sequence.
length.out
desired length of the sequence. A
non-negative number, which for seq and seq.int will be
rounded up if fractional.
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 "abIndex" vector, these
arguments will be concatenated to a new "abIndex" object.
Value
An abstract index vector, i.e., object of class
"abIndex".
See Also
the class abIndex documentation;
rep2abI() for another constructor;
rle (base).
Examples
stopifnot(identical(-3:20,
as(abIseq1(-3,20),"vector")))
try(## (arithmetic) not yet implemented
abIseq(1,50, by =3))
Detect or standardize a TsparseMatrix with
unsorted or duplicated (i,j) pairs.
Usage
anyDuplicatedT(x,...)
isUniqueT(x, byrow =FALSE, isT = is(x,"TsparseMatrix"))
asUniqueT(x, byrow =FALSE, isT = is(x,"TsparseMatrix"))
aggregateT(x)
Arguments
x
an R object. anyDuplicatedT and aggregateT
require x inheriting from TsparseMatrix.
asUniqueT requires x inheriting from
Matrix and coerces x
to TsparseMatrix if necessary.
...
optional arguments passed to the default method for
generic function anyDuplicated.
byrow
a logical indicating if x should be sorted
by row then by column.
isT
a logical indicating if x inherits from virtual
class TsparseMatrix.
Value
anyDuplicatedT(x) returns the index of the first duplicated
(i,j) pair in x (0 if there are no duplicated pairs).
isUniqueT(x) returns TRUE if x is a
TsparseMatrix with sorted, nonduplicated
(i,j) pairs and FALSE otherwise.
asUniqueT(x) returns the unique
TsparseMatrix representation of x with
sorted, nonduplicated (i,j) pairs. Values corresponding to
identical (i,j) pairs are aggregated by addition, where in the
logical case “addition” refers to logical OR.
Return the matrix obtained by setting to zero elements below a diagonal
(triu), above a diagonal (tril), or outside of a general
band (band).
Usage
band(x, k1, k2,...)
triu(x, k =0L,...)
tril(x, k =0L,...)
Arguments
x
a matrix-like object
k, k1, k2
integers specifying the diagonals that are not set to
zero, k1 <= k2. These are interpreted relative to the main
diagonal, which is k = 0.
Positive and negative values of k indicate
diagonals above and below the main diagonal, respectively.
...
optional arguments passed to methods, currently unused
by package Matrix.
Details
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).
Value
An object of a suitable matrix class, inheriting from
triangularMatrix where appropriate.
It inherits from sparseMatrix if
and only if x does.
Methods
x = "CsparseMatrix"
method for compressed, sparse,
column-oriented matrices.
x = "RsparseMatrix"
method for compressed, sparse,
row-oriented matrices.
x = "TsparseMatrix"
method for sparse matrices in
triplet format.
x = "diagonalMatrix"
method for diagonal matrices.
x = "denseMatrix"
method for dense matrices in
packed or unpacked format.
x = "matrix"
method for traditional matrices
of implicit class matrix.
See Also
bandSparse for the construction of a
banded sparse matrix directly from its non-zero diagonals.
Examples
## A random sparse matrix :
set.seed(7)
m <- matrix(0,5,5)
m[sample(length(m), size =14)]<- rep(1:9, length=14)(mm <- as(m,"CsparseMatrix"))
tril(mm)# lower triangle
tril(mm,-1)# strict lower triangle
triu(mm,1)# strict upper triangle
band(mm,-1,2)# general band(m5 <- Matrix(rnorm(25), ncol =5))
tril(m5)# lower triangle
tril(m5,-1)# strict lower triangle
triu(m5,1)# strict upper triangle
band(m5,-1,2)# general band(m65 <- Matrix(rnorm(30), ncol =5))# not square
triu(m65)# result not "dtrMatrix" unless square(sm5 <- crossprod(m65))# symmetric
band(sm5,-1,1)# "dsyMatrix": symmetric band preserves symmetry property
as(band(sm5,-1,1),"sparseMatrix")# often preferable(sm <- round(crossprod(triu(mm/2))))# sparse symmetric ("dsC*")
band(sm,-1,1)# remains "dsC", *however*
band(sm,-2,1)# -> "dgC"
Construct a sparse banded matrix by specifying its non-zero sup- and
super-diagonals.
Usage
bandSparse(n, m = n, k, diagonals, symmetric =FALSE,
repr ="C", giveCsparse =(repr =="C"))
Arguments
n, m
the matrix dimension (n,m)=(nrow,ncol).
k
integer vector of “diagonal numbers”, with identical
meaning as in band(*, k), i.e., relative to the main diagonal,
which is k=0.
diagonals
optional list of sub-/super- diagonals; if missing,
the result will be a pattern matrix, i.e., inheriting from
class nMatrix.
diagonals can also be n′×d matrix, where
d <- length(k) and n′>=min(n,m). In that case,
the sub-/super- diagonals are taken from the columns of
diagonals, where only the first several rows will be used
(typically) for off-diagonals.
symmetric
logical; if true the result will be symmetric
(inheriting from class symmetricMatrix) and
only the upper or lower triangle must be specified (via k and
diagonals).
(deprecated, replaced with repr):
logical indicating if the result should be a
CsparseMatrix or a
TsparseMatrix, where the default was TRUE,
and now is determined from repr; very often Csparse matrices are
more efficient subsequently, but not always.
Value
a sparse matrix (of classCsparseMatrix) of dimension n×m
with diagonal “bands” as specified.
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!)
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×k matrices,
the bdiag_m() function in the ‘Examples’ is an order of
magnitude faster.
Author(s)
Martin Maechler, built on a version posted by Berton Gunter to
R-help; earlier versions have been posted by other authors, notably
Scott Chasalow to S-news. Doug Bates's faster implementation builds
on TsparseMatrix objects.
bandSparse constructs a banded sparse matrix from
its non-zero sub-/super - diagonals.
Note that other CRAN R packages have own versions of bdiag()
which return traditional matrices.
Examples
bdiag(matrix(1:4,2), diag(3))## combine "Matrix" class and traditional matrices:
bdiag(Diagonal(2), matrix(1:3,3,4), diag(3:2))
mlist <- list(1,2:3, diag(x=5:3),27, cbind(1,3:6),100:101)
bdiag(mlist)
stopifnot(identical(bdiag(mlist),
bdiag(lapply(mlist, as.matrix))))
ml <- c(as(matrix((1:24)%%11==0,6,4),"nMatrix"),
rep(list(Diagonal(2, x=TRUE)),3))
mln <- c(ml, Diagonal(x =1:3))
stopifnot(is(bdiag(ml),"lsparseMatrix"),
is(bdiag(mln),"dsparseMatrix"))## random (diagonal-)block-triangular matrices:
rblockTri <-function(nb, max.ni, lambda =3){
.bdiag(replicate(nb,{
n <- sample.int(max.ni,1)
tril(Matrix(rpois(n * n, lambda = lambda), n, n))}))}(T4 <- rblockTri(4,10, lambda =1))
image(T1 <- rblockTri(12,20))##' Fast version of Matrix :: .bdiag() -- for the case of *many* (k x k) matrices:##' @param lmat list(<mat1>, <mat2>, ....., <mat_N>) where each mat_j is a k x k 'matrix'##' @return a sparse (N*k x N*k) matrix of class \code{"\linkS4class{dgCMatrix}"}.
bdiag_m <-function(lmat){## Copyright (C) 2016 Martin Maechler, ETH Zurichif(!length(lmat)) return(new("dgCMatrix"))
stopifnot(is.list(lmat), is.matrix(lmat[[1]]),(k <-(d <- dim(lmat[[1]]))[1])== d[2],# k x k
all(vapply(lmat, dim, integer(2))== k))# all of them
N <- length(lmat)if(N * k > .Machine$integer.max)
stop("resulting matrix too large; would be M x M, with M=", N*k)
M <- as.integer(N * k)## result: an M x M matrix
new("dgCMatrix", Dim = c(M,M),## 'i :' maybe there's a faster way (w/o matrix indexing), but elegant?
i = as.vector(matrix(0L:(M-1L), nrow=k)[, rep(seq_len(N), each=k)]),
p = k *0L:M,
x = as.double(unlist(lmat, recursive=FALSE, use.names=FALSE)))}
l12 <- replicate(12, matrix(rpois(16, lambda =6.4),4,4),
simplify=FALSE)
dim(T12 <- bdiag_m(l12))# 48 x 48
T12[1:20,1:20]
For 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).
Value
a pattern matrix, i.e., inheriting from "nMatrix",
or an "ldiMatrix" in case of a diagonal matrix.
Methods
We provide methods for both the “traditional” (R base) matrices
and numeric vectors and conceptually all matrices and
sparseVectors in package Matrix.
signature(x = "ANY", y = "ANY")
signature(x = "ANY", y = "Matrix")
signature(x = "Matrix", y = "ANY")
signature(x = "nMatrix", y = "nMatrix")
signature(x = "nMatrix", y = "nsparseMatrix")
signature(x = "nsparseMatrix", y = "nMatrix")
signature(x = "nsparseMatrix", y = "nsparseMatrix")
signature(x = "sparseVector", y = "sparseVector")
Note
These boolean arithmetic matrix products had been newly
introduced for Matrix 1.2.0 (March 2015). Its implementation
has still not been tested extensively.
Originally, it was left unspecified how non-structural zeros, i.e., 0's
as part of the M@x slot should be treated for numeric
("dMatrix") and logical ("lMatrix")
sparse matrices. We now specify that boolean matrix products should behave as if
applied to drop0(M), i.e., as if dropping such zeros from
the matrix before using it.
Equivalently, for all matrices M, boolean arithmetic should work as if
applied to M != 0 (or M != FALSE).
The current implementation ends up coercing both x and y to
(virtual) class nsparseMatrix which may be quite inefficient
for dense matrices. A future implementation may well return a matrix
with different class, but the “same” content, i.e., the
same matrix entries mij.
See Also
%*%, crossprod(), or tcrossprod(),
for (regular) matrix product methods.
Examples
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"
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.
for [cr]bind, vector- or matrix-like R objects
to be bound together; for [cr]bind2, further arguments
passed to or from methods; see cbind and
cbind2.
deparse.level
integer controlling the construction of labels
in the case of non-matrix-like arguments; see cbind.
x, y
vector- or matrix-like R objects to be bound together.
Value
typically a ‘matrix-like’ object of a similar
class as the first argument in ....
Note that sometimes by default, the result is a
sparseMatrix if one of the arguments is (even in
the case where this is not efficient). In other cases,
the result is chosen to be sparse when there are more zero entries is
than non-zero ones (as the default sparse in
Matrix()).
(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))
Computes the upper triangular Cholesky factor of an
n×n real, symmetric, positive semidefinite
matrix A, optionally after pivoting.
That is the factor L′ in
P1AP1′=LL′
or (equivalently)
A=P1′LL′P1
where
P1 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 P1 is an identity matrix.
Methods for sparseMatrix are built on
CHOLMOD routines cholmod_analyze and cholmod_factorize_p.
Usage
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",...)
Arguments
x
a finite, symmetric, positive
semidefinite matrix or Matrix to
be factorized. If x is square but not symmetric,
then it will be treated as symmetric; see uplo.
Methods for dense x require positive definiteness
when pivot = FALSE.
Methods for sparse (but not diagonal) x require
positive definiteness unconditionally.
pivot
a logical indicating if the rows and columns
of x should be pivoted. Methods for sparse x
employ the approximate minimum degree (AMD) algorithm
in order to reduce fill-in, i.e., without regard for
numerical stability.
tol
a finite numeric tolerance,
used only if pivot = TRUE.
The factorization algorithm stops if the pivot is less than
or equal to tol. Negative tol is equivalent
to nrow(x) * .Machine$double.eps * max(diag(x)).
uplo
a string, either "U" or "L",
indicating which triangle of x should be used
to compute the factorization. The default is "U",
even for lower triangular x, to be consistent with
chol from base.
...
further arguments passed to or from methods.
Details
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 P1 in addition
to the Cholesky factor L′, then call Cholesky
directly, as the result of chol(x, pivot = TRUE) specifies
L′ but not P1.
Value
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.
Chen, Y., Davis, T. A., Hager, W. W., & Rajamanickam, S. (2008).
Algorithm 887: CHOLMOD, supernodal sparse Cholesky factorization
and update/downdate.
ACM Transactions on Mathematical Software,
35(3), Article 22, 1-14.
doi:10.1145/1391989.1391995
Amestoy, P. R., Davis, T. A., & Duff, I. S. (2004).
Algorithm 837: AMD, an approximate minimum degree ordering algorithm.
ACM Transactions on Mathematical Software,
17(4), 886-905.
doi:10.1145/1024074.1024081
Golub, G. H., & Van Loan, C. F. (2013).
Matrix computations (4th ed.).
Johns Hopkins University Press.
doi:10.56021/9781421407944
See Also
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 factorL′.
Examples
showMethods("chol", inherited =FALSE)
set.seed(0)## ---- Dense ----------------------------------------------------------## chol(x, pivot = value) wrapping Cholesky(x, perm = value)
selectMethod("chol","dsyMatrix")## Except in packed cases where pivoting is not yet available
selectMethod("chol","dspMatrix")## .... Positive definite ..............................................(A1 <- new("dsyMatrix", Dim = c(2L,2L), x = c(1,2,2,5)))(R1.nopivot <- chol(A1))(R1 <- chol(A1, pivot =TRUE))## In 2-by-2 cases, we know that the permutation is 1:2 or 2:1,## even if in general 'chol' does not say ...
stopifnot(exprs ={
all.equal( A1 , as(crossprod(R1.nopivot),"dsyMatrix"))
all.equal(t(A1[2:1,2:1]), as(crossprod(R1 ),"dsyMatrix"))
identical(Cholesky(A1)@perm,2:1)# because 5 > 1})## .... Positive semidefinite but not positive definite ................(A2 <- new("dpoMatrix", Dim = c(2L,2L), x = c(1,2,2,4)))
try(R2.nopivot <- chol(A2))# fails as not positive definite(R2 <- chol(A2, pivot =TRUE))# returns, with a warning and ...
stopifnot(exprs ={
all.equal(t(A2[2:1,2:1]), as(crossprod(R2),"dsyMatrix"))
identical(Cholesky(A2)@perm,2:1)# because 4 > 1})## .... Not positive semidefinite ......................................(A3 <- new("dsyMatrix", Dim = c(2L,2L), x = c(1,2,2,3)))
try(R3.nopivot <- chol(A3))# fails as not positive definite(R3 <- chol(A3, pivot =TRUE))# returns, with a warning and ...## _Not_ equal: see details and examples in help("Cholesky")
all.equal(t(A3[2:1,2:1]), as(crossprod(R3),"dsyMatrix"))## ---- Sparse ---------------------------------------------------------## chol(x, pivot = value) wrapping## Cholesky(x, perm = value, LDL = FALSE, super = FALSE)
selectMethod("chol","dsCMatrix")## Except in diagonal cases which are handled "directly"
selectMethod("chol","ddiMatrix")(A4 <- toeplitz(as(c(10,0,1,0,3),"sparseVector")))(ch.A4.nopivot <- Cholesky(A4, perm =FALSE, LDL =FALSE, super =FALSE))(ch.A4 <- Cholesky(A4, perm =TRUE, LDL =FALSE, super =FALSE))(R4.nopivot <- chol(A4))(R4 <- chol(A4, pivot =TRUE))
det4 <- det(A4)
b4 <- rnorm(5L)
x4 <- solve(A4, b4)
stopifnot(exprs ={
identical(R4.nopivot, expand1(ch.A4.nopivot,"L."))
identical(R4, expand1(ch.A4,"L."))
all.equal(A4, crossprod(R4.nopivot))
all.equal(A4[ch.A4@perm +1L, ch.A4@perm +1L], crossprod(R4))
all.equal(diag(R4.nopivot), sqrt(diag(ch.A4.nopivot)))
all.equal(diag(R4), sqrt(diag(ch.A4)))
all.equal(sqrt(det4), det(R4.nopivot))
all.equal(sqrt(det4), det(R4))
all.equal(det4, det(ch.A4.nopivot, sqrt =FALSE))
all.equal(det4, det(ch.A4, sqrt =FALSE))
all.equal(x4, solve(R4.nopivot, solve(t(R4.nopivot), b4)))
all.equal(x4, solve(ch.A4.nopivot, b4))
all.equal(x4, solve(ch.A4, b4))})
Given formally upper and lower triangular matrices
U and L, compute (U′U)−1
and (LL′)−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, if R is constrained to have
positive diagonal entries.
Usage
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",...)
Arguments
x
a square matrix or Matrix,
typically the result of a call to chol.
If x is square but not (formally) triangular,
then only the upper or lower triangle is considered,
depending on optional argument uplo if x
is a Matrix.
uplo
a string, either "U" or "L",
indicating which triangle of x contains the
Cholesky factor. The default is "U", to be
consistent with chol2inv from base.
...
further arguments passed to or from methods.
Value
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.
See Also
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.
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,
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.
for graph2T(), an R object of class
"graph";
for T2graph(), a sparse matrix inheriting from
"TsparseMatrix".
use.weights
logical indicating if weights should be used, i.e.,
equivalently the result will be numeric, i.e. of class
dgTMatrix; otherwise the result will be
ngTMatrix or nsTMatrix,
the latter if the graph is undirected. The default looks if there
are weights in the graph, and if any differ from 1, weights
are used.
need.uniq
a logical indicating if from may need to be
internally “uniqified”; do not set this and hence rather use
the default, unless you know what you are doing!
edgemode
one of NULL, "directed", or
"undirected". The default NULL looks if the matrix is
symmetric and assumes "undirected" in that case.
Value
For graph2T(), a sparse matrix inheriting from
"TsparseMatrix".
For T2graph() an R object of class "graph".
See Also
Package igraph, which provides similar coercions
to and from its class igraph via functions
graph_from_adjacency_matrix and as_adjacency_matrix.
Examples
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"))}
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.
logical. Should missing values (including NaN)
be omitted from the calculations?
dims
completely ignored by the Matrix methods.
...
potentially further arguments, for method <->
generic compatibility.
sparseResult
logical indicating if the result should be sparse,
i.e., inheriting from class sparseVector. Only
applicable when x is inheriting from a
sparseMatrix class.
Value
returns a numeric vector if sparseResult is FALSE as per
default. Otherwise, returns a sparseVector.
(M <- bdiag(Diagonal(2), matrix(1:3,3,4), diag(3:2)))# 7 x 8
colSums(M)
d <- Diagonal(10, c(0,0,10,0,2,rep(0,5)))
MM <- kronecker(d, M)
dim(MM)# 70 80
length(MM@x)# 160, but many are '0' ; drop those:
MM <- drop0(MM)
length(MM@x)# 32
cm <- colSums(MM)(scm <- colSums(MM, sparseResult =TRUE))
stopifnot(is(scm,"sparseVector"),
identical(cm, as.numeric(scm)))
rowSums (MM, sparseResult =TRUE)# 14 of 70 are not zero
colMeans(MM, sparseResult =TRUE)# 16 of 80 are not zero## Since we have no 'NA's, these two are equivalent :
stopifnot(identical(rowMeans(MM, sparseResult =TRUE),
rowMeans(MM, sparseResult =TRUE, na.rm =TRUE)),
rowMeans(Diagonal(16))==1/16,
colSums(Diagonal(7))==1)## dimnames(x) --> names( <value> ) :
dimnames(M)<- list(paste0("r",1:7), paste0("V",1:8))
M
colSums(M)
rowMeans(M)## Assertions :
stopifnot(exprs ={
all.equal(colSums(M),
structure(c(1,1,6,6,6,6,3,2), names = colnames(M)))
all.equal(rowMeans(M),
structure(c(1,1,4,8,12,3,2)/8, names = paste0("r",1:7)))})
“Estimate”, i.e. compute approximately the CONDition number of
a (potentially large, often sparse) matrix A.
It works by apply a fast randomized approximation of the 1-norm,
norm(A,"1"), through onenormest(.).
Usage
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)
Arguments
A
a square matrix, optional for onenormest(), where
instead of A, A.x and At.x can be specified,
see there.
t
number of columns to use in the iterations.
normA
number; (an estimate of) the 1-norm of A, by
default norm(A, "1"); may be replaced by an estimate.
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 A is missing, these two must be given as
functions which compute A %% x, or t(A) %% x,
respectively.
n
== nrow(A), only needed when A is not specified.
iter.max
maximal number of iterations for the 1-norm estimator.
eps
the relative change that is deemed irrelevant.
Details
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().
Value
Both functions return a list;
condest() with components,
est
a number >0, the estimated (1-norm) condition number
κ^; when r:=rcond(A),
1/κ^≈r.
v
the maximal Ax column, scaled to norm(v) = 1.
Consequently, norm(Av)=norm(A)/est;
when est is large, v is an approximate null vector.
The function onenormest() returns a list with components,
est
a number >0, the estimated norm(A, "1").
v
0-1 integer vector length n, with an 1 at the index
j with maximal column A[,j] in A.
w
numeric vector, the largest Ax found.
iter
the number of iterations used.
Author(s)
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.
References
Nicholas J. Higham and Françoise Tisseur (2000).
A Block Algorithm for Matrix 1-Norm Estimation, with an Application
to 1-Norm Pseudospectra.
SIAM J. Matrix Anal. Appl.21, 4, 1185–1201.
William W. Hager (1984).
Condition Estimates.
SIAM J. Sci. Stat. Comput.5, 311–316.
data(KNex, package ="Matrix")
mtm <- with(KNex, crossprod(mm))
system.time(ce <- condest(mtm))
sum(abs(ce$v))## || v ||_1 == 1## Prove that || A v || = || A || / est (as ||v|| = 1):
stopifnot(all.equal(norm(mtm %*% ce$v),
norm(mtm)/ ce$est))## reciprocal1/ ce$est
system.time(rc <- rcond(mtm))# takes ca 3 x longer
rc
all.equal(rc,1/ce$est)# TRUE -- the approximation was good
one <- onenormest(mtm)
str(one)## est = 12.3## the maximal column:
which(one$v ==1)# mostly 4, rarely 1, depending on random seed
The 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.
Slots
Common to all matrix object in the package:
Dim:
Object of class "integer" - the dimensions
of the matrix - must be an integer vector with exactly two
non-negative values.
Dimnames:
list of length two; each component
containing NULL or a character vector length
equal the corresponding Dim element.
Methods
There are (relatively simple) group methods (see, e.g., Arith)
Arith
signature(e1 = "dMatrix", e2 = "dMatrix"): ...
Arith
signature(e1 = "dMatrix", e2 = "numeric"): ...
Arith
signature(e1 = "numeric", e2 = "dMatrix"): ...
Math
signature(x = "dMatrix"): ...
Math2
signature(x = "dMatrix", digits = "numeric"):
this group contains round() and signif().
Compare
signature(e1 = "numeric", e2 = "dMatrix"): ...
Compare
signature(e1 = "dMatrix", e2 = "numeric"): ...
Compare
signature(e1 = "dMatrix", e2 = "dMatrix"): ...
Summary
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:
expm
signature(x = "dMatrix"): computes the
“Matrix Exponential”, see expm.
The following methods are defined for all logical matrices:
which
signature(x = "lsparseMatrix") and many other
subclasses of "lMatrix": as the base function
which(x, arr.ind) returns the indices of the
TRUE entries in x; if arr.ind is true,
as a 2-column matrix of row and column indices. Since Matrix
version 1.2-9, if useNames is true, as by default, with
dimnames, the same as base::which.
See Also
The nonzero-pattern matrix class nMatrix, which
can be used to store non-NAlogical
matrices even more compactly.
The class "ddiMatrix" of numerical diagonal matrices.
Note that diagonal matrices now extend sparseMatrix, whereas
they did extend dense matrices earlier.
Objects from the Class
Objects can be created by calls of the form new("ddiMatrix", ...)
but typically rather via Diagonal.
Slots
x:
numeric vector. For an n×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 unit-diagonal, i.e., identity
matrices.
Dim,Dimnames:
matrix dimension and
dimnames, see the Matrix class
description.
denseLU is the class of dense, row-pivoted LU factorizations
of m×n real matrices A,
having the general form
P1A=LU
or (equivalently)
A=P1′LU
where
P1 is an m×m permutation matrix,
L is an m×min(m,n)
unit lower trapezoidal matrix, and
U is a min(m,n)×n
upper trapezoidal matrix. If m=n, then the factors
L and U are triangular.
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 P1 as a product of
transpositions. The corresponding permutation vector can
be obtained as asPerm(perm).
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).
Methods
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.
The dgCMatrix class is a class of sparse numeric
matrices in the compressed, sparse, column-oriented format. In this
implementation the non-zero elements in the columns are sorted into
increasing row order. dgCMatrix is the
“standard” class for sparse numeric matrices in the
Matrix package.
Objects from the Class
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().
Slots
x:
Object of class "numeric" - the non-zero
elements of the matrix.
...
all other slots are inherited from the superclass
"CsparseMatrix".
The dgRMatrix class is a class of sparse numeric
matrices in the compressed, sparse, row-oriented format. In this
implementation the non-zero elements in the rows are sorted into
increasing column order.
Note: The column-oriented sparse classes, e.g.,
dgCMatrix, are preferred and better supported in
the Matrix package.
Objects from the Class
Objects can be created by calls of the form new("dgRMatrix", ...).
Slots
j:
Object of class "integer" of length nnzero
(number of non-zero elements). These are the column numbers for
each non-zero element in the matrix.
p:
Object of class "integer" of pointers, one
for each row, to the initial (zero-based) index of elements in
the row.
x:
Object of class "numeric" - the non-zero
elements of the matrix.
Dim:
Object of class "integer" - the dimensions
of the matrix.
Methods
diag
signature(x = "dgRMatrix"): returns the diagonal
of x
dim
signature(x = "dgRMatrix"): returns the dimensions
of x
image
signature(x = "dgRMatrix"): plots an image of
x using the levelplot function
See Also
the RsparseMatrix class, the virtual class of all
sparse compressed row-oriented matrices, with its methods.
The dgCMatrix class (column compressed
sparse) is really preferred.
The "dgTMatrix" class is the class of sparse
matrices stored as (possibly redundant) triplets. The internal
representation is not at all unique, contrary to the one for class
dgCMatrix.
Objects from the Class
Objects can be created by calls of the form
new("dgTMatrix", ...), but more typically via
spMatrix() or sparseMatrix(*, repr = "T").
Slots
i:
integer row indices of non-zero
entries in 0-base, i.e., must be in 0:(nrow(.)-1).
j:
integer column indices of non-zero
entries. Must be the same length as slot i and
0-based as well, i.e., in 0:(ncol(.)-1).
x:
numeric vector - the (non-zero)
entry at position (i,j). Must be the same length as slot
i. If an index pair occurs more than once, the corresponding
values of slot x are added to form the element of the matrix.
Dim:
Object of class "integer" of length 2 -
the dimensions of the matrix.
Methods
+
signature(e1 = "dgTMatrix", e2 = "dgTMatrix")
image
signature(x = "dgTMatrix"): plots an image of
x using the levelplot function
t
signature(x = "dgTMatrix"): returns the transpose of
x
Note
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 xk's that belong to identical
(ik,jk) pairs.
However this means that a matrix typically can be stored in more than
one possible "TsparseMatrix" representations.
Use asUniqueT() in order to ensure uniqueness of the
internal representation of such a matrix.
m <- Matrix(0+1:28, nrow =4)
m[-3,c(2,4:5,7)]<- m[3,1:4]<- m[1:3,6]<-0(mT <- as(m,"TsparseMatrix"))
str(mT)
mT[1,]
mT[4, drop =FALSE]
stopifnot(identical(mT[lower.tri(mT)],
m [lower.tri(m)]))
mT[lower.tri(mT,diag=TRUE)]<-0
mT
## Triplet representation with repeated (i,j) entries## *adds* the corresponding x's:
T2 <- new("dgTMatrix",
i = as.integer(c(1,1,0,3,3)),
j = as.integer(c(2,2,4,0,0)), x=10*1:5, Dim=4:5)
str(T2)# contains (i,j,x) slots exactly as above, but
T2 ## has only three non-zero entries, as for repeated (i,j)'s,## the corresponding x's are "implicitly" added
stopifnot(nnzero(T2)==3)
Transform a triangular matrix x, i.e., of classtriangularMatrix,
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.
(optional, for speedup only:) class (definition) of x.
checkDense
logical indicating if dense (see
denseMatrix) matrices should be considered at
all; i.e., when false, as per default, the result will be sparse even
when x is dense.
Details
The concept of unit triangular matrices with a diag slot of
"U" stems from LAPACK.
Value
a triangular matrix of the same class but with a
different diag slot. For diagU2N (semantically) with
identical entries as x, whereas in diagN2U(x), the
off-diagonal entries are unchanged and the diagonal is set to all
1 even if it was not previously.
Note
Such internal storage details should rarely be of relevance to the
user. Hence, these functions really are rather internal
utilities.
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.
coerce
signature(from = "matrix", to = "diagonalMatrix"): ...
coerce
signature(from = "Matrix", to = "diagonalMatrix"): ...
coerce
signature(from = "diagonalMatrix", to = "generalMatrix"): ...
coerce
signature(from = "diagonalMatrix", to = "triangularMatrix"): ...
coerce
signature(from = "diagonalMatrix", to = "nMatrix"): ...
coerce
signature(from = "diagonalMatrix", to = "matrix"): ...
coerce
signature(from = "diagonalMatrix", to = "sparseVector"): ...
t
signature(x = "diagonalMatrix"): ...
and many more methods
solve
signature(a = "diagonalMatrix", b, ...): is
trivially implemented, of course; see also solve-methods.
which
signature(x = "nMatrix"), semantically
equivalent to base function which(x, arr.ind).
"Math"
signature(x = "diagonalMatrix"): all these
group methods return a "diagonalMatrix", apart from
cumsum() etc which return a vector also for
basematrix.
*
signature(e1 = "ddiMatrix", e2="denseMatrix"):
arithmetic and other operators from the Ops
group have a few dozen explicit method definitions, in order to
keep the results diagonal in many cases, including the following:
/
signature(e1 = "ddiMatrix", e2="denseMatrix"):
the result is from class ddiMatrix which is
typically very desirable. Note that when e2 contains
off-diagonal zeros or NAs, we implicitly use 0/x=0, hence
differing from traditional R arithmetic (where 0/0↦NaN), in order to preserve sparsity.
summary
(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.
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.
a matrix, possibly inheriting from virtual class
Matrix.
d1, d2, d
numeric vectors giving factors by which to scale
the rows or columns of x; they are recycled as necessary.
Details
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.
Value
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.
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)
For any n×m (typically) sparse matrix x
compute the Dulmage-Mendelsohn 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.
Usage
dmperm(x, nAns =6L, seed =0L)
Arguments
x
a typically sparse matrix; internally coerced to either
"dgCMatrix" or
"dtCMatrix".
nAns
an integer specifying the length of the
resulting list. Must be 2, 4, or 6.
seed
an integer code in -1,0,1; determining the (initial)
permutation; by default, seed = 0, no (or the identity) permutation;
seed = -1 uses the “reverse” permutation k:1; for
seed = 1, it is a random permutation (using R's RNG,
seed, etc).
Details
See the book section by Tim Davis; page 122–127, in the References.
integer vector with the permutation p, of length nrow(x).
q
integer vector with the permutation q, of length ncol(x).
r
integer vector of length nb+1, where block k is rows r[k] to r[k+1]-1 in A[p,q].
s
integer vector of length nb+1, where block k is cols s[k] to s[k+1]-1 in A[p,q].
rr5
integer vector of length 5, defining the coarse row
decomposition.
cc5
integer vector of length 5, defining the coarse column decomposition.
Author(s)
Martin Maechler, with a lot of “encouragement” by Mauricio
Vargas.
References
Section 7.4 Dulmage-Mendelsohn decomposition, pp. 122 ff of
Timothy A. Davis (2006)
Direct Methods for Sparse Linear Systems, SIAM Series
“Fundamentals of Algorithms”.
See Also
Schur, the class of permutation matrices; "pMatrix".
Examples
set.seed(17)(S9 <- rsparsematrix(9,9, nnz =10, symmetric=TRUE))# dsCMatrix
str( dm9 <- dmperm(S9))(S9p <- with(dm9, S9[p, q]))## looks good, but *not* quite upper triangular; these, too:
str( dm9.0<- dmperm(S9, seed=-1))# non-random too.
str( dm9_1 <- dmperm(S9, seed=1))# a random one## The last two permutations differ, but have the same effect!(S9p0 <- with(dm9.0, S9[p, q]))# .. hmm ..
stopifnot(all.equal(S9p0, S9p))# same as as default, but different from the random one
set.seed(11)(M <- triu(rsparsematrix(9,11,1/4)))
dM <- dmperm(M); with(dM, M[p, q])(Mp <- M[sample.int(nrow(M)), sample.int(ncol(M))])
dMp <- dmperm(Mp); with(dMp, Mp[p, q])
set.seed(7)(n7 <- rsparsematrix(5,12, nnz =10, rand.x =NULL))
str( dm.7<- dmperm(n7))
stopifnot(exprs ={
lengths(dm.7[1:2])== dim(n7)
identical(dm.7, dmperm(as(n7,"dMatrix")))
identical(dm.7[1:4], dmperm(n7, nAns=4))
identical(dm.7[1:2], dmperm(n7, nAns=2))})
The "dpoMatrix" class is the class of
positive-semidefinite symmetric matrices in nonpacked storage.
The "dppMatrix" class is the same except in packed
storage. Only the upper triangle or the lower triangle is
required to be available.
The "corMatrix" and "copMatrix" classes
represent correlation matrices. They extend "dpoMatrix"
and "dppMatrix", respectively, with an additional slot
sd allowing restoration of the original covariance matrix.
Objects from the Class
Objects can be created by calls of the
form new("dpoMatrix", ...) or from crossprod applied to
an "dgeMatrix" object.
Slots
uplo:
Object of class "character". Must be
either "U", for upper triangular, and "L", for lower triangular.
x:
Object of class "numeric". The numeric
values that constitute the matrix, stored in column-major order.
Dim:
Object of class "integer". The dimensions
of the matrix which must be a two-element vector of non-negative
integers.
Dimnames:
inherited from class "Matrix"
factors:
Object of class "list". A named
list of factorizations that have been computed for the matrix.
sd:
(for "corMatrix" and "copMatrix")
a numeric vector of length n containing the
(original) var(.) entries which allow
reconstruction of a covariance matrix from the correlation matrix.
Extends
Class "dsyMatrix", directly.
Classes "dgeMatrix", "symmetricMatrix", and many more
by class "dsyMatrix".
Methods
chol
signature(x = "dpoMatrix"):
Returns (and stores) the Cholesky decomposition of x, see
chol.
determinant
signature(x = "dpoMatrix"):
Returns the determinant of x, via
chol(x), see above.
rcond
signature(x = "dpoMatrix", norm = "character"):
Returns (and stores) the reciprocal of the condition number of
x. The norm can be "O" for the
one-norm (the default) or "I" for the infinity-norm. For
symmetric matrices the result does not depend on the norm.
solve
signature(a = "dpoMatrix", b = "...."), and
solve
signature(a = "dppMatrix", b = "....") work
via the Cholesky composition, see also the Matrix solve-methods.
Arith
signature(e1 = "dpoMatrix", e2 = "numeric") (and
quite a few other signatures): The result of (“elementwise”
defined) arithmetic operations is typically not
positive-definite anymore. The only exceptions, currently, are
multiplications, divisions or additions with positivelength(.) == 1 numbers (or logicals).
Note
Currently the validity methods for these classes such as
getValidity(getClass("dpoMatrix")) for efficiency reasons
only check the diagonal entries of the matrix – they may not be negative.
This is only necessary but not sufficient for a symmetric matrix to be
positive semi-definite.
A more reliable (but often more expensive) check for positive
semi-definiteness would look at the signs of diag(BunchKaufman(.))
(with some tolerance for very small negative values), and for (strict)
positive definiteness at something like
!inherits(tryCatch(chol(.), error=identity), "error") .
Indeed, when coercing to these classes, a version
of Cholesky() or chol() is
typically used, e.g., see selectMethod("coerce",
c(from="dsyMatrix", to="dpoMatrix")) .
Deletes “non-structural” zeros (i.e., zeros stored explicitly,
in memory) from a sparse matrix and returns the result.
Usage
drop0(x, tol =0, is.Csparse =NA, give.Csparse =TRUE)
Arguments
x
a Matrix, typically inheriting
from virtual class sparseMatrix.
denseMatrix and traditional vectors and
matrices are coerced to CsparseMatrix,
with zeros dropped automatically, hence users passing such
x should consider as(x, "CsparseMatrix") instead,
notably in the tol = 0 case.
tol
a non-negative number. If x is numeric,
then entries less than or equal to tol in absolute value
are deleted.
is.Csparse
a logical used only if give.Csparse
is TRUE, indicating if x already inherits from
virtual class CsparseMatrix, in which
case coercion is not attempted, permitting some (typically
small) speed-up.
give.Csparse
a logical indicating if the result must
inherit from virtual class CsparseMatrix.
If FALSE and x inherits from
RsparseMatrix,
TsparseMatrix, or
indMatrix,
then the result preserves the class of x.
The default value is TRUE only for backwards compatibility.
Value
A sparseMatrix, the result of deleting
non-structural zeros from x, possibly after coercion.
Note
drop0 is sometimes called in conjunction with
zapsmall, e.g., when dealing with sparse
matrix products; see the example.
(m <- sparseMatrix(i =1:8, j =2:9, x = c(0:2,3:-1),
dims = c(10L,20L)))
drop0(m)## A larger example:
t5 <- new("dtCMatrix", Dim = c(5L,5L), uplo ="L",
x = c(10,1,3,10,1,10,1,10,10),
i = c(0L,2L,4L,1L,3L,2L,4L,3L,4L),
p = c(0L,3L,5L,7:9))
TT <- kronecker(t5, kronecker(kronecker(t5, t5), t5))
IT <- solve(TT)
I. <- TT %*% IT ; nnzero(I.)# 697 ( == 625 + 72 )
I.0<- drop0(zapsmall(I.))## which actually can be more efficiently achieved by
I.. <- drop0(I., tol =1e-15)
stopifnot(all(I.0== Diagonal(625)), nnzero(I..)==625)
The dsCMatrix class is a class of symmetric, sparse
numeric matrices in the compressed, column-oriented format. In
this implementation the non-zero elements in the columns are sorted
into increasing row order.
The dsTMatrix class is the class of symmetric, sparse numeric
matrices in triplet format.
Objects from the Class
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).
Slots
uplo:
A character object indicating if the upper
triangle ("U") or the lower triangle ("L") is stored.
i:
Object of class "integer" of length nnZ
(half number of non-zero elements). These are the row
numbers for each non-zero element in the lower triangle of the matrix.
p:
(only in class "dsCMatrix":) an
integer vector for providing pointers, one for each
column, see the detailed description in CsparseMatrix.
j:
(only in class "dsTMatrix":) Object of
class "integer" of length nnZ (as i). These are the
column numbers for each non-zero element in the lower triangle of
the matrix.
x:
Object of class "numeric" of length nnZ –
the non-zero elements of the matrix (to be duplicated for full matrix).
factors:
Object of class "list" - a list
of factorizations of the matrix.
Dim:
Object of class "integer" - the dimensions
of the matrix - must be an integer vector with exactly two
non-negative values.
signature(a = "dsCMatrix", b = "...."): x
<- solve(a,b) solves Ax=b for x; see
solve-methods.
chol
signature(x = "dsCMatrix", pivot = "logical"):
Returns (and stores) the Cholesky decomposition of x, see
chol.
Cholesky
signature(A = "dsCMatrix",...):
Computes more flexibly Cholesky decompositions,
see Cholesky.
determinant
signature(x = "dsCMatrix", logarithm =
"missing"): Evaluate the determinant of x on the
logarithm scale. This creates and stores the Cholesky factorization.
determinant
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.
t
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.
t
signature(x = "dsTMatrix"): Transpose. The
uplo slot is swapped from "U" to "L" or vice
versa, as for a "dsCMatrix", see 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 dsRMatrix class is a class of symmetric, sparse
matrices in the compressed, row-oriented format. In this
implementation the non-zero elements in the rows are sorted into
increasing column order.
Objects from the Class
These "..RMatrix" classes are currently still mostly unimplemented!
Objects can be created by calls of the form new("dsRMatrix", ...).
Slots
uplo:
A character object indicating if the upper
triangle ("U") or the lower triangle ("L") is
stored. At present only the lower triangle form is allowed.
j:
Object of class "integer" of length
nnzero (number of non-zero elements). These are the row
numbers for each non-zero element in the matrix.
p:
Object of class "integer" of pointers, one
for each row, to the initial (zero-based) index of elements in
the row.
factors:
Object of class "list" - a list
of factorizations of the matrix.
x:
Object of class "numeric" - the non-zero
elements of the matrix.
Dim:
Object of class "integer" - the dimensions
of the matrix - must be an integer vector with exactly two
non-negative values.
The "dsyMatrix" class is the class of symmetric, dense matrices
in non-packed storage and
"dspMatrix" is the class of symmetric dense matrices in
packed storage, see pack(). Only the upper
triangle or the lower triangle is stored.
Objects from the Class
Objects can be created by calls of the form new("dsyMatrix",
...) or new("dspMatrix", ...), respectively.
Slots
uplo:
Object of class "character". Must be
either "U", for upper triangular, and "L", for lower triangular.
x:
Object of class "numeric". The numeric
values that constitute the matrix, stored in column-major order.
Dim,Dimnames:
The dimension (a length-2
"integer") and corresponding names (or NULL), see the
Matrix.
factors:
Object of class "list". A named
list of factorizations that have been computed for the matrix.
Extends
"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.
Methods
norm
signature(x = "dspMatrix", type = "character"), or
x = "dsyMatrix" or type = "missing": Computes the
matrix norm of the desired type, see, norm.
rcond
signature(x = "dspMatrix", type = "character"), or
x = "dsyMatrix" or type = "missing": Computes the
reciprocal condition number, rcond().
solve
signature(a = "dspMatrix", b = "...."), and
solve
signature(a = "dsyMatrix", b = "...."): x
<- solve(a,b) solves Ax=b for x; see
solve-methods.
t
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.
See Also
The positive (Semi-)definite dense (packed or non-packed
numeric matrix classes dpoMatrix,
dppMatrix and corMatrix,
## Only upper triangular part matters (when uplo == "U" as per default)(sy2 <- new("dsyMatrix", Dim = as.integer(c(2,2)), x = c(14,NA,32,77)))
str(t(sy2))# uplo = "L", and the lower tri. (i.e. NA is replaced).
chol(sy2)#-> "Cholesky" matrix(sp2 <- pack(sy2))# a "dspMatrix"## Coercing to dpoMatrix gives invalid object:
sy3 <- new("dsyMatrix", Dim = as.integer(c(2,2)), x = c(14,-1,2,-7))
try(as(sy3,"dpoMatrix"))# -> error: not positive definite## 4x4 example
m <- matrix(0,4,4); m[upper.tri(m)]<-1:6(sym <- m+t(m)+diag(11:14,4))(S1 <- pack(sym))(S2 <- t(S1))
stopifnot(all(S1 == S2))# equal "seen as matrix", but differ internally :
str(S1)
S2@x
The "dtCMatrix" class is a class of triangular, sparse
matrices in the compressed, column-oriented format. In this
implementation the non-zero elements in the columns are sorted into
increasing row order.
The "dtTMatrix" class is a class of triangular, sparse matrices
in triplet format.
Objects from the Class
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").
Slots
uplo:
Object of class "character". Must be
either "U", for upper triangular, and "L", for lower triangular.
diag:
Object of class "character". Must be
either "U", for unit triangular (diagonal is all ones), or
"N"; see triangularMatrix.
p:
(only present in "dtCMatrix":) an
integer vector for providing pointers, one for each
column, see the detailed description in CsparseMatrix.
i:
Object of class "integer" of length nnzero
(number of non-zero elements). These are the row numbers for
each non-zero element in the matrix.
j:
Object of class "integer" of length nnzero
(number of non-zero elements). These are the column numbers for
each non-zero element in the matrix. (Only present in the
dtTMatrix class.)
x:
Object of class "numeric" - the non-zero
elements of the matrix.
Dim,Dimnames:
The dimension (a length-2
"integer") and corresponding names (or NULL),
inherited from the Matrix, see there.
Extends
Class "dgCMatrix", directly.
Class "triangularMatrix", directly.
Class "dMatrix", "sparseMatrix", and more by class
"dgCMatrix" etc, see the examples.
Methods
solve
signature(a = "dtCMatrix", b = "...."):
sparse triangular solve (aka “backsolve” or
“forwardsolve”), see solve-methods.
t
signature(x = "dtCMatrix"): returns the transpose of
x
t
signature(x = "dtTMatrix"): returns the transpose of
x
The dtRMatrix class is a class of triangular, sparse
matrices in the compressed, row-oriented format. In this
implementation the non-zero elements in the rows are sorted into
increasing columnd order.
Objects from the Class
This class is currently still mostly unimplemented!
Objects can be created by calls of the form new("dtRMatrix", ...).
Slots
uplo:
Object of class "character". Must be
either "U", for upper triangular, and "L", for lower triangular.
At present only the lower triangle form is allowed.
diag:
Object of class "character". Must be
either "U", for unit triangular (diagonal is all ones), or
"N"; see triangularMatrix.
j:
Object of class "integer" of length
nnzero(.) (number of non-zero elements). These are
the row numbers for each non-zero element in the matrix.
p:
Object of class "integer" of pointers, one
for each row, to the initial (zero-based) index of elements in
the row. (Only present in the dsRMatrix class.)
x:
Object of class "numeric" - the non-zero
elements of the matrix.
Dim:
The dimension (a length-2 "integer")
Dimnames:
corresponding names (or NULL),
inherited from the Matrix, see there.
Extends
Class "dgRMatrix", directly.
Class "dsparseMatrix", by class "dgRMatrix".
Class "dMatrix", by class "dgRMatrix".
Class "sparseMatrix", by class "dgRMatrix".
Class "Matrix", by class "dgRMatrix".
Methods
No methods currently with class "dsRMatrix" in the signature.
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 from the Class
Objects can be created by calls of the form new("dtpMatrix",
...) or by coercion from other classes of matrices.
Slots
uplo:
Object of class "character". Must be
either "U", for upper triangular, and "L", for lower triangular.
diag:
Object of class "character". Must be
either "U", for unit triangular (diagonal is all ones), or
"N"; see triangularMatrix.
x:
Object of class "numeric". The numeric
values that constitute the matrix, stored in column-major order.
For a packed square matrix of dimension d×d,
length(x) is of length d(d+1)/2 (also when
diag == "U"!).
Dim,Dimnames:
The dimension (a length-2
"integer") and corresponding names (or NULL),
inherited from the Matrix, see there.
Extends
Class "ddenseMatrix", directly.
Class "triangularMatrix", directly.
Class "dMatrix" and more by class "ddenseMatrix" etc, see
the examples.
Methods
%*%
signature(x = "dtpMatrix", y = "dgeMatrix"):
Matrix multiplication; ditto for several other signature
combinations, see showMethods("%*%", class = "dtpMatrix").
determinant
signature(x = "dtpMatrix", logarithm = "logical"):
the determinant(x) trivially is
prod(diag(x)), but computed on log scale to prevent over-
and underflow.
diag
signature(x = "dtpMatrix"): ...
norm
signature(x = "dtpMatrix", type = "character"): ...
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 from the Class
Objects can be created by calls of the form new("dtrMatrix", ...).
Slots
uplo:
Object of class "character". Must be
either "U", for upper triangular, and "L", for lower triangular.
diag:
Object of class "character". Must be
either "U", for unit triangular (diagonal is all ones), or
"N"; see triangularMatrix.
x:
Object of class "numeric". The numeric
values that constitute the matrix, stored in column-major order.
Dim:
Object of class "integer". The dimensions
of the matrix which must be a two-element vector of non-negative
integers.
Extends
Class "ddenseMatrix", directly.
Class "triangularMatrix", directly.
Class "Matrix" and others, by class "ddenseMatrix".
signature(a = "dtrMatrix", b = "...."): efficiently
use a “forwardsolve” or backsolve for a lower or
upper triangular matrix, respectively, see also
solve-methods.
+, -, *, ..., ==, >=, ...
all the Ops group
methods are available. When applied to two triangular matrices,
these return a triangular matrix when easily possible.
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.
a matrix factorization, typically inheriting from
virtual class MatrixFactorization.
which
a character string indicating a matrix factor.
...
further arguments passed to or from methods.
Details
Methods for expand are retained only for backwards
compatibility with Matrix< 1.6-0. New code
should use expand1 and expand2, whose methods
provide more control and behave more consistently. Notably,
expand2 obeys the rule that the product of the matrix
factors in the returned list should reproduce
(within some tolerance) the factorized matrix,
including its dimnames.
Hence if x is a matrix and y is its factorization,
then
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.
Methods
The following table lists methods for expand1 together with
allowed values of argument which.
Methods for expand2 and expand are described
below. Factor names and classes apply also to expand1.
expand2
signature(x = "CHMsimpl"):
expands the factorization
A=P1′L1DL1′P1=P1′LL′P1
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=L1DL1′=LL′
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=UDUU′=LDLL′
where
U=∏k=1bUPkUk
and
L=∏k=1bLPkLk
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 2bU+1 or 2bL+1 matrix
factors is returned instead.
Pk are represented as pMatrix,
Uk and Lk are represented as
dtCMatrix, and
DU and DL are represented as
dsCMatrix.
expand2
signature(x = "Schur"):
expands the factorization
A=QTQ′
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=P1′LUP2′
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=P1′LU
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=P1′QRP2′=P1′Q1R1P2′
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 P1
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 P2′
but have opposite margin slots and inverted
perm slots. expand(x)[["P"]] represents
the permutation matrix P1 rather than its
transpose P1′; 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.
showMethods("expand1", inherited =FALSE)
showMethods("expand2", inherited =FALSE)
set.seed(0)(A <- Matrix(rnorm(9L,0,10),3L,3L))(lu.A <- lu(A))(e.lu.A <- expand2(lu.A))
stopifnot(exprs ={
is.list(e.lu.A)
identical(names(e.lu.A), c("P1.","L","U"))
all(sapply(e.lu.A, is,"Matrix"))
all.equal(as(A,"matrix"), as(Reduce(`%*%`, e.lu.A),"matrix"))})## 'expand1' and 'expand2' give equivalent results modulo## dimnames and representation of permutation matrices;## see also function 'alt' in example("Cholesky-methods")(a1 <- sapply(names(e.lu.A), expand1, x = lu.A, simplify =FALSE))
all.equal(a1, e.lu.A)## see help("denseLU-class") and others for more examples
a matrix, typically inheriting from the
dMatrix class.
Details
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...”.
Value
The matrix exponential of x.
Author(s)
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.
Cleve Moler and Charles Van Loan (2003)
Nineteen dubious ways to compute the exponential of a matrix,
twenty-five years later. SIAM Review45, 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 Review20, 4, 801–836.
doi:10.1137/1020098
Package expm, which provides newer (in some cases
faster, more accurate) algorithms for computing the matrix
exponential via its own (non-generic) function expm().
expm also implements logm(), sqrtm(), etc.
Read matrices stored in the Harwell-Boeing or MatrixMarket formats
or write sparseMatrix objects to one of these
formats.
Usage
readHB(file)
readMM(file)
writeMM(obj, file,...)
Arguments
obj
a real sparse matrix
file
for writeMM - the name of the file to be written.
For readHB and readMM the name of the file to read, as
a character scalar. The names of files storing matrices in the
Harwell-Boeing format usually end in ".rua" or ".rsa".
Those storing matrices in the MatrixMarket format usually end in
".mtx".
Alternatively, readHB and readMM accept connection objects.
...
optional additional arguments. Currently none are used in
any methods.
Value
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).
Note
The Harwell-Boeing format is older and less flexible than the
MatrixMarket format. The function writeHB was deprecated and
has now been removed. Please use writeMM instead.
Note that these formats do not know anything about
dimnames, hence these are dropped by writeMM().
A very simple way to export small sparse matrices S, is to use
summary(S) which returns a data.frame with
columns i, j, and possibly x, see summary in
sparseMatrix-class, and an example below.
str(pores <- readMM(system.file("external/pores_1.mtx", package ="Matrix")))
str(utm <- readHB(system.file("external/utm300.rua", package ="Matrix")))
str(lundA <- readMM(system.file("external/lund_a.mtx", package ="Matrix")))
str(lundA <- readHB(system.file("external/lund_a.rsa", package ="Matrix")))## https://math.nist.gov/MatrixMarket/data/Harwell-Boeing/counterx/counterx.htm
str(jgl <- readMM(system.file("external/jgl009.mtx", package ="Matrix")))## NOTE: The following examples take quite some time## ---- even on a fast internet connection:if(FALSE){## The URL has been corrected, but we need an untar step:
u. <- url("https://www.cise.ufl.edu/research/sparse/RB/Boeing/msc00726.tar.gz")
str(sm <- readHB(gzcon(u.)))}
data(KNex, package ="Matrix")## Store as MatrixMarket (".mtx") file, here inside temporary dir./folder:(MMfile <- file.path(tempdir(),"mmMM.mtx"))
writeMM(KNex$mm, file=MMfile)
file.info(MMfile)[,c("size","ctime")]# (some confirmation of the file's)## very simple export - in triplet format - to text file:
data(CAex, package ="Matrix")
s.CA <- summary(CAex)
s.CA # shows (i, j, x) [columns of a data frame]
message("writing to ", outf <- tempfile())
write.table(s.CA, file = outf, row.names=FALSE)## and read it back -- showing off sparseMatrix():
str(dd <- read.table(outf, header=TRUE))## has columns (i, j, x) -> we can use via do.call() as arguments to sparseMatrix():
mm <- do.call(sparseMatrix, dd)
stopifnot(all.equal(mm, CAex, tolerance=1e-15))
a character string indicating a factor in the
factorization represented by x, typically an element
of names(expand2(x, ...)).
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 y should be
multiplied on the left by the factor, rather than on the right.
...
further arguments passed to or from methods.
Details
facmul is experimental and currently no methods are
exported from Matrix.
Value
The value of op(M) %*% y or y %*% op(M),
depending on left, where M is the factor
(always withoutdimnames) and op(M)
is M or t(M), depending on trans.
Examples
## Conceptually, methods for 'facmul' _would_ behave as follows ...## Not run:
n <-3L
x <- lu(Matrix(rnorm(n * n), n, n))
y <- rnorm(n)
L <- unname(expand2(x)[[nm <-"L"]])
stopifnot(exprs ={
all.equal(facmul(x, nm, y, trans =FALSE, left =TRUE), L %*% y)
all.equal(facmul(x, nm, y, trans =FALSE, left =FALSE), y %*% L)
all.equal(facmul(x, nm, y, trans =TRUE, left =TRUE), crossprod(L, y))
all.equal(facmul(x, nm, y, trans =TRUE, left =FALSE), tcrossprod(y, L))})## End(Not run)
“Semi-API” functions used internally by Matrix,
often to bypass S4 dispatch and avoid the associated overhead.
These are exported to provide this capability to expert users.
Typical users should continue to rely on S4 generic functions
to dispatch suitable methods, by calling,
e.g., as(., <class>) for coercions.
a string (".", ",", "n", "l",
or "d") specifying the “kind” of the result.
"." indicates that the kind of from should be preserved.
"," is equivalent to "z" if from is complex
and to "d" otherwise.
"n" indicates that the result should inherit from
nMatrix or nsparseVector
(and so on).
shape
a string (".", "g", "s", or
"t") specifying the “shape” of the result. "."
indicates that the shape of from should be preserved.
"g" indicates that the result should inherit from
generalMatrix (and so on).
repr
a string (".", "C", "R", or
"T") specifying the sparse representation of the result.
"." is accepted only by .ind2sparse and indicates
the most efficient representation,
which is "C" ("R") for margin = 2 (1).
"C" indicates that the result should inherit from
CsparseMatrix (and so on).
a logical indicating if the result should inherit
from sparseMatrix rather than from
denseMatrix. If NA, then the result
will be formally sparse if and only if from is.
uplo
a string ("U" or "L") indicating whether
the result (if symmetric or triangular) should store the upper or
lower triangle of from. The elements of from in the
opposite triangle are ignored.
diag
a string ("N" or "U") indicating whether
the result (if triangular) should be formally nonunit or unit
triangular. In the unit triangular case, the diagonal elements
of from are ignored.
trans
a logical indicating if the result should be a 1-row
matrix rather than a 1-column matrix where from is a vector
but not a matrix.
class
a string whose first three characters specify the class
of the result. It should match the pattern
"^[.nld](ge|sy|tr|sp|tp)" for .m2dense and
"^[.nld][gst][CRT]" for .m2sparse,
where "." in the first position is equivalent to "l"
for logical arguments and "d" for numeric arguments.
a logical indicating if the first argument should be
tested for inheritance from dgCMatrix and
coerced if necessary. Set to FALSE for speed only if it
is known to already inherit from dgCMatrix.
object
a Cholesky factorization inheriting from virtual class
CHMfactor, almost always the result of a call to generic
function Cholesky.
a numeric vector of postive length. Only the first
element is used, and that must be finite.
Details
Functions with names of the form .<A>2<B> implement coercions
from virtual class A to the “nearest” non-virtual subclass of
virtual class B, where the virtual classes are abbreviated as follows:
Abbreviations should be seen as a guide, rather than as an
exact description of behaviour. Notably, .m2dense,
.m2sparse, and .m2V accept vectors that are
not matrices.
.tCRT(x)
If lazy = TRUE, then .tCRT constructs the transpose
of x using the most efficient representation,
which for ‘CRT’ is ‘RCT’. If lazy = FALSE,
then .tCRT preserves the representation of x,
behaving as the corresponding methods for generic function t.
.diag.dsC(x)
.diag.dsC computes (or uses if Chx is supplied)
the Cholesky factorization of x as LDL′ 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, ...).
Examples
D. <- diag(x = c(1,1,2,3,5,8))
D.0<- Diagonal(x = c(0,0,0,3,5,8))
S. <- toeplitz(as.double(1:6))
C. <- new("dgCMatrix", Dim = c(3L,4L),
p = c(0L,1L,1L,1L,3L), i = c(1L,0L,2L), x = c(-8,2,3))
stopifnot(exprs ={
identical(.M2tri (D.), as(D.,"triangularMatrix"))
identical(.M2sym (D.), as(D.,"symmetricMatrix"))
identical(.M2diag(D.), as(D.,"diagonalMatrix"))
identical(.M2kind(C.,"l"),
as(C.,"lMatrix"))
identical(.M2kind(.sparse2dense(C.),"l"),
as(as(C.,"denseMatrix"),"lMatrix"))
identical(.diag2sparse(D.0,".","t","C"),
.dense2sparse(.diag2dense(D.0,".","t",TRUE),"C"))
identical(.M2gen(.diag2dense(D.0,".","s",FALSE)),
.sparse2dense(.M2gen(.diag2sparse(D.0,".","s","T"))))
identical(S.,
.M2m(.m2sparse(S.,".sR")))
identical(S. * lower.tri(S.)+ diag(1,6L),
.M2m(.m2dense (S.,".tr","L","U")))
identical(.M2R(C.), .M2R(.M2T(C.)))
identical(.tCRT(C.), .M2R(t(C.)))})
A <- tcrossprod(C.)/6+ Diagonal(3,1/3); A[1,2]<-3; A
stopifnot(exprs ={
is.numeric( x. <- c(2.2,0,-1.2))
all.equal(x., .solve.dgC.lu(A, c(1,0,0), check=FALSE))
all.equal(x., .solve.dgC.qr(A, c(1,0,0), check=FALSE))})## Solving sparse least squares:
X <- rbind(A, Diagonal(3))# design matrix X (for L.S.)
Xt <- t(X)# *transposed* X (for L.S.)(y <- drop(crossprod(Xt,1:3))+ c(-1,1)/1000)# small rand.err.
str(solveCh <- .solve.dgC.chol(Xt, y, check=FALSE))# Xt *is* dgC..
stopifnot(exprs ={
all.equal(solveCh$coef,1:3, tol =1e-3)# rel.err ~ 1e-4
all.equal(solveCh$coef, drop(solve(tcrossprod(Xt), Xt %*% y)))
all.equal(solveCh$coef, .solve.dgC.qr(X, y, check=FALSE))})
Force a square matrix x to a symmetricMatrix,
without a symmetry check as it would be applied for as(x,
"symmetricMatrix").
Usage
forceSymmetric(x, uplo)
Arguments
x
any square matrix (of numbers), either “"traditional"”
(matrix) or inheriting from
Matrix.
uplo
optional string, "U" or "L" indicating which
“triangle” half of x should determine the result. The
default is "U" unless x already has a uplo slot
(i.e., when it is symmetricMatrix, or
triangularMatrix), where the default will be
x@uplo.
symmpart for the symmetric part of a matrix, or
the coercions as(x, <symmetricMatrix class>).
Examples
## 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.
character which should be used for
structural zeroes. The default "." may occasionally
be replaced by " " (blank); using "0" would look
almost like print()ing of non-sparse matrices.
align
a string specifying how the zero.print codes
should be aligned, see formatSpMatrix.
should the matrix be formatted as a logical matrix
(or rather as a numeric one); mostly for formatSparseM().
uniDiag
logical indicating if the diagonal entries of a sparse
unit triangular or unit-diagonal matrix should be formatted as
"I" instead of "1" (to emphasize that the 1's are
“structural”).
digits
significant digits to use for printing, see
print.default.
cx
(optional) character matrix; a formatted version of x, still
with strings such as "0.00" for the zeros.
iN0
(optional) integer vector, specifying the location of the
non-zeroes of x.
dimnames to be used; a list (of length two)
with row and column names (or NULL).
Value
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!
Author(s)
Martin Maechler
See Also
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.
Examples
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))
a list of
MatrixFactorization objects caching
factorizations of the matrix. Typically, it is initialized
as an empty list and updated “automagically” whenever
a factorization is computed.
Methods for function image in package
Matrix. An image of a matrix simply color codes all matrix
entries and draws the n×m matrix using an
n×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.
a Matrix object, i.e., fulfilling is(x, "Matrix").
xlim, ylim
x- and y-axis limits; may be used to “zoom
into” matrix. Note that x,y “feel reversed”:
ylim is for the rows (= 1st index) and xlim for the
columns (= 2nd index). For convenience, when the limits are integer
valued, they are both extended by 0.5; also, ylim is
always used decreasingly.
aspect
aspect ratio specified as number (y/x) or string;
see levelplot.
sub, xlab, ylab
axis annotation with sensible defaults;
see plot.default.
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,
panel.levelplot.raster (from lattice
package) is used, and the colorkey is also done via rasters, see
also levelplot and possibly
grid.raster.
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 abs(x) should be
shown; if TRUE, the former (implicit) default, the default
col.regions will be grey colors (and no
colorkey drawn). The default is FALSE unless the
matrix has no negative entries.
colorkey
logical indicating if a color key aka ‘legend’
should be produced. Default is to draw one, unless useAbs is
true. You can also specify a list, see
levelplot, such aslist(raster=TRUE) in
the case of rastering.
col.regions
vector of gradually varying colors; see
levelplot.
lwd
(only used when useRaster is false:) non-negative
number or NULL (default), specifying the line-width of the
rectangles of each non-zero matrix entry (drawn by
grid.rect). The default depends on the matrix
dimension and the device size.
border.col
color for the border of each rectangle. NA
means no border is drawn. When NULL as by default,
border.col <- if(lwd < .01) NA else NULL is used.
Consider using an opaque color instead of NULL which
corresponds to grid::get.gpar("col").
...
further arguments passed to methods and
levelplot, notably at for specifying
(possibly non equidistant) cut values for dividing the matrix
values (superseding cuts above).
Value
as all lattice graphics functions, image(<Matrix>)
returns a "trellis" object, effectively the result of
levelplot().
Methods
All methods currently end up calling the method for the
dgTMatrix class.
Use showMethods(image) to list them all.
showMethods(image)## And if you want to see the method definitions:
showMethods(image, includeDefs =TRUE, inherited =FALSE)
data(CAex, package ="Matrix")
image(CAex, main ="image(CAex)")-> imgC; imgC
stopifnot(!is.null(leg <- imgC$legend), is.list(leg$right))# failed for 2 days ..
image(CAex, useAbs=TRUE, main ="image(CAex, useAbs=TRUE)")
cCA <- Cholesky(crossprod(CAex), Imult =.01)## See ?print.trellis --- place two image() plots side by side:
print(image(cCA, main="Cholesky(crossprod(CAex), Imult = .01)"),
split=c(x=1,y=1,nx=2, ny=1), more=TRUE)
print(image(cCA, useAbs=TRUE),
split=c(x=2,y=1,nx=2,ny=1))
data(USCounties, package ="Matrix")
image(USCounties)# huge
image(sign(USCounties))## just the pattern# how the result looks, may depend heavily on# the device, screen resolution, antialiasing etc# e.g. x11(type="Xlib") may show very differently than cairo-based## Drawing borders around each rectangle;# again, viewing depends very much on the device:
image(USCounties[1:400,1:200], lwd=.1)## Using (xlim,ylim) has advantage : matrix dimension and (col/row) indices:
image(USCounties, c(1,200), c(1,400), lwd=.1)
image(USCounties, c(1,300), c(1,200), lwd=.5)
image(USCounties, c(1,300), c(1,200), lwd=.01)## These 3 are all equivalent :(I1 <- image(USCounties, c(1,100), c(1,100), useAbs=FALSE))
I2 <- image(USCounties, c(1,100), c(1,100), useAbs=FALSE, border.col=NA)
I3 <- image(USCounties, c(1,100), c(1,100), useAbs=FALSE, lwd=2, border.col=NA)
stopifnot(all.equal(I1, I2, check.environment=FALSE),
all.equal(I2, I3, check.environment=FALSE))## using an opaque border color
image(USCounties, c(1,100), c(1,100), useAbs=FALSE, lwd=3, border.col = adjustcolor("skyblue",1/2))if(interactive()|| nzchar(Sys.getenv("R_MATRIX_CHECK_EXTRA"))){## Using raster graphics: For PDF this would give a 77 MB file,## however, for such a large matrix, this is typically considerably## *slower* (than vector graphics rectangles) in most cases :if(doPNG <-!dev.interactive())
png("image-USCounties-raster.png", width=3200, height=3200)
image(USCounties, useRaster =TRUE)# should not suffer from anti-aliasingif(doPNG)
dev.off()## and now look at the *.png image in a viewer you can easily zoom in and out}#only if(doExtras)
The indMatrix class is the class of row and column
index matrices, stored as 1-based integer index vectors.
A row (column) index matrix is a matrix whose rows (columns)
are standard unit vectors. Such matrices are useful
when mapping observations to discrete sets of covariate values.
Multiplying a matrix on the left by a row index matrix is
equivalent to indexing its rows, i.e., sampling the rows
“with replacement”. Analogously, multiplying a matrix
on the right by a column index matrix is equivalent to
indexing its columns. Indeed, such products are implemented in
Matrix as indexing operations; see ‘Details’ below.
A matrix whose rows and columns are standard unit vectors
is called a permutation matrix. This special case is
designated by the pMatrix class, a direct
subclass of indMatrix.
Details
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
RM
=
R %*% M
=
M[p, ]
MC
=
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) column-indexing
of a column index matrix (when dimensions are not dropped).
Most other operations on indMatrix treat them as sparse
nonzero pattern matrices (i.e., inheriting from virtual class
nsparseMatrix). Hence vector-valued subsets
of indMatrix, such as those given by diag,
are always of type "logical".
Objects from the Class
Objects can be created explicitly with calls of the form
new("indMatrix", ...), but they are more commonly created
by coercing 1-based integer index vectors, with calls of the
form as(., "indMatrix"); see ‘Methods’ below.
Slots
margin
an integer, either 1 or 2, specifying
whether the matrix is a row (1) or column (2) index.
perm
a 1-based integer index vector, i.e.,
a vector of length Dim[margin] with elements
taken from 1:Dim[1+margin%%2].
signature(x = "indMatrix", y = "Matrix")
and others listed by showMethods("%*%", classes = "indMatrix"):
matrix products implemented where appropriate as indexing operations.
coerce
signature(from = "numeric", to = "indMatrix"):
supporting typical indMatrix construction from
a vector of positive integers. Row indexing is assumed.
coerce
signature(from = "list", to = "indMatrix"):
supporting indMatrix construction for row and
column indexing, including index vectors of length 0 and
index vectors whose maximum is less than the number of rows
or columns being indexed.
coerce
signature(from = "indMatrix", to = "matrix"):
coercion to a traditional matrix of logical type,
with FALSE and TRUE in place of 0 and 1.
t
signature(x = "indMatrix"):
the transpose, which is an indMatrix with identical
perm but opposite margin.
rowSums,rowMeans,colSums,colMeans
signature(x = "indMatrix"):
row and column sums and means.
rbind2,cbind2
signature(x = "indMatrix", y = "indMatrix"):
row-wise catenation of two row index matrices with equal numbers
of columns and column-wise catenation of two column index matrices
with equal numbers of rows.
kronecker
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”.
Author(s)
Fabian Scheipl at ‘uni-muenchen.de’, building on the existing class
pMatrix after a nice hike's conversation with
Martin Maechler. Methods for crossprod(x, y) and
kronecker(x, y) with both arguments inheriting from
indMatrix were made considerably faster thanks to a suggestion
by Boris Vaillant. Diverse tweaks by Martin Maechler and
Mikael Jagan, notably the latter's implementation of margin,
prior to which the indMatrix class was designated only for
row index matrices.
See Also
Subclass pMatrix of permutation matrices,
a special case of index matrices; virtual class
nMatrix of nonzero pattern matrices,
and its subclasses.
Examples
p1 <- as(c(2,3,1),"pMatrix")(sm1 <- as(rep(c(2,3,1), e=3),"indMatrix"))
stopifnot(all(sm1 == p1[rep(1:3, each=3),]))## row-indexing of a <pMatrix> turns it into an <indMatrix>:
class(p1[rep(1:3, each=3),])
set.seed(12)# so we know '10' is in sample## random index matrix for 30 observations and 10 unique values:(s10 <- as(sample(10,30, replace=TRUE),"indMatrix"))## Sample rows of a numeric matrix :(mm <- matrix(1:10, nrow=10, ncol=3))
s10 %*% mm
set.seed(27)
IM1 <- as(sample(1:20,100, replace=TRUE),"indMatrix")
IM2 <- as(sample(1:18,100, replace=TRUE),"indMatrix")(c12 <- crossprod(IM1,IM2))## same as cross-tabulation of the two index vectors:
stopifnot(all(c12 - unclass(table(IM1@perm, IM2@perm))==0))# 3 observations, 4 implied values, first does not occur in sample:
as(2:4,"indMatrix")# 3 observations, 5 values, first and last do not occur in sample:
as(list(2:4,5),"indMatrix")
as(sm1,"nMatrix")
s10[1:7,1:4]# gives an "ngTMatrix" (most economic!)
s10[1:4,]# preserves "indMatrix"-class
I1 <- as(c(5:1,6:4,7:3),"indMatrix")
I2 <- as(7:1,"pMatrix")(I12 <- rbind(I1, I2))
stopifnot(is(I12,"indMatrix"),
identical(I12, rbind(I1, I2)),
colSums(I12)== c(2L,2:4,4:2))
Class index is a virtual class designating index vectors,
or “subscripts”, for (possibly named) vectors and arrays.
It is typically used in signatures of methods for the subscript
and subassignment operators, namely [ and [<-.
It is implemented as a union of the atomic vector classes
numeric, logical,
and character.
invertPerm and signPerm compute the inverse and sign
of a length-n permutation vector. isPerm tests
if a length-n integer vector is a valid permutation vector.
asPerm coerces a length-m transposition vector to a
length-n permutation vector, where m <= n.
Usage
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)
Arguments
p
an integer vector of length n.
pivot
an integer vector of length m.
off
an integer offset, indicating that p is
a permutation of off+0:(n-1) or that pivot
contains m values sampled with replacement from
off+0:(n-1).
ioff
an integer offset, indicating that the result
should be a permutation of ioff+0:(n-1).
n
a integer greater than or equal to m,
indicating the length of the result. Transpositions
are applied to a permutation vector vector initialized
as seq_len(n).
zero.p
a logical. Equivalent to off=0 if TRUE
and off=1 if FALSE.
zero.res
a logical. Equivalent to ioff=0 if TRUE
and ioff=1 if FALSE.
Details
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.
Value
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).
