What & How & Why

这是本文档旧的修订版!


Matrix-Matrix Multiplication

LAFF Week 5 Notes

本章的重点内容:

  • 如何从分割的角度去看待矩阵乘法
  • 矩阵乘法的性质、转置矩阵、特殊矩阵与矩阵乘法
  • 从行、列、RANK 1 的观点看待矩阵乘法和算法

矩阵乘法的分割

假设我们有矩阵乘法运算 C=ABC = AB,那么该运算可以做如下的分割:

BB 的水平分割:

(C0C1)=A(B0B1)=(AB0AB1)\left(\begin{array}{c|c} C_0 & C_1 \end{array}\right) = A\left(\begin{array}{c|c} B_0 & B_1 \end{array}\right) =\left(\begin{array}{c|c} AB_0 & AB_1 \end{array}\right)

AA 的垂直分割: (C0C1)=(A0A1)B=(A0BA1B) \left (\frac{C_0}{C_1} \right) = \left (\frac{A_0}{A_1} \right)B = \left (\frac{A_0B}{A_1B} \right)

AA 的水平分割加上对 BB 的垂直分割 :

这种分割在运算中可以表现为子矩阵的点积

C=(A0A1)(B0B1)=A0B0+A1B1 C = \left(\begin{array}{c|c} A_0 & A_1 \end{array}\right)\left (\frac{B_0}{B_1} \right) = A_0B_0+A_1B_1

CC 的十字分割:

该分割实质上就是在做子矩阵的二维矩阵乘法(点积观点AA 乘以 BB

(C00C01C10C11)=(A00A01A10A11)(B00B01B10B11)=(A00B00+A01B10A00B01+A01B11A10B00+A11B10A10B01+A11B11)\left(\begin{array}{c|c} C_{00} & C_{01} \\ \hline C_{10} & C_{11} \end{array} \right) = \left(\begin{array}{c|c} A_{00} & A_{01} \\ \hline A_{10} & A_{11} \end{array} \right) \left(\begin{array}{c|c} B_{00} & B_{01} \\ \hline B_{10} & B_{11} \end{array} \right) \\ =\left(\begin{array}{c|c} A_{00}B_{00} + A_{01}B_{10} & A_{00}B_{01}+A_{01}B_{11} \\ \hline A_{10}B_{00} + A_{11}B_{10} & A_{10}B_{01}+A_{11}B_{11} \\ \end{array} \right)

需要注意的是,矩阵分割中必须保证对应子矩阵的维度匹配

矩阵乘法的属性

矩阵乘法的属性

矩阵乘法有如下的属性:

结合律:
(AB)C=A(BC)(AB)C=A(BC)
分配律:

分配律对于左乘和右乘都是适用的。
A(B+C)=AB+ACA(B+C)=AB+AC (A+B)C=AC+BC(A+B)C = AC + BC

矩阵乘法与转置

矩阵乘法的转置矩阵等于该乘法因子反转后的转置矩阵相乘:
(AB)T=BTAT(AB)^T = B^TA^T
该性质可以做如下推广:
(ABZ)T=ZTBTAT(AB \dots Z)^T = Z^T \dots B^TA^T

特殊矩阵的乘法性质


矩阵与单位矩阵

矩阵与单位矩阵的乘积还是其本身:
AI=IA=AA I = I A = A

矩阵与对角矩阵

矩阵与对角矩阵的乘积等于该矩阵的每列乘以对应对角矩阵的元素:
A=(a0a1an1).AD=(δ0a0δ1a1δn1an1). A= \left( \begin{array}{r|r|r|r} a_0 & a_1 & \dots & a_{n-1} \end{array} \right).\\ AD= \left ( \begin{array}{r|r|r|r} \delta_0 a_0 & \delta_1 a_1 & \dots & \delta_{n-1}a_{n-1} \end{array} \right).
该性质也适用于行:
A=(a~0Ta~1Ta~m1T)DA=(δ0a~0Tδ1a~1Tδm1a~m1T) A= \left( \begin{array}{c} \widetilde a_0^{T} \\ \hline \widetilde a_1^{T} \\ \hline \vdots \\ \hline \widetilde a_{m-1}^{T} \end{array} \right)\\ DA= \left ( \begin{array}{c} \delta_0 \widetilde a_0^{T} \\ \hline \delta_1 \widetilde a_1^{T} \\ \hline \vdots \\ \hline \delta_{m-1} \widetilde a_{m-1}^{T} \end{array} \right)

矩阵与三角矩阵

上(下)三角矩阵与上(下)三角矩阵的乘积也是上(下)三角矩阵。

矩阵与对称矩阵

矩阵与其转置矩阵的乘积为对称矩阵。
ATA  is symmetricA^TA \,\, \text {is symmetric}
转置矩阵与其原矩阵的乘积为对称矩阵。
AAT  is symmetricAA^T \,\, \text {is symmetric}
对称矩阵的和也是对称矩阵。
A+ATA  is symmetricA+A^TA \,\, \text {is symmetric}
对称矩阵的乘积也是对称矩阵。

矩阵乘法算法

从矩阵乘法的定义上看,我们可以通过三层循环来遍历参与运算的矩阵。而这三层循环的顺序是可以改变的;我们可以从行观点、列观点、rank1 观点来看待遍历。而这三种观点反映到算法中则意味着对应的循环在嵌套中的位置。但不管循环位置如何,这三层循环组合在一起实现的都是矩阵元素的遍历。因此不管怎么调整循环的位置,对最后的乘法的结果都是没有关系的。下面我们分别来看看这些不同观点的算法:

列观点算法

C=ABC = AB 为例,行观点的算法的遍历顺序:

  • 最外层循环为 BB 中的列循环,即 j=0n1j = 0 \to n-1
  • 中间层循环为 AA 中的行循环,即 i=0m1 i = 0 \to m-1;
  • 内层循环为 ABA、B 对应的行列点积,即 Ci,j=p=0k1αi,pβp,j+Ci,jC_{i,j} = \sum^{k-1}_{p=0}\alpha_{i,p}\beta_{p,j}+C_{i,j}



详情可以参考下图:

<html>

<img src=“/_media/math/linear_algebra/laff/m_m_col_first.svg” width=“500”>

</html>
明白了内部循环实际上是 AA 的行与 BB 的列的点积以后,那么实际上内部的两层循环(对点积的循环,对 AA 的行的循环)组合在一起,就可以表示为:ABjAB_j

因此,针对于最外层对 BB 的列 jj 的循环,我们可以把内部的循环写成: Cj=ABj+CjC_j = AB_j + C_jCC 的第 jj 列 等于 矩阵 AA 与 矩阵 BB 的第 jj 列的乘积结果。因此该算法可以写成:



行观点算法

相对于列观点算法,行观点算法的修改如下:

  • 最外层的循环为 AA 中的行循环(i=0n1i = 0 \to n-1
  • 中间层的循环为 BB 中的列循环 (j=0m1j = 0 \to m-1
  • 内部循环依然是 AA 中的行与 BB 中的列的点积

因为最外层是 AA 行循环,因此每一次总的循环都是 AA 的一行依次与 BB 中的每一列进行点积。过程参考下图:

<html>

<img src=“/_media/math/linear_algebra/laff/m_m_row_first.svg” width=“500”>

</html>

将内部两层循环合并在一起,就可以表示为:aiTBa_i^TB。因此,因此,针对于最外层对 AA 的列 ii 的循环,内部循环可以写成: Ci=aiTB+CiC_i = a_i^TB + C_iCC 的第 ii 行等于 矩阵 AA 的第 ii 行与矩阵 BB 的乘积结果。因此该算法可以写成:



Rank-1观点算法

比起之前的两种算法,Rank-1 观点下的矩阵乘法的结构如下:

  • 最外部的循环为 pp 的循环,也就是同时遍历 AA 的列 与 BB 的行。
  • 内部两层循环可以任意选择,无论选择 AA 列优先还是 BB 行优先,最终得到的结果都是 AA 列与 BB 行的外积

可以看到的是,Rank-1 观点下的算法改用外积得到一个与 CC 维度相同的矩阵,每一次内部循环对 CC 的更新实际上就是将外积得到的矩阵累加到 CC 上。来以 AA 行优先为例,看看示意图:

<html>

<img src=“/_media/math/linear_algebra/laff/m_m_p_first_r.svg” width=“500”>

</html>

根据上图,每一次内部循环的结果可以写成:AiBjTA_iB_j^T,因此该算法的内部循环为: C=AiBjT+CC = A_iB_j^T+C

具体算法如下:


矩阵乘法的优化

关键的思路:减少对内存的读写,尽量将矩阵的乘法分割到每次运算刚好可以在 CPUCPU 的一次操作内完成,(也就是最大限度的利用 CPU 的缓存,将需要反复使用的内容添加到寄存器,L1,L2 cache中)。

具体实现方法待完善。

参考资料