dyadFac {DyadiCarma} | R Documentation |
Efficient factorization of a positive definite symmetrically dyadic matrix.
Description
This function implement the efficient factorization of a
positive definite symmetrically dyadic matrix \boldsymbol \Sigma
.
It computes the vertically dyadic matrix \mathbf P
such that
\mathbf P^\top \boldsymbol \Sigma \mathbf P = \mathbf I
.
Usage
dyadFac(S, inv = FALSE, band = FALSE)
Arguments
S |
A |
inv |
The boolean value indicating whether the inverse of
|
band |
The boolean value indicating whether the input S is a band matrix. If TRUE, then a optimized band-focused algorithm is called. If band==TRUE, but the input matrix is not a band one, the function will return the corresponding result for the band part of the input matrix. |
Details
This function implement the efficient factorization of a positive definite symmetrically dyadic matrix.
Value
If inv == TRUE
, then the inverse of \boldsymbol \Sigma
,
which is a (2^(height)-1)*breadth x (2^(height)-1)*breadth
classic
matrix, is returned. Otherwise, the vertically Dyadic
object for
\mathbf P
is returned.
See Also
Dyadic-class
for a description of the class;
Examples
#-------------------------------------------------------------#
#---------Inverting a PD symmetrically dyadic matrix----------#
#-------------------------------------------------------------#
N <- 4
k <- 3
# A 45x45 vertically dyadic matrix
V <- construct(N, k, type = "vert", distr = "unif")
# A 45x45 symmetrically dyadic matrix
S <- t(V) %*% V
S@type <- "symm"
S@aentries <- list() # Convert S from "asymm" to "symm"
# Check what S looks like
matS <- as.matrix(S)
matS
# Find the vertically dyadic matrix that satisfies P^T S P = I
# using a dyadic factorization algorithm.
P <- dyadFac(S)
I1 <- as.matrix(t(P) %*% S %*% P)
I <- diag(dim(I1)[1])
max(abs(I1 - I)) # Should be trivially small
# Obtain the inverse of S via the dyadic algorithm
iS <- dyadFac(S, inv = TRUE)
I2 <- iS %*% matS
max(abs(I2 - I)) # Should be trivially small
iS_solve <- solve(matS)
I3 <- iS_solve %*% matS
max(abs(I3 - I)) # The result obtained using built-in method for inversion
#-------------------------------------------------------------#
#-----------------Inverting a PD band matrix------------------#
#-------------------------------------------------------------#
d <- k * (2^N - 1)
half_B <- matrix(0, nrow = d, ncol = d)
for (i in 1:d) {
half_B[i, i:min(d, (i + k - 1))] <- rnorm(min(d, (i + k - 1)) - i + 1, mean = N, sd = 1 / N)
}
matB <- t(half_B) %*% half_B # matB is a PD band matrix with half bandwidth 3.
# Convert matB into a dyadic object B
B <- as.dyadic(matB, "symm", N, k)
iB <- dyadFac(B, inv = TRUE)
I <- diag(dim(matB)[1])
max(abs(iB %*% matB - I)) # Should be trivially small
iB_band <- dyadFac(B, inv = TRUE, band = TRUE)
max(abs(iB_band %*% matB - I)) # Should be trivially small
iB <- dyadFac(B)
iB_band <- dyadFac(B, band = TRUE)
max(abs(as.matrix(iB) - as.matrix(iB_band))) # Should be trivially small