矩阵乘法优化:从高斯消元到Strassen与Winograd算法

📅 2026/8/24 20:38:48
矩阵乘法优化:从高斯消元到Strassen与Winograd算法
1. 从“昂贵”的乘法说起为什么我们要费尽心思减少它在计算机的世界里加法、减法和移位操作通常被认为是“廉价”的它们可以在一个时钟周期内快速完成。相比之下乘法操作就显得“昂贵”得多。这种“昂贵”体现在两个方面一是硬件实现更复杂需要更多的逻辑门和更长的数据通路导致其执行周期数通常是加法操作的数倍二是在高精度计算如大整数、浮点数场景下乘法的时间复杂度是O(n²)随着数字位数的增加其耗时呈平方级增长迅速成为性能瓶颈。最近网络上热议的“浮点数乘法速度”、“高精度乘法”等话题恰恰反映了工程师们在性能优化时对乘法操作的敏感。无论是SLAM同步定位与建图系统中的图优化算法还是各种昂贵的多模态优化算法其核心计算单元往往涉及大规模的矩阵运算而矩阵乘法的核心就是标量乘法的累积。每一次不必要的乘法都在消耗着宝贵的计算资源和时间。因此减少乘法次数尤其是减少核心算法中的乘法操作成为了算法优化中一个经典且至关重要的课题。这不仅仅是理论上的游戏它直接关系到程序的响应速度、服务器的计算成本乃至电池供电设备的续航能力。今天我们就来深入探讨三种在矩阵乘法领域里程碑式的优化算法经典的高斯消元法优化、颠覆性的Strassen算法以及在此基础上进一步精雕细琢的Winograd算法。我们将绕过枯燥的公式堆砌直接从“为什么有效”和“怎么用”的角度结合代码和性能对比把它们讲明白。2. 高斯消元法中的乘法优化朴素想法与实战技巧当我们提到“Gauss”时第一反应通常是高斯消元法用于求解线性方程组。但这里所说的“减少乘法次数的优化算法Gauss”更准确地是指在高斯消元法执行过程中通过巧妙的数学变换来减少不必要的乘法运算次数。这是一种针对特定计算过程的局部优化体现了“算法即优化”的思想。2.1 朴素高斯消元法的乘法开销分析假设我们有一个n阶线性方程组对应一个n x n的系数矩阵A。标准的无选主元高斯消元分为两个步骤前向消元和回代。前向消元对于第k列k从1到n-1我们需要将第k行以下所有行的第k个元素消为0。对于第i行i k计算乘数m A[i][k] / A[k][k]然后用第i行减去m乘以第k行。这个过程中对第i行的第j列j ≥ k元素进行更新A[i][j] A[i][j] - m * A[k][j]。回代从最后一行开始依次求解未知数。让我们粗略估算一下乘法次数。在前向消元中计算每个乘数m需要1次除法可视为一次特殊的乘法。对于每个m需要用它去乘第k行的 (n-k) 个元素然后做减法。所以对于第k步乘法次数约为(n-k) (n-k)*(n-k)。求和后前向消元的总乘法次数约为(2n³)/3 O(n²)。回代过程的乘法次数约为n²/2。因此朴素算法的主项是(2/3)n³量级。2.2 优化思路避免重复计算与缩放技巧优化的核心在于识别并消除重复计算。一个常见的技巧是“缩放行”Scaling。在消元开始前我们可以遍历每一行将其除以该行的第一个非零元素或绝对值最大的元素兼顾数值稳定性。这样做的目的是使每行用于消去其他行的“主元”系数变为1。为什么这能减少乘法在原始算法中每次消元都需要计算乘数m A[i][k] / A[k][k]这是一个除法。然后对于第i行的每一个元素A[i][j]都需要执行一次乘法m * A[k][j]。 如果我们预先将第k行除以了A[k][k]使得A[k][k] 1那么乘数m就简化为A[i][k]本身因为m A[i][k] / 1。更重要的是在更新A[i][j]时公式变为A[i][j] A[i][j] - A[i][k] * A[k][j]。这里A[i][k]是已知的A[k][j]是缩放后行中已知的。效果分析我们避免了每次计算乘数m的那一次除法。虽然这看起来只是一次操作但在整个O(n³)的过程中我们节省了大约n²/2次除法运算。更重要的是这种形式有时能让编译器或CPU更好地进行指令级优化。然而它并没有改变算法O(n³)的渐近复杂度。注意这种缩放操作本身需要O(n²)次除法这是一个一次性开销。对于只求解一次方程组的情况总节省可能不明显甚至可能因额外的缩放开销而变慢。但在需要反复对同一个矩阵A进行消元仅右侧常数向量b改变时缩放只需做一次后续每次消元都能受益这时优势才会凸显。这引出了矩阵的LU分解思想。2.3 实战心得与代码片段在实际编码中纯粹的“高斯乘法优化”往往融入更通用的优化策略内存访问优化确保按行连续访问内存C/C、Python NumPy数组是行优先。对于第i行的更新A[i][j]和A[k][j]的访问是连续的这对CPU缓存友好。循环展开手动或依靠编译器展开内层j循环可以减少循环开销提高指令流水线效率。使用BLAS库在性能关键的应用中没有人会手写原生循环。我们会直接调用高度优化的基础线性代数子程序库如OpenBLAS、Intel MKL。这些库使用了分块、向量化、多线程等技术将性能压榨到极致。高斯消元对应的就是LAPACK库中的gesv或getrfLU分解函数。import numpy as np import time def gauss_elimination_naive(A, b): 朴素高斯消元展示乘法点 n len(A) A A.astype(float).copy() b b.astype(float).copy() # 前向消元 for k in range(n-1): for i in range(k1, n): factor A[i, k] / A[k, k] # 一次除法昂贵操作 for j in range(k, n): A[i, j] - factor * A[k, j] # 一次乘法一次减法 b[i] - factor * b[k] # 回代 x np.zeros(n) for i in range(n-1, -1, -1): x[i] b[i] for j in range(i1, n): x[i] - A[i, j] * x[j] # 一次乘法 x[i] / A[i, i] # 一次除法 return x def gauss_elimination_with_scaling(A, b): 演示缩放思想实际常与列选主元结合 n len(A) A A.astype(float).copy() b b.astype(float).copy() # 可选的缩放步骤这里演示为每行除以该行最大元素数值稳定 # 但注意这改变了问题仅作为思想演示。 # 更常见的优化是进行LU分解并保存乘子。 for i in range(n): scale np.max(np.abs(A[i])) if scale ! 0: A[i] / scale b[i] / scale # 后续的消元过程乘数计算会略有不同 # ... (此处省略具体消元代码因缩放改变了算法语义) # 实际中我们直接调用优化库 # x np.linalg.solve(A, b) # 性能对比使用小矩阵演示大矩阵请用专业库 n 100 A np.random.randn(n, n) b np.random.randn(n) start time.time() x_naive gauss_elimination_naive(A, b) time_naive time.time() - start start time.time() x_lib np.linalg.solve(A, b) # 使用高度优化的LAPACK库 time_lib time.time() - start print(f朴素消元时间: {time_naive:.4f}秒) print(fNumPy solve时间: {time_lib:.4f}秒) print(f结果差异范数: {np.linalg.norm(x_naive - x_lib):.2e})个人体会在小型、一次性计算中手动优化高斯消元的乘法意义不大现代编译器和解释器如NumPy的优化远超普通人的手写代码。真正的价值在于理解其思想通过改变计算顺序或预处理数据将重复的、昂贵的运算如除以同一个主元提取出来从而减少总运算量。这个思想会延续到更复杂的算法中。3. Strassen算法打破传统用加法换乘法的魔法长久以来人们认为两个n×n矩阵相乘的复杂度下限就是O(n³)因为标准算法需要三层循环总共进行n³次乘法和大约n³次加法。直到1969年Volker Strassen发表了一个震惊学术界的算法证明了矩阵乘法可以在O(n^2.807)的时间内完成。这是如何做到的核心就是用大量的加法来替换少量的乘法。3.1 算法核心2x2矩阵乘法的七次乘法Strassen算法的基石是它将矩阵分块并递归地应用。我们从一个2x2的矩阵乘法看起 设C A * B其中A, B, C都是2x2矩阵。 标准算法需要8次乘法和4次加法C[0,0] A[0,0]*B[0,0] A[0,1]*B[1,0] // 2次乘1次加 C[0,1] A[0,0]*B[0,1] A[0,1]*B[1,1] // 2次乘1次加 C[1,0] A[1,0]*B[0,0] A[1,1]*B[1,0] // 2次乘1次加 C[1,1] A[1,0]*B[0,1] A[1,1]*B[1,1] // 2次乘1次加Strassen构造了7个辅助矩阵M1到M7每个M都是A和B的子矩阵的线性组合只做加法/减法后的乘积M1 (A[0,0] A[1,1]) * (B[0,0] B[1,1]) M2 (A[1,0] A[1,1]) * B[0,0] M3 A[0,0] * (B[0,1] - B[1,1]) M4 A[1,1] * (B[1,0] - B[0,0]) M5 (A[0,0] A[0,1]) * B[1,1] M6 (A[1,0] - A[0,0]) * (B[0,0] B[0,1]) M7 (A[0,1] - A[1,1]) * (B[1,0] B[1,1])然后C的四个子块可以通过M1到M7的加法组合得到C[0,0] M1 M4 - M5 M7 C[0,1] M3 M5 C[1,0] M2 M4 C[1,1] M1 - M2 M3 M6数一数乘法只有7次计算M1到M7而加法和减法则有18次。我们用了10次额外的加法/减法换来了1次乘法的减少。3.2 递归分治与复杂度分析对于更大的n×n矩阵假设n是2的幂次Strassen算法这样做将矩阵A、B和C都平均分成4个n/2 × n/2的子块。按照上述2x2的公式但把公式中的每个元素如A[0,0]替换为一个n/2 × n/2的子矩阵把标量乘法和加法替换为子矩阵的乘法和加法。计算M1到M7这7个子矩阵乘法。注意每个子矩阵乘法本身又是两个n/2 × n/2矩阵的相乘因此可以递归调用Strassen算法本身。通过M1到M7的加/减组合得到C的四个子块。复杂度推导设T(n)为计算n×n矩阵乘法的耗时。标准算法T(n) Θ(n³)。Strassen算法每次递归将问题划分为7个规模为n/2的子问题外加一些Θ(n²)的矩阵加/减法因为加/减需要操作n²个元素。因此递归式为T(n) 7 * T(n/2) Θ(n²)根据主定理Master Theorema7, b2, d2。由于log_b(a) log_2(7) ≈ 2.807 d2所以解为T(n) Θ(n^{log_2(7)}) ≈ Θ(n^{2.807})。这打破了O(n³)的壁垒。后续更复杂的算法如Coppersmith–Winograd算法及其优化将指数降到了更低但Strassen算法是第一个也是在实际中有时会被用到的。3.3 实战考量为什么我们不总是用Strassen理论上更优但Strassen算法在实践中有明显的开销巨大的常数因子算法中大量的矩阵加/减法Θ(n²)带来了显著的额外开销。这个常数因子比标准算法的常数因子大得多。递归开销函数调用、矩阵分块和合并都需要消耗时间。数值稳定性由于引入了更多的加/减步骤浮点数计算中的舍入误差可能会被放大导致结果精度略低于标准算法。对于病态矩阵这可能是个问题。内存占用递归过程中需要创建许多临时矩阵来存储M1到M7以及中间结果内存消耗更大。因此在实际的线性代数库如OpenBLAS、MKL中通常会设置一个阈值。当矩阵规模n小于这个阈值通常是几十到几百时使用经过极度优化的标准矩阵乘法可能已用汇编和SIMD指令优化只有当n大于阈值时才切换到Strassen算法或其变种以期用渐近复杂度的优势抵消常数因子的劣势。import numpy as np def strassen_multiply(A, B, threshold64): Strassen矩阵乘法递归实现 (简化版假设n是2的幂且大于threshold) n A.shape[0] # 递归基如果矩阵很小使用朴素算法 if n threshold: return np.dot(A, B) # 分块 mid n // 2 A11, A12, A21, A22 A[:mid, :mid], A[:mid, mid:], A[mid:, :mid], A[mid:, mid:] B11, B12, B21, B22 B[:mid, :mid], B[:mid, mid:], B[mid:, :mid], B[mid:, mid:] # 计算7个辅助矩阵 M1...M7 (递归调用) M1 strassen_multiply(A11 A22, B11 B22, threshold) M2 strassen_multiply(A21 A22, B11, threshold) M3 strassen_multiply(A11, B12 - B22, threshold) M4 strassen_multiply(A22, B21 - B11, threshold) M5 strassen_multiply(A11 A12, B22, threshold) M6 strassen_multiply(A21 - A11, B11 B12, threshold) M7 strassen_multiply(A12 - A22, B21 B22, threshold) # 组合结果矩阵C的子块 C11 M1 M4 - M5 M7 C12 M3 M5 C21 M2 M4 C22 M1 - M2 M3 M6 # 合并子块 C np.vstack((np.hstack((C11, C12)), np.hstack((C21, C22)))) return C # 测试与对比 n 256 # 必须是2的幂且大于阈值 A np.random.randn(n, n) B np.random.randn(n, n) print(开始计算...) import time start time.time() C_strassen strassen_multiply(A, B, threshold64) time_strassen time.time() - start start time.time() C_standard np.dot(A, B) # 使用高度优化的BLAS time_standard time.time() - start print(fStrassen (阈值64) 时间: {time_strassen:.4f}秒) print(fNumPy dot (BLAS) 时间: {time_standard:.4f}秒) print(f结果差异范数: {np.linalg.norm(C_strassen - C_standard):.2e}) # 注意由于递归和Python开销这个Python实现几乎不可能快过NumPy的CBLAS。个人踩坑点自己实现Strassen算法时最容易忽略的是递归基阈值的选择。阈值设得太高可能永远触发不了Strassen的优势设得太低巨大的递归和临时矩阵开销会让程序慢得无法忍受。需要在实际的硬件和问题上进行性能剖析Profiling来确定。此外处理矩阵规模不是2的幂次的情况需要额外的填充和裁剪逻辑这又会引入额外开销。4. Winograd算法在Strassen基础上的进一步精打细算Strassen算法展示了用加法换乘法的可能性而Winograd算法则可以看作是在这个方向上更极致的优化。它通常指的是由Shmuel Winograd提出的一类算法用于减少卷积或矩阵乘法中的乘法次数。在矩阵乘法语境下我们常讨论的是Winograd对Strassen算法的改进或者更广义的基于多项式插值或中国剩余定理的矩阵乘法算法。4.1 Winograd矩阵乘法的基本思想我们以计算两个2x2矩阵乘法为例展示经典的Winograd算法有时被称为Winograd’s variant of Strassen。目标同样是C A * B。 它定义如下6个辅助量其中包含5次乘法比Strassen的7次更少S1 A[1,0] A[1,1] S2 S1 - A[0,0] // 注意S2用了S1 S3 A[0,0] - A[1,0] S4 A[0,1] - S2 T1 B[0,1] - B[0,0] T2 B[1,1] - T1 T3 B[1,1] - B[0,1] T4 T2 - B[1,0] // 5次乘法 P1 A[0,0] * B[0,0] P2 A[0,1] * B[1,0] P3 S1 * T1 P4 S2 * T2 P5 S3 * T3 P6 S4 * B[1,1] // 这是第6个乘法等等我们数一下。 // 实际上标准的Winograd 2x2公式需要7次乘法但通过共享中间结果可以优化。 // 一个常见的Winograd变体Winograds identity用于卷积矩阵乘法有类似思想但具体公式不同。 // 更准确地说对于2x2矩阵已知的最少乘法次数是7次StrassenWinograd的工作是证明了对于某些特定规模如2x27次是乘法次数的下限并且他给出了另一种7次乘法的方案但加法次数可能更少或数值性质不同。实际上对于2x2矩阵乘法Winograd给出的一组公式也使用7次乘法但需要15次加法Strassen是18次加法。他的贡献在于系统地研究了乘法复杂度并给出了在某些规模下理论上最优的乘法次数。4.2 与Strassen的对比加法与乘法的权衡Winograd算法的核心价值在于它在保持乘法次数最小或次优的同时试图减少加法次数或改善计算的数值稳定性。对于2x2情况Strassen: 7次乘法18次加法/减法。Winograd (一种典型变体): 7次乘法15次加法/减法。虽然乘法次数相同但加法减少了3次。在递归过程中这些节省的加法会像滚雪球一样累积。递归方程仍然是T(n) 7T(n/2) O(n²)但O(n²)项的常数因子更小。因此在渐近复杂度相同O(n^2.807)的前提下Winograd变体的实际运行时间可能比原始Strassen算法略短因为它每一步的“额外工作”矩阵加减更少。4.3 实际应用与局限性Winograd算法及其思想在深度学习的卷积计算中得到了广泛应用例如著名的cuDNN库中就实现了基于Winograd变换的卷积算法。它将卷积核和输入进行变换使得核心计算变为更少的乘法操作从而极大加速了卷积神经网络的前向传播。然而在通用的稠密矩阵乘法库如BLAS的GEMM中Winograd算法并不像Strassen算法那样常见。主要原因有实现复杂性Winograd算法的推导和实现比Strassen更复杂分块和合并的逻辑更繁琐。数值稳定性一些Winograd变体可能比Strassen有更差的数值特性对于科学计算来说这是致命的。常数因子收益递减在高度优化的工业级BLAS实现中Strassen的常数因子已经很大Winograd带来的额外加法减少的收益可能无法抵消其增加的实现复杂性和潜在的数据搬运开销。更适合特定模式Winograd的思想在卷积这种固定小核、滑动窗口的操作中能发挥最大优势因为变换可以预先计算并复用。而通用矩阵乘法的两个矩阵都是任意的。# 以下是一个简化的Winograd风格2x2矩阵乘法示例基于一种常见形式 # 这主要用于展示思想并非最优或唯一实现。 def winograd_2x2(A, B): 计算2x2矩阵乘法的Winograd变体 (7次乘法15次加法) 注意此代码为演示原理未做优化。 # 提取元素 a11, a12, a21, a22 A[0,0], A[0,1], A[1,0], A[1,1] b11, b12, b21, b22 B[0,0], B[0,1], B[1,0], B[1,1] # 预计算S和T中间变量共8次加法 S1 a21 a22 S2 S1 - a11 S3 a11 - a21 S4 a12 - S2 T1 b12 - b11 T2 b22 - T1 T3 b22 - b12 T4 T2 - b21 # 7次乘法 P1 a11 * b11 P2 a12 * b21 P3 S1 * T1 P4 S2 * T2 P5 S3 * T3 P6 S4 * b22 # 注意这里需要第7个乘法。标准的Winograd公式中P6可能定义为S4*某物然后还有P7。 # 一个完整的7乘法版本需要更复杂的组合。此处为了演示我们假设P6是其中一个乘积。 # 实际上Winograd的完整公式需要精确构造7个乘积项。 # 让我们采用另一种已知的7乘法方案Winograds identity: # 实际上常见的展示是Strassen算法本身。Winograd的贡献是理论证明和给出另一种方案。 # 为简化我们回到Strassen的M1...M7来计算但强调Winograd的思想是优化加法。 # 下面我们计算一个已知的、加法更少的7乘法方案有时被称为Winograd变体: # 定义不同的中间变量 (以下为一组示例系数实际需严格推导) # 由于推导复杂此处我们改用一组演示性的计算来体现“用不同线性组合减少后续加法”的思想。 # 更实际的做法是直接比较Strassen和Winograd变体的加法次数。 # 结论是对于2x2Winograd变体用15次加法 vs Strassen的18次。 # 组合结果 (需要7次加法) # 这里我们用Strassen的结果组合公式但假设我们的P1..P7已经对应了Winograd的乘积。 # 实际上组合方式也会不同。 # 因此代码部分我们只给出概念说明 # Winograd算法的实现需要精确的数学推导其递归分治结构与Strassen类似 # 但每一步的线性组合系数和后续组合方式不同。 # 作为演示我们直接输出Strassen结果以示区别不大。 # 真正要实现Winograd矩阵乘法需要查阅其精确的矩阵分解公式。 print(Winograd算法具体实现涉及复杂的线性组合推导此处从略。) print(其核心价值在于在保持7次乘法的同时将加法次数从Strassen的18次降低到15次。) return None # 在实际应用中我们不会手写Winograd for GEMM而是使用库。 # 在卷积中调用如torch.nn.functional.conv2d后端可能会自动选择Winograd算法。经验之谈除非你是底层计算库的开发者或者正在为一个特定领域如深度学习推理框架寻找极致的卷积优化否则你通常不需要手动实现Winograd矩阵乘法。对于大多数应用理解其“通过线性变换预处理来减少核心运算次数”的思想更为重要。当你设计自己的算法时如果发现核心是密集的乘法运算可以思考能否通过预先计算一些加/减组合将多个乘法合并或简化这就是Strassen和Winograd带给我们的核心启示。5. 算法选择与工程实践理论复杂度不等于实际速度我们讨论了三种减少乘法次数的思路高斯消元中的局部预处理、Strassen的分治递归、Winograd的精细变换。面对一个实际项目我们该如何选择5.1 性能考量因素金字塔决定矩阵运算性能的因素从底层到高层如下硬件与指令集CPU的缓存层次L1/L2/L3、向量化指令如AVX-512、乘加融合指令FMA的支持。这是性能的物理基础。内存访问模式算法是否具有空间局部性是否能连续访问内存糟糕的访问模式如频繁跳转会让CPU缓存失效性能下降数十倍。Strassen/Winograd递归过程中的矩阵分块和合并会带来额外的数据拷贝开销。渐近复杂度即我们讨论的O(n³) vs O(n^2.807)。当n非常大时这项起主导作用。常数因子包括算法中的加法次数、数据搬运量、递归开销等。Strassen的常数因子很大。数值稳定性对于科学计算结果的精度至关重要。有些优化算法会引入更大的舍入误差。实现复杂度代码是否易于维护、调试和移植5.2 实战决策指南基于以上因素我的个人建议是对于中小规模矩阵n 100~500永远不要自己实现Strassen/Winograd来加速。直接调用高度优化的BLAS库如NumPy的dot、matmulPyTorch的运算符。这些库使用了循环展开、分块、向量化、多线程等技术其常数因子极小在中小规模上碾压任何手写的递归算法。高斯消元也同理用numpy.linalg.solve或scipy.linalg.solve。对于大规模矩阵n 几千通用计算仍然优先使用BLAS库如Intel MKL, OpenBLAS。这些库内部可能已经集成了Strassen算法并智能地在不同规模间切换。例如在矩阵规模极大时某些BLAS实现会递归调用Strassen算法。你应该信任这些经过数十年优化的工业级库。特殊场景如果你的问题有特殊结构如稀疏矩阵、带状矩阵、对称正定矩阵有更专用的算法和库如SuiteSparse、Eigen它们能带来数量级的提升。自研算法只有在你确信自己是计算性能的专家并且现有库无法满足需求例如你的矩阵乘法有极其特殊的模式或者你在为特定硬件设计算法时才考虑基于Strassen/Winograd思想进行定制化开发。关于数值稳定性如果你处理的是病态矩阵或要求高精度应使用经过严格数值验证的库如LAPACK并可能避免使用Strassen算法。在深度学习等容错性较高的领域数值稳定性要求相对宽松。5.3 一个简单的性能实验与误区警示让我们用Python做一个思想实验揭示理论复杂度和实际速度的差异import numpy as np import time import matplotlib.pyplot as plt def standard_matmul(A, B): 三层循环的标准矩阵乘法 O(n^3) n A.shape[0] C np.zeros((n, n)) for i in range(n): for j in range(n): s 0 for k in range(n): s A[i, k] * B[k, j] # 最内层是乘法累加 C[i, j] s return C sizes [10, 20, 40, 60, 80, 100] times_standard [] times_numpy [] for n in sizes: A np.random.randn(n, n) B np.random.randn(n, n) start time.time() C1 standard_matmul(A, B) times_standard.append(time.time() - start) start time.time() C2 np.dot(A, B) # 调用BLAS times_numpy.append(time.time() - start) # 验证结果正确性忽略浮点误差 # print(fn{n}, diff norm: {np.linalg.norm(C1 - C2):.2e}) plt.figure(figsize(10, 6)) plt.plot(sizes, times_standard, o-, labelNaive Triple Loop (O(n^3))) plt.plot(sizes, times_numpy, s-, labelNumPy dot (BLAS)) plt.xlabel(Matrix Size (n)) plt.ylabel(Time (seconds)) plt.yscale(log) # 对数坐标更清晰 plt.xscale(log) plt.title(Theoretical Complexity vs. Practical Performance) plt.legend() plt.grid(True, whichboth, ls--) plt.show() print(即使是最简单的O(n^3)算法在n100时也已经被高度优化的BLAS库彻底击败。) print(原因在于BLAS库利用了CPU的SIMD指令、缓存优化和多线程其常数因子极小。) print(而Strassen算法由于其较大的常数因子和递归开销在n1000时通常也跑不过优化的BLAS。)这个实验清楚地表明在工程实践中常数因子和硬件优化往往比渐近复杂度更重要。除非数据规模巨大否则选择经过充分优化的标准算法库通常是速度最快、最稳定、最省心的方案。6. 超越矩阵乘法优化思想在其他领域的映射减少乘法次数的思想并不局限于矩阵运算。它是一种重要的算法设计范式用廉价操作替换昂贵操作。我们可以将这种思想映射到其他领域卷积神经网络中的Winograd/FFT卷积这正是Winograd算法大放异彩的地方。将时域卷积通过变换转到频域用点乘代替卷积极大减少了计算量。大整数乘法Karatsuba算法类似于Strassen算法。计算两个n位大整数X和Y的乘积。标准方法需要O(n²)次单位乘法。Karatsuba算法发现可以将其分解用3次而不是4次递归乘法来计算复杂度降至O(n^1.585)。多项式乘法通过快速傅里叶变换FFT将多项式系数表示转化为点值表示在点值上进行乘法O(n)再变换回来总复杂度从O(n²)降至O(n log n)。这里的“昂贵操作”是卷积被FFT变换和点乘替代。稀疏矩阵运算当矩阵中绝大多数元素为零时存储和计算都只关注非零元乘法次数自然大幅减少。这不是减少单个乘法的代价而是减少了需要执行的乘法数量。查找表与预计算在图形学、信号处理中对于复杂函数如sin, cos, gamma校正如果输入范围有限且精度要求可接受可以预先计算好所有可能输入的输出值存入查找表。运行时只需一次内存访问廉价代替一次函数计算昂贵可能涉及多个乘加。核心思维模式当你面临一个计算瓶颈时问自己几个问题最耗时的操作是什么Profiling定位这个操作能否用更廉价的操作来近似或等价替换如加法代替乘法能否通过预处理预计算、变换将运行时开销转移到初始化阶段如Winograd的变换能否通过改变计算顺序或数据结构来避免重复计算如动态规划、高斯消元中的缩放掌握这种思维比死记硬背几个特定算法更为重要。它让你在面对新的性能问题时能够自己设计出创造性的优化方案。