5
\$\begingroup\$

I need to calculate a Bayesian likelihood in its negative logarithmic form:

$$-\ln L = N\ln M -\sum\limits_{i=1}^{N} \ln\left\{P_i \sum\limits_{j=1}^M \frac{\exp\left[ -\frac{1}{2}\left(\sum_{k=1}^D \frac{(f_{ik}-g_{jk})^2}{\sigma_{ik}^2 + \rho_{jk}^2}\right) \right]}{\prod_{k=1}^D \sqrt{\sigma_{ik}^2 + \rho_{jk}^2}} \right\}$$

where these $$P_i\;\;;\;\; f_{ik} \;\;;\;\; \sigma_{ik}^2\;\;;\;\; g_{jk} \;\;;\;\; \rho_{jk}^2$$ are all known floats. This is my code:

import numpy as np

# Define random data with correct shapes.
N, M = np.random.randint(10, 1000), np.random.randint(10, 1000)
D = np.random.randint(2, 5)
f_N, g_M = np.random.rand(N, D, 2), np.random.rand(M, D, 2)
P_i = np.random.uniform(0., 1., N)

# N summatory.    
# Store here each 'i' element from the M summatory.
sum_Mi = []
for f_i in f_N:

    # M summatory.
    sum_M_j = 0.
    for g_j in g_M:

        # Obtain ijk element inside the M summatory
        D_sum, sigma_ijk = 0., 1.
        for k, f_ik in enumerate(f_i):
            # Where: f_ik[0] = \f_{ik}       ; g_j[k][0] = \g_{jk}
            # and:   f_ik[1] = \sigma_{ik}^2 ; g_j[k][1] = \rho_{jk}^2
            D_sum += np.square(f_ik[0] - g_j[k][0]) / (f_ik[1] + g_j[k][1])
            sigma_ijk = sigma_ijk * (f_ik[1] + g_j[k][1])

        sum_M_j += np.exp(-0.5 * D_sum) / np.sqrt(sigma_ijk)

    sum_Mi.append(sum_M_j)

sum_Mi_Pi = P_i * sum_Mi
# Remove 0. elements before applying the logarithm below.
sum_N = sum_Mi_Pi[sum_Mi_Pi != 0.]

# Final negative logarithmic likelihood
lik = len(f_N) * np.log(len(g_M)) - sum(np.log(sum_N))
print(lik)

I'm after an optimization done at the code level, ie: not using multi-threading.

I've tried numpy broadcasting on the inner loop

# Obtain ijk element inside the M summatory
D_sum = sum(np.square((f_i[:, 0] - g_j[:, 0])) / (f_i[:, 1] + g_j[:, 1]))
sigma_ijk = np.prod(f_i[:, 1] + g_j[:, 1])

but that actually takes a longer time to complete than using the for loop!

I also tried using the logsumexp function from scipy combined with the broadcasted arrays:

# Store here each 'i' element from the M summatory.
sum_N = 0.
# N summatory.
for f_i in f_N:

    # M summatory.
    A, B = [], []
    for g_j in g_M:
        # Obtain ijk element inside the M summatory
        A.append(sum(
            np.square((f_i[:, 0] - g_j[:, 0])) / (f_i[:, 1] + g_j[:, 1])))
        B.append(np.prod(f_i[:, 1] + g_j[:, 1]))

    sum_N += logsumexp(-.5 * np.array(A), b=1. / np.sqrt(B))

# Final negative logarithmic likelihood
lik = len(f_N) * np.log(len(g_M)) - sum(np.log(P_i[P_i != 0.])) - sum_N
print(lik)

This resulted in a marginal improvement over just using broadcasting, but still slower than the original code.

\$\endgroup\$
9
  • \$\begingroup\$ There also seems to be a way to calculate the \$\sigma\$ from the \$f\$, which you do not show. \$\endgroup\$ – Graipher Mar 17 '17 at 18:19
  • \$\begingroup\$ $f_{ik}$ is not equal to $f_{jk}$, the notation implies it but this is not true, Ill fix that. The $\sigma$s are completely independent of the $f$s. (Why is my LateX not rendering?) \$\endgroup\$ – Gabriel Mar 17 '17 at 18:34
  • 1
    \$\begingroup\$ You need to add a \ in front of your $. \$\endgroup\$ – Graipher Mar 17 '17 at 18:42
  • \$\begingroup\$ Thank you @Graipher, I've fixed the notation above. Please let mw know if something doesn't make sense to you. (testing \$\alpha\$) \$\endgroup\$ – Gabriel Mar 17 '17 at 18:43
  • 2
    \$\begingroup\$ Vectorizing the inner loop doesn't help because D is so small compared to N and M, e.g. 3 v 500 and 500. It's just a few loops of a costly expression. \$\endgroup\$ – hpaulj Mar 17 '17 at 21:27
1
\$\begingroup\$

This is much faster. I added the seed to better compare results across runs.

import numpy as np

# Define random data with correct shapes.
np.random.seed(0)
N, M = np.random.randint(10, 1000), np.random.randint(10, 1000)
D = np.random.randint(2, 5)
f_N, g_M = np.random.rand(N, D, 2), np.random.rand(M, D, 2)
P_i = np.random.uniform(0., 1., N)

# replace loops with broadcasted array operation, followed by sums    
fgsum = f_N[:,None,:,1] + g_M[None,:,:,1]
fgdif = f_N[:,None,:,0] - g_M[None,:,:,0]
print(fgsum.shape)   # e.g. (694, 569, 3)

Dsum = (fgdif**2 / fgsum).sum(axis=-1)
sigma_ijk = np.prod(fgsum, axis=-1)
sum_M_j = np.exp(-0.5 * Dsum) / np.sqrt(sigma_ijk)
sum_Mi = np.sum(sum_M_j, axis=-1)

# as before
sum_Mi_Pi = P_i * sum_Mi
# Remove 0. elements before applying the logarithm below.
sum_N = sum_Mi_Pi[sum_Mi_Pi != 0.]

# Final negative logarithmic likelihood
lik = len(f_N) * np.log(len(g_M)) - sum(np.log(sum_N))
print(lik)
\$\endgroup\$
6
  • \$\begingroup\$ Excellent! This code is between 80-100 times faster in my machine, thank you very much hpaulj! \$\endgroup\$ – Gabriel Mar 17 '17 at 21:54
  • 1
    \$\begingroup\$ You can get another few percent faster by replacing sum(np.log(sum_N)) with np.log(sum_N).sum() or np.sum(np.log(sum_N)). \$\endgroup\$ – Graipher Mar 18 '17 at 9:59
  • 1
    \$\begingroup\$ okay ... nice that this is faster, but ... WHY??? \$\endgroup\$ – Vogel612 Mar 18 '17 at 16:41
  • \$\begingroup\$ Because it replaces Python loops with compiled numpy code. Sorry. I normally answer on SO where this kind of numpy 'vectorization' is an everyday question. \$\endgroup\$ – hpaulj Mar 18 '17 at 16:48
  • 1
    \$\begingroup\$ More about increasing performance with numpy: codereview.stackexchange.com/questions/17702/… \$\endgroup\$ – hpaulj Mar 18 '17 at 16:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.