Repository: pytorch
Issue: Wishart: problems with singular samples. 
Labels: module: distributions,triaged

### 🐛 Describe the bug

I have two issues with the sampling algorithm in `torch.distributions.Wishart`:
1. It generates a warning about singular samples, even when the samples are nowhere near singular.
2. I believe retrying additional samples to "correct" singular samples is statistically wrong.

Some related discussion can be found in an existing issue [here](https://github.com/pytorch/pytorch/issues/76160).   

### Issue 1: Unnecesary warning
This code demonstrates the unnecessary warning:
```
import torch
from torch import from_numpy as n2t
from torch import distributions as td
from scipy import stats as ss

torch.set_default_dtype(torch.float64)

dim = 4
nu = dim + 10.0
V = torch.eye(dim)
n = 100_000

Wt = td.Wishart(nu, V).sample((n,))
Ws = n2t(ss.wishart(nu, V.numpy()).rvs((n,)))
print(Wt.slogdet().logabsdet.min())
print(Ws.slogdet().logabsdet.min())
```
Which (for example), gives the output:
```
\torch\distributions\wishart.py:271: UserWarning: Singular sample detected.
  warnings.warn("Singular sample detected.")
tensor(5.6044)
tensor(5.4030)
```
As can be seen from the minimum log determinant over all samples, none of the samples are anywhere near singular---for which the log determinant would be at minus infinity.

### Issue 2: Correcting singular samples is incorrect!
Imagine sampling from the standard normal distribution, but you decide you don't like negative samples. Then you reject all negative samples and only use the positive ones. This is a form of rejection sampling and it returns samples from a [completely different disttribution](https://en.wikipedia.org/wiki/Half-normal_distribution). Rejecting singular Wishart samples is similar, you are then sampling from a different, _truncated_ distribution, _not_ the Wishart you wanted to sample from. 

#### Some comments
- The [scipy.stats Wishart documentation](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.wishart.html) advises:
"The algorithm implemented by the rvs method may produce numerically singular matrices with `dim-1 < nu < dim` ; the user may wish to check for this condition and generate replacement samples as necessary." I believe this advice is misleading, for the reason explained above. 

- For a small degree of freedom parameter (`nu` in my code), you can get numerically singular samples, but they are not mathematically singular. The [Bartlett sampling algorithm](https://en.wikipedia.org/wiki/Wishart_distribution#Bartlett_decomposition) generates chi-squared (equivalently [gamma](https://en.wikipedia.org/wiki/Chi-squared_distribution#Gamma,_exponential,_and_related_distributions)) variates and populates with their square roots the diagonal of the Cholesky factor of the Wishart sample matrix. With a small shape parameter, the gamma variates can underflow to zero. A range of gamma shapes is used along the diagonal, the smallest being `(nu+1-dim)/2`. If any of these diagonal elements are numerically zero, the determinant of the whole Wishart sample (the matrix) becomes zero.  Whether this will happen or not, is no mystery: with a bit of work, the probability that a Wishart sample will underflow to a singular matrix can be computed via the CDF of the gamma distribution evaluated at your floating point epsilon. An example of a similar calculation can be found [here](https://en.wikipedia.org/wiki/Gamma_distribution#Caveat_for_small_shape_parameter). 
- As also explained in the last link above, if you really need to work with such small gamma shape parameters (equivalently `dim-1 < nu < dim`), you can generate finite [log-gamma](https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.loggamma.html) variates instead. Using these you can compute a finite value for the log-determinant of your Wishart sample. For each log-gamma sample, `s`, you can do `(s/2).exp()` to recover the `sqrt` chi-squared diagonal elements and proceed with the Bartlett recipe. I suggest the log-determinant can be made an additional, or optional output of the Wishart sampler.  

#### Summary
- The so-called _correction_ of singular samples is statistically incorrect and it samples from the wrong distribution.
- I believe the correct thing to do is just to output the singular samples, without retries.
- A possible extra functionality in the code can compute the log-determinants (using log-gamma samples as explained above) and then optionally output those log-determinants alongside the Wishart samples. (These log determinants are available almost for free.)    
- The user should be forced to think carefully about what they are doing and avoid small `nu` unless they can deal with singular samples downstream.

#### A final suggestion
If the Wishart samples are generated using the Bartlett algorithm, it may be of considerable practical benefit for downstream use of the samples to have an option to instead output the Cholesky factors of the samples, for example, if these samples are to be used in a conjugate-prior recipe as precisions for multivariate normal distributions. 





  

### Versions

```
Collecting environment information...
PyTorch version: 2.2.1
Is debug build: False
CUDA used to build PyTorch: None
ROCM used to build PyTorch: N/A

OS: Microsoft Windows 10 Home Single Language
GCC version: Could not collect
Clang version: Could not collect
CMake version: Could not collect
Libc version: N/A

Python version: 3.10.9 | packaged by conda-forge | (main, Feb  2 2023, 20:14:58) [MSC v.1929 64 bit (AMD64)] (64-bit runtime)
Python platform: Windows-10-10.0.19045-SP0
Is CUDA available: False
CUDA runtime version: No CUDA
CUDA_MODULE_LOADING set to: N/A
GPU models and configuration: No CUDA
Nvidia driver version: No CUDA
cuDNN version: No CUDA
HIP runtime version: N/A
MIOpen runtime version: N/A
Is XNNPACK available: True

CPU:
Architecture=9
CurrentClockSpeed=2594
DeviceID=CPU0
Family=198
L2CacheSize=1024
L2CacheSpeed=
Manufacturer=GenuineIntel
MaxClockSpeed=2594
Name=Intel(R) Core(TM) i7-4720HQ CPU @ 2.60GHz
ProcessorType=3
Revision=15363

Versions of relevant libraries:
[pip3] flake8==6.0.0
[pip3] mypy-extensions==1.0.0
[pip3] numpy==1.24.2
[pip3] numpydoc==1.5.0
[pip3] torch==2.2.1
[conda] blas                      2.106                       mkl    conda-forge
[conda] cpuonly                   2.0                           0    pytorch
[conda] libblas                   3.9.0                     6_mkl    conda-forge
[conda] libcblas                  3.9.0                     6_mkl    conda-forge
[conda] liblapack                 3.9.0                     6_mkl    conda-forge
[conda] liblapacke                3.9.0                     6_mkl    conda-forge
[conda] mkl                       2020.4             hb70f87d_311    conda-forge
[conda] numpy                     1.24.2          py310hd02465a_0    conda-forge
[conda] numpydoc                  1.5.0              pyhd8ed1ab_0    conda-forge
[conda] pytorch                   2.2.1              py3.10_cpu_0    pytorch
[conda] pytorch-mutex             1.0                         cpu    pytorch
```

cc @fritzo @neerajprad @alicanb @nikitaved