dpoMatrix-class.RdThe "dpoMatrix" class is the class of
positive-semidefinite symmetric matrices in nonpacked storage.
The "dppMatrix" class is the same except in packed
storage. Only the upper triangle or the lower triangle is
required to be available.
The "corMatrix" and "copMatrix" classes
represent correlation matrices. They extend "dpoMatrix"
and "dppMatrix", respectively, with an additional slot
sd allowing restoration of the original covariance matrix.
Objects can be created by calls of the
form new("dpoMatrix", ...) or from crossprod applied to
an "dgeMatrix" object.
uplo:Object of class "character". Must be
either "U", for upper triangular, and "L", for lower triangular.
x:Object of class "numeric". The numeric
values that constitute the matrix, stored in column-major order.
Dim:Object of class "integer". The dimensions
of the matrix which must be a two-element vector of non-negative
integers.
Dimnames:inherited from class "Matrix"
factors:Object of class "list". A named
list of factorizations that have been computed for the matrix.
sd:(for "corMatrix" and "copMatrix")
a numeric vector of length n containing the
(original) \(\sqrt{var(.)}\) entries which allow
reconstruction of a covariance matrix from the correlation matrix.
Class "dsyMatrix", directly.
Classes "dgeMatrix", "symmetricMatrix", and many more
by class "dsyMatrix".
signature(x = "dpoMatrix"):
Returns (and stores) the Cholesky decomposition of x, see
chol.
signature(x = "dpoMatrix"):
Returns the determinant of x, via
chol(x), see above.
signature(x = "dpoMatrix", norm = "character"):
Returns (and stores) the reciprocal of the condition number of
x. The norm can be "O" for the
one-norm (the default) or "I" for the infinity-norm. For
symmetric matrices the result does not depend on the norm.
signature(a = "dpoMatrix", b = "...."), and
signature(a = "dppMatrix", b = "....") work
via the Cholesky composition, see also the Matrix solve-methods.
signature(e1 = "dpoMatrix", e2 = "numeric") (and
quite a few other signatures): The result of (“elementwise”
defined) arithmetic operations is typically not
positive-definite anymore. The only exceptions, currently, are
multiplications, divisions or additions with positive
length(.) == 1 numbers (or logicals).
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")) .
h6 <- Hilbert(6)
rcond(h6)
#> [1] 3.439939e-08
str(h6)
#> Formal class 'dpoMatrix' [package "Matrix"] with 5 slots
#> ..@ Dim : int [1:2] 6 6
#> ..@ Dimnames:List of 2
#> .. ..$ : NULL
#> .. ..$ : NULL
#> ..@ x : num [1:36] 1 0.5 0.333 0.25 0.2 ...
#> ..@ uplo : chr "U"
#> ..@ factors :List of 1
#> .. ..$ Cholesky:Formal class 'Cholesky' [package "Matrix"] with 5 slots
#> .. .. .. ..@ uplo : chr "U"
#> .. .. .. ..@ x : num [1:36] 1 0 0 0 0 ...
#> .. .. .. ..@ perm : int(0)
#> .. .. .. ..@ Dim : int [1:2] 6 6
#> .. .. .. ..@ Dimnames:List of 2
#> .. .. .. .. ..$ : NULL
#> .. .. .. .. ..$ : NULL
h6 * 27720 # is ``integer''
#> 6 x 6 Matrix of class "dpoMatrix"
#> [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,] 27720 13860 9240 6930 5544 4620
#> [2,] 13860 9240 6930 5544 4620 3960
#> [3,] 9240 6930 5544 4620 3960 3465
#> [4,] 6930 5544 4620 3960 3465 3080
#> [5,] 5544 4620 3960 3465 3080 2772
#> [6,] 4620 3960 3465 3080 2772 2520
solve(h6)
#> 6 x 6 Matrix of class "dpoMatrix"
#> [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,] 36 -630 3360 -7560 7560 -2772
#> [2,] -630 14700 -88200 211680 -220500 83160
#> [3,] 3360 -88200 564480 -1411200 1512000 -582120
#> [4,] -7560 211680 -1411200 3628800 -3969000 1552320
#> [5,] 7560 -220500 1512000 -3969000 4410000 -1746360
#> [6,] -2772 83160 -582120 1552320 -1746360 698544
str(hp6 <- pack(h6))
#> Formal class 'dppMatrix' [package "Matrix"] with 5 slots
#> ..@ uplo : chr "U"
#> ..@ Dim : int [1:2] 6 6
#> ..@ Dimnames:List of 2
#> .. ..$ : NULL
#> .. ..$ : NULL
#> ..@ x : num [1:21] 1 0.5 0.333 0.333 0.25 ...
#> ..@ factors :List of 1
#> .. ..$ Cholesky:Formal class 'Cholesky' [package "Matrix"] with 5 slots
#> .. .. .. ..@ uplo : chr "U"
#> .. .. .. ..@ x : num [1:36] 1 0 0 0 0 ...
#> .. .. .. ..@ perm : int(0)
#> .. .. .. ..@ Dim : int [1:2] 6 6
#> .. .. .. ..@ Dimnames:List of 2
#> .. .. .. .. ..$ : NULL
#> .. .. .. .. ..$ : NULL
### Note that as(*, "corMatrix") *scales* the matrix
(ch6 <- as(h6, "corMatrix"))
#> 6 x 6 Matrix of class "corMatrix"
#> [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,] 1.0000000 0.8660254 0.7453560 0.6614378 0.6000000 0.5527708
#> [2,] 0.8660254 1.0000000 0.9682458 0.9165151 0.8660254 0.8206518
#> [3,] 0.7453560 0.9682458 1.0000000 0.9860133 0.9583148 0.9270248
#> [4,] 0.6614378 0.9165151 0.9860133 1.0000000 0.9921567 0.9749960
#> [5,] 0.6000000 0.8660254 0.9583148 0.9921567 1.0000000 0.9949874
#> [6,] 0.5527708 0.8206518 0.9270248 0.9749960 0.9949874 1.0000000
stopifnot(all.equal(as(h6 * 27720, "dsyMatrix"), round(27720 * h6),
tolerance = 1e-14),
all.equal(ch6@sd^(-2), 2*(1:6)-1,
tolerance = 1e-12))
chch <- Cholesky(ch6, perm = FALSE)
stopifnot(identical(chch, ch6@factors$Cholesky),
all(abs(crossprod(as(chch, "dtrMatrix")) - ch6) < 1e-10))