Memory error in Quaternion.average

[enritoomey]

When calling average method with 37K quaternion list, the following error appears:

$len(q_yp2ym)
37743

$ q_yp2ym[0]
GeneralQuaternion(0.34208732518700213, 0.0030042303823605808, -0.71948547865943235, -0.60440705037616804)

$ Quaternion.average_and_std_naive(*q_yp2ym)
---------------------------------------------------------------------------
MemoryError                               Traceback (most recent call last)
<ipython-input-47-c7868c5f1e09> in <module>()
----> 1 Quaternion.average_and_std_naive(*q_yp2ym)

~/pyvenvs/calibration/lib/python3.5/site-packages/quaternions/quaternion.py in average_and_std_naive(weights, *quaternions)
    236         if weights is None:
    237             weights = np.ones(len(quaternions))
--> 238         q_average = Quaternion.average(*quaternions, weights=weights)
    239         diffs = np.array([q_average.distance(quat) for quat in quaternions])
    240         stddev = np.degrees(np.sqrt(sum((diffs ** 2) * weights) / np.sum(weights)))

~/pyvenvs/calibration/lib/python3.5/site-packages/quaternions/quaternion.py in average(cls, weights, *quaternions)
    110         if weights is None:
    111             weights = np.ones(len(quaternions))
--> 112         m = b.T.dot(np.diag(weights)).dot(b)
    113 
    114         return cls._first_eigenvector(m)

~/pyvenvs/calibration/lib/python3.5/site-packages/numpy/lib/twodim_base.py in diag(v, k)
    247     if len(s) == 1:
    248         n = s[0]+abs(k)
--> 249         res = zeros((n, n), v.dtype)
    250         if k >= 0:
    251             i = k

MemoryError: 

[matiasg]

The problem seems to be generating a huge diagonal matrix just to multiply. There are alternative ways (taken from https://stackoverflow.com/questions/18522216/multiplying-across-in-a-numpy-array):

In [39]: %timeit -n2 np.einsum('ij,i->ij', a, d)
217 µs ± 21.9 µs per loop (mean ± std. dev. of 7 runs, 2 loops each)

In [40]: %timeit -n2 np.array([d]).T * a
248 µs ± 45.3 µs per loop (mean ± std. dev. of 7 runs, 2 loops each)

In [41]: %timeit -n2 np.diag(d).dot(a)
2.36 s ± 26.6 ms per loop (mean ± std. dev. of 7 runs, 2 loops each)

In [42]: a.shape
Out[42]: (30102, 4)

In [43]: d.shape
Out[43]: (30102,)

I will implement the first version with Einstein sums. It also works better for smaller matrices:

In [49]: %timeit -n20 np.einsum('ij,i->ij', a, d)
1.88 µs ± 797 ns per loop (mean ± std. dev. of 7 runs, 20 loops each)

In [50]: %timeit -n20 np.array([d]).T * a
3.94 µs ± 317 ns per loop (mean ± std. dev. of 7 runs, 20 loops each)

In [51]: %timeit -n20 np.diag(d).dot(a)
5.69 µs ± 2.16 µs per loop (mean ± std. dev. of 7 runs, 20 loops each)

In [52]: a.shape
Out[52]: (30, 4)

[enritoomey]

Excelent!