数值计算技巧

本笔记用于记录在写代码过程中遇到的一些数据处理技巧与注意点。

对特定 axis 进行操作的结果：

axis	对应的维度	操作方向	结果	记忆技巧
axis=0	第一个维度 (行)	垂直向下 (跨行)	得到每列的计算结果	“压扁所有行”
axis=1	第二个维度 (列)	水平向右 (跨列)	得到每行的计算结果	“压扁所有列”

$\ell_2$ 范数计算： $\mid\mid x-y \mid\mid_2$ 的矩阵形式( $X$ 为 $(n, d)$ , $Y$ 为 $(k, d)$ )：

D = -2(X @ Y.T) + (X ** 2).sum(axis=1)[:None] + (Y ** 2).sum(axis=1)

X @ Y.T 返回的矩阵形状是 $(n,k)$ ，而对 $X$ 每一行进行求和返回的向量形状为 $n$ ，因此需要将它拓展为 $(n, 1)$ ，否则 $n$ 和 $k$ 不匹配、就无法使用 numpy 的广播机制了。而 $Y$ 就不需要这样，因为它的形状 $k$ 是匹配的。

筛选矩阵&向量：使用 numpy 的布尔矩阵&向量完成筛选：

np.sum((X - user_means[:,None]) * (X != 0))

I = (X[i,:]!=0) * (X[j,:]!=0)
# Use bool vector as index
xi = X[i,I] - user_means[i]

外积：当我们想用两个向量构建一个矩阵时，就可以使用外积。例如如果我们有向量 $x$ ，我们希望构建矩阵 $X$ ，其中 $X_{ij}=x_i \cdot x_j$ ，则：

X = np.outer(x, x)

布尔掩码的重要性质：当你用这个 mask 来索引 X 时（即 X[mask]），numpy 会将 X 中所有在 mask 里对应位置为 True 的元素全部提取出来，并把它们放进一个新的一维数组（向量）里。它不会保持原来的矩阵结构。
CNN实现逻辑：将对多个小窗口的循环操作转化为对一个高维大张量的单次批量操作。例如下面这个 CNN 的 Pool2D 的 forward 方法：

def forward(self, X: np.ndarray) -> np.ndarray:
    batch_size, in_rows, in_cols, channels = X.shape
    k_h, k_w = self.kernal_shape
    pad_h, pad_w = self.pad 
    
    out_rows = (in_rows + 2 * pad_h - k_h) // self.stride + 1 
    out_cols = (in_cols + 2 * pad_w - k_w) // self.stride + 1 

    X_pad = np.pad(
        X, 
        ((0, 0), (pad_h, pad_h), (pad_w, pad_w), (0, 0)),
        mode='constant'
    )

    X_pool = np.zeros((batch_size, out_rows, out_cols, channels))

    for b in range(batch_size):
        for c in range(channels):
            for h in range(out_rows):
                for w in range(out_cols):
                    h_start = h * self.stride 
                    h_end = h_start + k_h
                    w_start = w * self.stride 
                    w_end = w_start + k_w

                    window = X_pad[b, h_start:h_end, w_start:w_end, c]

                    X_pool[b, h, w, c] = self.pool_fn(window)

可以通过如下的 numpy 操作、一次性把这些滑动窗口提出来。我们可以使用 numpy.lib.stride_tricks.as_strided 方法，创建一个数组视图，这个视图的每个元素指向我们需要计算的小块。总体的流程如下：

一次性地、虚拟地提取出所有样本的所有小块 (patches)
用一个指令完成所有小块和卷积核 $W$ 的乘法与求和
得到最终输出

def forward(self, X: np.ndarray) -> np.ndarray:
    n_examples, in_rows, in_cols, channels = X.shape
    k_h, k_w = self.kernal_shape
    pad_h, pad_w = self.pad
    
    out_rows = (in_rows + pad_h * 2 - k_h) // self.stride + 1 
    out_cols = (in_cols + pad_w * 2 - k_w) // self.stride + 1 

    X_pad = np.pad(
        X, 
        ((0, 0), (pad_h, pad_h), (pad_w, pad_w), (0, 0)),
        mode='constant'
    )

    # Use as_stride to create slide windown view 
    view_shape = (n_examples, out_rows, out_cols, k_h, k_w, channels)
    s_n, s_r, s_c, s_ch = X_pad.strides
    view_strides = (s_n, s_r * self.stride, s_c * self.stride, s_r, s_c, s_ch)
    X_windows = np.lib.stride_tricks.as_strided(X_pad, shape=view_shape, strides=view_strides)

    # Rearrange dim to (N, out_H, out_W, C, k_H, k_W)
    X_windows_permuted = X_windows.transpose(0, 1, 2, 5, 3, 4)
    X_pool = self.pool_fn(X_windows_permuted, axis=(4, 5))

这里的 view_shape 定义了这块内存应该被看作几维，每维有多大，而 view_strides = (s_n, s_r * stride, s_c * stride, s_r, s_c, s_ch) 定义了如何在这块内存中移动：

要到下一个样本（新视图的维度0），移动 s_n 字节。
要到下一个输出行（新视图的维度1），在原数据里移动 stride 个行，即 s_r * stride 字节。
要到下一个输出列（新视图的维度2），在原数据里移动 stride 个列，即 s_c * stride 字节。
要到窗口内的下一行（新视图的维度3），在原数据里移动1个行，即 s_r 字节。
要到窗口内的下一列（新视图的维度4），在原数据里移动1个列，即 s_c 字节。
要到窗口内的下一个通道（新视图的维度5），在原数据里移动1个通道，即 s_ch 字节。

NumPy 对于你的数据代表什么（即“语义”）是完全无知的。它不知道什么是“通道”，什么是“样本”，什么是“像素”。在NumPy的眼里，一个数组仅仅是：

一块连续的内存。

一个 shape 元组，定义了这块内存应该被看作几维，每维有多大。

一个 strides 元组，定义了在内存中如何移动来遍历这些维度。

反向传播梯度计算时，**沿哪几个维度求和，axis就是多少！！！！并且返回的结果是剩下的那个维度！！！**例如下面的 einsum 使用就是在 nhw 维度求和：

dLdW = np.einsum('nhwijk,nhwo->ijko', X_windows, dLdZ)