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. Window下Ribbit MQ安裝
- 2. Linux下Redis安裝及集羣搭建
- 3. shiny搭建網站填坑戰略
- 4. Mysql8.0.22安裝與配置詳細教程
- 5. Hadoop安裝及配置
- 6. Python爬蟲初學筆記
- 7. 部署LVS-Keepalived高可用集羣
- 8. keepalived+mysql高可用集羣
- 9. jenkins 公鑰配置
- 10. HA實用詳解
- 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:左逆,右逆,僞逆