矩阵计算(Matrix Computations) 1.3~1.4

今天接着总结Matrix Computations的第一章. 这次总结一下1.3, 1.4节.

1.3 Block Matrices and Algorithms

这一节来看分块矩阵. 实际上, 分块矩阵就类同于把一般的数量矩阵的元素换成行列适当的矩阵. 简单的一些分块矩阵运算需要注意的就是矩阵作为元素, 它们间的乘法是没有交换律的, 并且要按照矩阵乘法规则来.

Kronecker

A=B⊗C=[b11b12b21b22b31b32]⊗[c11c12c13c21c22c23c31c32c33]A=B\otimes C=\left[ \begin{matrix} b_{11} & b_{12}\\ b_{21} & b_{22}\\ b_{31} & b_{32} \end{matrix}\right] \otimes \left[\begin{matrix} c_{11} & c_{12} & c_{13}\\ c_{21} & c_{22} & c_{23}\\ c_{31} & c_{32} & c_{33} \end{matrix}\right] A=B⊗C=⎣⎡b11b21b31b12b22b32⎦⎤⊗⎣⎡c11c21c31c12c22c32c13c23c33⎦⎤
意味着
A=[b11Cb12Cb21Cb22Cb31Cb32C]A=\left[ \begin{matrix} b_{11}C & b_{12}C\\ b_{21}C & b_{22}C\\ b_{31}C & b_{32}C \end{matrix} \right] A=⎣⎡b11Cb21Cb31Cb12Cb22Cb32C⎦⎤
这个计算的规模是瞬间由 Rn1×m1\mathbb{R}^{n_1\times m_1}Rn1×m1 和 Rn2×m2\mathbb{R}^{n_2\times m_2}Rn2×m2 增到 Rn1n2×m1m2\mathbb{R}^{n_1n_2\times m_1m_2}Rn1n2×m1m2, 在实际运算中处理的问题规模稍微大些就会很占空间. 怎么搞呢? 在一类带Kronecker积的运算中, 有一个用矩阵乘法来减小运算量的办法.

定义一个运算 vecvecvec. 对 X∈Rm×nX\in\mathbb{R}^{m\times n}X∈Rm×n,
vec(X)=[X(:,1)⋮X(:,n)]vec(X)=\left[ \begin{matrix} X(:,1)\\ \vdots\\ X(:,n) \end{matrix} \right] vec(X)=⎣⎢⎡X(:,1)⋮X(:,n)⎦⎥⎤
就是说把矩阵所有的列放在一起组成一个新的列向量.

给出B∈Rm1×n1B\in\mathbb{R}^{m_1\times n_1}B∈Rm1×n1, C∈Rm2×n2C\in \mathbb{R}^{m_2\times n_2}C∈Rm2×n2, X∈Rn2×n1X\in \mathbb{R}^{n_2\times n_1}X∈Rn2×n1我们有下面的关系
Y=CXBT⇔vec(Y)=(B⊗C)vec(X)Y=CXB^T\Leftrightarrow vec(Y) = (B\otimes C)vec(X) Y=CXBT⇔vec(Y)=(B⊗C)vec(X)
左边的运算量是 O(n3)O(n^3)O(n3), 右边是 O(n4)O(n^4)O(n4). 这个很牛逼.

上面定义的 vecvecvec 运算还可以用了定义 reshapereshapereshape. 如果 A∈Rm×nA\in \mathbb{R}^{m\times n}A∈Rm×n, 并且 m1n1=mnm_1n_1=mnm1n1=mn, 那么 B=reshape(A,m1,n1)B=reshape(A, m_1, n_1)B=reshape(A,m1,n1) 可以理解为
A→vec(A)→BA\to vec(A)\to B A→vec(A)→B
先把所有列向量排成一列, 再按长度 m1m_1m1 来切片成 n1n_1n1 个列向量组成矩阵 BBB.

用这两个运算来组合, 就可以把很多带Kronecker积的运算加速了.

Kronecker积的一些性质

(B⊗C)T=BT⊗CT(B\otimes C)^T=B^T\otimes C^T(B⊗C)T=BT⊗CT;
(B⊗C)(D⊗F)=BD⊗CF(B\otimes C)(D\otimes F)=BD\otimes CF(B⊗C)(D⊗F)=BD⊗CF;
(B⊗C)−1=B−1⊗C−1(B\otimes C)^{-1}=B^{-1}\otimes C^{-1}(B⊗C)−1=B−1⊗C−1;
B⊗(C⊗D)=(B⊗C)⊗DB\otimes (C\otimes D)=(B\otimes C)\otimes DB⊗(C⊗D)=(B⊗C)⊗D.

一般, 交换律是不存在的. 但是利用之前说的三类特殊交换矩阵之一: perfect shuffle, 对 B∈Rm1×n1B\in\R^{m_1\times n_1}B∈Rm1×n1, C∈Rm2×n2C\in\R^{m_2\times n_2}C∈Rm2×n2, 我们有下面的式子:
P(B⊗C)QT=C⊗BP(B\otimes C)Q^T=C\otimes B P(B⊗C)QT=C⊗B
这里的 P=Pm1,m2P=\mathcal{P}_{m_1,m_2}P=Pm1,m2, Q=Pn1,n2Q=\mathcal{P}_{n_1,n_2}Q=Pn1,n2.

复数矩阵的乘法

复数矩阵的运算实际上可以用实数矩阵运算来处理. 加减自不必说, 说下乘法(update).
C1+iC2=(C1+iC2)+(A1+iA2)(B1+iB2)C_1+iC_2=(C_1+iC_2)+(A_1+iA_2)(B_1+iB_2) C1+iC2=(C1+iC2)+(A1+iA2)(B1+iB2)
这个运算实际上可以用如下的实数矩阵运算来处理
[C1C2]=[C1C2]+[A1−A2A2A1][B1B2].\newcommand{\lmrm}[1]{\left[\begin{matrix}#1\end{matrix}\right]} \lmrm{C_1\\ C_2}=\lmrm{C_1\\ C_2}+\lmrm{A_1 & -A_2\\ A_2 & A_1}\lmrm{B_1\\ B_2}. [C1C2]=[C1C2]+[A1A2−A2A1][B1B2].

Strassen矩阵乘法

这是一个在大规模运算里非常好用的算法. 我们先用两个一般的二阶方阵间的乘法来说明.
B=[b11b12b21b22],C=[c11c12c21c22],A=[a11a12a21a22]=BC.B=\left[\begin{matrix} b_{11} & b_{12}\\ b_{21} & b_{22} \end{matrix}\right], C=\left[\begin{matrix} c_{11} & c_{12}\\ c_{21} & c_{22} \end{matrix}\right], A=\left[\begin{matrix} a_{11} & a_{12}\\ a_{21} & a_{22} \end{matrix}\right]=BC. B=[b11b21b12b22],C=[c11c21c12c22],A=[a11a21a12a22]=BC.
这个矩阵乘法对应了8次乘法和4次加法:
a11=b11c11+b12c21a12=b11c12+b12c22a21=b21c11+b22c21a22=b21c12+b22c22\begin{aligned} &a_{11}=b_{11}c_{11}+b_{12}c_{21}\\ &a_{12}=b_{11}c_{12}+b_{12}c_{22}\\ &a_{21}=b_{21}c_{11}+b_{22}c_{21}\\ &a_{22}=b_{21}c_{12}+b_{22}c_{22} \end{aligned} a11=b11c11+b12c21a12=b11c12+b12c22a21=b21c11+b22c21a22=b21c12+b22c22
可以用更多的加法(不到乘法的程度)减少一次乘法, 这就是Strassen算法. 首先, 创建10个中间量
S1=c12−c22,S2=b11+b12,S3=b21+b22,S4=c21−c11,S5=b11+b22,S6=c11+c22,S7=b12−b22,S8=c21+c22,S9=b11−b21,S10=c11+c12.\begin{aligned} & S_1=c_{12}-c_{22}, & S_2=b_{11}+b_{12},\\ & S_3=b_{21}+b_{22}, & S_4=c_{21}-c_{11},\\ & S_5=b_{11}+b_{22}, & S_6=c_{11}+c_{22},\\ & S_7=b_{12}-b_{22}, & S_8=c_{21}+c_{22},\\ & S_9=b_{11}-b_{21},&S_{10}=c_{11}+c_{12}. \end{aligned} S1=c12−c22,S3=b21+b22,S5=b11+b22,S7=b12−b22,S9=b11−b21,S2=b11+b12,S4=c21−c11,S6=c11+c22,S8=c21+c22,S10=c11+c12.
再用这十个量做7次乘法
P1=b11S1,P2=S2c22,P3=S3c11,P4=b22S4,P5=S5S6,P6=S7S8,P7=S9S10\begin{aligned} & P_1=b_{11}S_1, & P_2=S_2c_{22},\\ & P_3=S_3c_{11}, & P_4=b_{22}S_4,\\ & P_5=S_5S_6, & P_6=S_7S_8, \\ & P_7=S_9S_{10} \end{aligned} P1=b11S1,P3=S3c11,P5=S5S6,P7=S9S10P2=S2c22,P4=b22S4,P6=S7S8,
再用这7个量做加法得到 AAA:
a11=P1+P4−P5+P7;a12=P3+P5;a21=P2+P4;a22=P1+P3−P2+P6.\begin{aligned} & a_{11}=P_1+P_4-P_5+P_7; \\ & a_{12}=P_3+P_5;\\ & a_{21}=P_2+P_4;\\ & a_{22}=P_1+P_3-P_2+P_6. \end{aligned} a11=P1+P4−P5+P7;a12=P3+P5;a21=P2+P4;a22=P1+P3−P2+P6.
初看感觉费这么大功夫, 才减少一次乘法, 得不偿失啊. 但是如果把这个算法用到分块矩阵上, 它提升的效率就很客观了. 能把矩阵乘法复杂度降到 O(n2.807)O(n^{2.807})O(n2.807) 阶.

1.4 Fast Matrix-Vector Products

还没写完

快速傅里叶变换

快速三角变换