======Matrix-Vector to Matrix-Matrix Multiplication====== LAFF Week 4 Notes ---- 本章的核心概念主要介绍了如何对矩阵和向量的乘法进行象限分割,并利用这种象限分割优化特定矩阵与向量的乘法的算法。除此之外另外一个比较重要的内容是如何从矩阵向量乘法的角度去解释矩阵与矩阵的乘法。 ====矩阵向量乘法的象限分割==== 通过第三章,我们得知矩阵可以使用象限分割的方法来分割矩阵;同时,我们也知道通过对行或者列进行优先分割,我们可以将矩阵的乘法用点积或者 AXPY 的方式进行计算。那有没有办法将这两者联系到一起呢?来看看下面的修改: ===对 A, x, y 进行分割=== 我们知道在矩阵向量乘法中,每次迭代的最终结果都会加到对应的 $y$ 的分量上,无论是哪一种方式来计算。想象一下之前的象限分割下的各种算法,每一次迭代中 $y$ 的分量作为结果,都是有对应区域的 $A$ 与 $x$ 的分量产生的。因此,我们可以认为在分割向量乘法中,$A$、$x$、$y$ 中的元素是一一对应的。那么我们很容易得出一个结论: $$ \displaystyle \left( \begin{array}{ c | c | c | c } A_{0,0} & A_{0,1} & \cdots & A_{0,N-1} \\ \hline A_{1,0} & A_{1,1} & \cdots & A_{1,N-1} \\ \hline \vdots & \vdots & \ddots & \vdots \\ \hline A_{M-1,0} & A_{M-1,1} & \cdots & A_{M-1,N-1} \end{array} \right) \left( \begin{array}{ c } x_0 \\ \hline x_1 \\ \hline \vdots \\ \hline x_{N-1} \end{array} \right) \\ = \left( \begin{array}{c c c} A_{0,0} x_{0} + A_{0,1} x_{1} + \cdots + A_{0,N-1} x_{N-1} \\ \hline A_{1,0} x_{0} + A_{1,1} x_{1} + \cdots + A_{1,N-1} x_{N-1} \\ \hline \vdots \\ \hline A_{M-1,0} x_{0} + A_{M-1,1} x_{1} + \cdots + A_{M-1,N-1} x_{N-1} \end{array} \right) $$ ===象限分割下的矩阵向量乘法=== 根据上面的推论,我们可以改写一下使用象限分割后的矩阵向量乘法的算法。两种改写后的算法如下,左边是点积观点的,右边是 AXPY 观点的: \\ \\ {{ math:linear_algebra:laff:mvmult_n_un_var1-2b.png?600 |}} \\ \\ ==点积观点的算法改写== 点击观点的算法在一般情况下是行优先的算法,总的来说就是 $A$ 的每一行作为一个整体与 $x$ 做点积,得到的值加到 $y$ 对应分量上。经过象限分割后,我们的三个变量发生一定的变化: * $A$ 的行从整体到被分割成三部分(从左到右) * $x$ 从整体到被分割成三部分(从上到下) * $y$ 保持状态不变(之前的算法中就被分割了) 那么算法的变动实际上就很清楚了,请参考下图: \\ \\
y1 = a1t * x + y1 // M mult V via dot
y1 = a10t * x1 + a11 * chi1 + a12t * x2 + y1; //M mult V via dot in quadrants
==AXPY观点的算法改写==
AXPY 观点下的矩阵向量乘法是列优先的,也就是用 $x$ 的分量来乘以对应矩阵的列,再将结果加到整个 $y$ 上。我们也来看看经过象限分割后,三个变量发生了什么样的变化:
* $A$ 在行中也有分割,因此之前作为单位的列分割成了三部分 ((a01, a11,a21)t)
* $x$ 之前就有从上到下的分割,保持原样
* $y$ 不再是一个整体,而是从上到下分割成了三部分((y0, psi,y2)t)
算法的变动可以参考下图:
\\
\\
y = chi * a1 + y; //M mult V via AXPY
//--M mult V via AXPY in quadrants---//
y1 = chi * a01 + y1;
psi = chi * a11 + psi;
y2 = chi * a21 + y2;
===Transposing a Partitioned Matrix===
对分割后的矩阵转置分二步:
- 按区域为单位,对整个矩阵进行转置。区域内部的不转置
- 对区域内部的单位分别进行转置。
====特殊矩阵与矩阵向量乘法====
===转置矩阵的矩阵向量乘法===
假设我们要求 $y = A^Tx +y$。通常情况下,我们会先求 $A$ 的转置矩阵,再与 $x$ 相乘,最后累加到 $y$ 上。然而,在计算机中,如果我们有了转置的中间步骤,我们就需要找一块空间来存储 $A$ 的转置结果;更坏的是,为了转置 $A$,我们很可能需要遍历整个矩阵 $A$。因此,我们需要设计一种算法来避免这些情况的发生,该算法的核心就是避免转置操作的发生。
\\
\\
实际上,不管是通过点积观点还是AXPY观点,我们都会发现一个现象:那就是原矩阵的列就是转置矩阵的行,原矩阵的行就是转置矩阵的列。因此,我们只需要通过点积运算,就可以将当前的行(列)作为列(行)进行运算,即:
* 如果是点积观点的乘法,那么是**行优先**,那么需要按**列**划分 $A$, 然后做 $x$ 与 $A$ 中**列**的点积,即$A_1^Tx$;
* 如果是AXPY观点的乘法,那么是**列优先**,那么只需要按**行**划分 $A$,然后做 $x$ 与 $A$ 中**行**的AXPY 操作,即 $A_1x$(因为按行划分,因此原有的**行**为 $A_1^T$)。
算法参考如下:
\\
\\
{{ math:linear_algebra:laff:4311_algorithm.png?600 |}}
\\
\\
因此,只要通过替换分割的方向,我们就能在不修改原矩阵的情况下利用点积来求转置矩阵与向量的乘法。
===三角矩阵与矩阵向量乘法===
三角矩阵是一个很比较特殊的存在。在三角矩阵中,有一半的元素都是 0,而 0 意味着我们没有必要在这些 0 元素上浪费时间进行计算。因此在设计算法的时候,我们必须要考虑如何使用迭代将 0 的区域划分出来,并**舍弃**这一部分本应存在于一般矩阵与向量的乘法中的运算。
\\
\\
==从点积角度来看==
**上三角矩阵**:
\\
\\
以上三角矩阵为例,我们来看一看在正常的点积矩阵向量乘法中,属于上三角矩阵的 0 会在每一次迭代中处于什么样的位置:
\\
\\
y1 = a10t * x1 + a11 * chi1 + a12t * x2 + y1; //M mult V via dot in quadrants
y1 = u11 * chi1 + u12t * x2 +y1; //U mult V via dot in quadrants
**下三角矩阵**:
\\
\\
反过来看,如果是下三角矩阵,则是 $L_{12}^T$ 的部分在每次迭代中(图中绿色的区域)都为 0,因此下三角矩阵只需要将 $L_{10}^T$ 的运算部分去掉即可:
y1 = a10t * x1 + a11 * chi1 + a12t * x2 + y1; //M mult V via dot in quadrants
y1 = l10t * x1 + l11 * chi1 +y1; //L mult V via dot in quadrants
**严格三角矩阵**:
\\
\\
那么对角线为 0,只需要将迭代中的对角线部分($A_{11}$)去掉就好。
\\
\\
**x = Ux:**
\\
\\
以上是求 $y = Ax +y$ 的情况。再来看看 $x = Ux$ 的情况,这种情况下**输入和输出的都是** $x$。那么意味着:
* ''chi1'' 是每次迭代会被覆盖数据的那个分量。
* ''chi1'' 会参与其中的一次运算(点积的一部分)
因此我们需要注意的是在 ''chi1'' 被覆盖之前,先让其参与点积的运算,于是就有:
chi1 = u11 * chi1; // compute dot with original chi1, then overwrite it
chi1 = u12t * x2 + chi1;
**x = Lx:**
\\
\\
不过这里又会出现一个问题:上面的做法对于上三角矩阵是没有问题的;但对于下三角矩阵来说就会出问题了:对于下三角矩阵,舍去的是 $L_{12}^T$ 部分,$L_{10}^T$ 是需要先被计算的(遍历的顺序)。但我们不能将其存到 ''chi1'' 里,因为我们接下来还要对 $U_{11}\chi_1$ 进行计算,如果覆盖了 ''chi1'',那么接下来的结果就不正确了。那应该怎么办呢?
\\
\\
很简单,我们只需要将其作为一个 “上三角矩阵” 来处理就可以了;也就是说,我们可以从**右下角**开始我们的遍历,这样之前的 $L_{10}^T$ 就变成了 “$L_{12}^T$” 了,计算的第一位就变成了 $L_{11} \chi_1$ 了。
因此 $x = Lx$ 可以写成:
chi1 = u11 * chi1; // compute dot with original chi1, then overwrite it
chi1 = u10t * x0 + chi1; //这里的u10t / x0 都是从右下方开始计数的,因此遍历位置处于 u11 的后面
\\
\\
**y = U^Tx + y:**
\\
\\
除了以上情况,在象限分割下我们也能更好的处理转置矩阵与向量的乘法。我们知道在点积观点下进行转置矩阵的运算都是将行看成列然后与 $x$ 进行点积的;鉴于上三角矩阵对角线下方都是 0,那么只需要对对角线元素 $U_{11}$ 与对角线上方元素 $U_{01}$ 与对应 $x$ 分量进行点积就可以了:
psi1 = laff_dots(upsilon11, chi1, psi1);
psi1 = laff_dots(u01, x0, psi1);
**x = U^Tx :**
\\
\\
**注意:**该运算与 $x=Ux$ 不同。转置矩阵乘法使用列与 $x$ 做点积,如果从左上方开始遍历,那么县运算会是 $U_{01}$。由于 $U_{01}$ 的运算会覆盖 $\chi_1$ 的结果,因此会导致 $\chi_1$ 的值错误。所以该运算需要从**右下角开始遍历**(正好与 $x=Ux$ 相反)
chi1 = laff_dot(upsilon11, chi1);
chi1 = laff_dots(u01, x0, chi1); // 这里的 u01 实际上是 u01t,但laff_dots支持u01到u01t的点积转换。
**下三角矩阵的转置矩阵 :**:
\\
\\
该运算与向量的乘法与之前类似:只需要把 ''l01'' 替换成 ''l21'' 就可以了:
/*y = L^Tx + y*/
psi1 = laff_dots(lambda11, chi1, psi1);
psi1 = laff_dots(l21, x2, psi1); % l21t in the algorithm but OK for laff_dots()
/*x = L^Tx*/
chi1 = laff_dot(lambda11, chi1);
chi1 = laff_dots(l21, x2, chi1);
**注意:**相比起 $x=Lx$,$x = L^Tx$ 的算法中是从上到下分割(也就是先算 $L_{11}$,再算$L_{21}$)因此
自左上开始的遍历是满足该顺序的。这与 $x=Lx$ 也完全相反。
==从AXPY角度来看==
还是从上三角矩阵谈起。先来看看在 AXPY 算法的迭代中,0 都处于哪些区域:
\\
\\
y0 = u01 * chi1 +y0;
psi1 = u11 * chi1 +psi1;
**下三角矩阵:**
\\
\\
下三角矩阵需要去掉的部分为图中所有绿色的部分,也就是 $L_{01}$ 的部分。因此算法可以改写为:
psi1 = l11 * chi1 +psi1;
y2 = l21 * chi1 +y2;
**严格三角矩阵:**
\\
\\
只需去掉 ''u11/l11'' 部分即可。
\\
\\
**x = Ux**:
对于 AXPY 的迭代部分,''chi1'' 是需要保证在所有运算结束之前不变的。不过在该运算中,采用左上开始的遍历进行运算,并不会在运算结束之前改变 ''chi1'' 的值。
x0 = chi1 * u01 + x0;
chi1 = chi1 * u11;
**x = Lx**:
\\
\\
这个运算与点积观点下的运算类似。如果从左上开始,''chi1'' 与 ''l11'' 的运算会改变 ''chi1'' 本身的值,造成后面的运算出错。因此我们从左下开始遍历,将整个运算换成了上三角矩阵的运算,也就是下三角矩阵里的 ''l21'' 编程了 上三角矩阵里的 ''l01'':
x2 = chi1 * l21 +x2;
chi1 = l11 * chi1;
**转置矩阵:**
\\
\\
AXPY观点下的乘法只需要将**列视作行**进行 AXPY 操作,就可以得到转置矩阵与向量相乘一样的结果了。分别来看看算法的改动:
\\
\\
**y = U^Tx + y:**
\\
\\
之前参与运算的区域是 $U_{01}$,那么对应的转置区域就是 $U_{12}^T$。
psi1 = laff_axpy(chi1, upsilon11, psi1);
y2 = laff_axpy(chi1, u12t, y2);
**y = L^Tx + y:**
之前参与运算的区域是 $L_{21}$,那么对应的专职区域就应该是 $L_{10}^T$。
y0 = laff_axpy(chi1, l10t, y0); // l10 in the algorithm but OK for laff_axpy()
psi1 = laff_axpy(chi1, lambda11, psi1);
**x = U^Tx :**
\\
\\
与点积观点的乘法类似,因为需要先算 ''u12t'' 再算 ''u11'',因此必须从**右下角**开始遍历才能保证 ''u12t'' 在 ''u11'' 区域之前计算:
x2 = laff_axpy(chi1, u12t, x2); // u12 in the algorithm but OK for laff_dots()
chi1 = laff_axpy(chi1, upsilon11, 0);
**x = L^Tx :**
\\
\\
同样与点积观点类似,先算 ''u10t'' 再算 ''u11'',从**左上角**开始遍历正好可以满足该顺序。
y0 = laff_axpy(chi1, l10t, y0); // l10 in the algorithm but OK for laff_axpy()
psi1 = laff_axpy(chi1, lambda11, psi1);
===对称矩阵与矩阵向量乘法===
对称矩阵的元素相对于主对角线对称,因此在对阵矩阵与向量的乘法中,我们只需要利用**一半的矩阵空间**就可以表达整个对称矩阵,并提供足够的信息完成运算。
\\
\\
实现该算法实际上非常简单,只需要利用转置矩阵就可以。通过转置某一个区域,我们能用一个区域存储的内容表达两个区域的内容。下面来看看矩阵内容存储在下部和上部分别需要置换哪些区域:
\\
\\
==内容在下部==
**点积观点:**
\\
\\
内容在下部,意味着 $A_{12}^T$ 区域的内容在对角线以上,替换成 $A_{21}$ 即可:
psi1 = laff_dots( a10t, x0, psi1 );
psi1 = laff_dots( alpha11, chi1, psi1 );
psi1 = laff_dots( a21, x2, psi1 );
**AXPY观点:**
\\
\\
内容在下部,意味着 $A_{01}$ 区域的内容在对角线以上,替换成 $A_{10}^T$ 即可:
psi1 = laff_dots( a01, x0, psi1 );
psi1 = laff_dots( alpha11, chi1, psi1 );
psi1 = laff_dots( a12t, x2, psi1 );
**AXPY观点:**
\\
\\
内容在上部,意味着 $A_{21}$ 区域的内容在对角线以上,替换成 $A_{12}^T$ 即可: