Tukukkk

I need to run this code thousands of time, sequentially:

def update_posterior(y, x, prior_mu, prior_V, prior_a, prior_b):



    # Building blocks used to keep following calculation cleaner

    prior_cov_inverse = np.linalg.inv(prior_V)

    x_transpose = x.transpose()

    n = len(y)

    residuals = y - np.dot(x, prior_mu.transpose())



    # Calculation of posterior parameters

    V_posterior = np.linalg.inv((prior_cov_inverse + x_transpose * x))

    mu_posterior = (V_posterior * (prior_cov_inverse * prior_mu.transpose() + x_transpose * y)).transpose()

    a_posterior = prior_a + n/2

    b_posterior = np.asscalar(prior_b + (residuals.transpose() * np.linalg.inv((np.identity(n) + x * prior_V * x_transpose)) * residuals)/2)



    return mu_posterior, V_posterior, a_posterior, b_posterior

The way it works is that the output of the function is fed back into it, with mu_posterior becoming prior_mu, V_posterior becoming prior_V, a_posterior becoming prior_a, and b_posterior becoming prior_b for the next call. The y and x are different in each call.

This function is horrendously slow--it takes me about 8 seconds to run. This is because of the scale. I have ~5000 parameters, so prior_mu is (1, 5000), prior_V is (5000,5000) and symmetric positive-definite, and prior_a, and prior_b are scalars. y is a scalar, and x is (1, 5000).

Here is the time breakdown by line:

3.75s: prior_cov_inverse = np.linalg.inv(prior_V)

3.86s: V_posterior = np.linalg.inv((prior_cov_inverse + x_transpose * x))

0.13s: b_posterior = np.asscalar(prior_b + (residuals.transpose() * np.linalg.inv((np.identity(n) + x * prior_V * x_transpose)) * residuals)/2)

Any ideas how I can speed it up? I tried to use a Cholesky decomposition but it's even slower?! My guess is there is a more efficient way to implement Cholesky decomposition in Python.

prior_cov_inverse2 = np.linalg.inv(np.linalg.cholesky(prior_V))

prior_cov_inverse2 = np.dot(prior_cov_inverse2.transpose(), prior_cov_inverse2)

EDIT: Here are some sample data to illustrate the issue....

import numpy as np



prior_mu = np.asmatrix(np.full((1, 5040), 5))

prior_V = np.diagflat(np.asmatrix(np.full((1, 5040), 30))) #usually not diagonal, but always symmetric positive definitive

a = 2

b = 2

y = np.asmatrix([10])

x = np.asmatrix(np.concatenate(([1], np.zeros(5039))))



print(update_posterior(y, x, prior_mu, prior_V, a, b))

EDIT II:

I have been able to get this down from ~8s/run to ~1.4s/run by removing matrix inversions in favor of solves and also using the Sherman Morrison formula. Here is my current code. If anyone has any ideas on how to speed it up further, please share! :)

def update_posterior(y, x, prior_mu, prior_V, prior_a, prior_b, I):



    # Building blocks used to keep following calculation cleaner

    x_transpose = x.transpose()

    n = len(y)

    residuals = y - np.dot(x, prior_mu.transpose())



    # Calculation of posterior parameters

    # Below is equivalent to np.linalg.inv(prior_V_inverse + np.dot(x_transpose, x)) but significantly faster    

    V_posterior = prior_V - np.true_divide(np.linalg.multi_dot((prior_V, x_transpose, x, prior_V)), 1 + np.matmul(np.matmul(x, prior_V), x_transpose))



    # Below is equivalent to mu_posterior = np.dot(V_posterior, (np.matmul(prior_V_inverse, prior_mu.transpose()) + np.matmul(x_transpose, y))).transpose() but significantly faster

    mu_posterior = np.dot(V_posterior, np.linalg.solve(prior_V, prior_mu.transpose()) + np.matmul(x_transpose, y)).transpose()

    a_posterior = prior_a + n/2

    b_posterior = np.asscalar(prior_b + (np.matmul(np.matmul(residuals.transpose(), np.linalg.inv((np.identity(n) + np.matmul(np.matmul(x, prior_V), x_transpose)))), residuals))/2)



    return mu_posterior, V_posterior, a_posterior, b_posterior

Stability-wise, it's almost always better to write solve(A, unit_matrix) instead of inv(A). That won't help with performance though.

Performance of linear algebra here is almost surely fixed by the underlying LAPACK library. Stock ATLAS is likely slowest, OpenBLAS or MKL are better, sometimes much better.

However, I'm pretty sure the main improvement here is indeed algorithmic. First, for PSD matrices Cholecky (cholesky/cho_solve) should be better. Second, you seem to be doing a rank-one update (x.T @x), and that can typically be implemented in N**2 operations via some variant of the Shermann-Morrison formula, instead of N**3 for direct inversion.