有没有一种算法来乘以方阵原地？ [英] Is there an algorithm to multiply square matrices in-place?

问题描述

天真的算法乘以4X4矩阵看起来是这样的：

The naive algorithm for multiplying 4x4 matrices looks like this:

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];
            }
        }
    }
}

显然，这个算法给出虚假的结果，如果出== LHS 或出== RHS （这里 == 表示引用相等）。是否有一个版本，它允许一个并不简单地复制基体的情况下，或两者兼而有之？我很高兴，如果需要有不同的功能对每个案件。

Obviously, this algorithm gives bogus results if out == lhs or out == rhs (here == means reference equality). Is there a version that allows one or both of those cases that doesn't simply copy the matrix? I'm happy to have different functions for each case if necessary.

我发现这文章，但它讨论了施特拉森，威诺格拉德算法，这是矫枉过正的我的小矩阵。来这个问题的答案似乎表明，如果出== LHS和放大器;＆安培;出== RHS （即我们正在试图方形矩阵），那么它就不能要做到位，但即使在那里也没有令人信服的证据或证明。

I found this paper but it discusses the Strassen-Winograd algorithm which is overkill for my small matrices. The answers to this question seem to indicate that if out == lhs && out == rhs (i.e., we're attempting to square the matrix), then it can't be done in place, but even there there's no convincing evidence or proof.

推荐答案

我并不感到这个答案（我张贴这主要是沉默的，它显然不能这样做的人群），但我米怀疑，它可能做一个真正的原地算法好得多（O（1）储存多余的话乘以2 n阶矩阵）。让我们把两个矩阵进行multplied A和B.假设A和B没有别名。

I'm not thrilled with this answer (I'm posting it mainly to silence the "it obviously can't be done" crowd), but I'm skeptical that it's possible to do much better with a true in-place algorithm (O(1) extra words of storage for multiplying two n x n matrices). Let's call the two matrices to be multplied A and B. Assume that A and B are not aliased.

如果A的上三角，则乘法问题是这样的。

If A were upper-triangular, then the multiplication problem would look like this.

[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如下。通过 A11 乘B的第一行。添加的B A12 倍，第二排第一个。添加的B A13 倍，第三排第一位。 B的第四行次添加 A14 第一位。

We can compute the product into B as follows. Multiply the first row of B by a11. Add a12 times the second row of B to the first. Add a13 times the third row of B to the first. Add a14 times the fourth row of B to the first.

现在，我们已经覆盖B的第一行用正确的产品。幸运的是，我们并不需要它了。通过 A22 乘B的第二排。添加的B A23 倍，第三行第二。 （你的想法。）

Now, we've overwritten the first row of B with the correct product. Fortunately, we don't need it any more. Multiply the second row of B by a22. Add a23 times the third row of B to the second. (You get the idea.)

同样，如果是单位下三角，然后乘问题是这样的。

Likewise, if A were unit lower-triangular, then the multiplication problem would look like this.

[ 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的第四第三排。添加的B A42 倍，第二排第四位。添加的B A41 倍，第一排至第四。添加的B A32 倍，第二排第三。 （你的想法。）

Add a43 times to third row of B to the fourth. Add a42 times the second row of B to the fourth. Add a41 times the first row of B to the fourth. Add a32 times the second row of B to the third. (You get the idea.)

完整的算法是LU-à分解到位，到位乘UB到B，乘以LB为B，然后LU-undecompose A（我不知道是否有人曾经这样做，但它似乎很容易扭转的步骤）。大约有一万个理由不能在实践中实现这一点，两个是使A可能不LU分解，而A将不会被完全重建于一般以浮点运算。

The complete algorithm is to LU-decompose A in place, multiply U B into B, multiply L B into B, and then LU-undecompose A in place (I'm not sure if anyone ever does this, but it seems easy enough to reverse the steps). There are about a million reasons not to implement this in practice, two being that A may not be LU-decomposable and that A won't be reconstructed exactly in general with floating-point arithmetic.

这篇关于有没有一种算法来乘以方阵原地？的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！