是否有一种就地乘以平方矩阵的算法？

一尘不染

是否有一种就地乘以平方矩阵的算法？

algorithm

用于乘以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（此处==表示引用相等），此算法将给出假结果。是否有一个版本允许其中一种或两种情况不只是复制矩阵？如有必要，我很高兴为每种情况提供不同的功能。

我找到了这篇论文，但它讨论了Strassen-
Winograd算法，该算法对我的小型矩阵而言过于严格。这个问题的答案似乎表明，如果out == lhs && out == rhs（即，我们试图对矩阵求平方），那么就无法就位，但是即使没有令人信服的证据也无法证明。

阅读 261

2020-07-28

共1个答案

一尘不染

我对这个答案并不感到兴奋（我主要是为了使“显然不能完成的”人群沉默而发布它），但我怀疑使用真正的就地算法可以做得更好（
O（1）个用于存储两个nxn矩阵的额外存储字）。让我们将两个矩阵相乘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。将a12B的第二行与第一行相乘。将a13B的第三行与第一行相乘。将a14B的第四行与第一行相乘。

现在，我们用正确的乘积覆盖了B的第一行。幸运的是，我们不再需要它。将B的第二行乘以a22。将a23B的第三行与第二行相乘。（你明白了。）

同样，如果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]

将a43时间添加到B的第三行到第四行。将a42B的第二行与第四行相乘。将a41B的第一行与第四行相乘。将a32B的第二行与第三行相乘。（你明白了。）

完整的算法是先对A进行LU分解，再将UB乘以B，再将LB乘以B，然后对LU进行适当的分解（我不确定是否有人这样做，但是似乎很容易将A逆转脚步）。实际上有大约一百万个理由不执行此操作，其中两个原因是A可能无法进行LU分解，并且通常不会使用浮点算术精确地重构A。

2020-07-28