https://github.com/cran/Matrix
Tip revision: 23d9d53f0b3c5fcd18cc73551fe6e531e08c9f0e authored by Martin Maechler on 14 November 2023, 11:55:00 UTC
version 1.6-3
version 1.6-3
Tip revision: 23d9d53
posdef.R
## METHODS FOR CLASS: dpoMatrix, dppMatrix
## dense symmetric positive semidefinite matrices
## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Operations such as rounding can lose positive semidefiniteness
## but not symmetry, hence:
.indefinite <- function(x) {
cl <- .M.nonvirtual(x, TRUE)
if(any(cl == c("dpoMatrix", "corMatrix")))
as(x, "dsyMatrix")
else if(any(cl == c("dppMatrix", "pcorMatrix")))
as(x, "dspMatrix")
else x
}
.dsy2dpo <- .dsp2dpp <- function(from) {
if(is.null(tryCatch(Cholesky(from, perm = FALSE),
error = function(e) NULL)))
stop("not a positive definite matrix (and positive semidefiniteness is not checked)")
to <- new(.CLASS)
to@Dim <- from@Dim
to@Dimnames <- from@Dimnames
to@uplo <- from@uplo
to@x <- from@x
to@factors <- from@factors
to
}
body(.dsy2dpo)[[3L]][[3L]][[2L]] <- "dpoMatrix"
body(.dsp2dpp)[[3L]][[3L]][[2L]] <- "dppMatrix"
setAs("dppMatrix", "dpoMatrix", function(from) unpack(from))
setAs("dpoMatrix", "dppMatrix", function(from) pack(from))
setAs("dsyMatrix", "dpoMatrix", .dsy2dpo)
setAs("dspMatrix", "dppMatrix", .dsp2dpp)
setAs("Matrix", "dpoMatrix",
function(from) .dsy2dpo(.M2unpacked(.M2sym(.M2kind(from, "d")))))
setAs("Matrix", "dppMatrix",
function(from) .dsp2dpp(.M2packed (.M2sym(.M2kind(from, "d")))))
setAs("matrix", "dpoMatrix",
function(from) {
storage.mode(from) <- "double"
.dsy2dpo(.M2sym(from))
})
setAs("matrix", "dppMatrix",
function(from) {
storage.mode(from) <- "double"
.dsp2dpp(pack(from, symmetric = TRUE))
})
## METHODS FOR CLASS: corMatrix, pcorMatrix
## dense correlation matrices
## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
.dpo2cor <- function(from) {
if(!is.null(to <- from@factors$correlation))
return(to)
sd <- sqrt(diag(from, names = FALSE))
to <- new("corMatrix")
to@Dim <- d <- from@Dim
to@Dimnames <- from@Dimnames
to@uplo <- from@uplo
to@sd <- sd
n <- d[1L]
x <- from@x / sd / rep(sd, each = n)
x[indDiag(n)] <- 1
to@x <- x
.set.factor(from, "correlation", to)
}
.dpp2pcor <- function(from) {
if(!is.null(to <- from@factors$correlation))
return(to)
sd <- sqrt(diag(from, names = FALSE))
to <- new("pcorMatrix")
to@Dim <- d <- from@Dim
to@Dimnames <- from@Dimnames
to@uplo <- uplo <- from@uplo
to@sd <- sd
n <- d[1L]
u <- uplo == "U"
if(u) {
r <- seq_len(n)
s <- 1L
} else {
r <- seq.int(to = 1L, by = -1L, length.out = n)
s <- seq_len(n)
}
x <- from@x / rep.int(sd, r) / sd[sequence.default(r, s)]
x[indDiag(n, upper = u, packed = TRUE)] <- 1
to@x <- x
.set.factor(from, "correlation", to)
}
setAs("pcorMatrix", "corMatrix", function(from) unpack(from))
setAs( "corMatrix", "pcorMatrix", function(from) pack(from))
setAs("dpoMatrix", "corMatrix", .dpo2cor)
setAs("dppMatrix", "pcorMatrix", .dpp2pcor)
setAs("dsyMatrix", "corMatrix",
function(from) .dpo2cor (.dsy2dpo(from)))
setAs("dspMatrix", "pcorMatrix",
function(from) .dpp2pcor(.dsp2dpp(from)))
setAs("Matrix", "corMatrix",
function(from) .dpo2cor (.dsy2dpo(.M2unpacked(.M2sym(.M2kind(from, "d"))))))
setAs("Matrix", "pcorMatrix",
function(from) .dpp2pcor(.dsp2dpp(.M2packed (.M2sym(.M2kind(from, "d"))))))
setAs("matrix", "corMatrix",
function(from) {
storage.mode(from) <- "double"
.dpo2cor (.dsy2dpo(.M2sym(from)))
})
setAs("matrix", "pcorMatrix",
function(from) {
storage.mode(from) <- "double"
.dpp2pcor(.dsp2dpp(pack(from, symmetric = TRUE)))
})
## METHODS FOR GENERIC: cov2cor
## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
setMethod("cov2cor", signature(V = "unpackedMatrix"),
function(V) {
d <- V@Dim
if(d[1L] != d[2L] || .M.kind(V) == "z")
stop(gettextf("'%s' is not a square numeric matrix", "V"),
domain = NA)
as(forceSymmetric(V), "corMatrix")
})
setMethod("cov2cor", signature(V = "packedMatrix"),
function(V) {
d <- V@Dim
if(d[1L] != d[2L] || .M.kind(V) == "z")
stop(gettextf("'%s' is not a square numeric matrix", "V"),
domain = NA)
as(forceSymmetric(V), "pcorMatrix")
})
setMethod("cov2cor", signature(V = "sparseMatrix"),
function(V) {
d <- V@Dim
if(d[1L] != d[2L] || .M.kind(V) == "z")
stop(gettextf("'%s' is not a square numeric matrix", "V"),
domain = NA)
dn <- symDN(V@Dimnames)
V <- .M2kind(V, "d")
V.ii <- diag(V, names = FALSE)
if(length(V.ii) > 0L && is.na(m <- min(V.ii)) || m <= 0)
warning(gettextf("diag(%s) has non-positive or non-finite entries; finite result is doubtful",
"V"),
domain = NA)
D <- Diagonal(x = sqrt(1/V.ii))
r <- forceSymmetric(D %*% V %*% D)
diag(r) <- 1
r@Dimnames <- dn
r
})