algorithm - 是否有一种算法可以就地乘方阵？

Question

4x4 矩阵相乘的朴素算法如下所示:

void matrix_mul(double out[4][4], double lhs[4][4], double rhs[4][4]) {
    for (int i = 0; i < 4; ++i) {
        for (int j = 0; j < 4; ++j) {
            out[i][j] = 0.0;
            for (int k = 0; k < 4; ++k) {
                out[i][j] += lhs[i][k] * rhs[k][j];
            }
        }
    }
}

显然，如果 out == lhs 或 out == rhs(此处 == 表示引用相等)，此算法会给出虚假结果。是否有一个版本允许其中一种或两种情况不只是简单地复制矩阵？如有必要，我很乐意为每种情况提供不同的功能。

我找到了 this paper 但它讨论了 Strassen-Winograd 算法，这对我的小矩阵来说是过大的。 this的答案问题似乎表明，如果 out == lhs && out == rhs (即，我们试图对矩阵进行平方)，那么它就无法完成，但即使没有有说服力的证据或证明。

Answer 1

我对这个答案并不感到兴奋(我发布它主要是为了让“这显然无法完成”的人群安静下来)，但我怀疑是否有可能用真正的 in-放置算法(O(1) 额外的存储字用于将两个 n x n 矩阵相乘)。我们称这两个矩阵相乘为 A 和 B。假设 A 和 B 没有别名。

如果 A 是上三角矩阵，则乘法问题将如下所示。

[a11 a12 a13 a14] [b11 b12 b13 b14]
[ 0  a22 a23 a24] [b21 b22 b23 b24]
[ 0   0  a33 a34] [b31 b32 b33 b34]
[ 0   0   0  a44] [b41 b42 b43 b44]

我们可以按如下方式将乘积计算为 B。将 B 的第一行乘以 a11。将 B 的第二行添加到第一行的 a12 次。将 B 的第三行的 a13 次添加到第一行。将 a14 乘以 B 的第四行到第一行。

现在，我们已经用正确的产品覆盖了 B 的第一行。幸运的是，我们不再需要它了。将 B 的第二行乘以 a22。将 a23 乘以 B 的第三行到第二行。 (你明白了。)

同样，如果 A 是单位下三角矩阵，则乘法问题将如下所示。

[ 1   0   0   0 ] [b11 b12 b13 b14]
[a21  1   0   0 ] [b21 b22 b23 b24]
[a31 a32  1   0 ] [b31 b32 b33 b34]
[a41 a42 a43  1 ] [b41 b42 b43 b44]

在B的第三行到第四行添加a43次。将 B 的第二行添加到第四行的 a42 次。将 B 的第一行添加到第四行的 a41 次。将 a32 乘以 B 的第二行到第三行。 (你明白了。)

完整的算法是原地 LU 分解 A，将 U B 乘以 B，将 L B 乘以 B，然后原地 LU 分解 A(我不确定是否有人这样做过，但看起来很简单反转步骤)。有大约一百万个理由不在实践中实现这一点，其中两个原因是 A 可能不是 LU 可分解的，并且通常不会使用浮点运算准确地重构 A。