Use Cholesky decomposition in random.multivariate_normal()

In numpy.random.multivariate_normal(), decomposition of the covariance matrix is currently accomplished using SVD. As discussed in #1821 (without resolution), this is problematic since if a user passes a non-positive definite covariance matrix, the covariance of the resulting random variables may diverge considerably from what's expected, without the user receiving any warning. Currently the only mention of this issue is at the end of the docstring, which many users are unlikely to see. Other packages invariably seem to use a decomposition method that raises an exception if a non-positive definite matrix is provided, and this seems like sensible behavior.

This patch changes the method's default behavior to use Cholesky decomposition, which raises an exception when a non-positive definite matrix is passed. Unless the size of the requested array is very large, users are unlikely to perceive any difference in speed, and can always explicitly select svd with the newly added 'decomposition' argument.