Kann ich mit dask eine multivariate_normal-Matrix erstellen?

Question

Etwas im Zusammenhang mit diesem Beitrag versuche ich zu replizieren multivariate_normalin dask: Mit numpy kann ich eine multivariate normale Matrix mit einer bestimmten Kovarianz erstellen, indem ich:

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

Ich brauche jedoch eine deutlich große Matrix mit n_dim = 4_500_000und size = 100000. Die Berechnung sowohl für die CPU als auch für den Speicher ist teuer. Glücklicherweise habe ich Zugriff auf einen Cloudera DataScience Workbench-Cluster und habe versucht, dieses Problem zu lösen, indem ich 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

In der Dokumentation kann ich keine Funktion finden, die das zu tun scheint, was ich brauche. Kennt jemand eine Lösung / Problemumgehung, die möglicherweise xarrayein anderes Modul verwendet oder auf Clustern ausgeführt wird?

Answer 1

Eine Arbeit im Moment ist die Verwendung einer Cholesky-Zerlegung. Es ist zu beachten, dass jede Kovarianzmatrix C ausgedrückt werden kann als C = G * G '. Daraus folgt, dass x = G '* y wie in C angegeben korreliert, wenn y Standardnormal ist (siehe diesen ausgezeichneten Beitrag zu StackExchange Mathematic). In Code:

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()