有没有一种算法可以就地增加平方矩阵？

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(这里==指的是参考相等).是否有允许这些案例中的一个或两个不仅仅复制矩阵的版本？如果有必要,我很高兴为每个案例提供不同的功能.

我发现了这篇论文,但它讨论了Strassen-Winograd算法,这对我的小矩阵来说太过分了.这个问题的答案似乎表明,如果out == lhs && out == rhs(即,我们试图对矩阵进行平方),那么它就无法在适当的位置完成,但即使在那里也没有令人信服的证据或证据.

Answer 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.a12将第二行B 添加到第一行.a13将第三行B 添加到第一行.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]

将a43时间添加到第三行B到第四行.a42将第二行B 添加到第四行.a41将第一行B 添加到第四行.a32将第二行B 添加到第三行.(你明白了.)

完整的算法是LU分解A到位,将UB乘以B,将LB乘以B,然后LU-undecompose A到位(我不确定是否有人这样做,但似乎很容易扭转脚步).有大约一百万个理由不在实践中实现这一点,其中两个原因是A可能不是LU可分解的,并且A不会通过浮点算法精确地重建.

Answer 2

这个答案比我的另一个更明智,尽管它使用了一整列额外的存储空间,并且具有与原始复制算法相同的数据移动量.要将A与B相乘,将产品存储在B中(再次假设A和B分开存储):

For each column of B,
    Copy it into the auxiliary storage column
    Compute the product of A and the auxiliary storage column into that column of B

我切换伪代码首先进行复制,因为对于大型矩阵,缓存效果可能导致将A乘以连续的辅助列而不是B中的非连续条目更有效.