PYTHON替代MATLAB在线性代数学习中的应用(使用Python辅助MIT 18.06 Linear Algebra学习)
前言
MATLAB一贯是理工科学生的必备神器，但随着中美贸易冲突的一再升级，禁售与禁用的阴云也持续笼罩在高等学院的头顶。也许咱们都应当考虑更多的途径，来辅助咱们的学习和研究工做。
虽然PYTHON和众多模块也属于美国技术的范围，但开源软件的自由度毕竟不是商业软件可比拟的。python
本文是一篇入门性文章，以麻省理工学院（MIT) 18.06版本线性代数课程为例，按照学习顺序介绍PYTHON在代数运算中的基本应用。
介绍PYTHON代数计算的文章很是多，但一般都是按照模块做为划分顺序，在实际应用中仍然有较多困扰。
而按照代数课程学习的顺序，按部就班，集注在最经常使用和最实用的功能上，比较符合典型的应用逻辑。能够用较低的门槛逐步完成PYTHON在线性代数中各项功能的学习和应用。
MIT 2020版本的线性代数课程也已发布，但基本是在18.06版本上的修正。Gilbert教授的年龄已经很大，只录制了一个5节课的串讲。因此系统性仍是18.06版本更为完整。
很讽刺是吧，课程自己也是美国的-_-#。阿Q一下吧，就当是“师夷长技以制夷”。linux
首先给出几个相关连接：
MIT 18.06 Linear Algebra课程主页
B站完整版34讲Gilbert教授课程视频
配套第三版线性代数教材(百度网盘) 提取码：uhvc
最新发行的教材是第5版，建议听课时使用配套的第3版教材。课程完成后，把第5版教材做为辅助读物。否则在章节、内容方面会碰到不少困惑。git
版本选择
PYTHON版本的选择如今已经没有什么困惑了，PYTHON2中止了支持，PYTHON3如今是必选项。我是用Mac电脑，一般使用brew安装PYTHON3，每次有新版本的时候执行brew upgrade会自动升级。不使用内置的PYTHON3是为了防止安装不少扩展库的时候会有系统完整性检查致使的不兼容，不过只是跑一跑数学运算的话倒也不用担忧这些。
Linux各发行版各发行版是Python的开发环境，因此内置的PYTHON3就很好用。使用apt/yum等包管理工具升级的时候会自动完成版本维护。
PYTHON在Windows/Linux/Mac等各平台上兼容性很是好，特别是在数学计算方面基本不用担忧互相之间的通用问题。github
计算模块方面则有了不少的选择，常见的有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库：bash
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课程，如下简称课程）第1、二讲的内容。
在线性代数中，主要涉及3种数据类型，常量、标量（Scalar)、向量(Vector)、矩阵（Matrix)。 不管NumPy仍是SymPy，都直接使用了基本Python类型做为标量，好比：
c1 = 5
而对于向量和矩阵，处理方法则有很大区别，下面先讲述NumPy中的方法。
假设咱们有v一、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]])
第一、2种方法，虽然在不少状况下都能正常的使用，但都算不上规范化的矩阵表示方法。特别是对于向量的表示，向量原本是纵向的，表明矩阵中的一列。但在方法一、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格式输入。
求解线性方程
这也是课程第1、二讲中的内容。方程组是矩阵的起源，也是矩阵最初的目的。以课程中的方程组为例：
能够获得矩阵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.]])
方程组的最优解
内容一样来自课程第十四讲。
在实际的应用中，方程组的数据来源常常是测量的结果。在一组实验中，测到了多组结果，这表明方程有多行。但由于测量偏差以及干扰因素，这些不许确的测量值所造成的方程组，每每由于偏差致使的矛盾因素是无解的。
这时候，经过计算测量数据到方程组矩阵列空间的投影信息，造成新的方程组，能够获得最接近真实结果的答案，这就是最优解。
对于一个原始方程：
\(Ax=b\)
其最优解方程为：
\(A^TA\hat{x}=A^Tb\)
求得的\(\hat\)就是方程的最优解。它并非原来的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的投影矩阵：
\(P=A(A^TA)^{-1}A^T\)
下面以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^{-1}AS\)
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是吻合的，毋庸置疑。
矩阵乘幂、矩阵法求斐波那契数列、绘图
一样来自课程第二十二讲，对角化矩阵的一种典型应用就是简化矩阵的幂运算。
对于一个高维矩阵的高次幂来说，若是手工计算，其复杂程度是可想而知的。而经过对角化后的矩阵，矩阵幂的运算能够简化不少：
\(A^k = SΛ^kS^{-1}\)
使用计算机以后，这种简化手段意义显得再也不那么显著。但这种思想仍是很是有帮助。
好比对于计算序列值的时候，这成为了一种全新的思路，相似符合\(u_{k+1} = Au_k\)这样公式的数字序列，均可以用这个思路来计算。在\(u_0\)已知的状况下，公式能够变形为\(u_k=A^ku_0\)。
以电脑编程入门必学的斐波那契数列生成为例，一般是使用循环顺序生成（此处略去程序）。
使用矩阵运算的思想，联立方程（方程请参考原课程）能够获得矩阵：
\(A = \left[\begin{matrix}1&1\\1&0\\\end{matrix}\right]\)
以及初始向量为：
\(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指定到使用系统宋体，从而正常显示中文。在不一样的电脑上，要根据本身的电脑字体名称设置，选择一个替换。
对称矩阵、复矩阵
这部份内容来自课程第二十5、二十六讲。
对于实数矩阵来讲，对称矩阵就是转置与自身相同的矩阵，断定起来很容易。
实对称矩阵还有特征值所有为实数、特征向量相互垂直两个重要特性。
# 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 = \left[\begin{matrix}2&6\\6&c\\\end{matrix}\right]\)
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分别取值1九、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代入公式、展开：
\(x^TAx = 2x_1^2+12x_1x_2+cx_2^2\)
能够看出，第一个系数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的同窗带来帮助。 错误疏漏因水平所限难以免，欢迎指正。
- 1. 【线性代数公开课MIT Linear Algebra】 实际应用——python中的线性代数(1)
- 2. MIT学习笔记-Linear Algebra(二)
- 3. 【Linear Algebra】线性代数
- 4. MIT 18.06 线性代数总结（Part I）
- 5. MIT 18.06 linear algebra 第三十讲笔记
- 6. 线性代数《Linear Algebra and Its Application》学习总结
- 7. MIT 18.06 线性代数总结（Part II）
- 8. MIT 线性代数 Linear Algebra 20: 行列式的应用和物理意义
- 9. MIT 线性代数 Linear Algebra 15: 投影 projection
- 10. MIT 线性代数 Linear Algebra 32:左逆,右逆,伪逆
- 更多相关文章...
- • 在Spring中使用Redis - Redis教程
- • Redis在Java Web中的应用 - Redis教程
- • 适用于PHP初学者的学习线路和建议
- • TiDB 在摩拜单车在线数据业务的应用和实践
-
每一个你不满意的现在，都有一个你没有努力的曾经。
- 1. 融合阿里云，牛客助您找到心仪好工作
- 2. 解决jdbc(jdbctemplate)在测试类时不报错在TomCatb部署后报错
- 3. 解决PyCharm GoLand IntelliJ 等 JetBrains 系列 IDE无法输入中文
- 4. vue+ant design中关于图片请求不显示的问题。
- 5. insufficient memory && Native memory allocation (malloc) failed
- 6. 解决IDEA用Maven创建的Web工程不能创建Java Class文件的问题
- 7. [已解决] Error: Cannot download ‘https://start.spring.io/starter.zip?
- 8. 在idea让java文件夹正常使用
- 9. Eclipse启动提示“subversive connector discovery”
- 10. 帅某-技巧-快速转帖博主文章（article_content）
- 1. 【线性代数公开课MIT Linear Algebra】 实际应用——python中的线性代数(1)
- 2. MIT学习笔记-Linear Algebra(二)
- 3. 【Linear Algebra】线性代数
- 4. MIT 18.06 线性代数总结（Part I）
- 5. MIT 18.06 linear algebra 第三十讲笔记
- 6. 线性代数《Linear Algebra and Its Application》学习总结
- 7. MIT 18.06 线性代数总结（Part II）
- 8. MIT 线性代数 Linear Algebra 20: 行列式的应用和物理意义
- 9. MIT 线性代数 Linear Algebra 15: 投影 projection
- 10. MIT 线性代数 Linear Algebra 32:左逆,右逆,伪逆