Repository navigation
Deterministic (up to numerical error) multivariate_normal sampling #13386
Description
Currently, numpy.random.multivariate_normal makes use of SVD decomposition, numpy.dual.svd, on the provided covariance matrix to perform the sampling. This approach allows the use of positive semi-definite covariance matrix. However, it is not unique as the singular vectors are not guaranteed to be. In fact,
[[-1, 0],
[0, -1]]
and
[[1, 0],
[0, 1]]
can be the singular-vectors of a 2-by-2 identity matrix.
A Cholesky decomposition instead would result in an unique decomposition with the down side of only working on positive definite matrices. We could use instead a Cholesky decomposition with complete pivoting: http://www.netlib.org/lapack/explore-html/da/dba/group__double_o_t_h_e_rcomputational_ga31cdc13a7f4ad687f4aefebff870e1cc.html
It provides both advantages: works on positive semi-definite matrices and are unique (as far as I know).
Any thoughts?
Related issue: #13358
Activity
rkern commentedon Apr 22, 2019
See #3938 for previous discussion of this suggestion. There are several other discussions that you can find by searching for
multivariate_normal cholesky.rkern commentedon Apr 22, 2019
Though I see that the only remaining issue is stream compatibility, which will be manageable after the NEP 19 implementation PR gets in.
mattip commentedon May 30, 2019
@rkern can we close this as a duplicate?
mattip commentedon Jun 6, 2019
Closing, please reopen if I misunderstood.
horta commentedon Sep 13, 2019
It is not duplication because this thread tackles the problem of multiple decompositions (SVD) versus unique decomposition (pivoted Cholesky) to semi-definite positive matrices.
Rounding error is not within the scope of this discussion.