Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Math does not have the matrix (inverse) square root #1144

Open
bgoodri opened this issue Mar 7, 2019 · 2 comments
Open

Math does not have the matrix (inverse) square root #1144

bgoodri opened this issue Mar 7, 2019 · 2 comments
Assignees
@bgoodri
Labels
feature new function

Comments

@bgoodri
Copy link
Contributor

bgoodri commented Mar 7, 2019 •
edited by bob-carpenter
Loading

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
@bgoodri bgoodri self-assigned this Mar 7, 2019
@bgoodri bgoodri mentioned this issue Mar 7, 2019
Closed
5 tasks
bbbales2 added a commit that referenced this issue Mar 8, 2019
@bbbales2
Updated the input and output types of the sqrt_spd adj_jac_apply (Issue
f85bae9 
#1144)

Rearranged some includes in stan/math/rev/mat.hpp to get rid of some type traits errors. Not sure why these were moved around. Might need moved back.
@bob-carpenter
Copy link
Contributor

I'd like to boost the priority on this one and I suggest we name the function matrix_sqrt(). Sarah Heaps has an elegant unconstrained parameterization of stationary parameters for a vector autoregressive (VAR) model. It uses matrix square root. She sent along the code she was using for it in Stan, without the new tuple feature

matrix sqrtm(matrix A) {
    int m = rows(A);
    vector[m] root_root_evals = sqrt(sqrt(eigenvalues_sym(A)));
    matrix[m, m] evecs = eigenvectors_sym(A);
    matrix[m, m] eprod = diag_post_multiply(evecs, root_root_evals);
    return tcrossprod(eprod);
  }

and with:

matrix sqrtm(matrix A) {
    int m = rows(A);
    tuple(matrix[m, m], vector[m]) eigen = eigendecompose_sym(A);
    matrix[m, m] eprod = diag_post_multiply(eigen.1, sqrt(sqrt(eigen.2)));
    return tcrossprod(eprod);
  }

I'd also like to add the full constraining transform and the unconstraining transform, but that will be in a separate issue.

@bgoodri
Copy link
Contributor Author

bgoodri commented Jun 21, 2024

If we were to boost the priority of this, then https://github.com/KingJamesSong/FastDifferentiableMatSqrt is better than the closed PR I wrote a long time ago.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees

@bgoodri bgoodri

Labels
feature new function
Projects
None yet
Milestone
No milestone
Development

Successfully merging a pull request may close this issue.

[WIP] initial implementation of sqrt_spd and inv_sqrt_spd stan-dev/math
2 participants
@bob-carpenter @bgoodri