算法导论之矩阵运算

矩阵运算的重要性应该不亚于图算法。先重温下矩阵的相关概念和性质，为后续矩阵运算奠定数据理论基础。

矩阵A：数字的一个矩形阵列，形式化为A=(a_ij)，第i行j列元素为a_ij，如元素为实数的所有元素mXn矩阵组合的元素用R^mXn表示。

矩阵转置A^T：是矩阵A的行和列互相交换而产生的举证。

向量（Vector）：是数字的一维向量，列向量看成是nX1的矩阵，转置成行向量就是1Xn的矩阵。

单位向量e_i：矩阵第一个元素为1而其他元素均为0的常量。

零矩阵：所有元素都是0的矩阵。

对角矩阵：当i≠j时，a_ij=0，所有非对角线上的元素均为0。

nXn单位矩阵I_n：是对角线元素都是1的对角矩阵。

三对角矩阵T：满足|i-j|>1的元素t_ij=0的矩阵，非零元素仅出现在主对角线上、靠主对角线上面和对角线下面。

上三角矩阵U：满足i>j的元素U_ij=0的矩阵，对角线下面的所有元素均为0，若对角线元素为1，则是单位上三角矩阵。

下三角矩阵L：满足i<j的元素L_ij=0的矩阵，对角线上面的所有元素均为0，若对角线元素为1，则是单位下三角矩阵。

置换矩阵P：每一行或列中仅包含一个1，其他元素都为0，可以把一个向量x和一个置换矩阵相乘，结果是向量x中的元素的一种置换。

对称矩阵A：满足条件A= A^T。

矩阵加法：C(c_ij)=A(a_ij)+B(b_ij)，c_ij=a_ij+b_ij，都是mXn矩阵。

矩阵乘法：相容（A的列数等于B的行数）的两个矩阵才可以相乘，C(c_ij)=A(a_ik)+B(b_kj)，矩阵和单位矩阵等于自己、和零矩阵相乘等于零、和向量相乘得到向量，满足结合律和分配律，但不满足交换律。

逆矩阵A^-1：满足A A^-1= I_n= A^-1A，不可逆的矩阵称为奇异矩阵。

线性相关：设x1,x2,…,xn是n个向量，若存在不全为零的常系数c2,c2,…,cn，使c1x1+c2x2+…+cnxn=0成立，就称x1,x2,…,xn线性相关。

非零mXn矩阵A的秩：A的极大线性无关列向量组中向量的个数为列秩，A的极大线性无关行向量组中向量的个数为列秩。任意一个矩阵的行秩和列秩都是相等。mXn的秩是0和min(m,n)之间的一个整数。

非零mXn矩阵A的秩满足下列条件的最小的数r：存在mXr矩阵B和rXn矩阵C且有A=BC。如果一个nXn的矩阵秩为n，则是满秩，如果一个mXn矩阵的秩为n，则是列满秩。

一个方阵满秩当且仅当它为非奇异矩阵。

当前仅当A无空向量时，矩阵A为列满秩。

当且仅当A具有空向量时，方阵A是奇异的。

余子式：对于n>1，nXn矩阵A的第ij个余子式是把A的第i行和第j列元素去掉后所形成的一个(n-1)X(n-1)矩阵A_[ij]，表示为det(A_[ij])。

方阵A的行列式具有如下性质：如果A的任何行或列的元素为0，则det(A)=0；用常如a乘A的行列式任意一行（或任意一列）的诸元素，等于用a乘A的行列式；A的行列式中的一行（或一列）元素加上另一行（或另一列）中的相应元素，行列式的值不变；A的行列式值与其转置矩阵A^T的行列式的值相等；行列式的任意两行（或两列）互换，则其值异号。

正定矩阵：nXn矩阵A，如果对所有n维向量x≠0都有x^TAx>0，则矩阵A为正定矩阵。

对任意列满秩矩阵A，矩阵A^TA是正定的。

1）矩阵乘法的Strassen算法

Strassen算法运用分治法将两个nXn的矩阵乘积运行时间，从简易矩阵乘法算法⊙(n³)提升到⊙(n^lg7)= ⊙(n^2.81)。说明算法的主要思想和步骤。

一般情况下，矩阵乘积C=AB，其中A、B和C都是nXn方阵，假定n是2的幂，把A、B和C划分成四个n/2Xn/2矩阵：

对应四个等式：

r=ae+bg

s=af+bh

t=ce+dg

u=cf+dh

每个等式运行包含两次n/xXn/2矩阵乘法和两次乘积所得的n/2Xn/2矩阵的加法运算，推出nXn矩阵相乘所需时间T(n)的递归式：T(n)=8T(n/2)+⊙(n²)。

Strassen算法在这个一般性基础发现，发现只有进行7次n/xXn/2矩阵乘法即可解决，即T(n)=7T(n/2)+⊙(n²)= ⊙(n^lg7)=⊙(n^2.81)，具体算法步骤如下：

第一：输入矩阵A和B划分为n/2Xn/2的子矩阵；

第二：运行⊙(n²)标量加法和减法运算，计算出14个n/2Xn/2的矩阵A1,B1,A2,B2,…,A7,B7；

第三：递归计算出7个矩阵的乘积P_i=A_iXB_i，i=1,2,…,7；

第四：使用⊙(n²)次标量加法与减法运算，对P_i矩阵的各种组合进行求和或求差运算，从而获得结果矩阵C的四个子矩阵r\s\t\u。

Strassen算法的核心思想是什么呢？何以产生如此一个步骤的算法呢？依据是什么？

设想每个矩阵的积P_i可以写成如下形式：

P_i=A_iB_i=(α_i1a+α_i2b+α_i3c+α_i4d)( β_i1e+β_i2f+β_i3g+β_i4h)，其中系数α和β都属于集合{-1,0,1}，就是说，对矩阵A和B的子矩阵进行加减运算，并对所得的结果进行相乘法来计算出每个子矩阵的乘积。简单说，将子矩阵的乘法运算变成加减运算，从而将8次乘法减少到7次。这是Strassen算法的核心思想和依据。

算法导论中给出的描述很好理解算法的思路和步骤。正如算法导论中对该算法的讨论中所指出的问题，个人对于矩阵乘法的困难（多项式时间内完成的算法）是认同的。Strassen算法并不是最佳，但目前似乎没有更多的突破。

2）求解线性方程组（LUP分解思想）

求解一组同时成立的线性方程式在很多应用中都会出现，也是线性代数的基本问题。可以将一个线性系统表示为一个矩阵方程，其中每个矩阵或向量元素都属于一个域，如实数域R。

n个未知量的n个方程用矩阵表示如下：

设A=(a_ij)，x=(x_i)，b=(b_i)，有Ax=b。

如果A是非奇异矩阵（秩为n），存在逆矩阵A^-1，则x= A^-1b，且x是唯一解。方程数目少于未知量数目（或A的秩小于n）为欠定方程组，具有无穷多解，如果方程组不相容则可能无解；方程数目多于未知量数目为超定方程组，可能无解。

LUP分解在求解线性方程组具有数值稳定和速度快的优点。LUP分解的思想就是找出三个nXn矩阵L、U和P，满足PA=LU，其中L是一个单位下三角矩阵、U是一个上三角矩阵、P是一个置换矩阵。

定义y=Ux，其中x是要求解的未知向量。用正向替换（⊙(n²)）方法求解Ly=Pb，在用逆向替换（⊙(n²)）方法求解Ux=y，可得：

Ax=P^-1LUx==P^-1Ly==P^-1Pb=b

正向替换和逆向替换的方法通过导论中的说明很好理解，问题是如何有效找出LUP呢？只要找出LUP就可以通过正向替换和逆向替换求解x。

LU分解是采用高斯消元法。首先考虑A是nXn的非奇异矩阵，且P等价于I_n，找出A=LU，矩阵L和U称为A的LU分解。基于LU分解基础上开展LUP分解，数学基础一致，不失一般性而已，增加对置换矩阵P中值为1的元素的刻画。LU和LUP分解的算法要结合案例理解，这里不多描述。最重要还是要对矩阵运算和高斯消元法的数据基础有初步掌握。

基于LUP分解矩阵求解线性方程组的方法，同样适用于非奇异矩阵的求逆。矩阵乘法和计算逆矩阵问题具有相同难度的两个问题，在一定技术条件限制下可以使用一个算法在相同渐进时间内解决另外一个问题，即可用Strassen矩阵乘法算法来求一个矩阵的逆。算法导论中给出了证明以及求解过程描述。

3）最小二乘逼近

先理解对称正定矩阵的相关性质：任意对称正定矩阵都是非奇异矩阵；A是一个对称正定矩阵，则A的每一个主子式都是对称正定的，关于主子式的Schur补也是对称正定的；对一个对称正定矩阵进行LU分解不会出现除数为0的情形。

有了矩阵乘法和矩阵求逆的运算基础，加上对称正定矩阵的性质，来看看最小二乘逼近的应用。对给定一组数据的点进行曲线拟合是对称正定矩阵的重要应用。