理解 strassen 方法的递归算法

Question

所以，我正在尝试找出 strassen 的矩阵乘法方法，我使用的是 C++，但它可以是任何语言。现在，它看起来像：

 typedef vector<long int> ROW;
 typedef vector<ROW> MATRIX;

 void getQuad(const MATRIX& IN, MATRIX& OUT0, MATRIX& OUT1
                    MATRIX& OUT2, MATRIX& OUT3)
 { /*determine quadrants*/ }

 void strassen(const MATRIX& A, const MATRIX& B, MATRIX& C
 {
      if (A.size() == 2 && A[0] == 2) //know that its 2x2, stop
      {
           // Get M1-M7 vars and set MATRIX C with them
      }
      else
      {
           /*
             getQuad(...) returns the quadrants
             ___________
             | X0 | X1 |
             -----------
             | X2 | X3 |
             -----------
           */

        MATRIX A0,A1,A2,A3;
        getQuad(A,A0,A1,A2,A3);

        MATRIX B0,B1,B2,B3;
        getQuad(B,B0,B1,B2,B3);
      }
 }

我不知道下一步如何处理各个象限，即此时如何推导 M1-M7 矩阵。我想 M1-M7 矩阵（与基本情况下的原始数据类型相反）将以与基本情况相同的方式使用。我只是不确定这里的解散会是什么样子。

我知道阅读别人的代码有点困难，但希望它已经说清楚了。

我确信我的基本情况是正确的，并且我确信我正确地拆分了矩阵，我只是不清楚下一步该去哪里。也许我的算法写错了。

Answer 1

我相信你错过了 Strassen 算法的要点——它是recursive。在伪代码中，算法是这样的：

MATRIX strassen(const MATRIX&a, const MATRIX&b) {
    int aw = a.width();
    int ah = a.height();
    int bw = b.width();
    int bh = b.height();

    if (aw != bh)
        throw some_exception();

    // Strassen algorithm requires each size to be a power of 2
    int max_size = max(aw, ah, bw);
    int extended_size = next_pow_2(max_size);
    MATRIX aEx = a.extend(extended_size, extended_size);
    MATRIX bEx = a.extend(extended_size, extended_size);
    MATRIX cEx = strassenImpl(aEx, bEx);

    // truncate back from power of 2 to real size
    return cEx.truncate(ah, bw);
}


MATRIX strassenImpl(const MATRIX&A, const MATRIX&B) {
    // if matrix size is relatively small it is faster to do the usual straightforward multiplication
    if (A.size() <= threshold) {
        return usualMultiply(A, B);
    }
    // alternatively threshold is 1 so matrix multiplication is just multiplication of the single values
    //if (A.size() == 1) {
    //    return MATRIX(A[0][0]*B[0][0]);
    //} 
    else {
        MATRIX A11, A12, A21, A22;
        getQuad(A, A11, A12, A21, A22);

        MATRIX B11, B12, B21, B22;
        getQuad(B, B11, B12, B21, B22);

        // recursive calls, note that we don't need to go through the extension step
        // here because if the size is a power of 2, half of the size is also a power of 2
        MATRIX M1 = strassenImpl(A11 + A22, B11 + B22);
        MATRIX M2 = strassenImpl(A21 + A22, B11);
        MATRIX M3 = strassenImpl(A11, B12 - B22);
        MATRIX M4 = strassenImpl(A22, B21 - B11);
        MATRIX M5 = strassenImpl(A11 + A12, B22);
        MATRIX M6 = strassenImpl(A21 - A11, B11 + B12);
        MATRIX M2 = strassenImpl(A12 - A22, B21 + B22);

        MATRIX C11 = M1 + M4 - M5 + M7;
        MATRIX C12 = M3 + M5;           
        MATRIX C21 = M2 + M4;
        MATRIX C22 = M1 - M2 + M3 + M6;

        MATRIX C = buildFromQuads(C11, C12, C21, C22);
        return C;
    }
}