[Numpy-discussion] strange behavior of numpy.random.multivariate_normal, ticket:1842

[Charles R Harris]

On Thu, Feb 16, 2012 at 10:07 AM, <josef.pktd at gmail.com> wrote:

> On Thu, Feb 16, 2012 at 11:47 AM,  <josef.pktd at gmail.com> wrote:
> > On Thu, Feb 16, 2012 at 11:30 AM,  <josef.pktd at gmail.com> wrote:
> >> On Thu, Feb 16, 2012 at 11:20 AM, Warren Weckesser
> >> <warren.weckesser at enthought.com> wrote:
> >>>
> >>>
> >>> On Thu, Feb 16, 2012 at 10:12 AM, Pierre Haessig <
> pierre.haessig at crans.org>
> >>> wrote:
> >>>>
> >>>> Le 16/02/2012 16:20, josef.pktd at gmail.com a écrit :
> >>>>
> >>>>> I don't see any way to fix multivariate_normal for this case, except
> >>>>> for dropping svd or for random perturbing a covariance matrix with
> >>>>> multiplicity of singular values.
> >>>>
> >>>> Hi,
> >>>> I just made a quick search in what R guys are doing. It happens there
> are
> >>>> several codes (http://cran.r-project.org/web/views/Multivariate.html). For
> >>>> instance, mvtnorm
> >>>> (http://cran.r-project.org/web/packages/mvtnorm/index.html). I've
> attached
> >>>> the related function from the source code of this package.
> >>>>
> >>>> Interestingly enough, it seems they provide 3 different methods (svd,
> >>>> eigen values, and Cholesky).
> >>>> I don't have the time now to dive in the assessments of pros and cons
> of
> >>>> those three. Maybe one works for our problem, but I didn't check yet.
> >>>>
> >>>> Pierre
> >>>>
> >>>
> >>>
> >>> For some alternatives to numpy's multivariate_normal, see
> >>> http://www.scipy.org/Cookbook/CorrelatedRandomSamples.  Both versions
> >>> (Cholesky and eigh) are just a couple lines of code.
> >>
> >> Thanks both,
> >>
> >> The main point is that it is a "Needs decision"
> >>
> >> Robert argued several times on the mailing list why he chose svd.
> >> (with svd covariance can be closer to singular then with cholesky)
> >>
> >> In statsmodels we usually just use Cholesky for similar
> >> transformation, and I use occasionally an eigh version. (I need to
> >> look up the thread but I got puzzled about results with eig and
> >> multiplicity of eigenvalues before.)
> >>
> >> The R code is GPL, but the few lines of code look standard without any
> >> special provision for non-deterministic linear algebra.
> >>
> >> If multivariate_normal switches from svd to cholesky or eigh, we still
> >> need to check that we don't run into similar "determinacy" problems
> >> with numpy's linalg (I think in statsmodels we use mostly scipy, so I
> >> don't know.)
> >
> > np.linalg.eigh always produces the same eigenvectors, both running
> > repeatedly in the same session and running the script several times on
> > the command line.
> >
> > so eigh looks good as alternative to svd for this case, I don't know
> > if we buy numerical problems in other corner cases, but for near
> > singularity it's always possible to check the smallest eigenvalue
>
> cholesky is also deterministic in my runs
>
> What I would suggest is to use cholesky first, catch the singular
> exception and then use eigh. With eigh we would get perfectly
> correlated random variables.
>
> Again if my reading of linalg comments is correct, cholesky is the
> fastest way to detect singularity of a matrix, and is faster then eigh
> in the non-singular case.
>
> I have no idea if there is an almost singular case, where cholesky
> fails, but the current svd would produce a not perfectly correlated
> random sample (up to numerical precision).
>
> (Alternative, which I don't think I like so much, is to use a small
> Ridge correction (multiply diagonal by 1 + x*nulp ? This would bound
> it away from perfect correlation, I guess.)
>
>
Cholesky doesn't do any reordering of the matrix, but proceeds downward
factoring row by row, so to speak. It is Gauss elimination without row
pivoting and is only stable when the matrix is positive definite.
Fortunately, failure of positive definiteness shows up an attempt to take
the real square root of a negative number, so is detected.

The problem with svd is that the singular values are always non-negative,
hence the resulting factorization isn't always of the form R^T * R, which
is easily understood because that form is necessarily non-negative definite.

Chuck
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20120216/06d5d8f6/attachment.html>