线性代数及其应用（二）-矩阵代数

发表于 | 分类于 ，

1. 矩阵运算

1.1 矩阵乘法

若乘积AB有定义，AB的第i行第j列的元素是A的第i行与B的第j列对应元素乘积之和。若$(AB)_{ij}$表示AB的(i, j)元素，A为mxn矩阵，则$(AB)_{ij} = a_{i1}b_{1j} + a_{i2}b_{2j} + \dots + a_{in}b_{nj}$。

1.2 矩阵乘法的性质

设A为mxn矩阵，B、C的维数使下列各式的乘积有定义：

  • A(BC) = (AB)C
  • A(B + C) = AB + AC
  • (B + C)A = BA + CA
  • r(AB) = (rA)B = A(rB)
  • $I_{m}A = A = AI_{n}$

特别注意：

  • 一般情况下，$AB\neq BA$
  • 消去律对矩阵乘法不成立，即若$AB = AC$，一般情况下，B= C并不成立
  • 若乘积AB是零矩阵，一般情况下，不能断定A = 0或B = 0

1.3 矩阵的转置

设A与B表示矩阵，其维数使下列和与积有定义，则：

  • $(A^{T})^{T} = A$
  • $(A + B)^T = A^T + B^T$
  • 对任意数r，$(rA)^T = rA^T$
  • $(AB)^T = B^TA^T$

1.4 数值计算的注解

在计算机上求出AB的最快方法依赖于计算机存储矩阵的方法。标准的高性能算法（如LAPACK）中按列计算AB，正如我们所定义的那样。AB的定义使得我们可以在计算机上用并行算法计算，B的列可单独或分组分配给不同的处理器，因此可以同时计算AB的列。

代码：

  • Python版本
1
2
3
4
5
6
7
8
9
10
def Matrix_Mul(a,b):
    if a.shape[1] != b.shape[0]:
        print('这两个矩阵无法做乘法，请检查左边矩阵的列数是否与右边矩阵的行数相等！')
    else:
        c = np.zeros(a.shape[0]*b.shape[1]).reshape(a.shape[0],b.shape[1])
        for i in range(a.shape[0]):
            for j in range(b.shape[1]):
                for k in range(a.shape[1]):
                    c[i,j] = c[i,j] + a[i,k]*b[k,j]
    return c
  • C++版本
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
struct matrix{ ll a[M][M];
    ll* operator [](int x){ return a[x]; }
    matrix operator *(matrix& b){
        static matrix c;
        for(int i=1;i<=3;++i)
            for(int j=1;j<=3;++j)
                c[i][j]=0;
        for(int i=1;i<=3;++i)
            for(int j=1;j<=3;++j)
                for(int k=1;k<=3;++k)
                    c[i][j]+=a[i][k]*b[k][j],c[i][j]%=mod;
        return c;
    }
    inline void print(){
        for(int i=1;i<=n;++i,putchar('\n'))
            for(int j=1;j<=n;++j,putchar(' '))
                putchar('0'+a[i][j]);
    }
}F,T;

1.5 练习题

设A为4X4向量，x为$R^4$中的向量。计算$A^2x$的最快方法是什么？计算出乘法的次数。

答：首先计算Ax，再计算A(Ax)，总计需要32次乘法操作。因为Ax需要16次乘法运算，每个元素4次。A(Ax)也需要16次乘法操作。故总计需要32次乘法操作。

2. 矩阵的逆

类比于实数代数，即实数的倒数（乘法逆），矩阵的逆的定义如下：

  • 一个nXn矩阵A是可逆的，若存在一个nXn矩阵C使AC = I且CA = I，其中$I = I_n$是nXn的单位矩阵，这时称C是A的逆阵，且C是唯一的，一般将C记为$A^{-1}$，得到$AA^{-1} = I$或者$A^{-1}A = I$。
  • 不可逆矩阵有时称为奇异矩阵，而可逆矩阵也称为非奇异矩阵。

2.1 逆矩阵的求法

下面利用一个例子来说明逆矩阵的求法，并给出矩阵可逆的一些判断条件。

例：设$A = \left[ \begin{matrix} a & b \\ c & d \\ \end{matrix} \right]$，若$ad - bc \neq 0$，则A可逆且.

若$ad - bc = 0$，则A不可逆。其中，$det A = ad - bc$。

定理:

  • 若A是nXn可逆矩阵，则对每一$R^n$中的向量b，方程Ax = b有唯一解$x = A^{-1}b$
  • 若A是可逆矩阵，则$A^{-1}$也可逆且$(A^{-1})^{-1} = A$
  • 若A和B都是nXn可逆矩阵，AB也可逆，且其逆是A和B的逆矩阵按相反顺序的乘积，即$(AB)^{-1} = b^{-1}A^{-1}$
  • 若A可逆，则$A^{T}$也可逆，且其逆是$A^{-1}$的转置，即$(A^T)^{-1} = (A^{-1})^T$

代数余子式：代数余子式是从行列式的公式中提取出来的，它的作用是把n阶行列式化简为n – 1阶行列式。

2.2 初等矩阵

把单位矩阵进行一次行变换，就得到初等矩阵。

image-20191027003001123

定理：

  • nXn矩阵A是可逆的，当且仅当A行等价于$I_n$，这时，把A变为$I_n$的一系列初等行变换同时把$I_n$变成$A^{-1}$

2.3 $A^{-1}$的求法

把增广矩阵[A I]进行行化简。若A行等价于I，则[A I]行等价于$[I A^{-1}]$，否则A没有逆。

线性方程组解释：

用$e_1, \dots, e_n$表示$I_n$的各列，则把[A I]行变换成$[I A^{-1}]$的过程可看作是解n个方程组，，

其中这些方程组的“增广列”都放在A的右边，构成矩阵，。

方程$AA^{-1} = I$及矩阵乘法的定义说明$A^{-1}$的列正好是上述方程组的解。这一点是很有用的，因为在某些应用问题中，只需要$A^{-1}$的一列或两列，这时只需要解上述方程组的中的相应方程即可。

2.4 计算机中的矩阵求逆方法

(1) 数值计算的注解

在实际中，很少计算$A^{-1}$，除非需要$A^{-1}$的元素。计算$A^{-1}$和$A^{-1}b$总共需要的原酸次数大约是用行变换解方程组Ax = b的3倍，而且行变换可能更为精确。

(2) 高斯消元构造逆矩阵

  • 基本思想
    • 将原矩阵化成阶梯型矩阵
    • 朴素运算消成单位矩阵

  • 即上文2.2节内容

  • 代码：

    • C++版本
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
#include <iostream>
#include <stdlib.h>
using namespace std;
int n; static double a[50][50],b[50][50];
//交换当前行与下一行
void exchange(double a[][50],double b[][50],int current_line,int next_line,int all_line_number) {
    //交换两行
    int cl=current_line,nl=next_line,n=all_line_number;
    for(int i=0; i<n; i++)
        swap(a[cl][i],a[nl][i]),
        swap(b[cl][i],b[nl][i]); 
}
void the_data_change_to_1(double a[][50],double b[][50],int m,int n) { //将a[m][m]变成1
    if(a[m][m]) {
        for(int i=m+1; i<n; ++i) a[m][i]=a[m][i]/a[m][m];
        for(int i=0; i<n; ++i) b[m][i]=b[m][i]/a[m][m];
        return (void)(a[m][m]=1); //change_to_upper_angle_matrix(a,b,m,n);//将a[m][m]之下的元素全部变成0
    }
    while(m+1<n&&!a[m][m]) exchange(a,b,m,m+1,n);
    if(a[m][m]!=1) the_data_change_to_1(a,b,m,n);
}
//将a[m][m]之下的元素全部变成0
void change_to_upper_angle_matrix(double a[][50],double b[][50],int m,int n) {
    if(m+1>=n) return ;
    for(int i=m+1,t; i<n; ++i,a[i][m]=0) { t=a[i][m];
        for(int j=m; j<n; ++j) a[i][j]=a[i][j]-t*a[m][j];
        for(int k=0; k<n; ++k) b[i][k]=b[i][k]-t*b[m][k];
    }
}
/*将上三角矩阵变成单位阵*/
void change_to_unit_matrix(double a[][50],double b[][50],int l,int n) {
    if(l<=0) return ;
    for(int i=l-1; i>=0; a[i][l]=0,--i) //从a[l-1][l]开始向上，让每个元素都变成0
        for(int j=0; j<n; ++j)
            b[i][j]=b[i][j]-a[i][l]*b[l][j];
    --l,change_to_unit_matrix(a,b,l,n);
}
//打印结果
void print_result(double b[][50],int n) {
    for(int i=0; i<n; ++i,cout<<endl)
        for(int j=0; j<n; j++) cout<<b[i][j]<<" ";
}
int main() { cin>>n;
    for(int i=0; i<n; ++i)
        for(int j=0; j<n; ++j)
            scanf("%lf",a[i]+j);
    for(int i=0; i<n; ++i) b[i][i]=1;
    for(int i=0; i<n; ++i) {
        if(a[i][i]!=1) the_data_change_to_1(a,b,i,n);//将a[i][i]变成 1
        if(a[i][i]==1) change_to_upper_angle_matrix(a,b,i,n); //将a[i][i]之下的元素全部变成 0
    }
    change_to_unit_matrix(a,b,n-1,n);
    return print_result(b,n),0;
}

(3) 解方程组构造（非同阶）逆矩阵

  • 主要思想同2.3节
  • 最终构造的方程组也使用高斯消元法进行求解
  • 该解法时间复杂度为$O(n^4)$

3 可逆矩阵的特征

设A为nXn矩阵，则下列命题是等价的，即对某一特定的A，它们同时为真或同时为假：

  • A是可逆矩阵
  • A等价于nXn单位矩阵
  • A有n个主元位置
  • 方程Ax = 0仅有平凡解
  • A的各列线性无关
  • 线性变换$x\mapsto Ax$是一对一的
  • 对$R^n$中任意b，方程Ax = b至少有一个解
  • A的各列生成$R^n$
  • 线性变换$x\mapsto Ax$把$R^n$映射到$R^n$上
  • 存在nXn矩阵C使$CA = I$
  • 存在nXn矩阵D使$AD = I$
  • $A^T$是可逆矩阵

4. 分块矩阵

image-20191027155815026

4.1 分块矩阵的运算

  • 若矩阵A和B有相同维数且被同样的分块，则A+B的计算按照矩阵相加运算法则一样计算且被同样的分块
  • 分块矩阵与标量的乘法运算同矩阵与标量的乘法运算法则一致
  • 分块矩阵的乘法与矩阵的乘法运算法则一致

4.2 分块矩阵的逆

分块矩阵的求逆根据分块矩阵的维数和排列不同而计算结果不同，但计算过程一般如下例子所示。

例：设分块上三角矩阵$A = \left[ \begin{matrix} A_{11} & A_{22} \\ 0 & A_{22} \\ \end{matrix} \right]$，其中$A_{11}$是pXp矩阵，$A_{22}$是qXq矩阵，A为可逆矩阵，求$A^{-1}$。

解：设B为$A^{-1}$，那么得到，，

得到如下方程组：

  • $A_{11}B_{11} + A_{12}B_{21} I_{p}$ (1)
  • $A_{11}B_{12} + A_{12}B_{22} = 0$ (2)
  • $A_{22}B_{21} = 0$ (3)
  • $A_{22}B_{22} = I_{q}$ (4)

由(4)即$A_{22}$是方阵可知，$A_{22}$可逆。基于$A_{22}$可逆和公式 (4)即可求解上述方程组。解得，。

5. 矩阵因式分解

矩阵因式分解指的是将矩阵由多个矩阵进行表示。矩阵乘法是数据的综合，矩阵因式分解是数据的分解。矩阵的因式分解的主要意义在于便于计算。矩阵因式分解的方法有很多，例如LU分解、秩分解、QR分解、奇异值分解和谱分解等。

5.1 LU分解

给定一系列具有相同系数的方程组：$Ax = b_1, Ax = b_2, \dots, Ax = b_p$，如果A可逆，那么只需要求$A^{-1}b_1, \dots, A^{-1}b_{p}$等。接下来，我们首先对A进行LU分解。这样做的好处分析如下。

当A = LU时，方程Ax = b可写成LUx = b。另Ux = y，那么得到如下方程组：

  • Ly = b
  • Ux = y

这两个方程都比较容易解。

设A为mXn矩阵，另$A = LU$，其中L为mXm下三角矩阵，主对角元素全为1，U是A的等价的mXn的阶梯型矩阵。如下是矩阵LU分解的一个例子，其中L是可逆的，称为单位下三角矩阵。

image-20191027163408844

设A可以化为阶梯形U，化简过程仅用行倍加变换，即把一行的倍数加于它下面的另一行。这样，存在单位下三角初等矩阵$E_1, \dots, E_p$使。

之后，$A = (E_p\dots E_1)^{-1}U = LU$，其中$L = (E_p\dots E_1)^{-1}$。也就是说L就是使得A化为U的变换的逆。

  • 如果可能的话，用一系列的倍加变换把A化为阶梯形
  • 填充L的元素使相同的行变换把L变为I

例：

image-20191102225407207

image-20191102225430417

image-20191102225514477

计算复杂度分析，对于nXn稠密矩阵A，n相当大，例如$n\geq 30$：

  • 计算A的LU分解大约需要$2n^3/3$浮算，而求$A^{-1}$大约需要$2n^3$浮算
  • 解Ly = b和Ux = y大约需要$2n^2$浮算，因任意nXn三角方程组可以用大约$n^2$浮算解出
  • 把b乘以$A^{-1}$也需要$2n^2$浮算，但结果可能不如由L和U得出的精确（由于计算$A^{-1}$和$A^{-1}b$的舍入误差）
  • 若A是稀疏矩阵（大部分元素为0），则L和U可能也是稀疏的，然而$A^{-1}$很可能是稠密的。这时，用LU分解来解方程Ax = b很可能比用$A^{-1}$快很多。

5.2 QR分解

  • 是目前求一般矩阵全部特征值的最有效并广泛应用的方法
  • 将矩阵分解成一个正规正交矩阵Q与上三角形矩阵R
  • 分解过程中利用Gram−Schmidt正交化

其中，Gram−Schmidt正交化内容将在后续章节进行总结。

给定矩阵A，其列向量为$\alpha_1, \alpha_2, \dots, \alpha_n$，利用Gram−Schmidt正交化对矩阵A的列向量进行正交化，得到$\epsilon_1, \epsilon_2, \dots, \epsilon_j$，得到如下矩阵形式：,

其中，$Q = (\epsilon_1, \epsilon_2, \dots, \epsilon_j)$，$R = T^{-1}$。

求解步骤：

  • 通过Gram–Schmidt正交化求出正交矩阵Q
  • 再通过$R = Q^{T}A$得到矩阵R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
import numpy as np

def gram_schmidt(A):
    """Gram-schmidt正交化"""
    Q = np.zeros_like(A)
    cnt = 0
    for a in A.T:
        u = np.copy(a)
        for i in range(0, cnt):
            u -= np.dot(np.dot(Q[:, i].T, a), Q[:, i])
        e = u / np.linalg.norm(u)
        Q[:, cnt] = e
        cnt += 1
    R = np.dot(Q.T, A)
    return (Q, R)
  
 if __name__ == '__main__':
    np.set_printoptions(precision=4, suppress=True)
    A = np.array([[6, 5, 0], [5, -1, 4], [5, 1, -14], [0, 4, 3]], dtype=float)

    (Q, R) = gram_schmidt(A)
    print(Q)
    print(R)
    print(np.dot(Q, R))

5.3 SVD分解

在学习SVD分解前，先简单描述下如下定义，详细的内容将在后续进行总结。

  • 特征值
    • $Ax = \lambda x$
    • 矩阵是线性空间里的变换的描述。矩阵AA与向量相乘，本质上对向量x进行一次线性转换（旋转或拉伸），而该转换的效果为常数c乘以向量x（即只进行拉伸）。当我们求特征值与特征向量的时候，就是为了求矩阵A能使哪些向量（特征向量）只发生拉伸，而拉伸的程度，自然就是特征值λ了

对于一个矩阵A，有一组特征向量；再将这组向量进行正交化单位化，也就是我们学过的Schmidt正交化，就能得到一组正交单位向量。特征值分解，就是将矩阵A分解为如下方式：，其中，Q是矩阵A的特征向量组成的矩阵，Σ则是一个对角阵，对角线上的元素就是特征值。通过特征值分解得到的前N个特征向量，那么就对应了这个矩阵最主要的N个变化方向。利用这前N个变化方向，就可以近似这个矩阵（变换）。特征值分解可以得到特征值与特征向量，特征值表示的是这个特征到底有多重要，而特征向量表示这个特征是什么，可以将每一个特征向量理解为一个线性的子空间。

  • 特征值分解最大的问题是只能针对方阵
  • 为了解决该问题，SVD分解算法出现了

SVD分解的公式如下：。

SVD解决了特征值分解中只能针对方阵而没法对更一般矩阵进行分解的问题。

在大部分情况下，当我们把矩阵Σ里的奇异值按从大到小的顺序呢排列以后，很容易就会发现，奇异值减小的速度特别快。在很多时候，前10%甚至前1%的奇异值的和就占了全部奇异值和的99%以上。换句话说，大部分奇异值都很小，基本没什么用。因此，可以用前面r个奇异值来对这个矩阵做近似。于是，SVD也可以写为：，其中，$r\ll m$和$r\ll n$。实际在使用SVD的时候，需要我们根据不同的业务场景与需求还有资源情况，合理选择r的大小。本质是在计算精度与空间时间成本之间做个折中。

SVD算法流程如下：

  • 输入

    • 样本数据

  • 输出

    • 左奇异矩阵，奇异值矩阵，右奇异矩阵

  • 算法流程

    • 计算特征值

      • 特征值分解$AA^{T}$，其中$A\in R^{m\times n}$为原始样本数据，$AA^{T} = U\Sigma \Sigma^{T}U^{T}$
      • 得到左奇异矩阵$U\in R^{m\times m}$和奇异矩阵$\Sigma’ \in R^{m\times m}$

    • 间接求部分右奇异矩阵

      • 由$A = U\Sigma’V’$可得

    • 返回$U$，$\Sigma’$，$V’$

6. $R^{n}$子空间

  • $R^{n}$中的一个子空间是$R^n$中的集合H，具有以下三个性质
    • 零向量属于H
    • 对H中任意的向量u和v，u + v属于H
    • 对H中任意向量u和数c，cu属于H
  • 由$v_1, \dots, v_p$张成的子空间定义为$Span\{v_1, \dots, v_p\}$
  • 矩阵A的列空间是A的各列的线性组合的集合，记做Col A
  • 矩阵A的零空间是齐次方程Ax = 0的所有解的集合，记为Nul A
  • $R^n$中子空间H的一组基是H中一个线性无关集，它生成H
  • 矩阵A的主元列构成列空间的基

7. 维数和秩

  • 非零子空间H的维数，用dimH表示，是H的任意一个基的向量个数，零子空间的维数定义为零
  • 如果一个矩阵A有n列，则rank A +dim Nul A = n

8. 参考