¿Puedo crear una matriz multivariate_normal usando dask?

Question

Algo relacionado con esta publicación , estoy tratando de replicar multivariate_normalen dask: Usando numpy puedo crear una matriz normal multivariante con una covarianza específica usando:

import numpy as np
n_dim = 5
size = 300
A = np.random.randn(n_dim, n_dim) # a matrix
covm = A.dot(A.T) # A*A^T is positive semi-definite, as a covariance matrix
x = np.random.multivariate_normal(size=300, mean=np.zeros(len(covm)),cov=covm) # generate data

Sin embargo, necesito una matriz significativamente grande con n_dim = 4_500_000y size = 100000. Esto será costoso de calcular tanto con respecto a la CPU como a la memoria. Afortunadamente, tengo acceso a Cloudera DataScience Workbench Cluster y estaba tratando de resolver esto usando dask:

import dask.array as da
n_dim = 4_500_000
size = 100000
A = da.random.standard_normal((n_dim, n_dim))  
covm = A.dot(A.T)
#x = da.random.multivariate_normal(size=300, mean=np.zeros(len(covm)),cov=covm) # generate data

En la documentación , no puedo encontrar ninguna función que parezca hacer lo que necesito. ¿Alguien conoce una solución / solución alternativa, posiblemente usando xarrayo cualquier otro módulo que se ejecute en clústeres?

Answer 1

Un trabajo alrededor, por ahora, es usar una descomposición cholesky. Tenga en cuenta que cualquier matriz de covarianza C se puede expresar como C = G * G '. Luego se deduce que x = G '* y está correlacionado como se especifica en C si y es normal estándar (vea esta excelente publicación en StackExchange Mathematic). En codigo:

Numpy

n_dim =4
size = 100000
A = np.random.randn(n_dim, n_dim)
covm = A.dot(A.T)

x=  np.random.multivariate_normal(size=size, mean=np.zeros(len(covm)),cov=covm)
## verify numpys covariance is correct
np.cov(x, rowvar=False)
covm

Dask

## create covariance matrix
A = da.random.standard_normal(size=(n_dim, n_dim),chunks=(2,2))
covm = A.dot(A.T)

## get cholesky decomp
L = da.linalg.cholesky(covm, lower=True)

## drawn standard normal 
sn= da.random.standard_normal(size=(size, n_dim),chunks=(100,100))

## correct for correlation
x =L.dot(sn.T)
x.shape

## verify
covm.compute()
da.cov(x, rowvar=True).compute()