Posso criar uma matriz multivariate_normal usando dask?

Sep 15 2020

Um pouco relacionado a este post , estou tentando replicar multivariate_normalem dask: Usando numpy, posso criar uma matriz normal multivariada com uma covariância especificada 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

No entanto, preciso de uma matriz significativamente grande com n_dim = 4_500_000e size = 100000. Será caro computar tanto com relação à CPU quanto à memória. Felizmente, tenho acesso a um Cloudera DataScience Workbench Cluster e estava tentando resolver isso 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

Na documentação , não consigo encontrar nenhuma função que pareça fazer o que eu preciso. Alguém conhece uma solução / workarround, possivelmente usando xarrayou qualquer outro módulo que rode em clusters?

Respostas

1 SAFEX Sep 15 2020 at 21:15

Um arround de trabalho, por enquanto, é usar uma decomposição cholesky. Observe que qualquer matriz de covariância C pode ser expressa como C = G * G '. Segue-se então que x = G '* y é correlacionado conforme especificado em C se y for normal padrão (veja esta excelente postagem em StackExchange Mathematic). Em código:

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