Possible optimization in cython: numpy array

Question

Below is my Cython code for drawing from a multidimensional normal distribution. I use a loop because every time I have a different density. (conLSigma - Cholesky factor)

This takes a lot of time because I take the reverse and Cholesky decomposition for each cycle. This is faster than pure python code, but I was wondering if there is a way to increase speed.

from __future__ import division

import numpy as np 

cimport numpy as np 

ctypedef np.float64_t dtype_t

cimport cython
@cython.boundscheck(False)
@cython.wraparound(False)

def drawMetro(np.ndarray[dtype_t, ndim = 2] beta,
              np.ndarray[dtype_t, ndim = 3] H,
              np.ndarray[dtype_t, ndim = 2] Sigma,
              float s):

    cdef int ncons = betas.shape[0]
    cdef int nX = betas.shape[1]
    cdef int con

    cdef np.ndarray betas_cand = np.zeros([ncons, nX], dtype = np.float64)
    cdef np.ndarray conLSigma = np.zeros([nX, nX], dtype = np.float64)

    for con in xrange(ncons):
        conLSigma = np.linalg.cholesky(np.linalg.inv(H[con] + Sigma))
        betas_cand[con] = betas[con] + s * np.dot(conLSigma, np.random.standard_normal(size = nX))

    return(betas_cand)

Answer 1

Cholesky decomposition creates a lower triangular matrix. This means that almost half of the multiplications performed in np.dotare not required. If you change the line

betas_cand[con] = betas[con] + s * np.dot(conLSigma, np.random.standard_normal(size = nX))

at

tmp = np.random.standard_normal(size = nX)
for i in xrange(nX):
    for j in xrange(i+1):
        betas_cand[con,i] += s * conLSigma[i,j] * tmp[j]

However, you will also need to change

cdef np.ndarray betas_cand = np.zeros([ncons, nX], dtype = np.float64)

at

cdef np.ndarray betas_cand = np.array(betas)

, , , , . , , . , , .

Answer 2

. , linalg.cholesky(linalg.inv(S)).