每日算法 — 使用java实现Strassen矩阵乘法:分治策略与矩阵分块加速

算法概述

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算法适用于大规模稠密矩阵乘法场景,如科学计算中的数值模拟、机器学习中的大规模张量运算、计算机图形学中的坐标变换等。对于稀疏矩阵,应使用专门优化的稀疏矩阵乘法算法。