[科普中国]-施特拉森演算法

施特拉森演算法（Strassen算法）是一种矩阵乘法算法。 它比标准矩阵乘法算法更快，并且在大矩阵的实践中是有用的，但是比最大的矩阵的最快已知算法慢。

Strassen的算法适用于任何环，例如加/乘，但不是所有半环，例如min / plus或布尔代数，其中朴素算法仍然有效，所谓的组合矩阵乘法。

历史Volker Strassen于1969年首次发布该算法，证明了n3通用矩阵乘法算法不是最优的。 Strassen算法稍微好一些，但它的出版物导致了更多关于矩阵乘法的研究，这导致更快的方法，例如Coppersmith-Winograd算法1。

算法介绍矩阵相乘在进行3D变换的时候是经常用到的。在应用中常用矩阵相乘的定义算法对其进行计算。这个算法用到了大量的循环和相乘运算，这使得算法效率不高。而矩阵相乘的计算效率很大程度上的影响了整个程序的运行速度，所以对矩阵相乘算法进行一些改进是必要的。

这里要介绍的矩阵算法称为斯特拉森方法，它是由v.斯特拉森在1969年提出的一个方法。

我们先讨论二阶矩阵的计算方法。

先计算下面7个量(1)

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 = AB。根据矩阵相乘的规则，C的各元素为(2)

c11 = a11 * b11 + a12 * b21

c12 = a11 * b12 + a12 * b22

c21 = a21 * b11 + a22 * b21

c22 = a21 * b12 + a22 * b22

比较(1)(2)，C的各元素可以表示为(3)。

c11 = x1 + x4 - x5 + x7

c12 = x3 + x5

c21 = x2 + x4

c22 = x1 + x3 - x2 + x6

代码实现// 计算2X2矩阵void Multiply2X2(float& fOut_11, float& fOut_12, float& fOut_21, float& fOut_22, float f1_11, float f1_12, float f1_21, float f1_22, float f2_11, float f2_12, float f2_21, float f2_22){ const float x1((f1_11 + f1_22) * (f2_11 + f2_22)); const float x2((f1_21 + f1_22) * f2_11); const float x3(f1_11 * (f2_12 - f2_22)); const float x4(f1_22 * (f2_21 - f2_11)); const float x5((f1_11 + f1_12) * f2_22); const float x6((f1_21 - f1_11) * (f2_11 + f2_12)); const float x7((f1_12 - f1_22) * (f2_21 + f2_22)); fOut_11 = x1 + x4 - x5 + x7; fOut_12 = x3 + x5; fOut_21 = x2 + x4; fOut_22 = x1 - x2 + x3 + x6;}// 计算4X4矩阵void Multiply(CLAYMATRIX& mOut, const CLAYMATRIX& m1, const CLAYMATRIX& m2){ float fTmp[7][4]; // (ma11 + ma22) * (mb11 + mb22) Multiply2X2(fTmp[0][0], fTmp[0][1], fTmp[0][2], fTmp[0][3], m1._11 + m1._33, m1._12 + m1._34, m1._21 + m1._43, m1._22 + m1._44, m2._11 + m2._33, m2._12 + m2._34, m2._21 + m2._43, m2._22 + m2._44); // (ma21 + ma22) * mb11 Multiply2X2(fTmp[1][0], fTmp[1][1], fTmp[1][2], fTmp[1][3], m1._31 + m1._33, m1._32 + m1._34, m1._41 + m1._43, m1._42 + m1._44, m2._11, m2._12, m2._21, m2._22); // ma11 * (mb12 - mb22) Multiply2X2(fTmp[2][0], fTmp[2][1], fTmp[2][2], fTmp[2][3], m1._11, m1._12, m1._21, m1._22, m2._13 - m2._33, m2._14 - m2._34, m2._23 - m2._43, m2._24 - m2._44); // ma22 * (mb21 - mb11) Multiply2X2(fTmp[3][0], fTmp[3][1], fTmp[3][2], fTmp[3][3], m1._33, m1._34, m1._43, m1._44, m2._31 - m2._11, m2._32 - m2._12, m2._41 - m2._21, m2._42 - m2._22); // (ma11 + ma12) * mb22 Multiply2X2(fTmp[4][0], fTmp[4][1], fTmp[4][2], fTmp[4][3], m1._11 + m1._13, m1._12 + m1._14, m1._21 + m1._23, m1._22 + m1._24, m2._33, m2._34, m2._43, m2._44); // (ma21 - ma11) * (mb11 + mb12) Multiply2X2(fTmp[5][0], fTmp[5][1], fTmp[5][2], fTmp[5][3], m1._31 - m1._11, m1._32 - m1._12, m1._41 - m1._21, m1._42 - m1._22, m2._11 + m2._13, m2._12 + m2._14, m2._21 + m2._23, m2._22 + m2._24); // (ma12 - ma22) * (mb21 + mb22) Multiply2X2(fTmp[6][0], fTmp[6][1], fTmp[6][2], fTmp[6][3], m1._13 - m1._33, m1._14 - m1._34, m1._23 - m1._43, m1._24 - m1._44, m2._31 + m2._33, m2._32 + m2._34, m2._41 + m2._43, m2._42 + m2._44); // 第一块 mOut._11 = fTmp[0][0] + fTmp[3][0] - fTmp[4][0] + fTmp[6][0]; mOut._12 = fTmp[0][1] + fTmp[3][1] - fTmp[4][1] + fTmp[6][1]; mOut._21 = fTmp[0][2] + fTmp[3][2] - fTmp[4][2] + fTmp[6][2]; mOut._22 = fTmp[0][3] + fTmp[3][3] - fTmp[4][3] + fTmp[6][3]; // 第二块 mOut._13 = fTmp[2][0] + fTmp[4][0]; mOut._14 = fTmp[2][1] + fTmp[4][1]; mOut._23 = fTmp[2][2] + fTmp[4][2]; mOut._24 = fTmp[2][3] + fTmp[4][3]; // 第三块 mOut._31 = fTmp[1][0] + fTmp[3][0]; mOut._32 = fTmp[1][1] + fTmp[3][1]; mOut._41 = fTmp[1][2] + fTmp[3][2]; mOut._42 = fTmp[1][3] + fTmp[3][3]; // 第四块 mOut._33 = fTmp[0][0] - fTmp[1][0] + fTmp[2][0] + fTmp[5][0]; mOut._34 = fTmp[0][1] - fTmp[1][1] + fTmp[2][1] + fTmp[5][1]; mOut._43 = fTmp[0][2] - fTmp[1][2] + fTmp[2][2] + fTmp[5][2]; mOut._44 = fTmp[0][3] - fTmp[1][3] + fTmp[2][3] + fTmp[5][3];}本词条内容贡献者为:

李嘉骞 - 博士 - 同济大学