<p>@unutbu的答案对于将任何函数应用于数组的行非常有效。
在这种特殊情况下,如果使用大型数组,可以使用一些数学对称性来显著加快速度。</p>
<p>以下是函数的修改版本:</p>
<pre><code>def mahalanobis_sqdist3(x, mean, Sigma):
Sigma_inv = np.linalg.inv(Sigma)
xdiff = x - mean
return (xdiff.dot(Sigma_inv)*xdiff).sum(axis=-1)
</code></pre>
<p>如果您最终使用任何类型的大型<code>Sigma</code>,我建议您缓存<code>Sigma_inv</code>,并将其作为参数传递给函数。
因为在本例中是4x4,所以这无关紧要。
我将展示如何处理大型的<code>Sigma</code>不管是谁遇到这个。</p>
<p>如果您不打算重复使用同一个<code>Sigma</code>,您将无法缓存它,因此,您可以使用不同的方法来求解线性系统,而不是反转矩阵。
在这里,我将使用内置到SciPy中的LU分解。
这只会在<code>x</code>的列数相对于其行数较大时提高时间。</p>
<p>下面是一个函数,它显示了这种方法:</p>
<pre><code>from scipy.linalg import lu_factor, lu_solve
def mahalanobis_sqdist4(x, mean, Sigma):
xdiff = x - mean
Sigma_inv = lu_factor(Sigma)
return (xdiff.T*lu_solve(Sigma_inv, xdiff.T)).sum(axis=0)
</code></pre>
<p>这里有一些时间安排。
我将把另一个答案中提到的带有<code>einsum</code>的版本包括在内。</p>
<pre><code>import numpy as np
Sig1 = np.array([[ 0.16043333, 0.11808333, 0.02408333, 0.01943333],
[ 0.11808333, 0.13583333, 0.00625 , 0.02225 ],
[ 0.02408333, 0.00625 , 0.03916667, 0.00658333],
[ 0.01943333, 0.02225 , 0.00658333, 0.01093333]])
mean1 = np.array([ 5.028, 3.48 , 1.46 , 0.248])
x = np.random.rand(25, 4)
%timeit np.apply_along_axis(mahalanobis_sqdist, 1, x, mean1, Sig1)
%timeit mahalanobis_sqdist2(x, mean1, Sig1)
%timeit mahalanobis_sqdist3(x, mean1, Sig1)
%timeit mahalanobis_sqdist4(x, mean1, Sig1)
</code></pre>
<p>给予:</p>
<pre><code>1000 loops, best of 3: 973 µs per loop
10000 loops, best of 3: 36.2 µs per loop
10000 loops, best of 3: 40.8 µs per loop
10000 loops, best of 3: 83.2 µs per loop
</code></pre>
<p>但是,更改所涉及阵列的大小会更改计时结果。
例如,让<code>x = np.random.rand(2500, 4)</code>,计时如下:</p>
<pre><code>10 loops, best of 3: 95 ms per loop
1000 loops, best of 3: 355 µs per loop
10000 loops, best of 3: 131 µs per loop
1000 loops, best of 3: 337 µs per loop
</code></pre>
<p>让<code>x = np.random.rand(1000, 1000)</code>、<code>Sigma1 = np.random.rand(1000, 1000)</code>和<code>mean1 = np.random.rand(1000)</code>计时如下:</p>
<pre><code>1 loops, best of 3: 1min 24s per loop
1 loops, best of 3: 2.39 s per loop
10 loops, best of 3: 155 ms per loop
10 loops, best of 3: 99.9 ms per loop
</code></pre>
<p><strong>编辑</strong>:我注意到其他答案之一使用了Cholesky分解。
假设<code>Sigma</code>是对称的正定的,我们实际上可以做得比上面的结果更好。
通过SciPy,BLAS和LAPACK提供了一些很好的例程,可以处理对称正定矩阵。
这里有两个更快的版本。</p>
<pre><code>from scipy.linalg.fblas import dsymm
def mahalanobis_sqdist5(x, mean, Sigma_inv):
xdiff = x - mean
Sigma_inv = la.inv(Sigma)
return np.einsum('...i,...i->...',dsymm(1., Sigma_inv, xdiff.T).T, xdiff)
from scipy.linalg.flapack import dposv
def mahalanobis_sqdist6(x, mean, Sigma):
xdiff = x - mean
return np.einsum('...i,...i->...', xdiff, dposv(Sigma, xdiff.T)[1].T)
</code></pre>
<p>第一个仍然颠倒西格玛。
如果预先计算并重用逆运算,则速度会快得多(在我的机器上,使用预先计算的逆运算,1000x1000的情况需要35.6ms)。
我还使用了einsum来获取产品,然后沿着最后一个轴求和。
结果,这比做<code>(A * B).sum(axis=-1)</code>之类的事情快得多。
这两个功能提供以下计时:</p>
<p>第一个测试用例:</p>
<pre><code>10000 loops, best of 3: 55.3 µs per loop
100000 loops, best of 3: 14.2 µs per loop
</code></pre>
<p>第二个测试用例:</p>
<pre><code>10000 loops, best of 3: 121 µs per loop
10000 loops, best of 3: 79 µs per loop
</code></pre>
<p>第三个测试用例:</p>
<pre><code>10 loops, best of 3: 92.5 ms per loop
10 loops, best of 3: 48.2 ms per loop
</code></pre>