ENH: stats.gaussian_kde: replace use of inv_cov in pdf #16692
Merged
Conversation
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
|
I think all the failures are the same - |
|
Rather than all the permutations required to replace |
|
@steppi if you also like linear algebra, this may be interesting to you. |
Putting this next in my queue. |
mdhaber
commented
Aug 20, 2022
There was a problem hiding this comment.
Trying to fix series of errors like:
___________________ test_kde_output_dtype[float128-float128] ___________________
[gw1] darwin -- Python 3.10.5 /Users/runner/miniconda3/envs/scipy-dev/bin/python
scipy/stats/tests/test_kdeoth.py:328: in test_kde_output_dtype
result = k(points)
bw = 3.0
bw_type = <class 'numpy.float128'>
dataset = array([0., 1., 2., 3., 4.], dtype=float128)
dtype = <class 'numpy.float128'>
k = <scipy.stats._kde.gaussian_kde object at 0x140000160>
points = array([0., 1., 2., 3., 4.], dtype=float128)
weights = array([0., 1., 2., 3., 4.], dtype=float128)
scipy/stats/_kde.py:270: in evaluate
result = gaussian_kernel_estimate[spec](
d = 1
itemsize = 16
m = 5
output_dtype = <class 'numpy.float128'>
points = array([[0., 1., 2., 3., 4.]], dtype=float128)
self = <scipy.stats._kde.gaussian_kde object at 0x140000160>
spec = 'long double'
_stats.pyx:748: in scipy.stats._stats.gaussian_kernel_estimate
???
E ValueError: Buffer dtype mismatch, expected 'long double' but got 'double'
mdhaber
commented
Aug 20, 2022
steppi
reviewed
Aug 21, 2022
There was a problem hiding this comment.
I haven't had time to check everything, but things look OK mathematically so far. See the suggestion about adding some comments. I found the permutations a little inscrutable at first, so I think some comments to explain what's going on and link to more details would be helpful.
I should have time to complete my review next weekend.
Comment on lines
+585
to
+586
| self._data_cho_cov = linalg.cholesky( | ||
| self._data_covariance[::-1, ::-1]).T[::-1, ::-1] |
There was a problem hiding this comment.
I think there's enough going on here that there should be some comments to explain things. It took me a bit to figure out what's happening.
Just to make I'm on the same page, here are the details as I understand them:
Let $\Gamma$ be the covariance matrix for the Gaussian kernel and $LL^{\top} = \Gamma^{-1}$ be the Cholesky decomposition of the precision matrix $\Gamma^{-1}$ . Later on, we want to be able to transform $n \times d$ data matrices $X$ by multiplying on the right by the Cholesky factor $L$ for $\Gamma^{-1}$ . Equivalently, we want to be able to find $Z$ such that $Z = XL$ .
If we know
If
and thus
where we've used that
This means that if we know the Cholesky factor
I think we should have a comment explaining the types of equations we want to be able to solve
There was a problem hiding this comment.
It definitely does deserve an explanation, and I meant to include one before you got to this. Sorry to make you find it on your own. Yes, I think that is the original post I followed. I'll write a bit about the motivation and link to that.
There was a problem hiding this comment.
No worries. The explanation is very clear now. I think everything is in good shape now but still want to double check carefully.
There was a problem hiding this comment.
I didn't follow the whole argument but if applies here, use the lower=False keyword for starting with an upper triangular in calling cholesky. Might save a column swap or two.
There was a problem hiding this comment.
I didn't follow the whole argument but if applies here, use the
lower=Falsekeyword for starting with an upper triangular in callingcholesky. Might save a column swap or two.
It's a good thought but isn't quite what we want. According to the documentation for the underlying Lapack function, if
cholesky(Gamma, lower=True) will return L and cholesky(Gamma, lower=False) will return cholesky function can't do this so we're required to do the trick with reversing the rows and columns.
There was a problem hiding this comment.
Oh, but I guess it will save us one transpose when computing the Cholesky decomposition of
mdhaber
commented
Aug 22, 2022
mdhaber
commented
Aug 22, 2022
steppi
approved these changes
Aug 23, 2022
| @property | ||
| def inv_cov(self): | ||
| self.factor = self.covariance_factor() | ||
| self._data_covariance = atleast_2d(cov(self.dataset, rowvar=1, |
There was a problem hiding this comment.
Why do we recompute self._data_covariance here? Does it help to maintain backwards compatibility for existing subclasses? If this is needed, could probably use a comment to explain why but it’s fine if you think it isn’t necessary.
There was a problem hiding this comment.
I think I assumed that since _compute_covariance used to re-calculate the covariance every time, there must have been a reason. Otherwise, why not do it just once in the __init__ method? As it was, it was re-calculated every time set_bandwidth was called, so maybe people use set_bandwidth to recalculate everything after modifying the public attribute dataset? I don't really know, but figured it would be safer this way.
And maybe subconsciously I want use of inv_cov to be as slow as possible : ) See the discussion in the original incarnation of this issue - #5087 (comment).
mdhaber
commented
Aug 23, 2022
| @property | ||
| def inv_cov(self): | ||
| self.factor = self.covariance_factor() | ||
| self._data_covariance = atleast_2d(cov(self.dataset, rowvar=1, |
There was a problem hiding this comment.
I think I assumed that since _compute_covariance used to re-calculate the covariance every time, there must have been a reason. Otherwise, why not do it just once in the __init__ method? As it was, it was re-calculated every time set_bandwidth was called, so maybe people use set_bandwidth to recalculate everything after modifying the public attribute dataset? I don't really know, but figured it would be safer this way.
And maybe subconsciously I want use of inv_cov to be as slow as possible : ) See the discussion in the original incarnation of this issue - #5087 (comment).
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
Reference issue
Supersedes gh-5087
What does this implement/fix?
gh-5087 proposed replacing use of the inverse covariance matrix with the Cholesky decomposition of the covariance matrix throughout
gaussian_kdeto improve speed and avoid numerical instabilities associate with matrix inversion. There didn't seem to be disagreement from a technical standpoint; it looks like development just stopped.This PR implements the suggestion for
gaussian_kde.pdf. Sincelogpdfis being Cythonized in gh-15493, I'll leave that alone to avoid merge conflicts.Additional information
Here is the timing of the KDE benchmarks in main (results from CI of gh-16684):
stats.GaussianKDE.time_gaussian_kde_evaluate_few_points- 1.31±0msstats.GaussianKDE.time_gaussian_kde_evaluate_many_points- 1.61±0sIn this PR:
stats.GaussianKDE.time_gaussian_kde_evaluate_few_points- 645±0μsstats.GaussianKDE.time_gaussian_kde_evaluate_many_points- 905±0msWe might be able to do better with a closer look at new Cython/Python interactions.