附录A应用：矩阵乘积的快速算法 A.1矩阵乘积的普通方法 设A={al,B=[bl∈Rn×n,则C=[cl=AB,其中 c=） k=1 易知，总运算量为：n3（乘法)+n3(加法)=2n3 算法A1矩阵乘积：K顺序 1:for i=1 to n do 2 for j=1 to n do for k=1 to n do 5: end for 6: end for 7:end for 以上是计算矩阵乘积的常规实现方式，但并不是速度最快的实现方式.基于矩阵乘积在实际应用中 的重要性，人们在不断寻求计算矩阵乘积的快速方法.对于矩阵乘法的加速，有如下两种策略： ·基于算法的优化，如著名的Strassen算法 ·基于硬件的优化 A2基于硬件的加速方法 对于硬件优化，有循环展开、基于内存布局等的优化技巧.此外还可以多线程（比如C+I1有 库)并发执行.这里探讨一个最简单技巧：调换循环顺序.比如将原来的K顺序改成门顺序： 算法A2.矩阵乘积：K顺序 1:for i 1 to n do 2 for k 1 to n do 3: for j=1 to n do 4: cij =Cij+aikbkj 5: end for
附录 A 应用: 矩阵乘积的快速算法 A.1 矩阵乘积的普通方法 设 A = [aij ], B = [bij ] ∈ R n×n, 则 C = [cij ] = AB, 其中 cij = ∑n k=1 aikbkj . 易知, 总运算量为: n 3（乘法）+ n 3（加法）= 2n 3 . 算法 A.1. 矩阵乘积: IJK 顺序 1: for i = 1 to n do 2: for j = 1 to n do 3: for k = 1 to n do 4: cij = cij + aikbkj 5: end for 6: end for 7: end for 以上是计算矩阵乘积的常规实现方式, 但并不是速度最快的实现方式. 基于矩阵乘积在实际应用中 的重要性, 人们在不断寻求计算矩阵乘积的快速方法. 对于矩阵乘法的加速, 有如下两种策略: • 基于算法的优化, 如著名的 Strassen 算法. • 基于硬件的优化 A.2 基于硬件的加速方法 对于硬件优化, 有循环展开、基于内存布局等的优化技巧. 此外还可以多线程 (比如C++11有 库) 并发执行. 这里探讨一个最简单技巧: 调换循环顺序. 比如将原来的 IJK 顺序改成 IKJ 顺序: 算法 A.2. 矩阵乘积: IKJ 顺序 1: for i = 1 to n do 2: for k = 1 to n do 3: for j = 1 to n do 4: cij = cij + aikbkj 5: end for 1
2 附录A应用：矩阵乘积的快速算法 6:end for 7:end for 例A1在个人电脑(CPU9-10900K,内存128G)上用C语言测试 (test_dgemm.c) ·n=1824,gcc编译加-02选项，KN顺序大概比K顺序快2.5倍.如果加-03选项，则快3.5倍 ·n=2848,cc编译加-02选项，K冈顺序大概比K顺序快6倍.如果加-03选项，则快8倍. 以下分析摘自网络(https:/zhuanlan.zhihu.com/p/1462s0334),仅供参考. 下面分析为什么这两种顺序的计算效果差这么多. 首先从矩阵的存储方式说起.一般而言，矩阵有两种存储方式：一维数组和二维数组.对于矩 阵乘法而言，一维数组显然比二维数组好得多（下文的分析也能看出这一点）.但是如果用矩阵类实 现其他功能的话，可能其他功能用二维数组更方便一些.所以下面对于这两种实现方式都加以分 造成矩阵乘法慢的原因，除了算法本身的复杂度以外，还有内存访问的不连续，这会导致cah 命中率不高.所以为了加速，就要尽可能使内存访问连续，即不要跳来跳去.我们定义一个概念：跳 跃数，来衡量访问的不连续程度。 对于最普通的实现方式四K顺序，它是依次计算C中的每个元素，对应的计算公式是 因此计算c,需要将A的第i行与B的第j列依次相乘相加.按假设，矩阵是按行存储的，所以A 在第i行中不断向右移动时，内存访问是连续的.但B在第j列不断向下移动时，内存访问是不连 续的.计算完时，B已经间断地访问了n次，而A只间断1次（这一次就是算完后跳转到本行的 开头，故总共是n+1次.这样，计算完C中所有n2个元素，跳转了n3+n2次.但刚才没有统计 C的跳转次数，加上以后是n3+n2+n.(注意，在计算完C中每行的最后一个元素时，A是从相应 行末尾转到下一行开头.如果使用一维数组实现的话，这是连续地访问，要诚掉这n次.同时，C也 没有跳转次数了，还要减掉n次.因此对于一维数组，跳转数是n3+n2-n次) 而如果以N顺序实现，则对C一行一行计算，即 G=a41i1+a2b2+…+ambn,i=1,2,,n 其中亡和6分别表示C和B的第i行.当计算C的第i行时，先计算a1高1，即访问a1和B的 第i行（不间断往右移）.然后计算a2d2,此时A往右挪一个元素（不间断），B则跳转到下一行（如 果一维数组则间新一次，一维数组不间断.依次类推，算完C的第行后，恰好按顺序将B遍历一 遍，间断了n次（一维数组是1次），且恰好从左往右遍历了A的第i行，间断了1次（一维数组没有 间断)，加起来是n+1次（一维数组是1次）.故计算C的所有n行后，跳转了n2+n次（一维数组 是n次).刚才没有算C的跳转，算上后跳跃数是2n2+n次（一维数组是n2次）. 由此可见： http://math.ecnu.edu.cn/-jypan
· 2 · 附录 A 应用: 矩阵乘积的快速算法 6: end for 7: end for 例 A.1 在个人电脑 (CPU i910900K, 内存 128G) 上用 C 语言测试. (test_dgemm.c) • n=1024, gcc 编译加 ‐O2 选项, IKJ 顺序大概比 IJK 顺序快 2.5 倍. 如果加 ‐O3 选项, 则快 3.5 倍. • n=2048, gcc 编译加 ‐O2 选项, IKJ 顺序大概比 IJK 顺序快 6 倍. 如果加 ‐O3 选项, 则快 8 倍. 以下分析摘自网络 (https://zhuanlan.zhihu.com/p/146250334), 仅供参考. 下面分析为什么这两种顺序的计算效果差这么多. 首先从矩阵的存储方式说起. 一般而言, 矩阵有两种存储方式: 一维数组和二维数组. 对于矩 阵乘法而言, 一维数组显然比二维数组好得多 (下文的分析也能看出这一点). 但是如果用矩阵类实 现其他功能的话, 可能其他功能用二维数组更方便一些. 所以下面对于这两种实现方式都加以分 析。 造成矩阵乘法慢的原因, 除了算法本身的复杂度以外, 还有内存访问的不连续, 这会导致 cache 命中率不高. 所以为了加速, 就要尽可能使内存访问连续, 即不要跳来跳去. 我们定义一个概念: 跳 跃数, 来衡量访问的不连续程度. 对于最普通的实现方式 (IJK 顺序), 它是依次计算 C 中的每个元素, 对应的计算公式是 cij = ai1b1j + ai2b2j + · · · + ainbnj . 因此计算 cij , 需要将 A 的第 i 行与 B 的第 j 列依次相乘相加. 按假设, 矩阵是按行存储的, 所以 A 在第 i 行中不断向右移动时, 内存访问是连续的. 但 B 在第 j 列不断向下移动时, 内存访问是不连 续的. 计算完 cij 时, B 已经间断地访问了 n 次, 而 A 只间断 1 次 (这一次就是算完后跳转到本行的 开头), 故总共是 n + 1 次. 这样, 计算完 C 中所有 n 2 个元素, 跳转了 n 3 + n 2 次. 但刚才没有统计 C 的跳转次数, 加上以后是 n 3 + n 2 + n. (注意, 在计算完 C 中每行的最后一个元素时, A 是从相应 行末尾转到下一行开头. 如果使用一维数组实现的话, 这是连续地访问, 要减掉这 n 次. 同时, C 也 没有跳转次数了, 还要减掉 n 次. 因此对于一维数组, 跳转数是 n 3 + n 2 − n 次) 而如果以 IKJ 顺序实现, 则对 C 一行一行计算, 即 c˜i = ai1 ˜b1 + ai2 ˜b2 + · · · + ain˜bn, i = 1, 2, . . . , n, 其中 c˜i 和 ˜bi 分别表示 C 和 B 的第 i 行. 当计算 C 的第 i 行时, 先计算 ai1 ˜b1, 即访问 ai1 和 B 的 第 i 行 (不间断往右移). 然后计算 ai2 ˜b2, 此时 A 往右挪一个元素 (不间断), B 则跳转到下一行 (如 果二维数组则间断一次, 一维数组不间断). 依次类推, 算完 C 的第 i 行后, 恰好按顺序将 B 遍历一 遍, 间断了 n 次 (一维数组是 1 次), 且恰好从左往右遍历了 A 的第 i 行, 间断了 1 次 (一维数组没有 间断), 加起来是 n + 1 次 (一维数组是 1 次). 故计算 C 的所有 n 行后, 跳转了 n 2 + n 次 (一维数组 是 n 次). 刚才没有算 C 的跳转, 算上后跳跃数是 2n 2 + n 次 (一维数组是 n 2 次). 由此可见: http://math.ecnu.edu.cn/~jypan
A.2基于硬件的加速方法 .3 ·灯顺序的跳转数渐进地少于K顺序的跳转数 ·一维数组比二维数组好。 跳跃数总结 下面是各个循环顺序的跳跃数列表（写文章时现算的，可能会粗心犯错） ·I顺序：2n2+n(二维数组)n2(一维数组) ·K四顺序：3n2(二维数组)，2n2（一维数组) ·K顺序：n3+n2+n(二维数组)，n3+n2-n(一维数组) ·JK顺序：n3+2n2(仁维数组)，n3+n2+n(一维数组) ·K1顺序：2n3+n(一维数组).2n3(一维数组) ·JK顺序：2n3+n2(仁维数组)，2n3+n2(一维数组) 因此从跳跃数来看，不同顺序的执行效率排序为： KNKⅫ>JK~JK>~JK. 实测速度 测试环境：双CPUXeon4215R(8核16线程)，32G内存.带-02编译选项 NNK灯JK JIK KJI JKI 10240.5500.5703.4303.6107.9207.770 20184.3604.88025.05027.56097.10099.610 凸注记：需要指出的是，影响矩阵乘积的因素有许多，除了前面介绍的跳跃数外，还包括： ·向量化指令运算，如X86架构上的AVX指令集等。 ·循环展开(loop unrolling,编译器中加优化选项后，一般会自动进行循环展开 ·cache blocking(也称tiling.) ·多线程，如C/Ct+多线程，OpenMP 关于内存优化可参考 1]R.E.Bryant and D.R.O'Hallaron,Computer Systems:A Programmer Perspective,3rd edn,Pearson 2016.《深入理解计算机系统》，简称CSAPP.龚奕利，贺莲译 http://math.ecnu.edu.cn/-jypan
A.2 基于硬件的加速方法 · 3 · • IKJ 顺序的跳转数渐进地少于 IJK 顺序的跳转数; • 一维数组比二维数组好. 跳跃数总结 下面是各个循环顺序的跳跃数列表 (写文章时现算的, 可能会粗心犯错) • IKJ 顺序: 2n 2 + n (二维数组), n 2 (一维数组) • KIJ 顺序: 3n 2 (二维数组), 2n 2 (一维数组) • IJK 顺序: n 3 + n 2 + n (二维数组), n 3 + n 2 − n (一维数组) • JIK 顺序: n 3 + 2n 2 (二维数组), n 3 + n 2 + n (一维数组) • KJI 顺序: 2n 3 + n (二维数组), 2n 3 (一维数组) • JKI 顺序: 2n 3 + n 2 (二维数组), 2n 3 + n 2 (一维数组) 因此从跳跃数来看, 不同顺序的执行效率排序为: IKJ ∼ KIJ > IJK ∼ JIK > KJI ∼ JKI. 实测速度 测试环境: 双 CPU Xeon 4215R (8 核 16 线程), 32G 内存. 带 ‐O2 编译选项. N IKJ KIJ IJK JIK KJI JKI 1024 0.550 0.570 3.430 3.610 7.920 7.770 2018 4.360 4.880 25.050 27.560 97.100 99.610 b 注记：需要指出的是, 影响矩阵乘积的因素有许多, 除了前面介绍的跳跃数外, 还包括: • 向量化指令运算, 如 X86 架构上的 AVX 指令集等. • 循环展开 (loop unrolling), 编译器中加优化选项后, 一般会自动进行循环展开. • cache blocking (也称 tiling). • 多线程, 如 C/C++ 多线程, OpenMP. • ... ... 关于内存优化可参考 [1] R. E. Bryant and D. R. O’Hallaron, Computer Systems: A Programmer’s Perspective, 3rd edn, Pearson, 2016. 《深入理解计算机系统》, 简称 CSAPP, 龚奕利, 贺莲译. http://math.ecnu.edu.cn/~jypan
·4 附录A应用：矩阵乘积的快速算法 A.3 Strassen方法 德国数学家Strassen在1969年提出了计算矩阵乘积的快速算法，将运算量降为约O(2.81).下面以 二阶矩阵为例，描述Strassen方法的实现过程.设 A=an an] a21a22 B=bba 则C=AB的每个分量为 C11=a11b1+a12b21,c12=a11b12+a12b22, c21=a21b11+a22b21,c22=a21b12+a22b22. 在Strassen算法中，我们并不直接通过上面的公式来计算C,而是先计算下面7个量： x1=(a11+a22)(b11+b22), x2=(a21+a22)b1: x3=a11(b12-b22, x4=a22(b21-b11, t6=(a11+a12)b22 z6=(a21-a11)(b11+b12, x7=(a12-a22)(b21+b22). 于是，C的各元素可以表示为： C1=C1+x4-T5+r7, G2=x3+工6， c21=x2+T4, c22=x1+3-x2+x6 易知，总共需要做7次乘法和18次加法 下面考虑一般情形.我们采用分而治之的思想，先将矩阵A,B进行2×2分块，即 他a-岛剑 A=A2 Az 则C=AB也可以写成2×2分块形式，即 C1=X1+X4-X5+X7, C12=X3+X5, C21=X2+X, C22=X1+X3-X2+X6 http://math.ecnu.edu.cn/-jypan
· 4 · 附录 A 应用: 矩阵乘积的快速算法 A.3 Strassen 方法 德国数学家 Strassen 在 1969 年提出了计算矩阵乘积的快速算法, 将运算量降为约 O(n 2.81). 下面以 二阶矩阵为例, 描述 Strassen 方法的实现过程. 设 A = [ a11 a12 a21 a22] , B = [ b11 b12 b21 b22] . 则 C = AB 的每个分量为 c11 = a11b11 + a12b21, c12 = a11b12 + a12b22, c21 = a21b11 + a22b21, c22 = a21b12 + a22b22. 在 Strassen 算法中, 我们并不直接通过上面的公式来计算 C, 而是先计算下面 7 个量: x1 = (a11 + a22)(b11 + b22), x2 = (a21 + a22)b11, x3 = a11(b12 − b22), x4 = a22(b21 − b11), x5 = (a11 + a12)b22, x6 = (a21 − a11)(b11 + b12), x7 = (a12 − a22)(b21 + b22). 于是, C 的各元素可以表示为: c11 = x1 + x4 − x5 + x7, c12 = x3 + x5, c21 = x2 + x4, c22 = x1 + x3 − x2 + x6. 易知, 总共需要做 7 次乘法和 18 次加法. 下面考虑一般情形. 我们采用分而治之的思想, 先将矩阵 A, B 进行 2 × 2 分块, 即 A = [ A11 A12 A21 A22] , B = [ B11 B12 B21 B22] , 则 C = AB 也可以写成 2 × 2 分块形式, 即 C11 = X1 + X4 − X5 + X7, C12 = X3 + X5, C21 = X2 + X4, C22 = X1 + X3 − X2 + X6, http://math.ecnu.edu.cn/~jypan
A3 Strassen方法 .5 其中 X1=(411+A2)(B11+B22), X2=(421+A2)B11, X3=A11(B2-B2), X4=A22(B21-B11 X5=(A1+A12)B22, X6=(A21-A1(B1+B12), X7=(412-A2)(B21+B2). 需要7次子矩阵的乘积和18次子矩阵加法.假定采用普通方法计算子矩阵的乘积，即需要(/2)子乘法 和(n/2)3次加法，则采用Strassen方法计算A和B乘积的运算量为 7×(/22+a/2）+18xm/22-7n3+n2 大约是普通矩阵乘积运算量的。·在计算子矩阵的乘积时，我们仍然可以采用Strassen算法.依此类推 于是，由递归思想可知，则总运算量大约为（只考虑最高次项，并假定可以不断对分下去） 7oan=nga7≈n2.807… 课外阅读 记两个n阶矩阵乘积的运算量为O(n),在Strassen算法出现之前，大家一直认为a-3.但随着 Strassen算法的出现，大家意识到这个量级是可以小于3的.但a能不能进一步减少？如果能的话， 它的极限又是多少？这些问题引起了众多学者的极大兴趣.1978年，P证明了a<2.796.一年 后，Bini等又将这个上界改进到2.78.在Pan和Bini等的工作的基础上，Sch6 nhage于1981年证明 了a<2.522.而上界首次突破2.5是由Coppersmith和Winograd提出来的，他们证明了a<2.496 1990年，他们又在Strassen算法的基础上提出了一种新的矩阵乘积计算方法，将运算量级降至著名 的2.376.这个记录一直持续了近二十年，直到2010年前后，Stothers,Vassilevska-Williams和Le Gal 分别将这个上界降到2.37293,2.3728642和2.3728639.2021年，A1man和Williams将上界进一步降 到2.3728596.虽然有许多学者相信，a的极限应该是2，但至今无法证实. 相关参考文献 [1]V.Strassen,Gaussian elimination is not optimal,Numer.Math.,13(196),354-356 [2]V.Y.Pan,Strassen's algorithm is not optimal,In Proc.FOCS,19 (1978),166-176. 3]D.Bini,M.Capovani,FRomani and G.Lotti,O(n)complexity for napproximate matrix multiplication,Inf Proces.Lett,8(197).234-235. [4]A.Schonhage,Partial and total matrix multiplication,SIAM/.Comput.,10(1981),434-455. [5]D.Coppersmith and S.Winograd,On the asymptotic complexity of matrix multiplication,In Pro SFCS.82-90.1981. 6 D.Coppersmith and S.Winograd,Matrix multiplication via arithmetic progressions,.Symbolic Com putatio0n,9(1990),251-280. http://math.ecnu.edu.cn/-jypan
A.3 Strassen 方法 · 5 · 其中 X1 = (A11 + A22)(B11 + B22), X2 = (A21 + A22)B11, X3 = A11(B12 − B22), X4 = A22(B21 − B11), X5 = (A11 + A12)B22, X6 = (A21 − A11)(B11 + B12), X7 = (A12 − A22)(B21 + B22). 需要 7 次子矩阵的乘积和 18 次子矩阵加法. 假定采用普通方法计算子矩阵的乘积, 即需要 (n/2)3 乘法 和 (n/2)3 次加法, 则采用 Strassen 方法计算 A 和 B 乘积的运算量为 7 × ( (n/2)3 + (n/2)3 ) + 18 × (n/2)2 = 7 4 n 3 + 9 2 n 2 . 大约是普通矩阵乘积运算量的 7 8 . 在计算子矩阵的乘积时, 我们仍然可以采用 Strassen 算法. 依此类推, 于是, 由递归思想可知, 则总运算量大约为（只考虑最高次项, 并假定 n 可以不断对分下去） 7 log2 n = n log2 7 ≈ n 2.807... 课外阅读 记两个 n 阶矩阵乘积的运算量为 O(n α). 在 Strassen 算法出现之前, 大家一直认为 α = 3. 但随着 Strassen 算法的出现, 大家意识到这个量级是可以小于 3 的. 但 α 能不能进一步减少？如果能的话, 它的极限又是多少？这些问题引起了众多学者的极大兴趣. 1978 年, Pan 证明了 α < 2.796. 一年 后, Bini 等又将这个上界改进到 2.78. 在 Pan 和 Bini 等的工作的基础上, Schönhage 于 1981 年证明 了 α < 2.522. 而上界首次突破 2.5 是由 Coppersmith 和 Winograd 提出来的, 他们证明了 α < 2.496. 1990 年, 他们又在 Strassen 算法的基础上提出了一种新的矩阵乘积计算方法, 将运算量级降至著名 的 2.376. 这个记录一直持续了近二十年, 直到 2010 年前后, Stothers, VassilevskaWilliams 和 Le Gall 分别将这个上界降到 2.37293, 2.3728642 和 2.3728639. 2021 年, Alman 和 Williams 将上界进一步降 到 2.3728596. 虽然有许多学者相信, α 的极限应该是 2, 但至今无法证实. 相关参考文献: [1] V. Strassen, Gaussian elimination is not optimal, Numer. Math., 13 (1969), 354–356. [2] V. Y. Pan, Strassen’s algorithm is not optimal, In Proc. FOCS, 19 (1978), 166–176. [3] D. Bini, M. Capovani, F. Romani and G. Lotti, O(n 2.7799) complexity for n × n approximate matrix multiplication, Inf. Process. Lett., 8 (1979), 234–235. [4] A. Schönhage, Partial and total matrix multiplication, SIAM J. Comput., 10 (1981), 434–455. [5] D. Coppersmith and S. Winograd, On the asymptotic complexity of matrix multiplication, In Proc. SFCS, 82–90, 1981. [6] D. Coppersmith and S. Winograd, Matrix multiplication via arithmetic progressions, J. Symbolic Computation, 9 (1990), 251–280. http://math.ecnu.edu.cn/~jypan
6 附录A应用：矩阵乘积的快速算法 7A.)Stothers,On the complexity of matrix multiplication.PhD thesis.2010. [8]V.Vassilevska-Williams,Breaking the Coppersmith-Winograd barrier,In44th ACM Symposium on 7 eory of Computing(ST0C2012以，2012. 9]E.Le Gall,Powers of tensors and fast matrix multiplication,In39thInternational Symposimon Symbolic and Algebraic Computation (TSSAC 2014),2014,296-303.arXiv:1401.7714. [10]J.Alman and V.V.Williams,A refined laser method and faster matrix multiplication,Proceedings of the 2021 ACM-SIAM Symposium on Discrete Algorithms (SODA),2021. 2 here 24 23 22 https t1p11cat1on-20a1-20218323/ http://math.ecnu.edu.cn/-jypan
· 6 · 附录 A 应用: 矩阵乘积的快速算法 [7] A. J. Stothers, On the complexity of matrix multiplication, PhD thesis, 2010. [8] V. VassilevskaWilliams, Breaking the Coppersmith–Winograd barrier, In 44th ACM Symposium on Theory of Computing (STOC 2012), 2012. [9] F. Le Gall, Powers of tensors and fast matrix multiplication, In 39th International Symposium on Symbolic and Algebraic Computation (ISSAC 2014), 2014, 296–303. arXiv:1401.7714. [10] J. Alman and V. V. Williams, A refined laser method and faster matrix multiplication, Proceedings of the 2021 ACMSIAM Symposium on Discrete Algorithms (SODA), 2021. (https://www.quantamagazine.org/mathematicians‐inch‐closer‐to‐matrix‐multiplication‐goal‐20210323/) http://math.ecnu.edu.cn/~jypan