有没有办法减少这两个依赖矩阵之间的运算？答案

【问题标题】：Is there any way to reduce the operations between these two dependent matrices?有没有办法减少这两个依赖矩阵之间的运算？
【发布时间】：2017-04-22 17:44:30
【问题描述】：

我有两个矩阵，其中第一个矩阵是从第二个矩阵的数据创建的，用于转换第二个矩阵。我重复这个操作很多次。由于这两个矩阵的依赖关系，我一直没能找到一种方法来加快这里的操作。我会尝试用小矩阵大小向你展示我在说什么。

[ 1 0 0     0     ]     [ .11 .22 .33 .44 ]
[ 0 1 0     0     ]     [ .65 .42 .01 .92 ]
[ 0 0 .51^2 .85^2 ]  *  [ .31 .15 .51 .85 ]
[ 0 0 .44^2 .23^2 ]     [ .25 .78 .44 .23 ]
        A                        B

假设我执行此操作数百万次，计算值在 A 中的位置取决于我希望在 B 中发生旋转的位置。因此，在每次迭代中，矩阵 A 和 B 都不同，并且使用的值计算要放在 A 中的新值是不同的。

有谁知道加速这类代码的方法？考虑到矩阵乘法本质上是向量矩阵乘法，我希望创建一个组合 A（展开到 A 是一个完整矩阵的点，或者足够充分以利用 MMM 算法），但是新值的数据依赖性使看来我可能被卡住了。我得到这样的东西：

A * B = B'
A' * B' = B''
A'' * B'' = B'''

其中 A' 派生自 B'，A'' 派生自 B''，依此类推。

编辑：

为了澄清，第二轮可能是：

[ 1 0     0     0 ]     [ .11 .22 .33 .44 ]
[ 0 .65^2 .42^2 0 ]     [ .65 .42 .01 .92 ]
[ 0 .26^2 .60^2 0 ]  *  [ .26 .60 .45 .39 ]
[ 0 0     0     1 ]     [ .06 .04 .10 .17 ]
        A'                       B'

【问题讨论】：

Will A, A', A'', ... 的形式总是 A = [I, 0; 0，X]？如果是这样，您应该简化问题以删除 B 的上半部分以及 A 的下对角块以外的所有部分（例如，B' = [B_top; X*B_bottom]）。
不，A 并不总是如图所示。四个旋转元素可以在对角线上的任何位置。
矩阵的大小总是 4x4 还是更大？
此外，是否会有 2 行/列恰好在对角线上为 1，或者这个数字可以增加吗？
@NoseKnowsAll 两个矩阵都非常大。我只是使用一个较小的例子。 A 开始是一个非常大的单位矩阵，我写入了我想用来变换矩阵 B 的四个值。这四个值总是跨越对角线，如图所示。

标签： c matrix matrix-multiplication

【解决方案1】：

当矩阵A[rows][num]乘以矩阵B[num][cols]时，结果矩阵C[rows][cols]中的每个元素是

C[r][c] = A[r][0]*B[0][c] + ... + A[r][num-1]*B[num-1][c]
        = sum( A[r][k] * B[k][c], k = 0 .. num-1 )

这里，A 是一个单位矩阵，除了两行（r0 和 r1），每行在两列中有两个值（c0 和 c1）。因为行r0 和r1 否则为零，因此这些行上结果矩阵C 中的每个元素都是两个乘积的总和：

C[r0][i] = A[r0][c0] * B[c0][i] + A[r0][c1] * B[c1][i]
C[r1][i] = A[r1][c0] * B[c0][i] + A[r1][c1] * B[c1][i]

由于A 是一个单位矩阵，因此结果中的所有其他行都与B 矩阵中的相应行相同。

如果我们使用

a00 = A[r0][c0]        a01 = A[r0][c1]
a10 = A[r1][c0]        a11 = A[r1][c1]

那么我们可以将更新循环描述为

For i = 0 .. cols-1
    C[r0][i] = a00 * B[c0][i] + a01 * B[c1][i]
    C[r1][i] = a10 * B[c0][i] + a11 * B[c1][i]
End for

如果我们在更新循环中使用两个临时变量（以c0 == r0 或c1 == r0 为例，这样我们就不会在计算[r1][i] 处的元素之前覆盖[r0][i] 处的元素），整个操作可以就地完成，只使用B 矩阵。

如果左边矩阵A中的非恒等元素是矩阵B中的对应元素的平方，我们就地更新，就变得很简单了。在 C99 中：

#include <stdlib.h>

void update(const size_t n, double B[n][n], const size_t r, const size_t c)
{
    const double a00 = B[r  ][c  ] * B[r  ][c  ];
    const double a01 = B[r  ][c+1] * B[r  ][c+1];
    const double a10 = B[r+1][c  ] * B[r+1][c  ];
    const double a11 = B[r+1][c+1] * B[r+1][c+1];
    size_t       i;
    for (i = 0; i < n; i++) {
        const double b0i = B[c  ][i];
        const double b1i = B[c+1][i];
        B[r  ][i] = a00 * b0i + a01 * b1i;
        B[r+1][i] = a10 * b0i + a11 * b1i;
    }
}

上述操作的计算成本是4*n+4 乘法和2*n 加法（矩阵元素），即O(N)。

这种方法也适用于Givens 和Jacobi 旋转，只是不是使用给定的矩阵计算四个元素，而是通过传递给函数的参数计算它们；并且需要传递两个单独的行和两个列参数，而不是连续的行和列。

【讨论】：

【解决方案2】：

这更像是一个数学问题而不是编程问题，但是是的，有一种方法可以加快您的计算速度。此答案假定您正在乘以大矩阵，但A 是对单位矩阵的轻微修改，仅在位置 (i,i)、(i,i+1)、(i+1,i) 和 (@ 987654325@) 正如你在 cmets 中提到的那样。

请注意，将左侧的B 乘以矩阵A，其中A 的一行除对角线上的1 外全为零，相当于说AB 的第i 行是等于 B 的第 ith 行。这意味着您可以立即在迭代中节省大量计算。

鉴于A 的更新A 具有这些不错的属性，我们有以下算法：

将B' 初始化为精确的B。
执行矩阵乘法update = [A(i,:),A(i+1,:)]*B，其中左侧是(2xn) 矩阵。
将此更新的输出存储在正确的两行B':[B'(i,:),B'(i+1,:)]中。

每次更新都会失败 = O(2n^2)

如果您的大小 n 确实很大，您可能还希望通过利用 A 的两行的行稀疏性来节省这些小矩阵乘法中的更多计算。您可以执行(2x2) 与(2xn) 矩阵乘法，而不是执行(2xn) 与(nxn) 矩阵乘法。

这个新算法将再次更新正确的两行B'。但是，现在您的更新矩阵乘法只需：

update = [A(i  ,i), A(i  ,i+1)] * [B(i  ,:)]
         [A(i+1,i), A(i+1,i+1)]   [B(i+1,:)]

每次更新都会失败 = O(4n)

请注意，这两种算法都不会显式存储所有A，因此您将节省内存和计算量。此外，如果您以行优先格式存储矩阵，则上述两种算法的缓存效率都非常高。

【讨论】：