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個基，就是方程組的通解:
\[X_{Complete}= \left[\begin{matrix}-2\\0\\\frac32\\0\\\end{matrix}\right]+ k1\left[\begin{matrix}-2\\1\\0\\0\\\end{matrix}\right]+ k2\left[\begin{matrix}2\\0\\-2\\1\\\end{matrix}\right] \]
點積、獲取指定行向量和列向量、正交判定
點積也稱作點乘、內積，是向量、矩陣中最常用的一種運算。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}\)就是方程的最優解。它並不是原來的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指定到使用系統宋體，從而正常顯示中文。在不同的電腦上，要根據自己的電腦字型名稱設定，選擇一個替換。
對稱矩陣、復矩陣
這部分內容來自課程第二十五、二十六講。
對於實數矩陣來說，對稱矩陣就是轉置與自身相同的矩陣，判定起來很容易。
實對稱矩陣還有特徵值全部為實數、特徵向量相互垂直兩個重要特性。
# 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分別取值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代入公式、展開：
\[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的同學帶來幫助。
錯誤疏漏因水平所限難以避免，歡迎指正。