SPARK-5017 [MLlib] - Use SVD to compute determinant and inverse of covariance matrix #3871

tgaloppo · 2015-01-01T15:38:56Z

MultivariateGaussian was calling both pinv() and det() on the covariance matrix, effectively performing two matrix decompositions. Both values are now computed using the singular value decompositon. Both the pseudo-inverse and the pseudo-determinant are used to guard against singular matrices.

… matrix

SparkQA · 2015-01-01T15:42:32Z

Test build #24987 has started for PR 3871 at commit fd9784c.

This patch merges cleanly.

SparkQA · 2015-01-01T17:09:26Z

Test build #24987 has finished for PR 3871 at commit fd9784c.

This patch passes all tests.
This patch merges cleanly.
This patch adds no public classes.

AmplabJenkins · 2015-01-01T17:09:30Z

Test PASSed.
Refer to this link for build results (access rights to CI server needed):
https://amplab.cs.berkeley.edu/jenkins//job/SparkPullRequestBuilder/24987/
Test PASSed.

jkbradley · 2015-01-02T21:07:47Z

@tgaloppo Could you please add a description? It can be based off of the JIRA, just enough to cover the main points of the PR. Thanks!

… matrix. Code was calling both det() and pinv(), effectively performing two matrix decompositions. Futhermore, Breeze pinv() currently fails for singular matrices.

jkbradley · 2015-01-02T21:45:58Z

mllib/src/main/scala/org/apache/spark/mllib/stat/impl/MultivariateGaussian.scala

Perhaps you could add a note here about how this behaves when sigma is singular, plus a reference like [http://en.wikipedia.org/wiki/Multivariate_normal_distribution#Degenerate_case]

The doc could be a short version of what you have below for calculateCovarianceConstants

jkbradley · 2015-01-02T21:48:33Z

@tgaloppo @mengxr What are your thoughts about doing the computation in log space as much as possible, and then exponentiating at the end? I'm mainly thinking about numerical stability, but I could imagine wanting to provide pdf() and logpdf() methods eventually.

jkbradley · 2015-01-02T21:49:22Z

@tgaloppo The logic looks good; my comments are basically about clarity (except for the log space question). Thanks for the PR!

jkbradley · 2015-01-02T21:51:13Z

One more request: Could you please add a unit test with a singular matrix? Thank you! Perhaps in a new suite for MultivariateGaussian

…park-5017

SparkQA · 2015-01-02T21:57:33Z

Test build #24998 has started for PR 3871 at commit b4415ea.

This patch merges cleanly.

SparkQA · 2015-01-02T23:23:39Z

Test build #24998 has finished for PR 3871 at commit b4415ea.

This patch passes all tests.
This patch merges cleanly.
This patch adds no public classes.

AmplabJenkins · 2015-01-02T23:23:43Z

Test PASSed.
Refer to this link for build results (access rights to CI server needed):
https://amplab.cs.berkeley.edu/jenkins//job/SparkPullRequestBuilder/24998/
Test PASSed.

tgaloppo · 2015-01-03T00:24:07Z

@jkbradley I think performing the pdf calculation in log-space (and providing a logpdf() method) is a good idea. Perhaps we can make this part of transitioning MultivariateGaussian to public scope?

…ovariance matrix. Previous code called both pinv() and det(), effectively performing two matrix decompositions. Additionally, the pinv() implementation in Breeze is known to fail for singular matrices.

…ate case.

SparkQA · 2015-01-03T01:37:31Z

Test build #25001 has started for PR 3871 at commit d448137.

This patch merges cleanly.

jkbradley · 2015-01-03T03:03:13Z

mllib/src/main/scala/org/apache/spark/mllib/stat/impl/MultivariateGaussian.scala

"eigenvalues" --> "singular values" (here and in the next comment on line 79)

SparkQA · 2015-01-03T03:03:31Z

Test build #25001 has finished for PR 3871 at commit d448137.

This patch passes all tests.
This patch merges cleanly.
This patch adds no public classes.

SparkQA · 2015-01-03T05:49:03Z

Test build #25006 has finished for PR 3871 at commit dc3d0f7.

This patch passes all tests.
This patch merges cleanly.
This patch adds no public classes.

AmplabJenkins · 2015-01-03T05:49:06Z

Test PASSed.
Refer to this link for build results (access rights to CI server needed):
https://amplab.cs.berkeley.edu/jenkins//job/SparkPullRequestBuilder/25006/
Test PASSed.

jkbradley · 2015-01-03T19:29:05Z

mllib/src/test/scala/org/apache/spark/mllib/stat/impl/MultivariateGaussianSuite.scala

Good catch with the other var; could you please fix this and the test above too? var -> val

@jkbradley No problem; fixed.

jkbradley · 2015-01-03T19:30:02Z

@tgaloppo thanks for verifying about the test pdf values

SparkQA · 2015-01-03T20:02:32Z

Test build #25016 has started for PR 3871 at commit 629d9d0.

This patch merges cleanly.

SparkQA · 2015-01-03T21:29:15Z

Test build #25016 has finished for PR 3871 at commit 629d9d0.

This patch passes all tests.
This patch merges cleanly.
This patch adds no public classes.

AmplabJenkins · 2015-01-03T21:29:19Z

Test PASSed.
Refer to this link for build results (access rights to CI server needed):
https://amplab.cs.berkeley.edu/jenkins//job/SparkPullRequestBuilder/25016/
Test PASSed.

jkbradley · 2015-01-04T20:46:14Z

mllib/src/main/scala/org/apache/spark/mllib/stat/impl/MultivariateGaussian.scala

For link, use double-brackets:

(See [[http://en.wikipedia.org/wiki/Multivariate_normal_distribution#Degenerate_case]])

jkbradley · 2015-01-04T20:47:49Z

@tgaloppo 2 more small comments, but after those, I believe this will be ready. Thanks!

Fixed comment in MultivariateGaussian

tgaloppo · 2015-01-04T21:29:31Z

@jkbradley Good call on the test suite; I have added some non-center points to the tests. I also added the brackets to the in-comment link.

SparkQA · 2015-01-04T21:32:29Z

Test build #25038 has started for PR 3871 at commit a5b8bc5.

This patch merges cleanly.

jkbradley · 2015-01-04T21:45:01Z

@tgaloppo Thanks for the updates! LGTM after tests pass

CC: @mengxr

SparkQA · 2015-01-04T22:42:37Z

Test build #25038 has finished for PR 3871 at commit a5b8bc5.

This patch passes all tests.
This patch merges cleanly.
This patch adds no public classes.

AmplabJenkins · 2015-01-04T22:42:40Z

Test PASSed.
Refer to this link for build results (access rights to CI server needed):
https://amplab.cs.berkeley.edu/jenkins//job/SparkPullRequestBuilder/25038/
Test PASSed.

mengxr · 2015-01-05T19:57:37Z

mllib/src/test/scala/org/apache/spark/mllib/stat/impl/MultivariateGaussianSuite.scala

We can construct BDV and BDM directly.

MultivariateGaussianSuite - Create Breeze vectors and matrices directly instead of through MLlib vectors/matrices.

tgaloppo · 2015-01-06T00:21:48Z

@mengxr Changes made.

SparkQA · 2015-01-06T00:22:37Z

Test build #25069 has started for PR 3871 at commit 383b5b3.

This patch merges cleanly.

SparkQA · 2015-01-06T01:34:03Z

Test build #25069 has finished for PR 3871 at commit 383b5b3.

This patch passes all tests.
This patch merges cleanly.
This patch adds no public classes.

AmplabJenkins · 2015-01-06T01:34:07Z

Test PASSed.
Refer to this link for build results (access rights to CI server needed):
https://amplab.cs.berkeley.edu/jenkins//job/SparkPullRequestBuilder/25069/
Test PASSed.

tgaloppo · 2015-01-06T02:17:05Z

mllib/src/main/scala/org/apache/spark/mllib/stat/impl/MultivariateGaussian.scala

U D^-1 * U.t => U D^(-1/2) D^(-1/2) U.t => (U D^(-1/2)) (D^(-1/2) U.t)
...both are U and D are symmetric, so...
(U D^(-1/2)) = (U.t (D^(-1/2)).t) => (D^(-1/2) U).t
and
(D^(-1/2) U.t) = (D^(-1/2) U)
thus
U D^-1 U.t => (D^(-1/2) U).t (D^(-1/2) U)
... bringing in the delta we get
delta.t (D^(-1/2) U).t (D^(-1/2) U) delta
=> ((D^(-1/2) U) delta).t (D^(-1/2) U) delta = norm(D^(-1/2) U delta)^2
as indicated by @mengxr

(phew, hope I did that OK! :) )

@jkbradley Hmm. I am going to punt on this one. Perhaps @mengxr can point out my error. FWIW - Changing the code to U * D^(-1/2) causes the unit tests to fail.

I agree; I see no reason why U should always be symmetric (as you have demonstrated). Before I roll back the change, I just want to make sure that I did not misinterpret @mengxr , and/or that we are not missing something that he sees.

(u, d) is the eigendecomposition of sigma, so sigma = u * diag(d) * u^-1 ... but we have a special case since covariance matrices are always symmetric and positive semi-definite, in which case u * u.t = I, making it equivalent to the singular value decomposition... so sigma = u * diag(d) * u.t ... so in svd terms, v.t = u.t, then the inverse is v * inv(diag(d)) * u.t = u * inv(diag(d)) * u.t ...

Have I lost my bearings?

Ugh, no, I have. I was confused about which were inverses, and what you wrote looks perfectly fine. Sorry for the trouble! I'll remove the comments.d

jkbradley · 2015-01-06T21:13:23Z

@tgaloppo Now that my confusion is over...LGTM Thanks very much!
CC: @mengxr

mengxr · 2015-01-06T21:58:10Z

Merged into master. Thanks!

SPARK-5017 - Use SVD to compute determinant and inverse of covariance…

fd9784c

… matrix

tgaloppo changed the title ~~SPARK-5017 - Use SVD to compute determinant and inverse of covariance matrix~~ SPARK-5017 [MLlib] - Use SVD to compute determinant and inverse of covariance matrix Jan 1, 2015

SPARK-5017 - Use SVD to compute determinant and inverse of covariance…

6f11b6d

… matrix. Code was calling both det() and pinv(), effectively performing two matrix decompositions. Futhermore, Breeze pinv() currently fails for singular matrices.

jkbradley reviewed Jan 2, 2015
View reviewed changes

Merge branch 'spark-5017' of https://github.com/tgaloppo/spark into s…

b4415ea

…park-5017

tgaloppo added 2 commits January 2, 2015 20:29

SPARK-5017 - Fixed to use SVD to compute determinant and inverse of c…

1989be0

…ovariance matrix. Previous code called both pinv() and det(), effectively performing two matrix decompositions. Additionally, the pinv() implementation in Breeze is known to fail for singular matrices.

Added test suite for MultivariateGaussian, including test for degener…

d448137

…ate case.

jkbradley reviewed Jan 3, 2015
View reviewed changes

Moved some test values from var to val.

629d9d0

jkbradley reviewed Jan 4, 2015
View reviewed changes

Added additional points to tests in test suite.

a5b8bc5

Fixed comment in MultivariateGaussian

mengxr reviewed Jan 5, 2015
View reviewed changes

MultivariateGaussian - minor optimization in density calculation

383b5b3

MultivariateGaussianSuite - Create Breeze vectors and matrices directly instead of through MLlib vectors/matrices.

tgaloppo reviewed Jan 6, 2015
View reviewed changes

asfgit closed this in 4108e5f Jan 6, 2015

SPARK-5017 [MLlib] - Use SVD to compute determinant and inverse of covariance matrix #3871

SPARK-5017 [MLlib] - Use SVD to compute determinant and inverse of covariance matrix #3871

Uh oh!

Conversation

tgaloppo commented Jan 1, 2015

Uh oh!

SparkQA commented Jan 1, 2015

Uh oh!

SparkQA commented Jan 1, 2015

Uh oh!

AmplabJenkins commented Jan 1, 2015

Uh oh!

jkbradley commented Jan 2, 2015

Uh oh!

Choose a reason for hiding this comment

Uh oh!

jkbradley commented Jan 2, 2015

Uh oh!

jkbradley commented Jan 2, 2015

Uh oh!

jkbradley commented Jan 2, 2015

Uh oh!

SparkQA commented Jan 2, 2015

Uh oh!

SparkQA commented Jan 2, 2015

Uh oh!

AmplabJenkins commented Jan 2, 2015

Uh oh!

tgaloppo commented Jan 3, 2015

Uh oh!

SparkQA commented Jan 3, 2015

Uh oh!

Choose a reason for hiding this comment

Uh oh!

SparkQA commented Jan 3, 2015

Uh oh!

SparkQA commented Jan 3, 2015

Uh oh!

AmplabJenkins commented Jan 3, 2015

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

jkbradley commented Jan 3, 2015

Uh oh!

SparkQA commented Jan 3, 2015

Uh oh!

SparkQA commented Jan 3, 2015

Uh oh!

AmplabJenkins commented Jan 3, 2015

Uh oh!

Choose a reason for hiding this comment

Uh oh!

jkbradley commented Jan 4, 2015

Uh oh!

tgaloppo commented Jan 4, 2015

Uh oh!

SparkQA commented Jan 4, 2015

Uh oh!

jkbradley commented Jan 4, 2015

Uh oh!

SparkQA commented Jan 4, 2015

Uh oh!

AmplabJenkins commented Jan 4, 2015

Uh oh!

Choose a reason for hiding this comment

Uh oh!

tgaloppo commented Jan 6, 2015

Uh oh!

SparkQA commented Jan 6, 2015

Uh oh!

SparkQA commented Jan 6, 2015

Uh oh!

AmplabJenkins commented Jan 6, 2015

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!