[WIP] initial implementation of sqrt_spd and inv_sqrt_spd #1145

bgoodri · 2019-03-07T15:51:25Z

Summary

This implements the symmetric square root of a positive definite matrix (sqrt_spd) and its inverse (inv_sqrt_spd) and fixes #1144 .

These call the operatorSqrt and operatorInverseSqrt methods of Eigen's self-adjoint eigensolver http://eigen.tuxfamily.org/dox/classEigen_1_1SelfAdjointEigenSolver.html .

Tests

There are tests that the difference between an arbitrary positive definite matrix and its sqrt_spd multiplied by itself are zero, but there is an assert when evaluating adj_jac_apply that I am hoping @bbbales2 knows about. Same thing for the inverse square root.

Side Effects

No

Checklist

Math issue Math does not have the matrix (inverse) square root #1144
Copyright holder: Trustees of Columbia University

The copyright holder is typically you or your assignee, such as a university or company. By submitting this pull request, the copyright holder is agreeing to the license the submitted work under the following licenses:
- Code: BSD 3-clause (https://opensource.org/licenses/BSD-3-Clause)
- Documentation: CC-BY 4.0 (https://creativecommons.org/licenses/by/4.0/)
the basic tests are passing (Not for the new tests)
- unit tests pass (to run, use: ./runTests.py test/unit)
- header checks pass, (make test-headers)
- docs build, (make doxygen)
- code passes the built in C++ standards checks (make cpplint)
the code is written in idiomatic C++ and changes are documented in the doxygen
the new changes are tested

adj_jac_apply() and jacobian() seem to not work together

bbbales2

Is there an assert happening in the sqrt_spd_tests that are causing this to puke? Is that why they're commented out?

There's a specialization for sqrt_spd, but is there a specialization for inv_sqrt_spd?

stan/math/rev/mat/fun/sqrt_spd.hpp

bbbales2 · 2019-03-07T16:38:49Z

I reread the commit, looks like the answer is yes. I'll check this out in the afternoooooon.

bgoodri · 2019-03-07T18:21:15Z

Yeah, I used to be doing the sqrt_spd test like the inv_sqrt_spd test, which seemed fine. And then I added the adj_jac_apply stuff to rev mode sqrt_spd and jacobian() would not cooperate anymore.

…

On Thu, Mar 7, 2019 at 11:38 AM Ben Bales ***@***.***> wrote: I reread the commit, looks like the answer is yes. I'll check this out in the afternoooooon. — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#1145 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/ADOrqr1LtC1jVjK0TwuFBz2IS1WqZksDks5vUUCagaJpZM4bjg9y> .

bgoodri · 2019-03-07T22:54:07Z

Someone can kill the continuous integration tests for this. It is not ready to be merged and isn't necessary for 2.19.0.

#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.

…stable/2017-11-14)

bbbales2 · 2019-03-08T00:57:22Z

I wiggled around the types on the multiply_adjoint_jacobian. I dunno if I got the Kronecker stuff right. The sqrt_spd tests run now.

There was also some includes that got shuffled in stan/math/rev/mat.hpp. I moved them back. I was getting weird type traits errors.

bbbales2 · 2019-03-08T01:09:19Z

The type of the thing passed to multiply_adjoint_jacobian in one of these functors matches the output type of operator(). The return type of multiply_adjoint_jacobian is a tuple of the input types to operator() (excluding is_var), so that's where the type stuff came from.

stan/math/rev/mat/fun/sqrt_spd.hpp

bgoodri · 2019-03-08T04:04:49Z

stan/math/rev/mat/fun/sqrt_spd.hpp

    using Eigen::kroneckerProduct;
    using Eigen::MatrixXd;
+    Eigen::MatrixXd output(K_, K_);
+    Eigen::Map<const Eigen::VectorXd> map_input(adj.data(), adj.size());


And if adj is a MatrixXd, how can it be mapped to a VectorXd?

That mapping is doing this (https://en.wikipedia.org/wiki/Vectorization_(mathematics)). I was just pattern-matching the eqs in the Stackoverflow link. I'm not really familiar with that. I was just double checking myself and Abhishek makes a note at the end of his post that the Jacobian he wrote is in terms of the vectorization, which I assume is the Wikipedia definition.

Oh yeah fair warning, I was just kinda guessing on the math stuff to make things compile. Maybe add in a diagonal matrix test (with things that aren't 1 on the diagonal) so we know the answer and gradients and all just to make sure things are in the right direction.

bgoodri · 2019-03-08T04:07:12Z

stan/math/rev/mat/fun/sqrt_spd.hpp

    using Eigen::kroneckerProduct;
    using Eigen::MatrixXd;
+    Eigen::MatrixXd output(K_, K_);
+    Eigen::Map<const Eigen::VectorXd> map_input(adj.data(), adj.size());
+    Eigen::Map<Eigen::VectorXd> map_output(output.data(), output.size());


Similar, how can map_output be a Map<VectorXd> but output be a square MatrixXd?

And why is map_output necessary? Can we just do the linear algebra to form output and do std::make_tuple(output)?

The maps are just an eigen thing for mapping one data structure onto another. So the math happens in terms of the vectors like Abhishek had it but it just ends up in the right place in the matrices. I think it should be okay in terms of efficiency?

[CI skip]

…stable/2017-11-14)

bgoodri · 2019-03-08T04:52:05Z

Ah, OK. I had to print a bunch of stuff out but now I think I understand what it is doing.

…

On Thu, Mar 7, 2019 at 11:47 PM Ben Bales ***@***.***> wrote: ***@***.**** commented on this pull request. ------------------------------ In stan/math/rev/mat/fun/sqrt_spd.hpp <#1145 (comment)>: > multiply_adjoint_jacobian(const std::array<bool, size> &needs_adj, - const Eigen::VectorXd &adj) const { + const Eigen::MatrixXd &adj) const { adj was getting vectorized to be K^2 x 1. This is what the Map is doing. It's taking the KxK adj and then looking at it likes it's a K^2x1 vector. — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#1145 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/ADOrqg3bsdUbPX2_FejlIbZXbdH-QUHQks5vUetzgaJpZM4bjg9y> .

bgoodri · 2019-03-08T22:04:07Z

@bbbales2 Is there a way to automate the "second check" mentioned before the conclusions section of https://people.maths.ox.ac.uk/gilesm/files/NA-08-01.pdf when utilizing adj_jac_apply?

bbbales2 · 2019-03-09T02:23:04Z

Is there a way to automate the "second check" mentioned before the conclusions section of https://people.maths.ox.ac.uk/gilesm/files/NA-08-01.pdf when utilizing adj_jac_apply?

Yeah, looks like it'd work (in a hand-wavy I didn't actually think through the code kinda way). Looks like something that would be handy in other places too. @bob-carpenter, you seen this?

bgoodri · 2019-03-09T02:38:51Z

Would also be good in mix/

…

On Fri, Mar 8, 2019 at 9:23 PM Ben Bales ***@***.***> wrote: Is there a way to automate the "second check" mentioned before the conclusions section of https://people.maths.ox.ac.uk/gilesm/files/NA-08-01.pdf when utilizing adj_jac_apply? Yeah, looks like it'd work (in a hand-wavy I didn't actually think through the code kinda way). Looks like something that would be handy in other places too. @bob-carpenter <https://github.com/bob-carpenter>, you seen this? — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#1145 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/ADOrqjU4tjOulUffxCISSOFrCKL1ndCBks5vUxsJgaJpZM4bjg9y> .

…into feature/matrix_sqrt

…gs/RELEASE_500/final)

bob-carpenter · 2019-03-11T17:32:58Z

I have not seen that. My matrix fu is weak.

It says that if A, B lead to an output C, then the identity

Tr(transpose(bar(C)) * dot(C) 
= Tr(transpose(bar(A)) * dot(A))
  + Tr(transpose(bar(B)) * dot(B))

where I think the bar operation is gradient.

bgoodri · 2019-03-11T20:17:45Z

It is the same Giles paper cited in the Stan Math paper on archiv.org. The notation is given at the start of section 2.1

syclik · 2019-10-24T03:50:42Z

Closing due to inactivity. @bgoodri, I hope you get a chance to finish this.

bbbales2 · 2020-04-13T17:30:45Z

@bgoodri just pinging you in case you want to resurrect this.

SteveBronder · 2020-07-16T15:16:30Z

@bgoodri are you still working on this?

bgoodri · 2020-07-16T16:52:41Z

I want to finish it, but haven't been able to deal with anything except R packages and StanCon recently.

…

On Thu, Jul 16, 2020 at 11:16 AM Steve Bronder ***@***.***> wrote: @bgoodri <https://github.com/bgoodri> are you still working on this? — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#1145 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAZ2XKQFWENAOYW6QXVPBSTR34KWBANCNFSM4G4OB5ZA> .

syclik · 2020-07-31T04:54:27Z

Thanks, @bgoodri. I'm going to convert it to a draft so we know it's not ready yet.

syclik · 2020-11-10T17:38:40Z

@bgoodri, I'm going to close this with a WIP label so it can be found again. If you'd like to continue on it, please feel free to reopen.

initial implementation of sqrt_spd and inv_sqrt_spd

73bf8b3

adj_jac_apply() and jacobian() seem to not work together

bbbales2 requested changes Mar 7, 2019

View reviewed changes

stan/math/rev/mat/fun/sqrt_spd.hpp Show resolved Hide resolved

bbbales2 and others added 3 commits March 7, 2019 16:46

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.

Merge commit '9d0acc533a2c8656e14bc5db4f4da72ee99cf699' into HEAD

3f5807a

[Jenkins] auto-formatting by clang-format version 6.0.0 (tags/google/…

5480a40

…stable/2017-11-14)

bgoodri commented Mar 8, 2019

View reviewed changes

stan/math/rev/mat/fun/sqrt_spd.hpp Outdated Show resolved Hide resolved

bgoodri commented Mar 8, 2019

View reviewed changes

bgoodri and others added 3 commits March 7, 2019 23:08

fix order of headers

f89b6a5

bump size of test up

ca80605

[CI skip]

[Jenkins] auto-formatting by clang-format version 6.0.0 (tags/google/…

2fa0a27

…stable/2017-11-14)

add check for NaNs

2b1f2fe

bgoodri and others added 3 commits March 9, 2019 19:46

more usage of adj_jac_apply()

1751a76

Merge branch 'feature/matrix_sqrt' of https://github.com/stan-dev/math …

a6daa9e

…into feature/matrix_sqrt

[Jenkins] auto-formatting by clang-format version 5.0.0-3~16.04.1 (ta…

1adbee1

…gs/RELEASE_500/final)

serban-nicusor-toptal added this to the 2.20.0++ milestone Jul 18, 2019

serban-nicusor-toptal modified the milestones: 2.20.0++, v3.0.0 Oct 18, 2019

syclik closed this Oct 24, 2019

serban-nicusor-toptal modified the milestones: 3.0.0++, 3.1.0 Jan 24, 2020

mcol removed this from the 3.1.0 milestone Feb 3, 2020

bgoodri reopened this Apr 24, 2020

syclik marked this pull request as draft July 31, 2020 04:54

syclik added the WIP label Nov 10, 2020

syclik closed this Nov 10, 2020

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

[WIP] initial implementation of sqrt_spd and inv_sqrt_spd #1145

[WIP] initial implementation of sqrt_spd and inv_sqrt_spd #1145

bgoodri commented Mar 7, 2019

bbbales2 left a comment

bbbales2 commented Mar 7, 2019

bgoodri commented Mar 7, 2019 via email

bgoodri commented Mar 7, 2019

bbbales2 commented Mar 8, 2019

bbbales2 commented Mar 8, 2019

bgoodri Mar 8, 2019

bbbales2 Mar 8, 2019

bbbales2 Mar 8, 2019

bgoodri Mar 8, 2019

bgoodri Mar 8, 2019

bbbales2 Mar 8, 2019

bgoodri commented Mar 8, 2019 via email

bgoodri commented Mar 8, 2019

bbbales2 commented Mar 9, 2019

bgoodri commented Mar 9, 2019 via email

bob-carpenter commented Mar 11, 2019

bgoodri commented Mar 11, 2019

syclik commented Oct 24, 2019

bbbales2 commented Apr 13, 2020

SteveBronder commented Jul 16, 2020

bgoodri commented Jul 16, 2020 via email

syclik commented Jul 31, 2020

syclik commented Nov 10, 2020

[WIP] initial implementation of sqrt_spd and inv_sqrt_spd #1145

[WIP] initial implementation of sqrt_spd and inv_sqrt_spd #1145

Conversation

bgoodri commented Mar 7, 2019

Summary

Tests

Side Effects

Checklist

bbbales2 left a comment

Choose a reason for hiding this comment

bbbales2 commented Mar 7, 2019

bgoodri commented Mar 7, 2019 via email

bgoodri commented Mar 7, 2019

bbbales2 commented Mar 8, 2019

bbbales2 commented Mar 8, 2019

bgoodri Mar 8, 2019

Choose a reason for hiding this comment

bbbales2 Mar 8, 2019

Choose a reason for hiding this comment

bbbales2 Mar 8, 2019

Choose a reason for hiding this comment

bgoodri Mar 8, 2019

Choose a reason for hiding this comment

bgoodri Mar 8, 2019

Choose a reason for hiding this comment

bbbales2 Mar 8, 2019

Choose a reason for hiding this comment

bgoodri commented Mar 8, 2019 via email

bgoodri commented Mar 8, 2019

bbbales2 commented Mar 9, 2019

bgoodri commented Mar 9, 2019 via email

bob-carpenter commented Mar 11, 2019

bgoodri commented Mar 11, 2019

syclik commented Oct 24, 2019

bbbales2 commented Apr 13, 2020

SteveBronder commented Jul 16, 2020

bgoodri commented Jul 16, 2020 via email

syclik commented Jul 31, 2020

syclik commented Nov 10, 2020