Repository: scipy/scipy
Pull request: BUG: improved the stability of kde for highly (1000s) multi-dimension inputs
<!-- 
Thanks for contributing a pull request! Please ensure that
your PR satisfies the checklist before submitting:
http://scipy.github.io/devdocs/dev/contributor/development_workflow.html#checklist-before-submitting-a-pr

Also, please name and describe your PR as you would write a
commit message:
https://docs.scipy.org/doc/numpy/dev/development_workflow.html#writing-the-commit-message

Note that we are a team of volunteers; we appreciate your
patience during the review process.

Again, thanks for contributing!
-->

#### Reference issue
Closes: #12988 

#### What does this implement/fix?
It appeared that this problem was due to the instability in `linalg.det()` for very large matrices (as in line 572 of [kde.py](https://github.com/scipy/scipy/blob/d286f8525c16b2cd4e179dea2c77b6b09622aff9/scipy/stats/kde.py#L572)). 
This has been fixed by using the `np.linalg.slogdet()` function where possible (not possible for `float128` or `float96`objects which are not supported by `np.linalg`), and carrying the logarithmic result through in the `logpdf()` function.

#### Additional information
Simple `%%timeit` testing suggests that this implementation may actually be (slightly) faster. 

Problem statement:
BUG: Highly multi-dimensional `gaussian_kde` giving `-inf` logpdf
<!-- 
Thank you for taking the time to file a bug report. 
Please fill in the fields below, deleting the sections that 
don't apply to your issue. You can view the final output
by clicking the preview button above.
Note: This is a comment, and won't appear in the output.
-->

My issue is about the fact that a highly multidimensional `gaussian_kde` object (like >1500 dimensions) will regularly return `-inf` regardless of the input for `logpdf()`. 

#### Reproducing code example:
<!-- 
If you place your code between the triple backticks below, 
it will be rendered as a code block. 
-->

```python
import numpy as np
from scipy.stats import gaussian_kde
n = 2000
data = np.array([np.random.randn(5000) + (n*1) for n in range(0, n)])
gkde = gaussian_kde(data)
print(gkde.logpdf(np.arange(0, n)))
```

#### Error message:
<!-- If any, paste the *full* error message inside a code block
as above (starting from line Traceback)
-->

```
[-inf]
```

#### Scipy/Numpy/Python version information:

```
1.5.0 1.18.5 sys.version_info(major=3, minor=8, micro=3, releaselevel='final', serial=0)
```
<!-- You can simply run the following and paste the result in a code block
```
import sys, scipy, numpy; print(scipy.__version__, numpy.__version__, sys.version_info)
```
-->

Hints from issue discussion:

Commit message: account for lack of float96 setup in np.linalg

Patch to review:
diff --git a/scipy/stats/kde.py b/scipy/stats/kde.py
index 9015593bf61e..f0ff88c437c6 100644
--- a/scipy/stats/kde.py
+++ b/scipy/stats/kde.py
@@ -569,7 +569,13 @@ def _compute_covariance(self):
 
         self.covariance = self._data_covariance * self.factor**2
         self.inv_cov = self._data_inv_cov / self.factor**2
-        self._norm_factor = sqrt(linalg.det(2*pi*self.covariance))
+        itemsize = self.covariance.dtype.itemsize
+        if itemsize in (12, 16):
+            self.log_det = np.abs(np.log(linalg.det(2.0*pi*self.covariance)))
+            self.log_sign = np.log(np.sign(self.log_det))
+        else:
+            sign, self.log_det = np.linalg.slogdet(2.0*pi*self.covariance)
+            self.log_sign = np.log(sign)
 
     def pdf(self, x):
         """
@@ -607,18 +613,20 @@ def logpdf(self, x):
             for i in range(self.n):
                 diff = self.dataset[:, i, newaxis] - points
                 tdiff = dot(self.inv_cov, diff)
-                energy[i] = sum(diff*tdiff, axis=0) / 2.0
-            result = logsumexp(-energy.T,
-                               b=self.weights / self._norm_factor, axis=1)
+                energy[i] = sum(diff*tdiff, axis=0)
+            log_to_sum = 2.0 * np.log(
+                self.weights) - self.log_det - energy.T - self.log_sign
+            result = logsumexp(0.5 * log_to_sum, axis=1)
         else:
             # loop over points
             result = np.empty((m,), dtype=float)
             for i in range(m):
                 diff = self.dataset - points[:, i, newaxis]
                 tdiff = dot(self.inv_cov, diff)
-                energy = sum(diff * tdiff, axis=0) / 2.0
-                result[i] = logsumexp(-energy, b=self.weights /
-                                      self._norm_factor)
+                energy = sum(diff * tdiff, axis=0)
+                log_to_sum = 2.0 * np.log(
+                    self.weights) - self.log_det - energy - self.log_sign
+                result[i] = logsumexp(0.5 * log_to_sum)
 
         return result
 


Code review comments:
diff_hunk: @@ -569,7 +569,13 @@ def _compute_covariance(self):
 
         self.covariance = self._data_covariance * self.factor**2
         self.inv_cov = self._data_inv_cov / self.factor**2
-        self._norm_factor = sqrt(linalg.det(2*pi*self.covariance))
+        itemsize = self.covariance.dtype.itemsize
+        if itemsize in (12, 16):
+            self.log_det = np.abs(np.log(linalg.det(2.0*pi*self.covariance)))
+            self.log_sign = np.log(np.sign(self.log_det))
+        else:
+            sign, self.log_det = np.linalg.slogdet(2.0*pi*self.covariance)
+            self.log_sign = np.log(sign)
path: scipy/stats/kde.py
text: @user3:
For the linear algebra shims, @user1 or @user2 come to mind for a quick check (not my field of expertise)

@user1:
To avoid the typical overflow to inf or underflow to zero, log-det needs to be computed via first getting the cholesky factor of covariance matrix (assuming always pos-def) and then multiplying the diagonal elements of the factor and then taking the log (which is a sum by the log product rule) so:  

```python
L = cholesky(self.covariance*2*pi)
self.log_det = 2*np.log(np.diag(L)).sum()
```

@author:
Awesome, this is a point that I had overlooked, and it makes the code simpler with the same functionality. Added in https://github.com/scipy/scipy/pull/12989/commits/9af9cf3df18fa2e395f9c6b68228eb7468bac9fd.