bench_gpr_std_prediction.py

import numpy as np

n, p = 10000, 200
a = np.random.rand(n, p)
b = np.random.rand(p, p)

%timeit np.einsum('ki,kj,ij->k', a, a, b)
# 1 loop, best of 3: 728 ms per loop

% timeit np.einsum("ij,ij->i", np.dot(a, b), a)
# The slowest run took 6.79 times longer than the fastest.
# 100 loops, best of 3: 8.08 ms per loop

# Both call give the same result of course:
np.allclose(np.einsum('ki,kj,ij->k', a, a, b), np.einsum("ij,ij->i", np.dot(a, b), a))

n, p = 100000, 2000
a = np.random.rand(n, p)
b = np.random.rand(p, p)

# WARNING: very slow !
%timeit np.einsum('ki,kj,ij->k', a, a, b)
# 1 loop, best of 3: 13min 39s per loop

% timeit np.einsum("ij,ij->i", np.dot(a, b), a)
# 1 loop, best of 3: 9.57 s per loop