Repository navigation
Floating point rounding in Cholesky #13658
Description
Reproducing code example:
So I need correlated random variables for my Monte Carlo simulations. I use Cholesky method from numpy package to do so.
I generate them using this code. Since I am doing stock simulations, the specific matrix I am passing is important since it means that portion of the stocks is fully correlated as I want to test my strategy in that case.
import numpy as np
from scipy.linalg import eigh, cholesky
from scipy.stats import norm
# import numba
def cov2corr(cov):
cov = np.asanyarray(cov)
std_ = np.sqrt(np.diag(cov))
corr = cov / np.outer(std_, std_)
return corr
def correlated_random(correlation_matrix, method='choelsky'):
r = cov2corr(correlation_matrix)
sz = correlation_matrix.shape[0]
num_samples = 1
x = norm.rvs(size=(sz, num_samples))
if method == 'choelsky':
c = np.linalg.cholesky(r)
elif method == 'eigen':
evals, evecs = eigh(r)
c = np.dot(evecs, np.diag(np.sqrt(evals)))
y = np.dot(c, x)
return y.reshape((num_samples, sz))
correlated_random(
correlation_matrix=np.array([
[1.,1.,1.,1.,1.,1.,1.,0.,0.,0.],
[1.,1.,1.,1.,1.,1.,1.,0.,0.,0.],
[1.,1.,1.,1.,1.,1.,1.,0.,0.,0.],
[1.,1.,1.,1.,1.,1.,1.,0.,0.,0.],
[1.,1.,1.,1.,1.,1.,1.,0.,0.,0.],
[1.,1.,1.,1.,1.,1.,1.,0.,0.,0.],
[1.,1.,1.,1.,1.,1.,1.,0.,0.,0.],
[0.,0.,0.,0.,0.,0.,0.,1.,0.,0.],
[0.,0.,0.,0.,0.,0.,0.,0.,1.,0.],
[0.,0.,0.,0.,0.,0.,0.,0.,0.,1.]
])
)
When I run this I get error that matrix is not positive definite. If I run the following snippet to get eigenvalues, I get:
np.linalg.eigvals(np.array([
[1.,1.,1.,1.,1.,1.,1.,0.,0.,0.],
[1.,1.,1.,1.,1.,1.,1.,0.,0.,0.],
[1.,1.,1.,1.,1.,1.,1.,0.,0.,0.],
[1.,1.,1.,1.,1.,1.,1.,0.,0.,0.],
[1.,1.,1.,1.,1.,1.,1.,0.,0.,0.],
[1.,1.,1.,1.,1.,1.,1.,0.,0.,0.],
[1.,1.,1.,1.,1.,1.,1.,0.,0.,0.],
[0.,0.,0.,0.,0.,0.,0.,1.,0.,0.],
[0.,0.,0.,0.,0.,0.,0.,0.,1.,0.],
[0.,0.,0.,0.,0.,0.,0.,0.,0.,1.]
]))
Out[15]:
array([ 0.00000000e+00, 7.00000000e+00, 2.46519033e-32, -2.71947991e-16,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 1.00000000e+00,
1.00000000e+00, 1.00000000e+00])
We can see that almost all values are positive, except the 3rd one. Since it is really small I suspect that it is caused by floating point error. Is there any way I can circumvent this without writing custom procedure for Cholesky factorization.
Error message:
---------------------------------------------------------------------------
LinAlgError Traceback (most recent call last)
<ipython-input-194-cc640eb5f0a5> in <module>
99
100 np.random.seed(seed)
--> 101 series = simulate_portfolio(40, p)
102 if not np.any(np.isnan(series)):
103 if p == 0:
<ipython-input-194-cc640eb5f0a5> in simulate_portfolio(asset_N, correlation_percent, total_rtrn, sigma, N)
87 leverages=np.array([1 for _ in range(asset_N)]),
88 correlation_matrix=correlation_matrix,
---> 89 starting_values=np.array([1 for _ in range(asset_N)])
90 )
91
<ipython-input-194-cc640eb5f0a5> in simulate_multiple_series(total_rtrns, sigmas, N, leverages, correlation_matrix, starting_values)
50
51 for n in range(N):
---> 52 volatilities = correlated_random(correlation_matrix, method='choelsky')
53
54 prev_values = values[-1, :]
<ipython-input-194-cc640eb5f0a5> in correlated_random(correlation_matrix, method)
35
36 if method == 'choelsky':
---> 37 c = np.linalg.cholesky(r)
38 elif method == 'eigen':
39 evals, evecs = eigh(r)
~/workspace/jupyter/env/lib/python3.7/site-packages/numpy/linalg/linalg.py in cholesky(a)
757 t, result_t = _commonType(a)
758 signature = 'D->D' if isComplexType(t) else 'd->d'
--> 759 r = gufunc(a, signature=signature, extobj=extobj)
760 return wrap(r.astype(result_t, copy=False))
761
~/workspace/jupyter/env/lib/python3.7/site-packages/numpy/linalg/linalg.py in _raise_linalgerror_nonposdef(err, flag)
98
99 def _raise_linalgerror_nonposdef(err, flag):
--> 100 raise LinAlgError("Matrix is not positive definite")
101
102 def _raise_linalgerror_eigenvalues_nonconvergence(err, flag):
LinAlgError: Matrix is not positive definite
Numpy/Python version information:
1.16.2 3.7.2 (default, Feb 12 2019, 08:15:36)
[Clang 10.0.0 (clang-1000.11.45.5)]
Activity
tms1337 commentedon May 29, 2019
I solved the issue by setting
np.fill_diagonal(a, 1.01)so that diagonal is1.01instead of1, but still this is hacky perhaps there should be some way to control precision.seberg commentedon May 29, 2019
Without digging into it too much, these type of floating point errors need to be expected, there maybe symbolic libraries, such as sympy, that can give you an exact result if you really need it. Note that the errors are tiny with regards to the other values.
What I am not sure about is whether there is a recipe similar to an
rcondinpinvto decide a cutoff.tms1337 commentedon May 29, 2019
Yes would be nice to have something like
rcondewmoore commentedon May 29, 2019
choleskyuses lapack'sdpotrfto compute the decomposition. There isn't AFAIK anything likercondinpinvfor this. The issue is that within the numerical precision available, your matrix isn't positive definite. So there isn't a cholesky decomposition for it. Adding a little to the diagonal is probably a fine way to make this more robust. I don't think there is anything else without going to a different decomposition.tms1337 commentedon May 29, 2019
I think it is since all minor determinants are
D >= 0But I agree that the solution with the diagonal is not so bad, so if there is no way of adding something similar to
rcondinpinv, then I can close the issue.pmli commentedon May 30, 2019
@tms1337 As far as I see, you probably want to use
numpy.random.multivariate_normal, since it supports positive semidefinite covariance matrices.Maybe to make it clear, the covariance matrix you have is positive semidefinite and not positive definite, since it has zero eigenvalues. In particular, it is of rank 4, which you can see from
numpy.random.multivariate_normaluses SVD and there are suggestions to use pivoted Cholesky decomposition (see #13386, #3938).