check generalized inverse for full rank symmetric mat #2577

spinkney · 2021-09-08T14:08:02Z

Summary

Fixes #2576. Add a check on symmetric matrices for full rank. Uses the ~~FullPivLU~~ completeOrthogonalDecomposition to check the rank and returns Eigen's pseudoInverse from the complete ortho-decomp if rank < full rank (otherwise just do the regular inverse).

Tests

Added a test for a less than full rank symmetric matrix.

Side Effects

No.

Release notes

Replace this text with a short note on what will change if this pull request is merged in which case this will be included in the release notes.

Checklist

Math issue Generalized inverse should work for less than full rank symmetric matrices #2576
Copyright holder: Sean Pinkney
the basic tests are passing
the code is written in idiomatic C++ and changes are documented in the doxygen
the new changes are tested

…4.1 (tags/RELEASE_600/final)

stan/math/prim/fun/generalized_inverse.hpp

…/math into generalized_inverse_fix

…4.1 (tags/RELEASE_600/final)

This reverts commit cf36885.

…4.1 (tags/RELEASE_600/final)

…/math into generalized_inverse_fix

…4.1 (tags/RELEASE_600/final)

Instead of creating two identity matrices, create two vectors and add diag

This reverts commit 27c2647.

…4.1 (tags/RELEASE_600/final)

…/math into generalized_inverse_fix

…4.1 (tags/RELEASE_600/final)

stan/math/rev/fun/generalized_inverse.hpp

test/unit/math/mix/fun/generalized_inverse_test.cpp

stan-buildbot · 2021-09-10T14:01:28Z

Name	Old Result	New Result	Ratio	Performance change( 1 - new / old )
gp_pois_regr/gp_pois_regr.stan	3.47	3.51	0.99	-1.15% slower
low_dim_corr_gauss/low_dim_corr_gauss.stan	0.02	0.02	0.95	-4.8% slower
eight_schools/eight_schools.stan	0.09	0.09	0.98	-1.71% slower
gp_regr/gp_regr.stan	0.14	0.14	0.99	-1.21% slower
irt_2pl/irt_2pl.stan	5.13	5.1	1.0	0.47% faster
performance.compilation	91.12	89.1	1.02	2.22% faster
low_dim_gauss_mix_collapse/low_dim_gauss_mix_collapse.stan	8.27	8.21	1.01	0.75% faster
pkpd/one_comp_mm_elim_abs.stan	29.69	29.66	1.0	0.09% faster
sir/sir.stan	122.13	120.58	1.01	1.26% faster
gp_regr/gen_gp_data.stan	0.03	0.03	0.99	-0.93% slower
low_dim_gauss_mix/low_dim_gauss_mix.stan	2.98	2.98	1.0	-0.05% slower
pkpd/sim_one_comp_mm_elim_abs.stan	0.37	0.37	1.01	0.98% faster
arK/arK.stan	2.05	2.04	1.01	0.61% faster
arma/arma.stan	0.25	0.25	1.0	0.27% faster
garch/garch.stan	0.61	0.6	1.02	1.65% faster
Mean result: 0.999224096962

Jenkins Console Log
Blue Ocean
Commit hash: 2c686b4

Machine information

ProductName: Mac OS X ProductVersion: 10.11.6 BuildVersion: 15G22010

CPU:
Intel(R) Xeon(R) CPU E5-1680 v2 @ 3.00GHz

G++:
Configured with: --prefix=/Applications/Xcode.app/Contents/Developer/usr --with-gxx-include-dir=/usr/include/c++/4.2.1
Apple LLVM version 7.0.2 (clang-700.1.81)
Target: x86_64-apple-darwin15.6.0
Thread model: posix

Clang:
Apple LLVM version 7.0.2 (clang-700.1.81)
Target: x86_64-apple-darwin15.6.0
Thread model: posix

SteveBronder · 2021-09-10T15:55:31Z

stan/math/prim/fun/generalized_inverse.hpp

+  if (G_ref.rows() == G_ref.cols()) {
+    Eigen::CompleteOrthogonalDecomposition<
+        Eigen::Matrix<value_type_t<EigMat>, EigMat::RowsAtCompileTime,
+                      EigMat::ColsAtCompileTime>>
+        complete_ortho_decomp_G = G_ref.completeOrthogonalDecomposition();
+    if (!(complete_ortho_decomp_G.rank() < G_ref.rows()))
+      return inverse(G_ref);


So my only Q is how much this costs / accuracy of computing the completeOrthogonalDecomposition() and then using it for the check and then sometimes using it for the pseudoinverse. If someone passes an NxM matrix here i sort of feels like they either

Would know ahead of time that they will or will not have a square matrix

Know their matrix is low rank and they can't use inverse()

Are those not good assumptions? If they are then should we always just do the pseudoinverse()? Is the accuracy of pseudoinverse pretty low compared to inverse()? I almost feel like even if so someone using this would expect a pseudoinverse instead of the inverse

These are great questions!

On

If someone passes an NxM matrix here i sort of feels like they either

Would know ahead of time that they will or will not have a square matrix

Know their matrix is low rank and they can't use inverse()

I'm interested in the case where you'd specifically pass a square matrix but not know ahead of time if the matrix is low rank or not. It's probably low-rank but it may not be (say in a factor model).

This got me to revisit the paper that the code is based on. The author does something that I missed because he notates L as a cholesky factor when it is actually a modified cholesky factor that works even for low rank matrices. This is how he is able to not do a rank find (which one usually does either as I've done here or using SVD). I've coded up the prim version at https://github.com/spinkney/math/tree/generalized_inverse_low_rank_chol. The point is that some rank issue must be dealt with or there's numerical instability issues.

The code now follows the paper explicitly. The major change involves adding a cholesky_low_rank_decomposition lambda. I've coded it using Eigen block notation. I believe there's further optimizations that could be done, if someone with a background in high performance computing helps. It could be added to Stan as a separate function, which is cool in it's own right.

The prim tests pass. I haven't benchmarked though to see how it compares to the completeOrthogonalDecomposition() version here.

I've put together an example how low rank chol can transform a low rank matrix of "factors" into a valid correlation matrix. I'm not sure if this works all the time but it's promising. low_rank_chol_eg.R.zip

Oh nice! So with the low rank chol version should we close this PR and open one up based on the low rank chol? Or should we merge this then make a separate PR for that?

I'm thinking that the low rank chol should be the final version. Though I think it may be worth just getting this update to fix the symmetric low-rank case (for this release cycle) and then in the next release having 2 PRs: 1) add low rank chol as a function and 2) update gen inverse using it in a new PR.

stan/math/rev/fun/generalized_inverse.hpp

SteveBronder

lgtm!

…4.1 (tags/RELEASE_600/final)

SteveBronder

lgtm!

spinkney · 2021-09-18T09:37:16Z

@rok-cesnovar @SteveBronder this passed before and that commit shouldn't have failed. I'm not sure what that strange CI error is, there's no info, is this a server issue?

rok-cesnovar · 2021-09-18T09:41:01Z

Yeah seems like a machine went down. Restarted only the upstream tests, should be done in an hour: https://jenkins.mc-stan.org/blue/organizations/jenkins/Math%20Pipeline/detail/PR-2577/26/pipeline

stan-buildbot · 2021-09-18T10:36:07Z

Name	Old Result	New Result	Ratio	Performance change( 1 - new / old )
gp_pois_regr/gp_pois_regr.stan	3.49	3.51	0.99	-0.56% slower
low_dim_corr_gauss/low_dim_corr_gauss.stan	0.02	0.02	0.97	-2.64% slower
eight_schools/eight_schools.stan	0.09	0.09	1.02	1.49% faster
gp_regr/gp_regr.stan	0.14	0.14	1.01	1.47% faster
irt_2pl/irt_2pl.stan	5.07	5.09	1.0	-0.22% slower
performance.compilation	90.58	88.82	1.02	1.94% faster
low_dim_gauss_mix_collapse/low_dim_gauss_mix_collapse.stan	8.15	8.2	0.99	-0.6% slower
pkpd/one_comp_mm_elim_abs.stan	31.33	30.25	1.04	3.45% faster
sir/sir.stan	121.09	119.27	1.02	1.5% faster
gp_regr/gen_gp_data.stan	0.03	0.03	1.0	-0.26% slower
low_dim_gauss_mix/low_dim_gauss_mix.stan	3.03	3.05	0.99	-0.68% slower
pkpd/sim_one_comp_mm_elim_abs.stan	0.37	0.37	0.99	-0.66% slower
arK/arK.stan	2.08	2.06	1.01	1.07% faster
arma/arma.stan	0.25	0.25	1.0	0.01% faster
garch/garch.stan	0.6	0.6	1.0	0.0% slower
Mean result: 1.00376021383

Jenkins Console Log
Blue Ocean
Commit hash: 3644737

Machine information

ProductName: Mac OS X ProductVersion: 10.11.6 BuildVersion: 15G22010

CPU:
Intel(R) Xeon(R) CPU E5-1680 v2 @ 3.00GHz

G++:
Configured with: --prefix=/Applications/Xcode.app/Contents/Developer/usr --with-gxx-include-dir=/usr/include/c++/4.2.1
Apple LLVM version 7.0.2 (clang-700.1.81)
Target: x86_64-apple-darwin15.6.0
Thread model: posix

Clang:
Apple LLVM version 7.0.2 (clang-700.1.81)
Target: x86_64-apple-darwin15.6.0
Thread model: posix

spinkney · 2021-09-18T12:32:35Z

@SteveBronder this is good to go!

rok-cesnovar · 2021-09-18T12:36:51Z

@spinkney go ahead and merge. You should be able to.

check symmetric mat for full rank

e1a375b

spinkney changed the title ~~check symmetric mat for full rank~~ check generalized inverse for full rank symmetric mat Sep 8, 2021

[Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.0…

88e3ae2

…4.1 (tags/RELEASE_600/final)

spinkney commented Sep 8, 2021

View reviewed changes

stan/math/prim/fun/generalized_inverse.hpp Outdated Show resolved Hide resolved

spinkney and others added 20 commits September 8, 2021 14:54

remove fullpivLU

294fbb4

Merge branch 'generalized_inverse_fix' of https://github.com/spinkney…

9ef5618

…/math into generalized_inverse_fix

[Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.0…

d459414

…4.1 (tags/RELEASE_600/final)

Update generalized_inverse.hpp

f176c35

fix for expression tests

cf36885

Revert "fix for expression tests"

8a1786d

This reverts commit cf36885.

[Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.0…

cadfa12

…4.1 (tags/RELEASE_600/final)

fix expression tests

e20d39e

Merge branch 'generalized_inverse_fix' of https://github.com/spinkney…

3df6a8d

…/math into generalized_inverse_fix

[Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.0…

0ebe908

…4.1 (tags/RELEASE_600/final)

rank revealing decomp for symmetric matrix

456bd74

[Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.0…

c78fc21

…4.1 (tags/RELEASE_600/final)

small optimization in derivative

27c2647

Instead of creating two identity matrices, create two vectors and add diag

Revert "small optimization in derivative"

9bae32e

This reverts commit 27c2647.

[Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.0…

52c1be3

…4.1 (tags/RELEASE_600/final)

Merge branch 'generalized_inverse_fix' of https://github.com/spinkney…

9c357d7

…/math into generalized_inverse_fix

[Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.0…

3f1ea94

…4.1 (tags/RELEASE_600/final)

Update generalized_inverse.hpp

8d60e00

add braces for if

1670e27

[Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.0…

70a21d0

…4.1 (tags/RELEASE_600/final)

spinkney mentioned this pull request Sep 9, 2021

Release 2.28 checklist stan-dev/cmdstan#1037

Closed

23 tasks

spinkney commented Sep 9, 2021

View reviewed changes

stan/math/rev/fun/generalized_inverse.hpp Outdated Show resolved Hide resolved

SteveBronder reviewed Sep 9, 2021

View reviewed changes

stan/math/rev/fun/generalized_inverse.hpp Show resolved Hide resolved

SteveBronder reviewed Sep 9, 2021

View reviewed changes

stan/math/rev/fun/generalized_inverse.hpp Outdated Show resolved Hide resolved

update for review comments

9a3e82a

SteveBronder reviewed Sep 9, 2021

View reviewed changes

test/unit/math/mix/fun/generalized_inverse_test.cpp Outdated Show resolved Hide resolved

update test names

2c686b4

SteveBronder reviewed Sep 10, 2021

View reviewed changes

SteveBronder requested changes Sep 10, 2021

View reviewed changes

stan/math/rev/fun/generalized_inverse.hpp Outdated Show resolved Hide resolved

SteveBronder previously approved these changes Sep 15, 2021

View reviewed changes

revert add_diag

1de51d0

spinkney dismissed SteveBronder’s stale review via 1de51d0 September 17, 2021 15:49

yashikno and others added 2 commits September 17, 2021 15:55

Merge commit '8920da4297b72d50f3fa456ed8f2cd2bd714265c' into HEAD

7722ddf

[Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.0…

3644737

…4.1 (tags/RELEASE_600/final)

SteveBronder approved these changes Sep 17, 2021

View reviewed changes

spinkney merged commit f8765cf into stan-dev:develop Sep 18, 2021

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

check generalized inverse for full rank symmetric mat #2577

check generalized inverse for full rank symmetric mat #2577

spinkney commented Sep 8, 2021 •

edited

Loading

stan-buildbot commented Sep 10, 2021

SteveBronder Sep 10, 2021

spinkney Sep 11, 2021

spinkney Sep 11, 2021

SteveBronder Sep 15, 2021

spinkney Sep 16, 2021 •

edited

Loading

SteveBronder left a comment

SteveBronder left a comment

spinkney commented Sep 18, 2021

rok-cesnovar commented Sep 18, 2021

stan-buildbot commented Sep 18, 2021

spinkney commented Sep 18, 2021

rok-cesnovar commented Sep 18, 2021

check generalized inverse for full rank symmetric mat #2577

check generalized inverse for full rank symmetric mat #2577

Conversation

spinkney commented Sep 8, 2021 • edited Loading

Summary

Tests

Side Effects

Release notes

Checklist

stan-buildbot commented Sep 10, 2021

SteveBronder Sep 10, 2021

Choose a reason for hiding this comment

spinkney Sep 11, 2021

Choose a reason for hiding this comment

spinkney Sep 11, 2021

Choose a reason for hiding this comment

SteveBronder Sep 15, 2021

Choose a reason for hiding this comment

spinkney Sep 16, 2021 • edited Loading

Choose a reason for hiding this comment

SteveBronder left a comment

Choose a reason for hiding this comment

SteveBronder left a comment

Choose a reason for hiding this comment

spinkney commented Sep 18, 2021

rok-cesnovar commented Sep 18, 2021

stan-buildbot commented Sep 18, 2021

spinkney commented Sep 18, 2021

rok-cesnovar commented Sep 18, 2021

spinkney commented Sep 8, 2021 •

edited

Loading

spinkney Sep 16, 2021 •

edited

Loading