GitModel-动手学线性代数（python）

1 动手学线性代数（GitModel学习笔记）

目录

• 1 动手学线性代数（GitModel学习笔记）
• • 1.1 理论知识
  • • 1.1.1 线性方程组与向量
    • 1.1.2 向量空间、矩阵、行列式及范数
    • • 向量的运算法则
      • 向量的线性相关与线性无关
      • 矩阵：矩阵不仅仅是用来解方程的
      • 特殊矩阵
      • 矩阵的逆
      • 非奇异阵的伪逆
    • 1.1.3 对角化、矩阵的特征值与特征向量、正交化
  • 1.2 实战：基于矩阵变换的图像变换

线性代数知识主要用于解决离散型矩阵形状的数据。python中用于线性代数矩阵处理的库主要是numpy。 numpy官网还介绍了numpy包中各种函数与MATLAB相同操作的一一对应关系。

1.1 理论知识 1.1.1 线性方程组与向量

• 引入：鸡兔同笼问题，今有头共10只，脚共28只。问鸡兔各有几何？

• 分析：利用高中知识建立方程组，设有鸡、兔各 x x x、 y y y只。由题意得：
  { x + y = 10 2 x + 4 y = 28 left { begin{aligned} x+y = 10 \ 2x+4y = 28 end{aligned} right. {x+y=102x+4y=28​

• 求解：由方程组得增广矩阵：
  [ 1 1 10 2 4 28 ] begin{bmatrix} 1 & 1 & 10 \ 2 & 4 & 28 end{bmatrix} [12​14​1028​]
  经初等变换：
  [ 1 1 10 2 4 28 ] → 第 二 行 − 第 一 行 ∗ 2 [ 1 1 10 0 2 8 ] → 第 二 行 / 2 [ 1 1 10 0 1 4 ] → 第 一 行 − 第 二 行 [ 1 0 6 0 1 4 ] begin{bmatrix} 1 & 1 & 10 \ 2 & 4 & 28 end{bmatrix} overset{第二行-第一行*2}{rightarrow} begin{bmatrix} 1 & 1 & 10 \ 0 & 2 & 8 end{bmatrix} overset{第二行/2}{rightarrow} begin{bmatrix} 1 & 1 & 10 \ 0 & 1 & 4 end{bmatrix} overset{第一行-第二行}{rightarrow} begin{bmatrix} 1 & 0 & 6 \ 0 & 1 & 4 end{bmatrix} [12​14​1028​]→第二行−第一行∗2[10​12​108​]→第二行/2[10​11​104​]→第一行−第二行[10​01​64​]

根据变换结果可知， x = 6 , y = 4 x=6, y =4 x=6,y=4。

python代码

import numpy as np

A = np.array([[1, 1], [2, 4]])
b = np.array([[10], [28]])
x = np.linalg.solve(A, b)
print("方程组的解为：n{}".format(x))

方程组的解为：
[[6.]
 [4.]]

1.1.2 向量空间、矩阵、行列式及范数

• 引入：对于刚才的题目，现又往里加入鸡鸭若干只，头变为14个，脚变为40只。问鸡、兔、鸭各几何？

• 分析：依然是建立方程组, 设有鸡、兔、鸭各 x x x、 y y y、 z z z只。由题意得：
  { x + y + z = 14 2 x + 4 y + 2 z = 40 left { begin{aligned} x+y+z=14 \ 2x+4y+2z=40 end{aligned} right. {x+y+z=142x+4y+2z=40​

• 求解：由方程得增广矩阵：
  [ 1 1 1 14 2 4 2 40 ] begin{bmatrix} 1 & 1 & 1 & 14\ 2 & 4 & 2 & 40 end{bmatrix} [12​14​12​1440​]

经初等变换：
[ 1 1 1 14 2 4 2 40 ] → 第 二 行 / 2 [ 1 1 1 14 1 2 1 20 ] → 第 二 行 − 第 一 行 [ 1 1 1 14 0 1 0 6 ] → 第 一 行 − 第 二 行 [ 1 0 1 8 0 1 0 6 ] begin{bmatrix} 1 & 1 & 1 & 14\ 2 & 4 & 2 & 40 end{bmatrix} overset{第二行/2}{rightarrow} begin{bmatrix} 1 & 1 & 1 & 14\ 1 & 2 & 1 & 20 end{bmatrix} overset{第二行-第一行}{rightarrow} begin{bmatrix} 1 & 1 & 1 & 14\ 0 & 1 & 0 & 6 end{bmatrix} overset{第一行-第二行}{rightarrow} begin{bmatrix} 1 & 0 & 1 & 8\ 0 & 1 & 0 & 6 end{bmatrix} [12​14​12​1440​]→第二行/2[11​12​11​1420​]→第二行−第一行[10​11​10​146​]→第一行−第二行[10​01​10​86​]

根据变换结果知： x + z = 8 , y = 6 x+z=8, y=6 x+z=8,y=6。这说明 x , z x, z x,z具有多个解。

python代码

import numpy as np

A = np.array([[1, 1, 1], [2, 4, 2]])
b = np.array([[14], [40]])
# x = np.linalg.solve(A, b)   # 不能用solve方法，因为solve方法要求 A为方阵（满秩方阵）
x = np.linalg.lstsq(A, b, rcond=None)     # lstsq方法返回的解是利用最小二乘法求解出来的，返回的解在在所有的解中范数最小
print("方程组的解为：n{}".format(x))

方程组的解为：
(array([[4.],
       [6.],
       [4.]]), array([], dtype=float64), 2, array([5.16724093, 0.54737667]))

向量的运算法则

给定 n n n维向量：
x = [ x 1 x 2 x 3 ⋯ x n ] y = [ y 1 y 2 y 3 ⋯ y n ] bf{x} = begin{bmatrix} x_{1} \ x_{2} \ x_{3} \ cdots \ x_{n} end{bmatrix} bf{y} = begin{bmatrix} y_{1} \ y_{2} \ y_{3} \ cdots \ y_{n} end{bmatrix} x=⎣⎢⎢⎢⎢⎡​x1​x2​x3​⋯xn​​⎦⎥⎥⎥⎥⎤​y=⎣⎢⎢⎢⎢⎡​y1​y2​y3​⋯yn​​⎦⎥⎥⎥⎥⎤​

• 向量相加:两个相同维度得向量相加。
  x + y = [ x 1 + y 1 x 2 + y 2 x 3 + y 3 ⋯ x n + y n ] bf{x+y} = begin{bmatrix} x_{1} + y_{1} \ x_{2} + y_{2}\ x_{3} + y_{3}\ cdots \ x_{n} + y_{n} end{bmatrix} x+y=⎣⎢⎢⎢⎢⎡​x1​+y1​x2​+y2​x3​+y3​⋯xn​+yn​​⎦⎥⎥⎥⎥⎤​

• 向量数乘：一个常数 k k k乘以一个向量。
  k x = [ k x 1 k x 2 k x 3 ⋯ k x n ] kbf{x}= begin{bmatrix} kx_{1} \ kx_{2} \ kx_{3} \ cdots \ kx_{n} end{bmatrix} kx=⎣⎢⎢⎢⎢⎡​kx1​kx2​kx3​⋯kxn​​⎦⎥⎥⎥⎥⎤​

注（习惯）：默认所有的向量都采用列向量的形式。

import numpy as np

x = np.array([1, 2, 3])
y = np.array([4, 5, 6])
print("x={}, y={}".format(x, y))

# 向量加法
print("x+y={}".format(x+y))

# 向量数乘
k = 3
print("k*x={}".format(k*x))

x=[1 2 3], y=[4 5 6]
x+y=[5 7 9]
k*x=[3 6 9]

向量的线性相关与线性无关

• 定义：给定一组向量 ( α 1 , α 2 , ⋯   , α k ) (alpha_1, alpha_2, cdots, alpha_k) (α1​,α2​,⋯,αk​)，对于向量 β beta β，如果能被存在一组不全为0的常数 m 1 , m 2 , ⋯   , m k m_1, m_2, cdots, m_k m1​,m2​,⋯,mk​，使得
  β = m 1 α 1 + m 2 α 2 + ⋯ + m k α k beta = m_1alpha_1 + m_2alpha_2 + cdots + m_kalpha_k β=m1​α1​+m2​α2​+⋯+mk​αk​
  则称向量 β beta β与向量组 ( α 1 , α 2 , ⋯   , α k ) (alpha_1, alpha_2, cdots, alpha_k) (α1​,α2​,⋯,αk​)是线性相关的，或称 β beta β可以被向量组 ( α 1 , α 2 , ⋯   , α k ) (alpha_1, alpha_2, cdots, alpha_k) (α1​,α2​,⋯,αk​)线性表出。否则，则成为非线性相关。一旦向量是线性相关的，也就说明 β beta β是一个多余的向量，因为它可以由其他的向量去表示。

对于之前的矩阵 [ 1 1 1 14 2 4 2 40 ] begin{bmatrix} 1 & 1 & 1 & 14\ 2 & 4 & 2 & 40 end{bmatrix} [12​14​12​1440​]，可以将未知数系数组成的矩阵 [ 1 1 1 2 4 2 ] begin{bmatrix} 1 & 1 & 1 \ 2 & 4 & 2 end{bmatrix} [12​14​12​]看成一组向量，由于第一列和第三列线性相关，故方程又多个解，因此，通过判断方程系数矩阵各列（或各行）是否线性相关可以推断出方程是否有唯一解或多解。

• 引出：然而，如果当一个方程组未知数的量很大之后，你需要去判断哪些方程是“有效的”也是一件非常花时间的工作，有没有一个更好的方法呢，答案是有的，那就是行列式。

注：只有 n n n维矩阵才有行列式。

import numpy as np

A = np.array([[1, 1, 1],
              [2, 4, 2],
              [2, 2, 2]])

np.linalg.det(A) # 计算方阵A的行列式
print("A的行列式的值为：",np.linalg.det(A))

B = np.array([[1,1,1,1],
              [1,2,0,0],
              [1,0,3,0],
              [1,0,0,4]])
B_det = np.linalg.det(B)
print("B的行列式的值为：",B_det)

A的行列式的值为： 0.0
B的行列式的值为： -2.0

判断方程是否有唯一解：

• n n n个未知数要有 n n n个方程
• 使用线性无关去判断“有效的方程”
• 系数行列式 ∣ A ∣ ≠ 0 |A| ne 0 ∣A∣​=0,则这个方程组是有唯一解的(克莱姆法则)

• 对于方程的解，克拉姆法则提出了方程的解的结构：
  设线性方程组的表达式为： { a 11 x 1 + a 12 x 2 + ⋯ + a 1 n x n = b 1 a 21 x 1 + a 22 x 2 + ⋯ + a 2 n x n = b 2 ⋯ ⋯ a n 1 x 1 + a n 2 x 2 + ⋯ + a n n x n = b n left{ begin{array}{c}a_{11} x_{1}+a_{12} x_{2}+cdots+a_{1 n} x_{n}=b_{1} \\ a_{21} x_{1}+a_{22} x_{2}+cdots+a_{2 n} x_{n}=b_{2} \\ cdots cdots \\ a_{n 1} x_{1}+a_{n 2} x_{2}+cdots+a_{n n} x_{n}=b_{n}end{array}right. ⎩⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎧​a11​x1​+a12​x2​+⋯+a1n​xn​=b1​a21​x1​+a22​x2​+⋯+a2n​xn​=b2​⋯⋯an1​x1​+an2​x2​+⋯+ann​xn​=bn​​，系数行列式为： D = ∣ a 11 a 12 ⋯ a 1 n a 21 a 22 ⋯ a 2 n ⋯ ⋯ ⋯ ⋯ a n 1 a n 2 ⋯ a m n ∣ ≠ 0 D = left|begin{array}{cccc}a_{11} & a_{12} & cdots & a_{1 n} \\ a_{21} & a_{22} & cdots & a_{2 n} \\ cdots & cdots & cdots & cdots \\ a_{n 1} & a_{n 2} & cdots & a_{m n}end{array}right| neq 0 D=∣∣∣∣∣∣∣∣∣∣∣∣∣∣​a11​a21​⋯an1​​a12​a22​⋯an2​​⋯⋯⋯⋯​a1n​a2n​⋯amn​​∣∣∣∣∣∣∣∣∣∣∣∣∣∣​​=0，则该线性方程组有且仅有唯一解:
  x 1 = D 1 D , x 2 = D 2 D , ⋯   , x n = D n D x_{1}=frac{D_{1}}{D}, x_{2}=frac{D_{2}}{D}, cdots, x_{n}=frac{D_{n}}{D} x1​=DD1​​,x2​=DD2​​,⋯,xn​=DDn​​
  其中， D j = ∣ a 11 ⋯ a 1 , j − 1 b 1 a 1 , j + 1 ⋯ a 1 n a 21 ⋯ a 2 , j − 1 b 2 a 2 , j + 1 ⋯ a 2 n ⋯ ⋯ ⋯ ⋯ ⋯ ⋯ ⋯ a n 1 ⋯ a n , j − 1 b n a n , j + 1 ⋯ a n n ∣ D_{j}=left|begin{array}{ccccccc}a_{11} & cdots & a_{1, j-1} & b_{1} & a_{1, j+1} & cdots & a_{1 n} \\ a_{21} & cdots & a_{2, j-1} & b_{2} & a_{2, j+1} & cdots & a_{2 n} \\ cdots & cdots & cdots & cdots & cdots & cdots & cdots \\ a_{n 1} & cdots & a_{n, j-1} & b_{n} & a_{n, j+1} & cdots & a_{n n}end{array}right| Dj​=∣∣∣∣∣∣∣∣∣∣∣∣∣∣​a11​a21​⋯an1​​⋯⋯⋯⋯​a1,j−1​a2,j−1​⋯an,j−1​​b1​b2​⋯bn​​a1,j+1​a2,j+1​⋯an,j+1​​⋯⋯⋯⋯​a1n​a2n​⋯ann​​∣∣∣∣∣∣∣∣∣∣∣∣∣∣​

例子：解线性方程组 { 2 x 1 + x 2 − 5 x 3 + x 4 = 8 x 1 − 3 x 2 − 6 x 4 = 9 2 x 2 − x 3 + 2 x 4 = − 5 x 1 + 4 x 2 − 7 x 3 + 6 x 4 = 0 left{begin{array}{l}2 x_{1}+x_{2}-5 x_{3}+x_{4}=8 \ x_{1}-3 x_{2}-6 x_{4}=9 \ 2 x_{2}-x_{3}+2 x_{4}=-5 \ x_{1}+4 x_{2}-7 x_{3}+6 x_{4}=0end{array}right. ⎩⎪⎪⎨⎪⎪⎧​2x1​+x2​−5x3​+x4​=8x1​−3x2​−6x4​=92x2​−x3​+2x4​=−5x1​+4x2​−7x3​+6x4​=0​

**解：**方程组的系数行列式
D = ∣ 2 1 − 5 1 1 − 3 0 − 6 0 2 − 1 2 1 4 − 7 6 ∣ = 27 ≠ 0 D=left|begin{array}{cccc} 2 & 1 & -5 & 1 \ 1 & -3 & 0 & -6 \ 0 & 2 & -1 & 2 \ 1 & 4 & -7 & 6 end{array}right|=27 neq 0 D=∣∣∣∣∣∣∣∣​2101​1−324​−50−1−7​1−626​∣∣∣∣∣∣∣∣​=27​=0
由克莱姆法则知：方程组有唯一解.

D 1 = ∣ 8 1 − 5 1 9 − 3 0 − 6 − 5 2 − 1 2 0 4 − 7 6 ∣ = 81 ⇒ x 1 = D 1 D = 81 27 = 3 D_{1}=left|begin{array}{cccc}8 & 1 & -5 & 1 \\ 9 & -3 & 0 & -6 \\ -5 & 2 & -1 & 2 \\ 0 & 4 & -7 & 6end{array}right|=81 Rightarrow x_{1}=frac{D_{1}}{D}=frac{81}{27} = 3 D1​=∣∣∣∣∣∣∣∣∣∣∣∣∣∣​89−50​1−324​−50−1−7​1−626​∣∣∣∣∣∣∣∣∣∣∣∣∣∣​=81⇒x1​=DD1​​=2781​=3

D 2 = ∣ 2 8 − 5 1 1 9 0 − 6 0 − 5 − 1 2 1 0 − 7 6 ∣ = − 108 ⇒ x 2 = D 2 D = − 108 27 = 4 D_{2}=left|begin{array}{cccc}2 & 8 & -5 & 1 \\ 1 & 9 & 0 & -6 \\ 0 & -5 & -1 & 2 \\ 1 & 0 & -7 & 6end{array}right|=-108Rightarrow x_{2}=frac{D_{2}}{D} =frac{-108}{27}= 4 D2​=∣∣∣∣∣∣∣∣∣∣∣∣∣∣​2101​89−50​−50−1−7​1−626​∣∣∣∣∣∣∣∣∣∣∣∣∣∣​=−108⇒x2​=DD2​​=27−108​=4

D 3 = ∣ 2 1 8 1 1 − 3 9 − 6 0 2 − 5 2 1 4 0 6 ∣ = − 27 ⇒ x 3 = D 3 D = = − 27 27 = − 1 D_{3}=left|begin{array}{cccc}2 & 1 & 8 & 1 \\ 1 & -3 & 9 & -6 \\ 0 & 2 & -5 & 2 \\ 1 & 4 & 0 & 6end{array}right|=-27 Rightarrow x_{3}=frac{D_{3}}{D} = =frac{-27}{27}=-1 D3​=∣∣∣∣∣∣∣∣∣∣∣∣∣∣​2101​1−324​89−50​1−626​∣∣∣∣∣∣∣∣∣∣∣∣∣∣​=−27⇒x3​=DD3​​==27−27​=−1

D 4 = ∣ 2 1 − 5 8 1 − 3 0 9 0 2 − 1 − 5 1 4 − 7 0 ∣ = 27 ⇒ x 4 = D 4 D = 27 27 = 1 D_{4}=left|begin{array}{cccc}2 & 1 & -5 & 8 \\ 1 & -3 & 0 & 9 \\ 0 & 2 & -1 & -5 \\ 1 & 4 & -7 & 0end{array}right|=27 Rightarrow x_{4}=frac{D_{4}}{D} = frac{27}{27} = 1 D4​=∣∣∣∣∣∣∣∣∣∣∣∣∣∣​2101​1−324​−50−1−7​89−50​∣∣∣∣∣∣∣∣∣∣∣∣∣∣​=27⇒x4​=DD4​​=2727​=1

python 代码

import numpy as np

D = np.array([[2.,1,-5,1],[1,-3,0,-6],[0,2,-1,2],[1,4,-7,6]])
D_det = np.linalg.det(D)

D1 = np.array([[8.,1,-5,1],[9,-3,0,-6],[-5,2,-1,2],[0,4,-7,6]])
D1_det = np.linalg.det(D1)

D2 = np.array([[2.,8,-5,1],[1,9,0,-6],[0,-5,-1,2],[1,0,-7,6]])
D2_det = np.linalg.det(D2)

D3 = np.array([[2.,1,8,1],[1,-3,9,-6],[0,2,-5,2],[1,4,0,6]])
D3_det = np.linalg.det(D3)

D4 = np.array([[2.,1,-5,8],[1,-3,0,9],[0,2,-1,-5],[1,4,-7,0]])
D4_det = np.linalg.det(D4)

x1 = D1_det / D_det
x2 = D2_det / D_det
x3 = D3_det / D_det
x4 = D4_det / D_det
print("克拉默法则解线性方程组的解为：n x1={:.2f},n x2={:.2f},n x3={:.2f},n x4={:.2f}".format(x1,x2,x3,x4))

克拉默法则解线性方程组的解为：
 x1=3.00,
 x2=-4.00,
 x3=-1.00,
 x4=1.00

矩阵：矩阵不仅仅是用来解方程的

• 矩阵乘法： A n × m B m × k = C n × k A_{ntimes m}B_{mtimes k}=C_{ntimes k} An×m​Bm×k​=Cn×k​

注：两矩阵相乘一定要注意两矩阵维度之间的对应关系。

import numpy as np

A = np.array([[1, 2, 3], [2, 3, 4]])
B = np.array([[2, 3, 4, 5], [4, 5, 6, 7], [2, 3, 4, 5]])
print("矩阵A的形状：{}".format(A.shape))
print(f"矩阵B的形状：{B.shape}")

print("AB=n {}".format(np.matmul(A, B)))
# BA 会报错，因为维度不一致

矩阵A的形状：(2, 3)
矩阵B的形状：(3, 4)
AB=
 [[16 22 28 34]
 [24 33 42 51]]

• 矩阵加法：与向量加类似，两个维度相同的矩阵才可以进行相加操作。

import numpy as np

A = np.array([[1, 2, 3], [2, 3, 4]])
B = np.array([[2, 3, 4], [4, 5, 6]])

print("A + B=n {}".format(np.add(A, B)))

A + B=
 [[ 3  5  7]
 [ 6  8 10]]

特殊矩阵

1. 单位矩阵：主对角线元素全为1， E = [ 1 0 ⋯ 0 0 1 ⋯ 0 ⋯ ⋯ ⋯ ⋯ 0 0 ⋯ 1 ] E=begin{bmatrix} 1 & 0 & cdots & 0 \ 0 & 1 & cdots & 0 \ cdots & cdots & cdots & cdots \ 0 & 0 & cdots & 1 end{bmatrix} E=⎣⎢⎢⎡​10⋯0​01⋯0​⋯⋯⋯⋯​00⋯1​⎦⎥⎥⎤​。

2. 初等矩阵：经初等变换后的矩阵称为初等矩阵。

初等变换：

• 对换变换：互换矩阵第 i i i行（列）和第 j j j行（列）的位置，记为 r i ↔ r j （ c i ↔ c j ） r_{i}leftrightarrow r_{j}（c_{i}leftrightarrow c_{j}） ri​↔rj​（ci​↔cj​）。
• 数乘变换：用一个非0常数 k k k乘矩阵的第 i i i行（列），记为 k r i ( k c i ) kr_{i}(kc_{i}) kri​(kci​)。
• 倍加变换：将矩阵第 j j j行（列）的元素的 k k k倍加到第 i i i行（列），记为 r i + k r j ( c i + k c j ) r_{i}+kr_{j}(c_{i}+kc_{j}) ri​+krj​(ci​+kcj​)。

• 初等矩阵性质
  • 1. 交换第 i i i个与第 j j j个方程↔️系数矩阵第 i i i行与第 j j j行互换↔️左乘一个( P n × n P_{n times n} Pn×n​)， P P P为单位阵 E n E_{n} En​第 i i i行与第 j j j行互换
  • 2. 第 i i i个方程左右乘非零常数 k k k倍↔️系数矩阵第 i i i行乘 k k k↔️左乘一个( P n × n P_{n times n} Pn×n​)， P P P为单位阵 E n E_{n} En​第 i i i行乘 k k k
  • 3. 第 i i i个方程加到第 j j j个方程中↔️系数矩阵第 i i i行加到第 j j j行↔️左乘一个( P n t i m e s n P_{n \times n} Pntimesn​)， P P P为单位阵 E n E_{n} En​第 i i i行加到第 j j j行

注：矩阵 A A A左乘一个初等矩阵，表示对矩阵 A A A做行变换；矩阵 A A A右乘一个初等矩阵，表示对矩阵 A A A做列变换

import numpy as np

print("单位矩阵：n{}".format(np.eye(3)))

print('初等矩阵的性质：n',"#"*50)

A = np.array([[1, 2, 3], [3, 4, 5], [3, 1, 2]])
print("矩阵A:n{}".format(A))

P = np.array([[0, 1, 0], [1, 0, 0], [0, 0, 1]])  
print("交换矩阵A第1行和第2行的位置：n{}".format(np.matmul(P, A)))

P = np.array([[1, 0, 0], [0, 2, 0], [0, 0, 1]])
print("矩阵A的第2行元素乘以2：n{}".format(np.matmul(P, A)))

P = np.array([[1, 0, 0], [0, 1, 0], [1, 0, 1]])
print("将矩阵A的第1行加到第3行：n{}".format(np.matmul(P, A)))

单位矩阵：
[[1. 0. 0.]
 [0. 1. 0.]
 [0. 0. 1.]]
初等矩阵的性质：
 ##################################################
矩阵A:
[[1 2 3]
 [3 4 5]
 [3 1 2]]
交换矩阵A第1行和第2行的位置：
[[3 4 5]
 [1 2 3]
 [3 1 2]]
矩阵A的第2行元素乘以2：
[[ 1  2  3]
 [ 6  8 10]
 [ 3  1  2]]
将矩阵A的第1行加到第3行：
[[1 2 3]
 [3 4 5]
 [4 3 5]]

矩阵的逆

• 定义：对于 n n n阶方阵 A A A，如果存在一个 n n n阶方阵 B B B，使得
  A B = B A = E n AB=BA=E_n AB=BA=En​
  则称 A A A是可逆的， B B B是 A A A的逆矩阵，记做 B = A − 1 B=A^{-1} B=A−1。

矩阵 A A A可逆的充分必要条件： ∣ A ∣ ≠ 0 |A|ne 0 ∣A∣​=0。

import numpy as np

A = np.array([[1, 2], [3, 4]])
print(np.linalg.det(A),"行列式不为0，非奇异阵") # 检验是否奇异
print("A的逆矩阵：n", np.linalg.inv(A)) # 矩阵求逆

A_inv = np.linalg.inv(A)

print("验证AA_inv = E n", np.matmul(A, A_inv))

-2.0000000000000004 行列式不为0，非奇异阵
A的逆矩阵：
 [[-2.   1. ]
 [ 1.5 -0.5]]
验证AA_inv = E 
 [[1.0000000e+00 0.0000000e+00]
 [8.8817842e-16 1.0000000e+00]]

非奇异阵的伪逆

• 定义：对于任意一个矩阵 A ∈ R n × m A in mathbb{R}^{ntimes m} A∈Rn×m，存在一个矩阵 A g ∈ R m × n A^g in mathbb{R}^{mtimes n} Ag∈Rm×n，使得 A A g A = A AA^gA=A AAgA=A，则称 A g A^g Ag为 A A A的伪逆(广义逆)。

import numpy as np

B = np.array([[1, 2], [2, 4]])
print("det(B)：{}".format(np.linalg.det(B)))
if  np.linalg.det(B) == 0:
    print("矩阵B的逆不存在。")

pB_inv = np.linalg.pinv(B)
print("B的伪逆为：n{}".format(pB_inv))
print("验证B*pB_inv*B?=B:n")
print("B*pB_inv*B=n{}.nB=n{}".format(np.matmul(B.dot(pB_inv), B), B))

det(B)：0.0
矩阵B的逆不存在。
B的伪逆为：
[[0.04 0.08]
 [0.08 0.16]]
验证B*pB_inv*B?=B:

B*pB_inv*B=
[[1. 2.]
 [2. 4.]].
B=
[[1 2]
 [2 4]]

1.1.3 对角化、矩阵的特征值与特征向量、正交化

• 引入：假设一个向量在坐标系 1 mathbb{1} 1下表示的坐标为 x x x，当这个向量 x x x经过一个线性变换形成一个新的向量 y y y，用矩阵表示这个变换就是： y = A x y = Ax y=Ax，矩阵 A A A对应着 x → y x rightarrow y x→y的线性变换. 同时，向量也可以在坐标系 2 mathbb{2} 2下表示，其坐标为 x ′ x^{'} x′，那么 x ′ = P x x^{'} = Px x′=Px. 同理， x ′ x^{'} x′也可以经过同一个线性变换变成 y ′ y^{'} y′，即： y ′ = B x ′ = B P x y^{'} = Bx^{'}=BPx y′=Bx′=BPx. 最后我们把 y ′ y^{'} y′转化为同一个坐标系下表达，即 y = P − 1 y ′ = P − 1 B P x y=P^{-1}y^{'}=P^{-1}BPx y=P−1y′=P−1BPx. 因此，我们可以得到： A x = P − 1 B P x Ax = P^{-1}BPx Ax=P−1BPx，即：
  A = P − 1 B P A = P^{-1}BP A=P−1BP
  称满足上式的矩阵A、B称为相似矩阵. 总结：一个向量在空间位置里，选取不同的坐标系，其坐标值是不同的. 对于空间中同一个线性变换，在不同的坐标系下，用于描述这个变换的矩阵也是不同的, 而这些不同矩阵所描述的线性变换是相似的，因此称他们为相似矩阵.

• 相似矩阵的作用：一个矩阵代表着一个线性变换，而不同的坐标系又会得到不同的相似矩阵，那能不能选用一个最佳的坐标系，使得我们描述的这个线性变换的矩阵是最佳的呢？什么矩阵才能称得上是最佳矩阵呢？答案就是对角矩阵！因为当我们同时需要经历很多次线性变换的时候，对角矩阵能极大的减少我们的计算量，即：
  A n = [ a 1 a 2 a 3 ] n = [ a 1 n a 2 n a 3 n ] A^{n}=left[begin{array}{lll} a_{1} & & \ & a_{2} & \ & & a_{3} end{array}right]^{n}=left[begin{array}{lll} a_{1}^{n} & & \ & a_{2}^{n} & \ & & a_{3}^{n} end{array}right] An=⎣⎡​a1​​a2​​a3​​⎦⎤​n=⎣⎡​a1n​​a2n​​a3n​​⎦⎤​
  那给定一组基，怎么将这组基转化成最优的基，即对角矩阵对应的基呢？答案是利用矩阵的特征值和特征向量。

• 特征值和特征向量：设 A A A是 n n n阶矩阵，如果存在数 λ lambda λ和非0的 n n n维向量 p p p，使得：
  A p = λ p Ap=lambda p Ap=λp
  则称数 λ lambda λ为矩阵 A A A的一个特征值，非0向量 p p p为 A A A的对应于（或属于）特征值 λ lambda λ的特征向量。

注：特征向量可以在某个矩阵的变换下保持在同一直线上，但没有发生角度的偏转。

import numpy as np

A = np.array([[1, 2, 3], [4, 4, 5], [4, 2, 1]])
print("矩阵A：n{}".format(A))

lamb, p = np.linalg.eig(A)
print("矩阵A的特征值为：n{}".format(lamb))
print("矩阵A的特征向量为：n{}".format(p))
print("矩阵A的对角化矩阵为：n{}".format(np.matmul(np.linalg.inv(p), np.matmul(A, p))))
# 将特别小的值过滤
res = np.matmul(np.linalg.inv(p), np.matmul(A, p))
res[np.abs(res) < 1e-6] = 0
print("矩阵A的对角化矩阵（过滤后）为：n{}".format(res))

矩阵A：
[[1 2 3]
 [4 4 5]
 [4 2 1]]
矩阵A的特征值为：
[ 8.4986663  -2.40063774 -0.09802855]
矩阵A的特征向量为：
[[ 0.38817063  0.5319676   0.2964301 ]
 [ 0.81766508  0.28923706 -0.84117071]
 [ 0.42514396 -0.7958344   0.45228423]]
矩阵A的对角化矩阵为：
[[ 8.49866630e+00 -3.51961651e-15  1.37857939e-15]
 [ 4.39211080e-16 -2.40063774e+00  1.93632180e-16]
 [-8.35826809e-16 -1.55490979e-15 -9.80285522e-02]]
矩阵A的对角化矩阵（过滤后）为：
[[ 8.4986663   0.          0.        ]
 [ 0.         -2.40063774  0.        ]
 [ 0.          0.         -0.09802855]]

• 正交矩阵：如果满足以下等式的矩阵 A A A称为正交矩阵：
  A T A = A A T = E A^TA = AA^T = E ATA=AAT=E
  通过上式可知：正交矩阵的逆矩阵与转置矩阵相同，即： A − 1 = A T A^{-1} = A^T A−1=AT，因为： A − 1 A = E A^{-1}A=E A−1A=E。

正交矩阵的作用：令一个向量进行旋转变换或者镜像变换。

正交矩阵规律：

• 正交矩阵的各行是单位向量且两两正交（垂直），或者说正交矩阵的各列是单位向量且两两正交（垂直）。
• α i ⋅ α i T = 1 alpha_{i} cdot alpha_{i}^{T}=1 αi​⋅αiT​=1，也就是说$alpha_{i} 是 单 位 向 量 （ 是单位向量（ 是单位向量（alpha_{i}$ 为正交矩阵的行向量）。
• α i ⋅ α j T = 0 ( i ≠ j ) alpha_{i} cdot alpha_{j}^{T}=0 quad(i neq j) αi​⋅αjT​=0(i​=j)，向量之间两两正交。
• 对于正交矩阵 A A A， ∣ A ∣ = ± 1 |A|=pm 1 ∣A∣=±1

给定一组基，怎么将其变成标准正交基？答案是利用施密特正交化。

原理：随便找一个原向量，把它固定住，接着找来第二个向量，往上面做投影，剩下的垂直分量就是第二个正交基，·····重复的做下去，直到所有基都变换过一遍后，得到的这组新的基就是正交基. 更进一步，如果我们把所有基进行标准化，使得每个基自己与自己做内积都为1，即模长为1，则得到的是一组标准正交基。

python代码

import numpy as np
import scipy

A = np.array([[2, 4, 5],
              [4, 5, 6],
              [5, 8 , 9]])

B = scipy.linalg.orth(A)    # 正交化
print("B*B.T=n{}".format(np.matmul(B, np.transpose(B))))

res = np.matmul(B, np.transpose(B))
res[np.abs(res) < 1e-6] = 0
print(res)

B*B.T=
[[ 1.00000000e+00 -5.39825374e-16  4.27489209e-16]
 [-5.39825374e-16  1.00000000e+00 -1.17069301e-16]
 [ 4.27489209e-16 -1.17069301e-16  1.00000000e+00]]
[[1. 0. 0.]
 [0. 1. 0.]
 [0. 0. 1.]]

1.2 实战：基于矩阵变换的图像变换

• 任务：将下面的图片旋转为便于正常人观看的的位置。

• 原理：设图像中某一点 ( x , y ) (x, y) (x,y)经旋转之后变为点 ( x ′ , y ′ ) (x^{'}, y^{'}) (x′,y′), 只需要找出原点 ( x , y ) (x, y) (x,y)和变换之后点 ( x ′ , y ′ ) (x^{'}, y^{'}) (x′,y′)之间的关系即可实现图像的旋转。

建立如上图所示的坐标轴，其中点 B = ( x , y ) B=(x, y) B=(x,y)为原点，点 C = ( x ′ , y ′ ) C=(x^{'}, y^{'}) C=(x′,y′)为点 B B B经旋转角度 α alpha α之后的点，点 B B B与横坐标轴的夹角为 θ theta θ，点 B B B到原点的距离为 r r r。

由上图可知,对于点 B B B:
{ x = r c o s θ y = r s i n θ left { begin{aligned} x = rcostheta \ y = rsintheta end{aligned} right. {x=rcosθy=rsinθ​

点 B B B经旋转角度 α alpha α到点 C C C之后，对于点 C C C：
{ x ′ = r c o s ( θ + α ) = r [ c o s θ c o s α − s i n θ s i n α ] = x c o s α − y s i n α y ′ = r s i n ( θ + α ) = r [ s i n θ c o s α + c o s θ s i n α ] = y c o s α + x s i n α left { begin{aligned} x^{'} = rcos(theta + alpha)=r[costheta cosalpha - sintheta sinalpha]=xcosalpha - ysinalpha \ y^{'} = rsin(theta + alpha)=r[sintheta cosalpha + costheta sinalpha]=ycosalpha + xsinalpha end{aligned} right. {x′=rcos(θ+α)=r[cosθcosα−sinθsinα]=xcosα−ysinαy′=rsin(θ+α)=r[sinθcosα+cosθsinα]=ycosα+xsinα​

三角函数公式

• c o s ( α + β ) = c o s α c o s β − s i n α s i n β cos(alpha + beta) = cosalpha cosbeta - sinalpha sinbeta cos(α+β)=cosαcosβ−sinαsinβ
• s i n ( α + β ) = s i n α c o s β + c o s α s i n β sin(alpha + beta) = sinalpha cosbeta + cosalpha sinbeta sin(α+β)=sinαcosβ+cosαsinβ

因此，可得变换矩阵：
[ x ′ , y ′ ] = [ x , y ] [ c o s α s i n α − s i n α c o s α ] [x^{'}, y^{'}] = [x, y]begin{bmatrix} cosalpha & sinalpha \ -sinalpha & cosalpha end{bmatrix} [x′,y′]=[x,y][cosα−sinα​sinαcosα​]

python实现

import numpy as np
from math import *

def point_rotate(x, y, alpha):
    Trans = np.array([[cos(alpha), sin(alpha)],
                      [-sin(alpha), cos(alpha)]])
    temp = np.array([x, y])
    res = np.matmul(temp, Trans)
    return res

res = point_rotate(3, 2, pi/2)
print(res)

[-2.  3.]

注意：但如果这样的话，会出现一个问题，对于一张图片而言，旋转中心在左上角，导致整张图片旋转不是中心旋转的. 下面我们需要对坐标轴进行平移，完善我们的变换公式。

假设图片宽度为 W W W，高度为 H H H，则在第一个坐标系下(左图)的坐标 ( x ′ , y ′ ) (x',y') (x′,y′)，变换之后的坐标为 ( x ′ ′ , y ′ ′ ) (x'',y'') (x′′,y′′)，则
x ′ ′ = x ′ − 1 2 W y ′ ′ = − y ′ + 1 2 H begin{aligned} x'' &= x'- frac{1}{2}W \ y'' &= -y'+ frac{1}{2}H end{aligned} x′′y′′​=x′−21​W=−y′+21​H​

则转换矩阵变为：
[ x ′ ′ , y ′ ′ , 1 ] = [ x ′ , y ′ ] [ 1 0 0 0 − 1 0 − 1 2 W 1 2 H 1 ] [x^{''}, y^{''}, 1] = [x^{'}, y^{'}]begin{bmatrix} 1 & 0 & 0 \ 0 & -1 & 0\ -frac{1}{2} W & frac{1}{2} H & 1\ end{bmatrix} [x′′,y′′,1]=[x′,y′]⎣⎡​10−21​W​0−121​H​001​⎦⎤​

记 B = [ 1 0 0 0 − 1 0 − 1 2 W 1 2 H 1 ] B = begin{bmatrix} 1 & 0 & 0 \ 0 & -1 & 0\ -frac{1}{2} W & frac{1}{2} H & 1\ end{bmatrix} B=⎣⎡​10−21​W​0−121​H​001​⎦⎤​

则：
[ x ′ , y ′ ] = [ x ′ ′ , y ′ ′ , 1 ] B − 1 [x^{'}, y^{'}] = [x^{''}, y^{''}, 1]B^{-1} [x′,y′]=[x′′,y′′,1]B−1

对于之前的变换 [ x ′ , y ′ ] = [ x , y ] [ c o s α s i n α − s i n α c o s α ] [x^{'}, y^{'}] = [x, y]begin{bmatrix} cosalpha & sinalpha \ -sinalpha & cosalpha end{bmatrix} [x′,y′]=[x,y][cosα−sinα​sinαcosα​]需要做进一步的修改，以使其满足矩阵运算的维度要求。具体变换如下：
[ x ′ , y ′ , 1 ] = [ x , y , 1 ] [ c o s α s i n α 0 − s i n α c o s α 0 0 0 1 ] [x^{'}, y^{'}, 1] = [x, y, 1]begin{bmatrix} cosalpha & sinalpha & 0 \ -sinalpha & cosalpha & 0 \ 0 & 0 & 1 end{bmatrix} [x′,y′,1]=[x,y,1]⎣⎡​cosα−sinα0​sinαcosα0​001​⎦⎤​

记 C = [ c o s α s i n α 0 − s i n α c o s α 0 0 0 1 ] C = begin{bmatrix} cosalpha & sinalpha & 0 \ -sinalpha & cosalpha & 0 \ 0 & 0 & 1 end{bmatrix} C=⎣⎡​cosα−sinα0​sinαcosα0​001​⎦⎤​

最后将图像上的一点旋转到一定位置需要经过以下步骤：

• 将图像坐标 ( x ′ , y ′ ) (x^{'}, y^{'}) (x′,y′)转换为标准坐标: [ x ′ , y ′ , 1 ] B [x^{'}, y^{'}, 1]B [x′,y′,1]B
• 将转换的标准坐标进行旋转： [ x ′ , y ′ , 1 ] B C [x^{'}, y^{'}, 1]BC [x′,y′,1]BC
• 将旋转后的标准坐标转换回图像坐标： [ x ′ , y ′ , 1 ] B C B − 1 [x^{'}, y^{'}, 1]BCB^{-1} [x′,y′,1]BCB−1

得到转换图像坐标的矩阵变换：
[ x ′ ′ , y ′ ′ , 1 ] = [ x ′ , y ′ , 1 ] B C B − 1 = [ x ′ , y ′ , 1 ] [ 1 0 0 0 − 1 0 − 1 2 W 1 2 H 1 ] [ c o s α s i n α 0 − s i n α c o s α 0 0 0 1 ] [ 1 0 0 0 − 1 0 − 1 2 W 1 2 H 1 ] = [ x ′ , y ′ , 1 ] [ c o s α − s i n α 0 s i n α c o s α 0 − 0.5 W c o s α − 0.5 H s i n α − 0.5 W 0.5 W s i n α − 0.5 H c o s α + 0.5 H 1 ] begin{aligned} [x^{''}, y^{''}, 1] &=[x^{'}, y^{'}, 1]BCB^{-1} \ &= [x^{'}, y^{'}, 1] begin{bmatrix} 1 & 0 & 0 \ 0 & -1 & 0 \ -frac{1}{2} W & frac{1}{2} H & 1\ end{bmatrix} begin{bmatrix} cosalpha & sinalpha & 0 \ -sinalpha & cosalpha & 0 \ 0 & 0 & 1 end{bmatrix} begin{bmatrix} 1 & 0 & 0 \ 0 & -1 & 0 \ -frac{1}{2} W & frac{1}{2} H & 1\ end{bmatrix} \ &= [x^{'}, y^{'}, 1] begin{bmatrix} cosalpha & -sinalpha & 0 \ sinalpha & cosalpha & 0 \ -0.5Wcosalpha-0.5Hsinalpha-0.5W & 0.5Wsinalpha-0.5Hcosalpha+0.5H & 1\ end{bmatrix} end{aligned} [x′′,y′′,1]​=[x′,y′,1]BCB−1=[x′,y′,1]⎣⎡​10−21​W​0−121​H​001​⎦⎤​⎣⎡​cosα−sinα0​sinαcosα0​001​⎦⎤​⎣⎡​10−21​W​0−121​H​001​⎦⎤​=[x′,y′,1]⎣⎡​cosαsinα−0.5Wcosα−0.5Hsinα−0.5W​−sinαcosα0.5Wsinα−0.5Hcosα+0.5H​001​⎦⎤​​

python代码实现

import numpy as np
from math import *
from PIL import Image
img = Image.open('figures/chiken.jpg')
print(img.size)
img

(250, 376)

def point_rotate(x0, y0, W, H, alpha):
    trans = np.array([[cos(alpha), -sin(alpha), 0],
                       [sin(alpha), cos(alpha), 0],
                       [-0.5*W*cos(alpha)-0.5*H*sin(alpha)-0.5*W, 0.5*W*sin(alpha)-0.5*H*cos(alpha)+0.5*H, 1]])
    temp = np.array([x0, y0, 1])
    res = np.matmul(temp, trans)
    return (int(res[0]), int(res[1]))


def image_rotate(image, alpha):
    (H, W, C) = image.shape
    img = np.zeros((max(W, H), max(W, H), C))
    for k in range(C):
        for i in range(W):
            for j in range(H):
                (x1, y1) = point_rotate(i, j, W, H, alpha)
                img[y1, x1, k] = image[j, i, k]
    return img

a = np.array(img)
res = image_rotate(a, pi/2)
b = Image.fromarray(np.uint8(res))
b