算法概述
Strassen矩阵乘法算法是Volker Strassen于1969年提出的革命性分治算法。传统矩阵乘法的时间复杂度为O(n³),而Strassen算法通过巧妙的代数变换,将两个2×2矩阵的乘法从8次标量乘法降至7次,从而将整体复杂度降低到O(n^log₂7)≈O(n^2.807)。这一突破首次证明矩阵乘法可以在低于O(n³)的时间内完成,是算法史上重要的里程碑。
核心思想
传统分治将n×n矩阵划分为4个n/2×n/2的子矩阵:
A = | A11 A12 | B = | B11 B12 |
| A21 A22 | | B21 B22 |
常规方法需要8次子矩阵乘法和4次加法。Strassen的关键在于通过7个中间矩阵M₁~M₇的构造,仅用7次乘法和18次加法/减法完成计算:
- M₁ = (A11 + A22) × (B11 + B22)
- M₂ = (A21 + A22) × B11
- M₃ = A11 × (B12 – B22)
- M₄ = A22 × (B21 – B11)
- M₅ = (A11 + A12) × B22
- M₆ = (A21 – A11) × (B11 + B12)
- M₇ = (A12 – A22) × (B21 + B22)
结果子矩阵:
– C11 = M₁ + M₄ – M₅ + M₇
– C12 = M₃ + M₅
– C21 = M₂ + M₄
– C22 = M₁ – M₂ + M₃ + M₆
Java实现
import java.util.Arrays;
/**
* Strassen矩阵乘法算法实现
* 时间复杂度:O(n^2.807)
* 空间复杂度:O(n^2)
*/
public class StrassenMatrixMultiplication {
// 当矩阵规模小于此阈值时,切换为普通矩阵乘法
private static final int THRESHOLD = 64;
/**
* 标准矩阵乘法(三层循环)
*/
public static int[][] standardMultiply(int[][] A, int[][] B) {
int n = A.length;
int[][] C = new int[n][n];
for (int i = 0; i < n; i++) {
for (int k = 0; k < n; k++) {
for (int j = 0; j < n; j++) {
C[i][j] += A[i][k] * B[k][j];
}
}
}
return C;
}
/**
* Strassen矩阵乘法入口
*/
public static int[][] strassenMultiply(int[][] A, int[][] B) {
int n = A.length;
// 确保矩阵维度为2的幂次,不足则补零
int size = 1;
while (size < n) size <<= 1;
int[][] paddedA = padMatrix(A, size);
int[][] paddedB = padMatrix(B, size);
int[][] paddedC = strassenRecursive(paddedA, paddedB);
// 截取原始维度
return unpadMatrix(paddedC, n);
}
/**
* 递归实现Strassen算法
*/
private static int[][] strassenRecursive(int[][] A, int[][] B) {
int n = A.length;
// 小矩阵使用标准乘法,避免递归开销
if (n <= THRESHOLD) {
return standardMultiply(A, B);
}
int half = n / 2;
// 划分子矩阵
int[][] A11 = new int[half][half];
int[][] A12 = new int[half][half];
int[][] A21 = new int[half][half];
int[][] A22 = new int[half][half];
int[][] B11 = new int[half][half];
int[][] B12 = new int[half][half];
int[][] B21 = new int[half][half];
int[][] B22 = new int[half][half];
splitMatrix(A, A11, A12, A21, A22);
splitMatrix(B, B11, B12, B21, B22);
// 计算7个M矩阵
int[][] M1 = strassenRecursive(addMatrix(A11, A22), addMatrix(B11, B22));
int[][] M2 = strassenRecursive(addMatrix(A21, A22), B11);
int[][] M3 = strassenRecursive(A11, subtractMatrix(B12, B22));
int[][] M4 = strassenRecursive(A22, subtractMatrix(B21, B11));
int[][] M5 = strassenRecursive(addMatrix(A11, A12), B22);
int[][] M6 = strassenRecursive(subtractMatrix(A21, A11), addMatrix(B11, B12));
int[][] M7 = strassenRecursive(subtractMatrix(A12, A22), addMatrix(B21, B22));
// 计算结果子矩阵
int[][] C11 = addMatrix(subtractMatrix(addMatrix(M1, M4), M5), M7);
int[][] C12 = addMatrix(M3, M5);
int[][] C21 = addMatrix(M2, M4);
int[][] C22 = addMatrix(addMatrix(subtractMatrix(M1, M2), M3), M6);
// 合并结果
return joinMatrix(C11, C12, C21, C22);
}
/**
* 矩阵加法
*/
private static int[][] addMatrix(int[][] A, int[][] B) {
int n = A.length;
int[][] C = new int[n][n];
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
C[i][j] = A[i][j] + B[i][j];
}
}
return C;
}
/**
* 矩阵减法
*/
private static int[][] subtractMatrix(int[][] A, int[][] B) {
int n = A.length;
int[][] C = new int[n][n];
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
C[i][j] = A[i][j] - B[i][j];
}
}
return C;
}
/**
* 将矩阵划分为4个子矩阵
*/
private static void splitMatrix(int[][] parent, int[][] c11, int[][] c12,
int[][] c21, int[][] c22) {
int half = parent.length / 2;
for (int i = 0; i < half; i++) {
System.arraycopy(parent[i], 0, c11[i], 0, half);
System.arraycopy(parent[i], half, c12[i], 0, half);
System.arraycopy(parent[i + half], 0, c21[i], 0, half);
System.arraycopy(parent[i + half], half, c22[i], 0, half);
}
}
/**
* 将4个子矩阵合并为一个大矩阵
*/
private static int[][] joinMatrix(int[][] c11, int[][] c12,
int[][] c21, int[][] c22) {
int half = c11.length;
int n = half * 2;
int[][] result = new int[n][n];
for (int i = 0; i < half; i++) {
System.arraycopy(c11[i], 0, result[i], 0, half);
System.arraycopy(c12[i], 0, result[i], half, half);
System.arraycopy(c21[i], 0, result[i + half], 0, half);
System.arraycopy(c22[i], 0, result[i + half], half, half);
}
return result;
}
/**
* 将矩阵补零扩展至指定大小
*/
private static int[][] padMatrix(int[][] matrix, int size) {
int n = matrix.length;
int[][] padded = new int[size][size];
for (int i = 0; i < n; i++) {
System.arraycopy(matrix[i], 0, padded[i], 0, n);
}
return padded;
}
/**
* 截取矩阵至指定大小
*/
private static int[][] unpadMatrix(int[][] matrix, int size) {
int[][] result = new int[size][size];
for (int i = 0; i < size; i++) {
System.arraycopy(matrix[i], 0, result[i], 0, size);
}
return result;
}
public static void main(String[] args) {
// 测试用例
int[][] A = {
{1, 2, 3, 4},
{5, 6, 7, 8},
{9, 10, 11, 12},
{13, 14, 15, 16}
};
int[][] B = {
{17, 18, 19, 20},
{21, 22, 23, 24},
{25, 26, 27, 28},
{29, 30, 31, 32}
};
int[][] standard = standardMultiply(A, B);
int[][] strassen = strassenMultiply(A, B);
System.out.println("标准乘法结果:");
printMatrix(standard);
System.out.println("Strassen结果:");
printMatrix(strassen);
System.out.println("结果一致:" + Arrays.deepEquals(standard, strassen));
}
private static void printMatrix(int[][] matrix) {
for (int[] row : matrix) {
System.out.println(Arrays.toString(row));
}
}
}
复杂度分析
| 算法 | 时间复杂度 | 空间复杂度 |
|---|---|---|
| 标准矩阵乘法 | O(n³) | O(n²) |
| Strassen算法 | O(n^2.807) | O(n²) |
Strassen算法的递归关系式为T(n) = 7T(n/2) + O(n²),根据主定理可得时间复杂度为O(n^log₂7)≈O(n^2.807)。实际应用中,当矩阵规模小于64~128时,应切换为标准乘法以避免递归和额外加减法带来的开销。
应用场景
Strassen算法适用于大规模稠密矩阵乘法场景,如科学计算中的数值模拟、机器学习中的大规模张量运算、计算机图形学中的坐标变换等。对于稀疏矩阵,应使用专门优化的稀疏矩阵乘法算法。