为什么np.linalg.norm(..，axis=1)比写出向量范数的公式慢

Question

要将矩阵X的行标准化为单位长度，我通常使用：

X /= np.linalg.norm(X, axis=1, keepdims=True)

试图为算法优化此操作时，我非常惊讶地发现，在我的机器上写出规范化大约快40%：

X /= np.sqrt(X[:,0]**2+X[:,1]**2+X[:,2]**2)[:,np.newaxis]
X /= np.sqrt(sum(X[:,i]**2 for i in range(X.shape[1])))[:,np.newaxis]

怎么回事？np.linalg.norm()的性能损失在哪里？

import numpy as np
X = np.random.randn(10000,3)
%timeit X/np.linalg.norm(X,axis=1, keepdims=True)
# 276 µs ± 4.55 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
%timeit X/np.sqrt(X[:,0]**2+X[:,1]**2+X[:,2]**2)[:,np.newaxis]
# 169 µs ± 1.38 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
%timeit X/np.sqrt(sum(X[:,i]**2 for i in range(X.shape[1])))[:,np.newaxis]
# 185 µs ± 4.17 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

我在支持OpenBLAS的MacbookPro 2015上观察到了(1) python3.6 + numpy v1.17.2和(2) python3.9 + numpy v1.19.3的情况。

我不认为这是这篇文章的重复，这篇文章涉及矩阵范数，而这篇文章是关于向量的L2范数。

Answer 1

行L2规范的源代码可以归结为以下几行代码：

def norm(x, keepdims=False):
x = np.asarray(x)
s = x**2
return np.sqrt(s.sum(axis=(1,), keepdims=keepdims))

简化代码假定x为实值，并利用了np.add.reduce(s, ...)等价于s.sum(...)的事实。

因此，OP问题与询问为什么np.sum(x,axis=1)比sum(x[:,i] for i in range(x.shape[1]))慢相同：

%timeit X.sum(axis=1, keepdims=False)
# 131 µs ± 1.6 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
%timeit sum(X[:,i] for i in range(X.shape[1]))
# 36.7 µs ± 91.2 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

这个问题已经在这里得到了回答。简言之，缩减(.sum(axis=1)(带来了开销成本，这些开销通常在浮点精度和速度(例如缓存机制、并行性(方面得到了回报，但在仅缩减三列的特殊情况下却没有。在这种情况下，与实际计算相比，开销相对较大。

如果X有更多的列，情况就会改变。numpy增强的规范化现在比使用python for loop:的减少要快得多

X = np.random.randn(10000,100)
%timeit X/np.linalg.norm(X,axis=1, keepdims=True)
# 3.36 ms ± 132 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%timeit X/np.sqrt(sum(X[:,i]**2 for i in range(X.shape[1])))[:,np.newaxis]
# 5.92 ms ± 168 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

另一个相关的SO线程可以在这里找到：numpy-ufuncsvs.forloop。

问题仍然存在，为什么numpy不显式处理归约的常见特殊情况(例如低轴维度矩阵的列或行上的求和(。也许是因为这种优化的效果往往强烈依赖于目标机器，并大大增加了代码的复杂性。