Positive Semi-definite Dense (Packed | Non-packed) Numeric Matrices
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 slotsdallowing 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") anumericvector of lengthncontaining the (original) \(\sqrt{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 ofx, seechol.- determinant
signature(x = "dpoMatrix"): Returns thedeterminantofx, viachol(x), see above.- rcond
signature(x = "dpoMatrix", norm = "character"): Returns (and stores) the reciprocal of the condition number ofx. Thenormcan 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 Matrixsolve-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(.) == 1numbers (orlogicals).
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")) .
Examples
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))