lgspline (version 0.2.0)

matinvsqrt: Calculate Matrix Inverse Square Root

Description

Calculate Matrix Inverse Square Root

Usage

matinvsqrt(mat)

Value

A matrix \(\textbf{B}\) such that \(\textbf{B}\textbf{B} = \textbf{M}^{-1}\)

Arguments

mat

A symmetric, positive-definite matrix \(\textbf{M}\)

Details

For matrix \(\textbf{M}\), computes \(\textbf{B}\) where \(\textbf{B}\textbf{B} = \textbf{M}^{-1}\) using eigenvalue decomposition:

1. Compute eigendecomposition \(\textbf{M} = \textbf{V}\textbf{D}\textbf{V}^T\)

2. Set eigenvalues below sqrt(.Machine$double.eps) to 0

3. Take elementwise reciprocal square root: \(\textbf{D}^{-1/2}\)

4. Reconstruct as \(\textbf{B} = \textbf{V} \textbf{D}^{-1/2} \textbf{V}^T\)

For nearly singular matrices, eigenvalues below the numerical threshold are set to 0, and their reciprocals in \(\textbf{D}^{-1/2}\) are also set to 0.

This implementation is particularly useful for whitening procedures in GLMs with correlation structures and for computing variance-covariance matrices under constraints.

You can use this to help construct a custom VhalfInv_fxn for fitting correlation structures, see lgspline.

Examples

## Identity matrix
m1 <- diag(2)
matinvsqrt(m1)  # Returns identity matrix

## Compound symmetry correlation matrix
rho <- 0.5
m2 <- matrix(rho, 3, 3) + diag(1-rho, 3)
B <- matinvsqrt(m2)
# Verify: B %**% B approximately equals solve(m2)
all.equal(B %**% B, solve(m2))

## Example for GLM correlation structure
n_blocks <- 2  # Number of subjects
block_size <- 3  # Measurements per subject
rho <- 0.7  # Within-subject correlation
# Correlation matrix for one subject
R <- matrix(rho, block_size, block_size) +
     diag(1-rho, block_size)
## Full correlation matrix for all subjects
V <- kronecker(diag(n_blocks), R)
## Create whitening matrix
VhalfInv <- matinvsqrt(V)

# Example construction of VhalfInv_fxn for lgspline
VhalfInv_fxn <- function(par) {
  rho <- tanh(par)  # Transform parameter to (-1, 1)
  R <- matrix(rho, block_size, block_size) +
       diag(1-rho, block_size)
  kronecker(diag(n_blocks), matinvsqrt(R))
}