PYTHON替代MATLAB在线性代数学习中的应用
使用Python辅助MIT 18.06 Linear Algebra学习
前言
MATLAB一向是理工科学生的必备神器，但随着中美贸易冲突的一再升级，禁售与禁用的阴云也持续笼罩在高等学院的头顶。也许我们都应当考虑更多的途径，来辅助我们的学习和研究工作。
虽然PYTHON和众多模块也属于美国技术的范围，但开源软件的自由度毕竟不是商业软件可比拟的。
本文是一篇入门性文章，以麻省理工学院（MIT) 18.06版本线性代数课程为例，按照学习顺序介绍PYTHON在代数运算中的基本应用。
介绍PYTHON代数计算的文章非常多，但通常都是按照模块作为划分顺序，在实际应用中仍然有较多困扰。
而按照代数课程学习的顺序，循序渐进，集注在最常用和最实用的功能上，比较符合典型的应用逻辑。可以用较低的门槛逐步完成PYTHON在线性代数中各项功能的学习和应用。
MIT 2020版本的线性代数课程也已发布，但基本是在18.06版本上的修正。Gilbert教授的年龄已经很大，只录制了一个5节课的串讲。所以系统性还是18.06版本更为完整。
很讽刺是吧，课程本身也是美国的-_-#。阿Q一下吧，就当是“师夷长技以制夷”。
首先给出几个相关链接：
MIT 18.06 Linear Algebra课程主页
B站完整版34讲Gilbert教授课程视频
配套第三版线性代数教材(百度网盘) 提取码：uhvc
最新发行的教材是第5版，建议听课时使用配套的第3版教材。课程完成后，把第5版教材作为辅助读物。不然在章节、内容方面会碰到很多困惑。
版本选择
PYTHON版本的选择现在已经没有什么困惑了，PYTHON2停止了支持，PYTHON3现在是必选项。我是用Mac电脑，通常使用brew安装PYTHON3，每次有新版本的时候执行brew upgrade会自动升级。不使用内置的PYTHON3是为了防止安装很多扩展库的时候会有系统完整性检查导致的不兼容，不过只是跑一跑数学运算的话倒也不用担心这些。
Linux各发行版各发行版是Python的开发环境，所以内置的PYTHON3就很好用。使用apt/yum等包管理工具升级的时候会自动完成版本维护。
PYTHON在Windows/Linux/Mac等各平台上兼容性非常好，特别是在数学计算方面基本不用担心互相之间的通用问题。
计算模块方面则有了很多的选择，常见的有NumPy/SciPy/SymPy。
其中在数值计算方面，NumPy应用非常广泛，特别是TensorFlow/PyTorch等机器学习平台也把NumPy当做默认支持之后。所以本文会把NumPy当做一个选择。
在课程学习和理论研究方面，符号计算更为重要。SymPy在这方面有比较多的优势，所以本文会把SymPy当做另外一个选择。
SciPy以及还有一些小众计算模块同样非常优秀，但限于篇幅，本文中只好做一些取舍。
在PYTHON3环境下安装NumPy/SymPy模块的方法很简单：
pip3 install numpy sympy
如果碰到麻烦，一般是因为网络速度造成的。特别是默认的国外软件源。修改软件源到国内的服务器会提高不少下载速度，方法是修改文件~/.pip/pip.conf，默认是没有这个文件的，要自己建立~/.pip目录和新建对应的文本文件，内容为：
[global]
timeout = 6000
index-url = https://pypi.tuna.tsinghua.edu.cn/simple
trusted-host = pypi.tuna.tsinghua.edu.cn
这里使用了清华大学的镜像服务器。
以上是在Linux/Mac之上的操作方法。Windows用户，虽然PYTHON3本身没有兼容问题，但还是建议你使用Windows10内置的Linux子系统来学习。因为Mac/Linux下Python的资料更为丰富，能让你节省很多时间。
矩阵的表达
在Pyhton中使用扩展库，首先要做引用，比如引入NumPy库：
import numpy as np
意思是引用numpy计算库，并重新命名为np。使用numpy中的方法时，首先要以“np.”开头。
SymPy库的引用，通常会直接从中将所有资源直接引用到当前作用域，像使用原生方法一样使用SymPy中定义的方法，这也是SymPy官方推荐的：
from sympy import *
出于个人习惯，我还是更喜欢同使用NumPy一样使用SymPy:
import sympy as sp
虽然因此所有的SymPy的方法都要冠以“sp.”前缀，但这样不容易造成混淆从而导致错误。
以下内容大致对应课程（MIT 18.06 Linear Algebra课程，以下简称课程）第一、二讲的内容。
在线性代数中，主要涉及3种数据类型，常量、标量（Scalar)、向量(Vector)、矩阵（Matrix)。
无论NumPy还是SymPy，都直接使用了基本Python类型作为标量，比如：
c1 = 5
而对于向量和矩阵，处理方法则有很大区别，下面先讲述NumPy中的方法。
假设我们有v1、v2两个向量，及A、B两个矩阵：
\[v1 = \left[\begin{matrix}1\\2\\\end{matrix}\right]\] \[v2 = \left[\begin{matrix}3\\4\\\end{matrix}\right]\] \[A = \left[\begin{matrix}1&2\\3&4\\\end{matrix}\right]\] \[B = \left[\begin{matrix}5&6\\7&8\\\end{matrix}\right]\]- 首先，NumPy接受Python原生的数组当做向量和矩阵
# 除非特别注明，我们的示例都在交互方式使用Python # 每一行开始的“>>>”就是交互方式下Python给出的提示符 >>> v1c = [1,2] >>> v2c = [3,4] >>> Ac = [[1,2],[3,4]] >>> Bc = [[5,6],[7,8]] >>> v1c #交互方式直接给出变量名是显示变量内容的意思 [1, 2] >>> v2c [3, 4] >>> Ac [[1, 2], [3, 4]] >>> Bc [[5, 6], [7, 8]] - 其次，NumPy内置的数组类型（Array)也可以用来表示向量和矩阵
>>> import numpy as np #别忘记引用numpy,以后不再特别注明 >>> v1n=np.array([1,2]) >>> v2n=np.array(v2c) #直接使用前面定义好的内部数组类型初始化向量是一样的 >>> An=np.array([[1,2],[3,4]]) >>> Bn=np.array(Bc) #直接使用前面定义好的内部数组类型初始化矩阵 >>> v1n array([1, 2]) >>> v2n array([3, 4]) >>> An array([[1, 2], [3, 4]]) >>> Bn array([[5, 6], [7, 8]]) - 正式的矩阵表示方法，使用NumPy内置的Matrix类型
>>> v1 = np.mat([[1],[2]]) >>> v2 = np.mat([[3],[4]]) >>> A = np.mat([[1,2],[3,4]]) >>> B = np.mat(Bc) #使用以前定义过的内部数组类型来定义矩阵 >>> v1 matrix([[1], [2]]) >>> v2 matrix([[3], [4]]) >>> A matrix([[1, 2], [3, 4]]) >>> B matrix([[5, 6], [7, 8]])第1、2种方法，虽然在很多情况下都能正常的使用，但都算不上规范化的矩阵表示方法。特别是对于向量的表示，向量本来是纵向的，代表矩阵中的一列。但在方法1、2中，都成为了横向的。这很容易造成概念的混淆，和计算中的错误。
当然Python内置的列表类型以及NumPy内置的列表类型并非不能使用，实际上它们在计算速度上有非常明显的优势。简单说就是功能少的类型，往往有高的速度。所以渐渐的我们能看到很多需要速度的运算，还是会使用np.array类型操作。实际上对各种类型熟悉了之后，有了自己的习惯和原则，什么时候用什么类型，你自然会有自己的标准。
为了让大家对这种差异有更清晰的认识，这里举几个例子，也顺便看一看最基本的矩阵计算：# 计算 矩阵*常量 >>> Ac*3 #Python内部列表类型得到完全错误的结果，不可用 [[1, 2], [3, 4], [1, 2], [3, 4], [1, 2], [3, 4]] >>> An*3 array([[ 3, 6], [ 9, 12]]) >>> A*3 matrix([[ 3, 6], [ 9, 12]]) # 计算 向量*常量 >>> v1c*3 #Python内部列表类型得到完全错误的结果，不可用 [1, 2, 1, 2, 1, 2] >>> v1n*3 array([3, 6]) >>> v1*3 matrix([[3], [6]]) # 计算向量加法 >>> v1c+v2c #Python内部列表类型得到完全错误的结果，不可用 [1, 2, 3, 4] >>> v1n+v2n array([4, 6]) >>> v1+v2 matrix([[4], [6]]) # 计算矩阵乘法 >>> Ac*Bc #Python内部没有定义对应操作，不可用 Traceback (most recent call last): File "<stdin>", line 1, in <module> TypeError: can‘t multiply sequence by non-int of type 'list' >>> An*Bn #NumPy列表类型只是对应行列数字相乘，不是真正的矩阵乘法 array([[ 5, 12], [21, 32]]) >>> A*B matrix([[19, 22], [43, 50]])以上的示例可以明显看出，对于Python内置数组类型，并没有定义对应的矩阵操作，所以不能直接用于线性代数的计算。NumPy的很多方法都接受使用Python内部数组作为参数来表达向量和矩阵，所以给人的印象，这些类型之间没有什么区别。
NumPy内置的数组类型和矩阵类型，在简单运算中都能得到正确的结果，可以用于常用的计算。但实际上很多高级函数及算法，对两种类型的处理仍然存在很大区别，就类似示例中出现的矩阵乘法。所以在彻底了解之前，不建议使用np.array类型当做矩阵类型来使用。否则在复杂的项目中，很多莫名其妙的计算错误会让你排错工作异常复杂。
NumPy还提供了一种更方便的方法来定义向量和矩阵，是我当前最常用的方法：
>>> v1=np.mat("1;2")
>>> v2=np.mat("3;4")
>>> A=np.mat("1 2;3 4")
>>> B=np.mat("5 6;7 8")
>>> v1
matrix([[1],
[2]])
>>> v2
matrix([[3],
[4]])
>>> A
matrix([[1, 2],
[3, 4]])
>>> B
matrix([[5, 6],
[7, 8]])
>>> As=sp.Matrix(A) #sympy的Matrix定义方法太麻烦了，有的时候你会喜欢用np.mat转换过来
>>> As
Matrix([
[1, 2],
[3, 4]])
熟悉MATLAB的同学应当开心了，这跟MATLAB中定义矩阵的方法完全一样，算是Python环境中最方便的方法。
在线性代数课程中，经常会需要选取一个典型矩阵，做一些计算的练习。课堂上Gilbert教授经常随手就可以举出一个矩阵的例子，并且各行列线性无关，而我们往往很难做到这一点。这时候可以使用随机数初始化矩阵的方法：
>>> C=np.random.rand(3,3) #使用随机数初始化一个3x3矩阵
>>> C
array([[0.58896997, 0.45879526, 0.34384609],
[0.78480365, 0.19043321, 0.69538183],
[0.66016878, 0.81037627, 0.75616191]])
#也可以选择整数的版本，第一个参数100的意思是产生1-100的随机数
>>> np.random.randint(100,size=(2,3))
array([[51, 8, 75],
[64, 67, 20]])
与此类似的，还有初始化一个全0矩阵、全1矩阵、单位方阵I的方法，如果你打算用程序逻辑建立一个矩阵，那这些往往是开始的第一步：
#定义并初始化一个全0矩阵，参数为维度形状(m,n)，所以是两个括号嵌套
>>> np.zeros((4,3))
array([[0., 0., 0.],
[0., 0., 0.],
[0., 0., 0.],
[0., 0., 0.]])
#定义并初始化一个全1矩阵，第一个参数为维度形状(m,n)
>>> np.ones((3,3),dtype=float)
array([[1., 1., 1.],
[1., 1., 1.],
[1., 1., 1.]])
#定义并初始化一个单位矩阵I
>>> np.eye(3,3)
array([[1., 0., 0.],
[0., 1., 0.],
[0., 0., 1.]])
下面看看SymPy定义向量、矩阵的方法。
>>> import sympy as sp #别忘记引入函数库，以后将不再提醒
#喜欢使用from sympy import *方式的自行修改对应代码
>>> v1s=sp.Matrix([[1],[2]])
>>> v2s=sp.Matrix([[3],[4]])
>>> As=sp.Matrix([[1,2],[3,4]])
>>> Bs=sp.Matrix(Bc) #使用Python内置数组当做参数初始化矩阵
>>> v1s
Matrix([
[1],
[2]])
>>> v2s
Matrix([
[3],
[4]])
>>> As
Matrix([
[1, 2],
[3, 4]])
>>> Bs
Matrix([
[5, 6],
[7, 8]])
# 基本运算示例
>>> v1s+v2s
Matrix([
[4],
[6]])
>>> As*Bs
Matrix([
[19, 22],
[43, 50]])
>>> As*3
Matrix([
[3, 6],
[9, 12]])
# sympy定义并初始化一个随机矩阵
>>> sp.randMatrix(2,3)
Matrix([
[44, 34, 34],
[33, 15, 96]])
# 定义并初始化一个全0矩阵
>>> sp.zeros(2,3)
Matrix([
[0, 0, 0],
[0, 0, 0]])
# 定义并初始化一个全1矩阵
>>> sp.ones(2,3)
Matrix([
[1, 1, 1],
[1, 1, 1]])
# 定义并初始化一个单位矩阵I
>>> sp.eye(3,3)
Matrix([
[1, 0, 0],
[0, 1, 0],
[0, 0, 1]])
作为符号计算的代表，SymPy的计算结果通常都是公式形式，所以SymPy专门提供了LaTeX的输出方式：
>>> sp.latex(As)
'\\left[\\begin{matrix}1 & 2\\\\3 & 4\\end{matrix}\\right]'
这种输出格式对通常的程序没有什么意义。但如果是用于论文写作的话，可以直接拷贝到LaTex编辑器，成为一个精致的公式输出。就类似本文中的公式，通常也是采用LaTeX格式输入。
求解线性方程
这也是课程第一、二讲中的内容。方程组是矩阵的起源，也是矩阵最初的目的。以课程中的方程组为例：
\[\left\{ \begin{array}{l} 2x-y=0\\ -x+2y=3 \end{array} \right.\]可以得到矩阵A及向量b:
\[A = \left[\begin{matrix}2&-1\\-1&2\\\end{matrix}\right]\] \[b = \left[\begin{matrix}0\\3\\\end{matrix}\right]\] \[Ax=b\]这里的x实际代表两个未知数组成的向量：
\[x = \left[\begin{matrix}x\\y\\\end{matrix}\right]\]使用NumPy解方程组示例：
>>> A=np.mat("2 -1;-1 2")
>>> b=np.mat("0;3")
>>> np.linalg.solve(A,b)
matrix([[1.], #未知数x
[2.]]) #未知数y
使用SymPy解方程组示例：
>>> A=sp.Matrix([[2,-1],[-1,2]])
>>> b=sp.Matrix([[0],[3]])
>>> A.LDLsolve(b)
Matrix([
[1],
[2]])
作为符号计算的优势，SymPy中可以定义未知数符号之后，再使用跟NumPy中同名的方法solve()来直接对一个方程组求解，但那个不属于本文的主题范畴，所以不做介绍。有兴趣的话也可以参考这篇老博文《从零开始学习PYTHON3讲义（十一）》。
SymPy跟NumPy语法差异还是比较大的，使用中需要特别注意。两个软件包，虽然都是Python中的实现，但并不是由同一支开发团队完成的。所以这种差异感始终是存在的。
矩阵乘法和逆矩阵
这是课程第三讲的内容，其中矩阵同矩阵的乘法运算在前面一开始就演示过了，对于手工计算来讲，这是最繁琐的部分。而对于Python,这是最简单的部分。
矩阵的逆在线性代数中会频繁出现，经常用到，两个软件包中都已经有了内置的方法。
下面在同一个代码块中分别演示NumPy和SymPy的方法：
#numpy矩阵定义及矩阵乘法
>>> A=np.mat("1 2;3 4")
>>> B=np.mat("5 6;7 8")
>>> A*B
matrix([[19, 22],
[43, 50]])
#sympy矩阵定义及矩阵乘法
>>> As=sp.Matrix([[1,2],[3,4]])
>>> Bs=sp.Matrix([[5,6],[7,8]])
>>> As*Bs
Matrix([
[19, 22],
[43, 50]])
>>> A*Bs #numpy的矩阵*sympy矩阵，两个软件包的变量是可以相互通用的，但通常尽量不这样做
Matrix([
[19, 22],
[43, 50]])
# numpy求逆
>>> np.linalg.inv(A)
matrix([[-2. , 1. ],
[ 1.5, -0.5]]) #数值计算会尽量求得精确小数
>>> A ** -1 #使用求幂的方式获得逆矩阵，**是Python内置的求幂运算符，numpy做了重载
matrix([[-2. , 1. ],
[ 1.5, -0.5]])
>>> np.linalg.matrix_power(A,-1) #numpy中的求幂函数
matrix([[-2. , 1. ],
[ 1.5, -0.5]])
>>> A.I #矩阵类型求逆的快捷方式
matrix([[-2. , 1. ],
[ 1.5, -0.5]])
# sympy求逆
>>> As.inv()
Matrix([
[ -2, 1],
[3/2, -1/2]]) #符号计算会保持分数形式
#numpy也可以从sympy的计算结果中，获取计算数值，通常，这能提供更高的精度
#当然，sympy并不以速度见长
#后面的参数是将结果转换为浮点数，否则sympy数据会当做对象存储在numpy矩阵
>>> np.mat(As.inv(),dtype=float)
matrix([[-2. , 1. ],
[ 1.5, -0.5]])
#sympy中使用求幂的方式获得逆矩阵
>>> As ** -1 #sympy所重载的求幂运算符
Matrix([
[ -2, 1],
[3/2, -1/2]])
>>> As.pow(-1) #sympy标准的求幂函数
Matrix([
[ -2, 1],
[3/2, -1/2]])
# 分别证明A的逆*A=单位矩阵I
>>> np.linalg.inv(A)*A
matrix([[1.00000000e+00, 0.00000000e+00],
[1.11022302e-16, 1.00000000e+00]]) #注意左边的值e-16,是一个很接近0的小数
>>> As*As.inv()
Matrix([
[1, 0],
[0, 1]]) #符号计算通常能精确的还原应有的整数
上面代码非常明显的体现出了NumPy数值计算和SymPy符号计算的不同。前者会因精度所限有极小的误差，而后者通常能保持美观的整数数字。但前者的数字可以直接应用到机器学习等业务系统。而后者是对人的理解更有益，归根结底还是符号，不能当做数值使用。
好在Python之中，如果不考虑转换速度，不同模块之间共享数据非常容易。前面的演示中已经有了将NumPy矩阵转换为SymPy矩阵，以及将SymPy的计算结果转换到NumPy的实例。这对用户来说，是非常方便的。
矩阵的LU分解
课程第四讲重点讲解了矩阵的LU分解。对于一个给定矩阵A,可以表现为一个下三角矩阵和一个上三角矩阵乘积的形式：
\[A=LU\]其中上三角矩阵U是求解方程组的初步中间产物。由这一步开始，逐步求解靠后的主元，再回代至方程，以求解更多的未知数主元。重复这个步骤，直到完成所有未知数的求解。
NumPy中，并没有提供矩阵的LU分解功能。可能是因为觉得L、U矩阵用途并不是那么广泛，并且可以直接用方程求解来替代。
如果需要用到的话，通常方式是使用其它软件包替代，比如SciPy。
这里也提供一个架构于NumPy之上的子程序，来完成LU分解的功能。子程序内部是将矩阵类型转换为数组类型，从而方便遍历。接着是使用手工消元相同的方式循环完成LU分解。
需要说明的是，这类附带了子程序的Python片段，建议还是保存到一个文本文件中，以脚本方式执行。在交互式方式下很容易出现各种错误。
import numpy as np
def LUdecomposition(mat):
A=np.array(mat)
n=len(A[0])
L = np.zeros([n,n])
U = np.zeros([n, n])
for i in range(n):
L[i][i]=1
if i==0:
U[0][0] = A[0][0]
for j in range(1,n):
U[0][j]=A[0][j]
L[j][0]=A[j][0]/U[0][0]
else:
for j in range(i, n):#U
temp=0
for k in range(0, i):
temp = temp+L[i][k] * U[k][j]
U[i][j]=A[i][j]-temp
for j in range(i+1, n):#L
temp = 0
for k in range(0, i ):
temp = temp + L[j][k] * U[k][i]
L[j][i] = (A[j][i] - temp)/U[i][i]
return np.mat(L),np.mat(U)
A=np.mat("1 2;3 4")
print(LUdecomposition(A))
程序执行可以获得类似这样的输出：
(matrix([[1., 0.],
[3., 1.]]), matrix([[ 1., 2.],
[ 0., -2.]]))
偏重于计算分析的SymPy则直接内置了LU分解功能，对速度不敏感的场合，使用SymPy做LU分解，再转换到NumPy矩阵也是一个好办法：
>>> As=sp.Matrix(np.mat("1 2;3 4"))
>>> L,U,_=As.LUdecomposition()
>>> L
Matrix([
[1, 0],
[3, 1]])
>>> U
Matrix([
[1, 2],
[0, -2]])
LU分解那一行的代码，使用下划线忽略的部分是函数返回的行交换矩阵。
在消元过程中，对应主元位置如果为0的话会导致消元失败，此时会产生行交换。这种情况下，会由单位矩阵I变换而来的行交换矩阵先同矩阵A相乘，从而将主元为0的行交换到其它行，保证消元的顺利进行。
使用Python辅助解方程，这些步骤都是很少需要手工操作了，如果有必要，就自行赋值给矩阵变量保留吧。
顺便提一句，讲到置换矩阵的时候，教授还提到了对于一个n*n的方阵，置换矩阵可能有多少种呢？答案是n!，也就是n的阶乘。
在Python内置的数学库、NumPy、SymPy中，都有求阶乘的函数：
>>> math.factorial(4) #Python内置数学库求阶乘
24
>>> np.math.factorial(4) #numpy求阶乘
24
>>> sp.factorial(4) #sympy求阶乘
24
第四讲还介绍了矩阵的转置，这是线性代数中使用极为高频的功能：
>>> A=np.mat("1 2;3 4") #定义一个numpy矩阵
>>> As=sp.Matrix(A) #定义一个相同的sympy矩阵
>>> A.T #numpy求转置
matrix([[1, 3],
[2, 4]])
>>> As.T #sympy求转置
Matrix([
[1, 3],
[2, 4]])
简化行阶梯矩阵、0空间基、特解、通解
课程第五至第十讲围绕着矩阵的四个基本空间展开。推导和计算很多，但都是基础线性组合，用Python当成计算器就够用了。
在空间维度判断方面，我们倒是能帮上一些小忙，就是计算矩阵的轶。
矩阵的行空间、列空间轶都是相同的。0空间维度是n-r,左0空间维度是m-r。
>>> A=np.mat("1 2 3 1;1 1 2 1;1 2 3 1") #numpy在这里只是帮助简化输入
>>> As=sp.Matrix(A)
>>> np.linalg.matrix_rank(A) #numpy求矩阵的轶
2
>>> As.rank() #sympy求矩阵的轶
2
如果方程组满轶，也就是方程组有解的情况下，开始一节介绍的解线性方程组很不错。
非满轶的情况，求方程组的特解和通解。将矩阵化简为“简化行阶梯矩阵（Reduced Row Echelon Form）”会非常有用。可惜的是，NumPy依然没有提供内置的支持。自己实现的话，代码还挺长的，远不如使用现有的软件包更方便。所以在这里我们主要看SymPy的实现：
>>> A=np.mat("1 2 3 1;1 1 2 1;1 2 3 1")
>>> As=sp.Matrix(A)
>>> As.rref()
(Matrix([
[1, 0, 1, 1],
[0, 1, 1, 0],
[0, 0, 0, 0]]), (0, 1))
#另一个例子
>>> Bs=sp.Matrix(np.mat("1 2 2 2;2 4 6 8; 3 6 8 10"))
>>> Bs.rref()
(Matrix([
[1, 2, 0, -2],
[0, 0, 1, 2],
[0, 0, 0, 0]]), (0, 2))
函数返回一个RREF矩阵还有一个元组，元组指示RREF矩阵中主元所在的列。这个元组是非常必要的，在第二个例子中就能明显看出，主列并不一定是从左到右相邻排列的。
此时，可以通过RREF最下面的全0行跟方程组b向量的情况判断函数可解性。以及根据自由变量F子矩阵的情况获得方程的0空间解。
当然，如同前面的解方程一样，SymPy中直接提供了函数获取0空间解。RREF这些中间过程主要用于分析学习，对于使用工具自动解方程意义并不大：
>>> As.nullspace()
[Matrix([
[-1],
[-1],
[ 1],
[ 0]]), Matrix([
[-1],
[ 0],
[ 0],
[ 1]])]
>>> Bs.nullspace()
[Matrix([
[-2],
[ 1],
[ 0],
[ 0]]), Matrix([
[ 2],
[ 0],
[-2],
[ 1]])]
方程组的通解包括特解和0空间基两部分。前面获得的是0空间的基。特解则需要方程式右侧向量的配合：
#设置b值,这代表方程组右侧的常数
>>> b=sp.Matrix(np.mat("1;5;6"))
>>> sp.linsolve((As,b)) #(As,b)这种写法是将As矩阵跟b矩阵组合在一起，以增广矩阵的方式求解
EmptySet #参考前面rref矩阵，因为有全0行，b不符合可解性要求，所以方程组使用b向量不可解
>>> sp.linsolve((Bs,b))
FiniteSet((-2*tau0 + 2*tau1 - 2, tau0, 3/2 - 2*tau1, tau1))
Bs矩阵同b向量的组合获得一个有限集的解，那么这个解中的tau0/tau1是什么意思呢？
参考前面的rank计算或者rref矩阵，我们知道Bs矩阵有两个自由变量（由n-r得来），tau0/tau1就是这两个自由变量。这也是因为我们没有定义未知数符号所导致的自动命名。如果需要，我们可以定义x1/x2…这样的未知数。不过这不是我们的重点，请忽略这个命名。
方程的特解是当自由变量为0的时候，方程的解。因此将tau0/tau1都设为0，化简一下，很容易得到方程的特解为：
(-2,0,3/2,0)。
再结合上面计算的Bs矩阵在0空间的2个基，就是方程组的通解:
点积、获取指定行向量和列向量、正交判定
点积也称作点乘、内积，是向量、矩阵中最常用的一种运算。NumPy和SymPy中都有支持。
#numpy中计算点积
>>> a=np.mat("1;2;3")
>>> b=np.mat("4;5;6")
>>> a.T.dot(b)
matrix([[32]])
#sympy中计算点积
>>> a1=sp.Matrix(a)
>>> b1=sp.Matrix(b)
>>> a1.T.dot(b1)
32
矩阵运算中，使用一个向量的转至乘另外一个向量，或者乘自己用于求平方，都是非常常用的使用方法。在使用NumPy做运算的时候要特别注意一点，这样点积的结果仍然是一个矩阵，只是1维*1维。
在线性代数课程上，都会直接把这个点积结果继续用于计算，但在使用NumPy的时候，要特别注意应当将其转换为浮点数，然后再用于计算。不然会出现矩阵维度不符的错误。示例如下：
>>> a.T.dot(b) * a
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
File "/usr/local/lib/python3.8/site-packages/numpy/matrixlib/defmatrix.py", line 218, in __mul__
return N.dot(self, asmatrix(other))
File "<__array_function__ internals>", line 5, in dot
ValueError: shapes (1,1) and (3,1) not aligned: 1 (dim 1) != 3 (dim 0)
>>> float(a.T.dot(b))
32.0
>>> float(a.T.dot(b)) * a
matrix([[32.],
[64.],
[96.]])
SymPy作为主要用于实验分析的符号计算工具，点积运算的结果直接就是可继续用于计算的数字，不需要另行转换。
获取矩阵的特定行向量和列向量，在NumPy/SymPy中都是重载了Python语言的列表（数组）操作符，所以方法都是相同的。
需要注意的是在数学中，矩阵行列的计数通常从1开始，第1行、第2行…第1列、第2列。而在Python中，遵循了计算机语言一贯的习俗，是从0开始计数的。Python中矩阵的第0行，就相当于通常数学课程上所说的第1行。
先来看获取矩阵中特定元素的方法：
#一下方法由numpy演示，在sympy中是相同的，不再另外演示
>>> a=np.mat("1 2 3;4 5 6;7 8 9")
>>> a
matrix([[1, 2, 3],
[4, 5, 6],
[7, 8, 9]])
#获取a矩阵第0行、第2列的元素，也既通常所说的第一行、第三列
>>> a[0,2]
3
#修改矩阵中某一特定元素的值
>>> a[0,2]=24
>>> a
matrix([[ 1, 2, 24],
[ 4, 5, 6],
[ 7, 8, 9]])
获取行向量、列向量，相当于获取矩阵某一行或者某一列所有的数据。在Python中，使用’:’字符放置在行、列参数的位置，就代表获取完整行或者列的数据：
#获取第1列的列向量，也就是通常数学课上所说的第二列(后略)
#在行参数位置使用':'字符，表示任意一行的数据都要，从而组成了列
>>> a[:,1]
matrix([[2],
[5],
[8]])
#获取第0行的行向量
>>> a[0,:]
matrix([[ 1, 2, 24]])
#判断第1列同第0列是否正交
>>> a[:,1].T.dot(a[:,0])
matrix([[78]]) #结果不为0，表示没有正交
点积和正交判断是在课程第十四讲中引入的。
判断两个向量是否正交，就是用一个向量的转置，点积另外一个向量。相互正交的向量，点积结果为0。上面的例子说明，我们随意定义的矩阵，前两列并不正交。
单位矩阵I的每一行、每一列都是正交的，我们测试一下：
#定义一个5x5的单位矩阵，eye方法默认返回是多维列表，在本实验中可以直接使用，
#但为了良好的习惯，还是转换为mat类型。
>>> I=np.mat(np.eye(5))
>>> I
matrix([[1., 0., 0., 0., 0.],
[0., 1., 0., 0., 0.],
[0., 0., 1., 0., 0.],
[0., 0., 0., 1., 0.],
[0., 0., 0., 0., 1.]])
#判断第0行跟第1行的向量是否正交
>>> I[0,:].dot(I[1,:].T)
matrix([[0.]]) #说明两行是正交的
此外在NumPy和SymPy的运算符重载中，乘法运算符’*‘直接就定义为了点积运算，是可以直接使用的：
>>> I[0,:] * I[1,:].T
matrix([[0.]])
方程组的最优解
内容同样来自课程第十四讲。
在实际的应用中，方程组的数据来源经常是测量的结果。在一组实验中，测到了多组结果，这代表方程有多行。但因为测量误差以及干扰因素，这些不准确的测量值所形成的方程组，往往因为误差导致的矛盾因素是无解的。
这时候，通过计算测量数据到方程组矩阵列空间的投影信息，形成新的方程组，可以得到最接近真实结果的答案，这就是最优解。
对于一个原始方程：
其最优解方程为：
\[A^TA\hat{x}=A^Tb\]求得的\(\hat{x}\)就是方程的最优解。它并不是原来的x,而是最接近合理值的近似解，所以称为最优解。
下面使用SymPy为例演示方程组求解最优解,NumPy可以使用同样的方法：
>>> a=sp.Matrix(np.mat("1 1; 1 2; 1 5")) #定义A矩阵
>>> b=sp.Matrix(np.mat("1;2;2")) #定义向量b
#先尝试求解Ax=b
>>> a.solve(b) #报错信息提示A矩阵不可逆，无法求解
Traceback (most recent call last):
File "/usr/local/lib/python3.8/site-packages/sympy/matrices/solvers.py", line 727, in _solve
soln, param = M.gauss_jordan_solve(rhs)
File "/usr/local/lib/python3.8/site-packages/sympy/matrices/matrices.py", line 2212, in gauss_jordan_solve
return _gauss_jordan_solve(self, B, freevar=freevar)
File "/usr/local/lib/python3.8/site-packages/sympy/matrices/solvers.py", line 553, in _gauss_jordan_solve
raise ValueError("Linear system has no solution")
ValueError: Linear system has no solution
During handling of the above exception, another exception occurred:
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
File "/usr/local/lib/python3.8/site-packages/sympy/matrices/matrices.py", line 2218, in solve
return _solve(self, rhs, method=method)
File "/usr/local/lib/python3.8/site-packages/sympy/matrices/solvers.py", line 734, in _solve
raise NonInvertibleMatrixError("Matrix det == 0; not invertible.")
sympy.matrices.common.NonInvertibleMatrixError: Matrix det == 0; not invertible.
#使用映射的方式将b投影到A的列空间
>>> a1=a.T*a
>>> b1=a.T*b
#求解最优解
>>> a1.solve(b1)
Matrix([ #得到的最优解
[15/13],
[ 5/26]])
投影矩阵
投影矩阵的概念来自课程第十五讲。
使用投影矩阵公式可以求得矩阵A的投影矩阵：
下面以NumPy为例，演示计算投影矩阵：
#定义一个求投影矩阵的子程序
def project_matrix(mat):
return mat*np.linalg.inv(mat.T*mat)*mat.T
#定义一个矩阵A,注意A矩阵需要是列满轶的
>>> a=np.mat("3 7;1 5;2 4")
#求A的映射矩阵
>>> p=project_matrix(a)
>>> p
matrix([[ 0.65384615, 0.11538462, 0.46153846],
[ 0.11538462, 0.96153846, -0.15384615],
[ 0.46153846, -0.15384615, 0.38461538]])
#下面验证几个投影矩阵的性质
#1.投影矩阵是对称的
#因为numpy是数值计算，小数值很长，所以提供了专门方法np.allclose()用于比较两个矩阵
>>> np.allclose(p.T,p)
True #返回True,表示两个矩阵相等
#2.投影矩阵的平方等于投影矩阵自身,表示多次投影也是同一个垂足本身
>>> np.allclose(p**2,p)
True
#3.一个可逆矩阵，其投影矩阵为单位矩阵I。这代表对于可逆矩阵，b直接就在其列空间，所以投影为自身
>>> a=np.mat(np.random.randint(100,size=(3,3))) #随机生成的矩阵为多维列表模式，要转换为矩阵
>>> project_matrix(a)
matrix([[ 1.00000000e+00, 1.02695630e-15, -8.25728375e-16],
[ 5.20417043e-16, 1.00000000e+00, 3.69496100e-16],
[-1.66533454e-16, 2.77555756e-17, 1.00000000e+00]])
>>> np.allclose(project_matrix(a),np.eye(3))
True
正交矩阵和正交化法
这部分内容来自课程第十七讲。
按照教授的说法，标准正交矩阵是能得到的最好的矩阵，有很多优良的性质，便于计算和分析。
标准正交矩阵每一列都互相垂直，模长为1。通常把标准正交矩阵记为Q。
但很可惜，通常的矩阵都不是标准正交矩阵。课程中介绍了格拉姆-施密特(Graham-Schmidt)正交化法，将一个列满轶的矩阵A，转换为一个由标准正交向量组构成的矩阵Q。
SymPy内置了这个算法，用于将一组线性无关的向量正交化，来看看示例：
import sympy as sp
vlist=[] #定义一个列表用于保存多个希望进行正交化的列向量
Q=sp.zeros(5,5) #定义一个空白的5*5矩阵，用于保存最终生成的标准正交矩阵
for i in range(0,5): #循环5次，随机生成5个向量
vlist.append(sp.randMatrix(5,1))
#格拉姆-施密特正交化法，orthonormal参数表示希望进行标准化
qlist=sp.GramSchmidt(vlist,orthonormal=True)
for i in range(0,5): #循环5次，使用我们前面介绍过读写矩阵特定一列的方法，
Q[:,i]=qlist[i] #将标准正交向量组成标准正交矩阵
print(Q) #输出标准正交矩阵
print(Q.T*Q) #测试标准正交矩阵的特性，转置*自身=单位矩阵I
这个小程序段需要单独保存为一个脚本来执行，输出因为SymPy符号计算的特点，会变得极为复杂。这种复杂主要来自于标准化除以模长所导致的分数化。
Matrix([[41*sqrt(16898)/8449, -94603*sqrt(269690238118)/134845119059, -15748588*sqrt(1302577778973635969)/186082539853376567, -425130641780*sqrt(23692552673045367710006107)/3384650381863623958572301, 5002949*sqrt(290294034089473)/290294034089473], [45*sqrt(16898)/16898, 317381*sqrt(269690238118)/269690238118, 781784552*sqrt(1302577778973635969)/1302577778973635969, -111786279859*sqrt(23692552673045367710006107)/23692552673045367710006107, 3273660*sqrt(290294034089473)/290294034089473], [33*sqrt(16898)/8449, 129105*sqrt(269690238118)/134845119059, -718289251*sqrt(1302577778973635969)/1302577778973635969, 1062149555036*sqrt(23692552673045367710006107)/23692552673045367710006107, -3858716*sqrt(290294034089473)/290294034089473], [33*sqrt(16898)/16898, -158161*sqrt(269690238118)/269690238118, 7790318*sqrt(1302577778973635969)/1302577778973635969, 3492057859131*sqrt(23692552673045367710006107)/23692552673045367710006107, 9758746*sqrt(290294034089473)/290294034089473], [26*sqrt(16898)/8449, -101825*sqrt(269690238118)/134845119059, 404026822*sqrt(1302577778973635969)/1302577778973635969, 1225299826763*sqrt(23692552673045367710006107)/23692552673045367710006107, -12017690*sqrt(290294034089473)/290294034089473]])
Matrix([[1, 0, 0, 0, 0], [0, 1, 0, 0, 0], [0, 0, 1, 0, 0], [0, 0, 0, 1, 0], [0, 0, 0, 0, 1]])
毕竟不是编程的课程，所以虽然是很短一个小程序，非IT专业的同学看起来可能也会觉得晕。这是由于SymPy中内置的格拉姆-施密特算法主要用于处理向量所导致的。我们不得不把矩阵变为向量，完成正交化后，再转换回矩阵。
实际上有更好的办法，就是使用QR分解。QR分解计算起来更麻烦，在课程中并没有介绍，不过还是老话，计算机最不怕的就是清晰的计算。
QR分解的大意是，任何一个列满轶的矩阵A,都可以分解为一个标准正交向量Q和一个上三角矩阵R的乘积形式。上三角矩阵前面见过，就是我们使用高斯消元的中间步骤产物U。
SymPy和NumPy中都内置了QR分解算法，请看示例：
#先是sympy的操作
>>> a1=sp.randMatrix(3,3) #随机生成一个3*3的矩阵，这次用小一点的维度，容易看清楚
>>> q,r=a1.QRdecomposition() #QR分解
>>> q #标准正交矩阵
Matrix([
[33*sqrt(4166)/4166, 1257*sqrt(736017635)/43295155, 17*sqrt(706690)/41570],
[31*sqrt(4166)/4166, 379*sqrt(736017635)/147203527, -147*sqrt(706690)/141338],
[23*sqrt(4166)/2083, -16607*sqrt(736017635)/736017635, 144*sqrt(706690)/353345]])
>>> r #上三角矩阵
Matrix([
[2*sqrt(4166), 2034*sqrt(4166)/2083, 4415*sqrt(4166)/4166],
[ 0, 3*sqrt(736017635)/2083, 219945*sqrt(736017635)/147203527],
[ 0, 0, 4939*sqrt(706690)/141338]])
>>> q[:,0].T.dot(q[:,1]) #验证第0列跟第1列垂直
0
>>> q[:,0].T.dot(q[:,0]) #验证列模长为1
1
>>> q.T * q #标准正交矩阵，逆*自身=I
Matrix([
[1, 0, 0],
[0, 1, 0],
[0, 0, 1]])
>>> q.T == q**-1 #验证标准正交矩阵重要特征：逆=转置
True
#下面是numpy操作部分，生成一个numpy随机矩阵
>>> a=np.random.randint(100,size=(3,3))
>>> a
array([[29, 47, 45],
[48, 0, 76],
[30, 84, 64]])
>>> q,r=np.linalg.qr(a) #使用numpy的QR分解
>>> q.T * q #此时q是多重列表类型，进行矩阵操作会得到错误的结果
array([[ 0.207911 , -0.19433572, 0.40185179],
[-0.19433572, 0.38341184, 0.16081412],
[ 0.40185179, 0.16081412, 0.22721918]])
>>> q=np.mat(q) #将q转换为矩阵，这也是我们前面一再强调的，一定要用矩阵类型做矩阵运算
>>> q.T * q #验证转置*自身=I，输出结果请注意看e幂小数点的位置
matrix([[ 1.00000000e+00, 2.37714715e-17, -1.52168377e-16],
[ 2.37714715e-17, 1.00000000e+00, 1.37111751e-16],
[-1.52168377e-16, 1.37111751e-16, 1.00000000e+00]])
行列式、伴随矩阵、特征值、特征向量
这几个概念可以说是线性代数的核心，因为计算太复杂，贯穿了多讲内容，从第十八讲一直延续到了第二十一讲。
其中为了降低行列式的计算量，还穿插了代数余子式。但计算机的发展让这些复杂计算都成为了一行函数的事情，所以很多基本的加法、乘法的运算，我们就忽略掉了。
这部分没有太多可说的，直接用示例来说明吧：
#使用numpy随机生成一个3*3矩阵
>>> a=np.random.randint(10,size=(3,3))
>>> a #先试一下生成的矩阵
array([[4, 0, 0],
[2, 3, 5],
[1, 7, 1]])
>>> np.linalg.det(a) #numpy计算行列式值
-127.99999999999997 #numpy可怜的精度误差
>>> a1=sp.Matrix(a) #同样的矩阵，生成一个sympy矩阵
>>> a1.det() #sympy计算行列式
-128
# 伴随矩阵本身是为了降低求逆的难度而出现的
# 但这种中间结果numpy/sympy均没有提供，需要使用的话只能用逆反求
>>> np.linalg.det(a)*np.linalg.inv(a) #numpy求伴随矩阵
array([[-32., -0., -0.],
[ 3., 4., -20.],
[ 11., -28., 12.]])
>>> a1.det()*a1.inv() #sympy求伴随矩阵
Matrix([
[-32, 0, 0],
[ 3, 4, -20],
[ 11, -28, 12]])
#numpy求特征值及特征向量
>>> l,v = np.linalg.eig(a)
>>> v=np.mat(v) #别忘了把列表类型转换成矩阵类型
>>> l
array([-4., 8., 4.]) #3*3矩阵，得到3个特征值
#验证三个特征值的和=矩阵对角线主元之和(迹，trace)
>>> np.allclose(l[0]+l[1]+l[2],a[0,0]+a[1,1]+a[2,2])
True
#验证三个特征值的乘积=行列式值
>>> np.allclose(l[0]*l[1]*l[2],np.linalg.det(a))
True
>>> v #三列分别代表三个与上面特征值对应的特征向量
matrix([[ 0. , 0. , 0.86454916],
[ 0.58123819, -0.70710678, -0.29718877],
[-0.81373347, -0.70710678, -0.40525742]])
#验证公式 A * x = lambda * x
>>> np.allclose(l[0]*v[:,0],a*v[:,0])
True
>>> a1.eigenvects() #使用sympy求矩阵特征值、特征向量
[(-4, 1, [Matrix([
[ 0],
[-5/7],
[ 1]])]), (4, 1, [Matrix([
[-32/15],
[ 11/15],
[ 1]])]), (8, 1, [Matrix([
[0],
[1],
[1]])])]
#结果比较复杂，需要解释一下，首先整体是一个列表类型，列表的每一项包含有3个元素：
# 元素1：特征值
# 元素2：本特征值对应特征向量的数量
# 元素3：一个特征向量组成的数组，数组长度跟元素2的数量相同
# 本例中的特征值3个，没有重复，所以特征值对应特征向量数量都是1，后面的数组也都只有一个特征向量
对角矩阵和对角化
这部分内容来自课程第二十二讲。
矩阵A对角化为Λ的公式为：
S是矩阵A特征向量所组成的矩阵。
下面用numpy示例一下：
#生成一个3*3矩阵
>>> a=np.mat("5 2 0;4 4 0;7 0 0")
>>> a
array([[5, 2, 0],
[4, 4, 0],
[7, 0, 0]])
#求特征值特征向量
>>> l,v=np.linalg.eig(a)
>>> v=np.mat(v) #转换成矩阵类型，虽然并不一定必要，但习惯很重要
>>> np.linalg.inv(v).dot(a).dot(v) #S^-1 * A * S
matrix([[ 0.00000000e+00, 6.33735745e-16, 1.15838209e-15],
[ 0.00000000e+00, 1.62771868e+00, -7.40457886e-16],
[ 0.00000000e+00, 1.76816846e-16, 7.37228132e+00]])
#毕竟不是I，得到的对角矩阵看上去不直观，做一个取整来便于观察，
#但请记住np.round在最终结果之前一定不要使用，否则很影响精度
>>> np.round(np.linalg.inv(v).dot(a).dot(v))
matrix([[ 0., 0., 0.], #获得的对角阵
[ 0., 2., -0.],
[ 0., 0., 7.]])
嗯，为了验证课程中的公式，故意搞复杂了点。这样的计算其实完全没有必要，对角化矩阵实际就是矩阵特征值排列在对角线所组成的矩阵。所以实际上把特征值乘单位矩阵I，转化到对角线就好了：
>>> l,v=np.linalg.eig(a)
>>> l*np.eye(3)
array([[0. , 0. , 0. ],
[0. , 1.62771868, 0. ],
[0. , 0. , 7.37228132]])
SymPy也可以使用对角化公式计算，但SymPy计算的特征向量需要自己解析、组合成矩阵S，有点麻烦。
幸运的是，SymPy直接提供了矩阵对角化的函数：
#直接使用上面numpy的矩阵
>>> a1=sp.Matrix(a)
>>> S,D=a1.diagonalize()
>>> S #特征向量组成的矩阵
Matrix([
[0, 9/14 - sqrt(33)/14, sqrt(33)/14 + 9/14],
[0, 3/7 - sqrt(33)/7, 3/7 + sqrt(33)/7],
[1, 1, 1]])
>>> D #得到的对角矩阵
Matrix([
[0, 0, 0],
[0, 9/2 - sqrt(33)/2, 0],
[0, 0, sqrt(33)/2 + 9/2]])
SymPy的计算结果看上去总是那么工整。既然有了S和Λ，是不是想用SymPy算回到矩阵A验证一下？不不不，你不会想那么做的，不信我给你练练：
>>> S*D*S.inv()
Matrix([
[24*(9/14 - sqrt(33)/14)*(9/2 - sqrt(33)/2)/(49*(198/343 - 18*sqrt(33)/343)) + (7/6 - 7*sqrt(33)/66)*(sqrt(33)/14 + 9/14)*(sqrt(33)/2 + 9/2), (9/14 - sqrt(33)/14)*(9/2 - sqrt(33)/2)*(7 + 7*sqrt(33))/(-66 + 6*sqrt(33)) + (-7/12 + 7*sqrt(33)/44)*(sqrt(33)/14 + 9/14)*(sqrt(33)/2 + 9/2), 0],
[ 24*(3/7 - sqrt(33)/7)*(9/2 - sqrt(33)/2)/(49*(198/343 - 18*sqrt(33)/343)) + (3/7 + sqrt(33)/7)*(7/6 - 7*sqrt(33)/66)*(sqrt(33)/2 + 9/2), (3/7 - sqrt(33)/7)*(9/2 - sqrt(33)/2)*(7 + 7*sqrt(33))/(-66 + 6*sqrt(33)) + (-7/12 + 7*sqrt(33)/44)*(3/7 + sqrt(33)/7)*(sqrt(33)/2 + 9/2), 0],
[ 24*(9/2 - sqrt(33)/2)/(49*(198/343 - 18*sqrt(33)/343)) + (7/6 - 7*sqrt(33)/66)*(sqrt(33)/2 + 9/2), (9/2 - sqrt(33)/2)*(7 + 7*sqrt(33))/(-66 + 6*sqrt(33)) + (-7/12 + 7*sqrt(33)/44)*(sqrt(33)/2 + 9/2), 0]])
看起来很晕是吧？
符号计算有符号计算的优点，但缺点也那么显而易见。速度慢就不说了，复杂计算的时候，表达式化简能力，特别是灵活程度毕竟不能同人相比，这就是一个很典型的例子。这样的结果，肯定是还不如用NumPy计算的近似值。
怀疑计算出了错？那倒肯定没有，我们把上面矩阵用NumPy计算出最终数值看一下：
>>> np.mat(P*D*P.inv(),dtype=float)
matrix([[5.00000000e+000, 2.00000000e+000, 0.00000000e+000],
[4.00000000e+000, 4.00000000e+000, 0.00000000e+000],
[7.00000000e+000, 4.72728505e-125, 0.00000000e+000]])
跟最初的矩阵A是吻合的，毋庸置疑。
矩阵乘幂、矩阵法求斐波那契数列、绘图
同样来自课程第二十二讲，对角化矩阵的一种典型应用就是简化矩阵的幂运算。
对于一个高维矩阵的高次幂来讲，如果手工计算，其复杂程度是可想而知的。而通过对角化后的矩阵，矩阵幂的运算可以简化很多：
使用计算机之后，这种简化手段意义显得不再那么显著。但这种思想还是非常有帮助。
比如对于计算序列值的时候，这成为了一种全新的思路，类似符合\(u_{k+1} = Au_k\)这样公式的数字序列，都可以用这个思路来计算。在\(u_0\)已知的情况下，公式可以变形为\(u_k=A^ku_0\)。
以电脑编程入门必学的斐波那契数列生成为例，通常是使用循环顺序生成（此处略去程序）。
使用矩阵运算的思想，联立方程（方程请参考原课程）可以得到矩阵：
以及初始向量为：
\[u_0 = \left[\begin{matrix}1\\1\\\end{matrix}\right]\]想得到序列中任意一个数值，无需循环，一次计算就可以直接得到。来看一下获取斐波那契数列的代码片段：
import numpy as np
#获取指定位置斐波那契数列值的子程序
def Fibonacci_Matrix_tool(n):
Matrix = np.matrix("1 1;1 0", dtype='int64') #定义矩阵，使用整数类型，速度可以更快
# 因为u0向量刚好是[1;1],省略了，所以计算完成后仍是一个矩阵，但我们只需要左上角的值
return np.linalg.matrix_power(Matrix, n)[0,0] #矩阵求幂，[0,0]是获取矩阵左上角元素值
print(Fibonacci_Matrix_tool(0)) #序列中第0个元素
print(Fibonacci_Matrix_tool(1)) #序列中第1个元素
print(Fibonacci_Matrix_tool(2)) #序列中第2个元素
print(Fibonacci_Matrix_tool(100)) #序列中第100个元素
把上面代码保存为脚本文件，执行后获得输出为：
1
1
2
1298777728820984005
线性代数是研究向量和空间的学科，绘图能够在很大程度上帮助问题的思考。前面一直没有合适的例子，在这里牵强的引入一下，就是将计算结果用绘图的方式展示出来。
接着上面程序代码继续在后面加入：
import matplotlib.pyplot as plt
x = np.arange(0,20,1)
y = []
for i in x:
y.append(Fibonacci_Matrix_tool(i))
plt.plot(x,y)
plt.show()
程序执行前首先要安装绘图软件包pip3 install matplotlib。安装完成后再次执行程序，除了得到上面相同的4行输出外，还会绘出如下的曲线，表现了斐波那契数列序号同数值之间的关系：
在Win10的Linux子系统使用绘图功能的话，还需要配置X11协议的远程显示，请参考这篇文章：win10配置linux子系统使用python绘图并显示
马尔科夫矩阵
内容来自课程第二十四讲。
马尔科夫矩阵也叫状态转移矩阵，用于表现在总数不变的情况下，随着时间的延续，不同节点之间某项指标的动态情况。
课程中的例子是讲伴随着硅谷的崛起，加州人口逐渐增加。而相邻麻省人口受加州经济吸引，移居至加州，从而导致本省人口减少。这种情况持续一段时间之后，最终形成一个稳态。例子中的数字显然就是示意性。
马尔科夫矩阵代表了使用矩阵思想，在具体事例中的应用。
下面代码使用课程中的数据做了一个计算推演，并且绘制了曲线图供参考。代码需要保存为脚本文件执行：
import matplotlib.pyplot as plt
import numpy as np
current_status = np.array([0, 1000]) #初始状态
transfer_matrix = np.array([[0.9, 0.2], [0.1, 0.8]]) #马尔科夫转移矩阵
pc=[] #加州人口变动列表
pm=[] #麻省人口变动列表
for i in range(100): #假设延续100年
current_status = np.dot(transfer_matrix,current_status) #计算下一个状态
pc.append(current_status[0]) #拆分数据到两个两个列表，以便分别加图例注释
pm.append(current_status[1]) #不拆分直接绘制也可以，但无法单独加图例注释
plt.rcParams['font.sans-serif']=['Songti SC'] #本行在不同电脑可能需要修改，这是在mac电脑上选择中文字体
plt.plot(pc,label="加州人口") #绘图及图例注释
plt.plot(pm,label="麻省人口")
plt.legend() #显示图例
plt.show()
程序执行绘制的曲线如下：
程序在绘图中，使用了中文的图例注释。如果不做任何处理，出现的会是乱码。plt.rcParams['font.sans-serif']=['Songti SC']这一行代码，就是将默认绘图字体名sans-serif指定到使用系统宋体，从而正常显示中文。在不同的电脑上，要根据自己的电脑字体名称设置，选择一个替换。
对称矩阵、复矩阵
这部分内容来自课程第二十五、二十六讲。
对于实数矩阵来说，对称矩阵就是转置与自身相同的矩阵，判定起来很容易。
实对称矩阵还有特征值全部为实数、特征向量相互垂直两个重要特性。
# numpy定义一个对称矩阵
>>> a=np.mat("1 2 3;2 2 6;3 6 4")
>>> a
matrix([[1, 2, 3],
[2, 2, 6],
[3, 6, 4]])
>>> a.T == a #判定转置是否等于本身
matrix([[ True, True, True],
[ True, True, True],
[ True, True, True]])
>>> np.allclose(a.T,a) #使用np.allclose方法判定
True
>>> e,v = np.linalg.eig(a) #获取特征值特征向量
>>> e #显示3个特征值，全部是实数
array([10.44316671, -0.30514935, -3.13801736])
>>> v = np.mat(v) #特征向量转换为np.mat类型
>>> v[:,0].T.dot(v[:,1]) #验证向量0、向量1垂直
matrix([[-5.55111512e-17]])
>>> v[:,0].T.dot(v[:,2]) #验证向量0、向量2垂直
matrix([[-1.66533454e-16]])
SymPy的相关操作前面都已经出现过，此处不再重复。
复矩阵就是元素中存在复数的矩阵。关键是复数如何表达，NumPy中延续了Python中对复数的定义方式；SymPy中定义了自己的虚数符号类。两种方式都离我们日常数学中的习惯区别很大。
# python及numpy中表示复数
>>> x = 3+2j
>>> x
(3+2j)
# sympy中表示复数
>>> xs= 3+2*sp.I
>>> xs
3 + 2*I
# numpy定义复矩阵
>>> a=np.mat("3 2+7j 4+6j;2-7j 4 8+1j;4-6j 8-1j 5")
>>> a
matrix([[3.+0.j, 2.+7.j, 4.+6.j],
[2.-7.j, 4.+0.j, 8.+1.j],
[4.-6.j, 8.-1.j, 5.+0.j]])
# sympy 定义复矩阵
>>> bs=sp.Matrix([[1,2+3*sp.I],[4+sp.I, 8]])
>>> bs
Matrix([
[ 1, 2 + 3*I],
[4 + I, 8]])
#sympy使用numpy中定义的复矩阵
>>> a1=sp.Matrix(a)
>>> a1
Matrix([
[ 3.0, 2.0 + 7.0*I, 4.0 + 6.0*I],
[2.0 - 7.0*I, 4.0, 8.0 + 1.0*I],
[4.0 - 6.0*I, 8.0 - 1.0*I, 5.0]])
对称复矩阵(埃尔米特矩阵)的定义跟实数矩阵有所区别，在复矩阵中，对称是指矩阵做完共轭、转置操作后，同本身相等。这也意味着，在对称复矩阵对角线上的元素必须都是实数。否则不可能做到共轭后与自身相同。
复矩阵组成的正交矩阵称为酉矩阵。
#numpy中对矩阵取共轭
>>> a.conjugate()
matrix([[3.-0.j, 2.-7.j, 4.-6.j],
[2.+7.j, 4.-0.j, 8.-1.j],
[4.+6.j, 8.+1.j, 5.-0.j]])
#sympy中对矩阵取共轭，跟numpy完全相同
>>> a1.conjugate()
Matrix([
[ 3.0, 2.0 - 7.0*I, 4.0 - 6.0*I],
[2.0 + 7.0*I, 4.0, 8.0 - 1.0*I],
[4.0 + 6.0*I, 8.0 + 1.0*I, 5.0]])
>>> a.H #numpy中对矩阵同时取共轭、转置
matrix([[3.-0.j, 2.+7.j, 4.+6.j],
[2.-7.j, 4.-0.j, 8.+1.j],
[4.-6.j, 8.-1.j, 5.-0.j]])
>>> a1.H #sympy中对矩阵同时取共轭、转置
Matrix([
[ 3.0, 2.0 + 7.0*I, 4.0 + 6.0*I],
[2.0 - 7.0*I, 4.0, 8.0 + 1.0*I],
[4.0 - 6.0*I, 8.0 - 1.0*I, 5.0]])
>>> a.H == a #numpy中判断复矩阵对称
matrix([[ True, True, True],
[ True, True, True],
[ True, True, True]])
>>> a1.H == a1 #sympy中判断复矩阵对称
True
>>> e,v=np.linalg.eig(a) #numpy获取复矩阵的特征向量
>>> np.round(v.H*v) #复对称矩阵的特征向量组成的矩阵是酉矩阵，共轭转置*自身为I
matrix([[ 1.+0.j, -0.+0.j, -0.-0.j],
[-0.-0.j, 1.+0.j, -0.+0.j],
[-0.+0.j, -0.-0.j, 1.+0.j]])
正定矩阵
正定矩阵是课程第二十七讲的内容。
首先经常用到的是正定矩阵的判定。课程中教授提了一个问题：
A矩阵定义如上，右下角c取值是多少，使得A矩阵成为正定矩阵。
老师给了几个人工判定的标准：
- 矩阵为对称方阵。
- 所有特征值为正。
- 所有主元为正。
- 从左上角开始的子对称矩阵行列式为正。
- 对于任意非零向量x,xᵀAx的结果为正。
对于SymPy来讲比较容易，内置提供了正定矩阵判定的方法。NumPy没有内置此种功能，但可以根据上面的标准，用一小段程序来判断，难度也不大。不过NumPy还有一个取巧的办法，NumPy中有矩阵的霍尔斯基分解函数，霍尔斯基分解是要求矩阵为正定矩阵的。如果提供的矩阵参数不是正定矩阵，函数会报错。此时截获这个错误，也可以准确、方便的判定矩阵的正定性。
（霍尔斯基分解不属于本课程内容，这里只是利用它来判断矩阵的正定性，所以不用管其具体功能。）
下面用代码来看看具体方法：
# numpy定义一个函数，使用霍尔斯基分解来判定矩阵的正定性
def is_pos_def(A):
if np.array_equal(A, A.T): #首先需要对称
try:
np.linalg.cholesky(A) #霍尔斯基分解、平方根分解，要求参数必须是正定矩阵
return True
except np.linalg.LinAlgError: #报错就是非正定性
return False
else:
return False
接着以上面的矩阵为例，来实际判定测试一下：
# 定义两个SymPy矩阵，c分别取值19、7
>>> a=sp.Matrix([[2,6],[6,19]])
>>> b=sp.Matrix([[2,6],[6,7]])
# 这次反过来，numpy使用sympy定义的矩阵
>>> a1=np.mat(a,dtype=float)
>>> b1=np.mat(b,dtype=float)
# numpy中使用自定义的函数来判断a1/b1是否为正定矩阵
>>> is_pos_def(a1)
True
>>> is_pos_def(b1)
False
# 直接使用sympy内置属性判定矩阵是否为正定矩阵
>>> a.is_positive_definite
True
>>> b.is_positive_definite
False
回到老师出的问题，课程中讲过，通过xᵀAx>0来判断矩阵的正定性是最全面、可靠的。
以题中的两维矩阵为例，向量也就是两维，假设为：[x1;x2]，把矩阵A代入公式、展开：
可以看出，第一个系数2，就是矩阵A左上角的值；第二个系数12是A第1行第2列及第2行第1列的和；第三个系数就是c了。实际上这来自于行列式的计算。由此可判断c>18可以保证矩阵A是正定矩阵。
对于课程的内容我没有什么要补充的。但是在Python的帮助下，如果将上面公式图示出来，肯定可以帮助我们更深入的理解矩阵A中c取值对于矩阵正定性的影响。
因为上面公式有x1/x2两个变量，加上最终整体公式的取值算作一个维度，我们需要绘制的是三维图。
下面程序中，我们分别使用c=7以及c=20，绘制两幅三维图片。程序使用了matplotlib绘图软件包的mpl_toolkits三维绘图模块。这个模块是matplotlib新版本才引入的，所以如果找不到这个模块的话，请升级matplotlib模块。
# 正定矩阵的正定性分析
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np
F = plt.figure()
ax = Axes3D(F)
x=np.arange(-10,10,1)
y=np.arange(-10,10,1)
X, Y = np.meshgrid(x, y) #网格化，或者叫2D化,实际上是形成X/Y的平面，因为我们是3D图
Z1=2*X**2+12*X*Y+7*Y**2 #使用c=7参数计算此时公式在Z轴的值
Z2=2*X**2+12*X*Y+20*Y**2 #使用c=20参数计算此时公式在Z轴的值
plt.xlabel('x')
plt.ylabel('y')
ax.text(0,0,1800,"c=7",None) #图例文字
ax.plot_surface(X, Y, Z1, rstride=1, cstride=1, cmap='rainbow') #绘制值平面
plt.show()
# 在3D方式无法使用子图，所以先绘制一幅，关闭后再绘制第二幅
F = plt.figure()
ax = Axes3D(F)
plt.xlabel('x')
plt.ylabel('y')
ax.text(0,0,1800,"c=20",None)
ax.plot_surface(X, Y, Z2, rstride=1, cstride=1, cmap='rainbow')
plt.show()
绘制结果请看图：
在第一张图片中，可以看到当c取值7时，有相当一部分图都位于0（Z轴)之下。
而第二张图片中，c取值20，所有曲线都会在0之上了，代表xᵀAx>0，矩阵是正定矩阵。
绘制的三维图片，可以使用鼠标拖动，从各个角度观察。还可以旋转、缩放、保存为图片文件。Python实在是数学学习和研究的利器。
奇异值分解
这是课程第二十九讲的内容。奇异值分解的公式如下：
\[A = U∑V^T\]其中，U是AAᵀ矩阵的特征向量形成的标准正交矩阵；V是AᵀA矩阵的特征向量形成的标准正交矩阵；∑则是两个矩阵特征值开根号后形成的对角矩阵。
SVD分解几乎串起了整个线性代数课程的知识点，手工计算的话，还是相当的麻烦。
NumPy中已经内置了奇异值分解的函数：
>>> a=np.mat("4 4;-3 3")
>>> u, s, vt = np.linalg.svd(a) # 这里vt为V的转置
>>> u
matrix([[-1., 0.],
[ 0., 1.]])
>>> s
array([5.65685425, 4.24264069])
>>> vt
matrix([[-0.70710678, -0.70710678],
[-0.70710678, 0.70710678]])
这次终于碰到了SymPy没有内置对应功能的情况。实际你也看出来了，SVD是典型的数值计算应用，对于符号运算分析的话作用并不大。而且因为运算步骤多，符号计算保留过多的符号操作很容易造成计算溢出，可读性更是没有了保障。
所以在SymPy的官方推荐中，也是使用mpmath运算包完成SVD分解。在新版本的SymPy中，这个包已经分离并且需要单独安装，所以你还不如直接使用NumPy计算了。
上面的计算中，变量s代表了SVD分解之后的∑对角矩阵，实际是AAᵀ矩阵或者AᵀA矩阵特征值再开方的值。使用NumPy做完SVD分解后，直接保存为列表类型。对角阵除了对角线的数据，其它元素都是0，这样保存非常合理。
把列表还原回对角阵，前面介绍了乘单位矩阵I的方法，这里再介绍一个NumPy内置的diag函数，用起来更方便：
>>> np.diag(s)
array([[5.65685425, 0. ],
[0. , 4.24264069]])
最后
本文基本涵盖了MIT18.06版线性代数课程中所需要用到的计算方法。很多章节中，计算量虽然不小，但已经在前面讲过的，就被省略掉了。
希望能为学习线性代数的同学，或者希望从MATLAB迁移至Python的同学带来帮助。
错误疏漏因水平所限难以避免，欢迎指正。
