Math does not have the matrix (inverse) square root

Description

You can apply the sqrt function elementwise to a matrix, but Stan Math does not have a function to take the (symmetric) square root of a positive definite matrix (or its inverse), defined as V * diag_matrix(sqrt(v)) * V' where V is a matrix of eigenvectors and v is a vector of eigenvalues. In the case of the inverse square root of a matrix, there is a diag_matrix(inv_sqrt(v)) in the middle. Eigen already calculates these, so it is just a matter of exposing them and working through the Jacobians.

Example

data {
  int<lower = 1> K;
  matrix[K, K] X;
}
transformed data {
  matrix[K, K] sqrt_X = sqrt_spd(X);
  matrix[K, K] inv_sqrt_X = inv_sqrt_spd(X);
  real<lower = 0, upper = 1e-15> check = fabs(X - sqrt_X * sqrt_X);
}

[Note: The sqrt_spd and inv_sqrt_spd functions do not exist in Stan yet.]

Expected Output

Satisfies inequality restrictions on check.

Current Version:

v2.33