1.5 使用矩陣和線性代數

NumPy數組也可以用作矩陣，這是數學和計算編程的基礎。簡單地說，矩陣就是一個二維數組。矩陣在許多應用中都是核心，如幾何變換和聯立方程，而在統計學等其他領域也作為有用的工具出現。與其他數組相比，矩陣本身只有在我們為其賦予矩陣算術運算時才是獨特的。矩陣具有元素級的加法和減法運算，就像NumPy數組一樣。還有一種稱為標量乘法的第三種運算，就是將矩陣的每個元素乘以一個常數，它是一種不同的矩陣乘法概念。我們將在后面看到，矩陣乘法與其他乘法概念有著根本的不同。

矩陣最重要的屬性之一是其形狀，其定義與NumPy數組完全相同。通常具有m行和n列的矩陣被描述為m×n矩陣。行數與列數相同的矩陣稱為方陣，這些矩陣在向量和矩陣理論中起著特殊的作用。

單位矩陣（大小為n）是一個n×n矩陣，其中（i,i）位置上的元素為1，當i≠j時，（i,j）位置上的元素為零。下面的數組創建例程可以通過指定n的值來生成n×n單位矩陣：

顧名思義，單位矩陣是一個特殊的矩陣，具有如下性質：任何矩陣與單位矩陣相乘等于其本身、兩個逆矩陣的乘積為單位矩陣等。

1.5.1 基本方法和性質

與單位矩陣相關的術語和量有很多。我們在這里只介紹兩個性質，因為稍后會用到它們。這兩個性質分別是矩陣的轉置（transpose）和方陣的跡（trace）。其中矩陣轉置表示行和列互換；跡是方陣主對角線的元素之和，主對角線由元素ai,i構成，即沿矩陣左上角到右下角的對角線上的元素。

NumPy數組可以通過在array對象上調用transpose方法輕松地進行轉置。事實上，由于這是一種常見的操作，數組還有一個便捷的屬性T，它也返回矩陣的轉置。轉置將矩陣（數組）的形狀順序顛倒，使得行變為列，列變為行。例如，如果我們有一個3×2矩陣（3行，2列），那么它的轉置將是一個2×3矩陣，如下例所示：

與矩陣相關的另一個有用的量是跡。如前面的代碼所示，方陣A的跡被定義為沿矩陣左上角到右下角主對角線上的元素之和。跡的公式如下所示：

NumPy數組有一個trace方法，可以返回矩陣的跡：

也可以使用np.trace函數獲取跡，該函數并不是數組的綁定操作。

1.5.2 矩陣乘法

矩陣乘法是對兩個矩陣執行的運算，它保留了兩個矩陣的一些結構和特性。形式上，假設我們有如下l×m矩陣A和m×n矩陣B兩個矩陣：

矩陣A和B的矩陣積C是一個l×n矩陣，其（p,q）項由以下方程給出：

注意，為了定義兩個矩陣的乘法，第一個矩陣的列數必須與第二個矩陣的行數相匹配。通常，我們用AB表示矩陣A和B的乘積（如果根據定義它們能夠相乘的話）。矩陣乘法是一種特殊的運算，它不像大多數其他算術運算那樣滿足交換律：即使AB和BA都可以計算，它們也不一定相等。在實踐中，這意味著矩陣乘法的順序很重要。這源于矩陣代數作為線性映射表示的起源，其中乘法對應于函數的組合。

Python有一個為矩陣乘法保留的運算符@，它是在Python 3.5中新添加的。NumPy數組應用該運算符能夠實現矩陣乘法。請注意，@與數組的逐元素乘法運算符*有本質區別：

單位矩陣是矩陣乘法的中性元素。也就是說，如果A是任意k×m矩陣，而I是m×m單位矩陣，那么AI=A；同樣，如果B是m×k矩陣，那么IB=B。使用NumPy數組可以很容易地檢驗特定示例：

可以看出，輸出的結果矩陣等于原始矩陣。如果我們顛倒A和I的順序并執行乘法IA，結果也是相同的。在1.5.3節中，我們將討論矩陣的逆：若矩陣B與矩陣A相乘得到單位矩陣，則稱矩陣B是A的逆矩陣。

1.5.3 行列式和逆

方陣的行列式在很多應用中都很重要，因為它與求解矩陣的逆有著緊密的聯系。如果矩陣的行數和列數相等，則它是一個方陣。特別地，具有非零行列式的矩陣有一個（唯一的）逆矩陣，這可以轉化為某些方程組的唯一解。矩陣的行列式是通過遞歸定義的。假設我們有一個通用的2×2矩陣，如下所示：

這個通用矩陣A的行列式由以下公式定義：

detA=a1,1a2,2-a1,2a2,1

對于一般的n×n矩陣，其中n＞2，我們通過遞歸的形式定義行列式。對于1≤i，j≤n，A的第i個0子矩陣Ai,j是從矩陣A中刪除第i行和第j列后的結果。子矩陣Ai,j是（n-1）×（n-1）矩陣，因此我們可以計算它的行列式。然后，我們將以下量定義為A的行列式：

實際上，上述方程中出現的索引1可以替換為任意1≤i≤n，結果是相同的。

在NumPy中，用于計算行列式的例程位于名為linalg的獨立模塊中。該模塊包含許多線性代數（涵蓋向量和矩陣代數的數學分支）的常用函數。用于計算方陣行列式的具體例程為det：

請注意，在計算行列式時出現了浮點舍入誤差。

如果安裝了SciPy包，它也提供了一個linalg模塊，該模塊擴展了NumPy的linalg。SciPy版本的linalg模塊不僅包括額外的例程，而且它在編譯時始終支持基礎線性代數子程序包（BLAS）和線性代數包（LAPACK），而它們在NumPy版本中則是可選的。因此，如果更注重速度，根據NumPy的編譯方式，使用SciPy版本可能更好。

n×n矩陣A的逆矩陣是另一個（必然唯一的）n×n矩陣B，使得AB=BA=I，其中I表示n×n單位矩陣，這里的乘法是矩陣乘法。并不是每個方陣都有逆矩陣，那些沒有逆矩陣的矩陣有時被稱為奇異矩陣。事實上，當且僅當一個矩陣的行列式不為0時，該矩陣才是非奇異的（即具有逆矩陣）。當A有逆矩陣時，它的逆矩陣通常用A^-1表示。

用linalg模塊中的inv例程可以計算矩陣的逆（如果存在）：

將矩陣A乘以其逆矩陣（在任一側）得到的結果是2×2單位矩陣，依據這一事實我們可以驗證inv例程給出的矩陣確實是A的逆矩陣：

由于計算逆矩陣的方式不同，在這些計算中會出現浮點誤差，這些誤差已經隱藏在以“Approximately”（約等于）開頭的注釋后面。

linalg包還包含許多其他方法，如用于計算矩陣各種范數的norm方法。它還包括以各種方式分解矩陣和求解方程組的函數。

除此之外，它還包含用于矩陣運算的指數函數expm、對數函數logm、正弦函數sinm、余弦函數cosm和正切函數tanm。請注意，這些函數與NumPy基本包中的標準exp、log、sin、cos和tan函數不同，后者是按逐元素方式執行相應函數的。相反，矩陣指數函數是使用矩陣的冪級數定義的：

這是對任意n×n矩陣A定義的，其中A^k表示A的k次矩陣冪，即A矩陣連乘自己k次。請注意，這個“冪級數”總是在適當的意義上收斂。按照慣例，我們取A0=I，其中I是n×n單位矩陣。這完全類似于實數或復數的指數函數的常規冪級數定義，但用矩陣和矩陣乘法代替了數字和常規乘法。其他函數以類似的方式定義，但我們將略過這些詳細信息。

在1.5.4節中，我們將看到矩陣及其理論在求解方程組領域的應用。

1.5.4 方程組

求解線性方程組是數學中研究矩陣的主要動機之一，這類問題在各種應用中頻繁出現。我們從如下線性方程組開始：

這里，n至少為2，ai,j和bi是已知量，xi是我們希望求解的未知數。

在求解這樣的方程組之前，我們需要將問題轉化為矩陣方程。這是通過將方程中的系數ai,j收集到一個n×n矩陣中，并利用矩陣乘法的性質將這個矩陣與方程組關聯起來來實現的。因此，我們構造以下矩陣，矩陣元素為方程的系數：

然后，如果我們將x定義為包含xi值的未知（列）向量，將b定義為包含bi值的已知（列）向量，那么我們可以將方程組重寫為以下單矩陣方程：

Ax=b

我們現在可以使用矩陣技術來求解這個矩陣方程。在這種情況下，我們將列向量視為n×1矩陣，因此上述方程中的乘法是矩陣乘法。為了求解該矩陣方程，我們使用linalg模塊中的solve例程。為了說明該技術，我們將以求解以下方程組為例：

這個方程組有三個未知量：x1、x2和x3。首先，我們創建系數矩陣和向量b。由于我們使用NumPy來處理矩陣和向量，我們為矩陣A創建一個二維NumPy數組，并為b創建一個一維數組：

現在，可以使用linalg.solve例程求出方程組的解：

這確實是方程組的解，可以通過計算A@x并將結果與b數組進行比較來輕松驗證。在這個計算中可能存在浮點舍入誤差。

solve例程需要兩個輸入，即系數矩陣A和等式右側的向量b。它使用的例程將矩陣A分解為更簡單的矩陣，從而快速將問題簡化為可通過簡單替換求解的形式。這種求解矩陣方程的技術非常強大和高效，而且不太容易出現困擾其他方法的浮點舍入誤差。例如，如果矩陣的逆是已知的，那么可以通過左乘A的逆矩陣來計算方程組的解。然而，這通常不如使用solve例程效果好，因為它可能更慢或者會導致更大的數值誤差。

在我們使用的例子中，系數矩陣A是方陣。也就是說，方程的數量與未知數的數量相同。在這種情況下，當且僅當矩陣A的行列式不為0時，方程組才有唯一解。當A的行列式為0時，可能發生兩種情況：一是方程組沒有解，這時我們說方程組是不相容的；二是方程組可能有無窮多個解。方程組是相容還是不相容通常由向量b決定。例如，考慮以下方程組：

左邊的方程組是相容的，并且有無窮多個解：例如，取x=1和y=1，或者取x=0和y=2都是解。右邊的方程組是不相容的，沒有解。在這兩個方程組中，solve例程都會失敗，因為系數矩陣是奇異的。

求解方程組時，不一定要求系數矩陣為方陣——例如，如果方程的數量多于未知數的數量（系數矩陣的行數多于列數）。這樣的方程組稱為“超定方程組”，只要它是相容的，就會有解。如果方程的數量少于未知數的數量，則稱該方程組為“欠定方程組”。欠定方程組沒有足夠的信息來唯一確定所有的未知數，因此如果它是相容的，它通常有無窮多個解。遺憾的是，對于系數矩陣不是方陣的方程組，即使方程組確實有解，solve例程也無法找到解。

在1.5.5節中，我們將討論特征值和特征向量，與你在前面看到的類似，它們是通過觀察一種非常特殊的矩陣方程產生的。

1.5.5 特征值和特征向量

考慮矩陣方程Ax=λx，其中A是一個n×n方陣，x是一個向量，λ是一個系數。若存在x是此方程的解，則稱λ為A的特征值，對應的向量x稱為特征向量。特征值和對應的特征向量對矩陣A的信息進行編碼，因此在許多涉及矩陣的應用中非常重要。

我們使用以下矩陣演示如何計算特征值和特征向量：

我們首先需要將其定義為NumPy數組：

linalg模塊中的eig例程可以用來找到方陣的特征值和特征向量。該函數返回一個二元組（v,B），其中v是包含特征值的一維數組，B是一個二維數組，其列是相應的特征向量：

只有實數元素的矩陣完全可能具有復特征值和特征向量。因此，eig例程的返回類型有時會是復數類型，如complex32或complex64。在某些應用中，復特征值具有特殊含義，而在其他情況下，我們只考慮實特征值。

我們可以使用以下序列從eig例程的輸出結果中提取特征值和特征向量：

eig例程返回的特征向量已經被歸一化，使得它們的范數（長度）為1。（歐氏范數被定義為數組元素平方和的平方根。）我們可以通過使用linalg模塊中的norm例程計算特征向量的范數來驗證這一點：

最后，我們可以通過計算A @ x0，并檢查它是否在浮點精度范圍內等于lamb-da0*x0，來驗證這些值確實滿足特征值和特征向量的定義：

這里計算得到的范數表示方程Ax=λx的左側（lhs）和右側（rhs）之間的距離。由于這個距離非常?。ㄐ迭c后14位為零），我們可以相當確信它們實際上是相同的。事實上，這種距離不為零可能是由于浮點精度誤差引起的。

求特征值和特征向量的理論過程是先通過解以下方程來求特征值λ：

det (A-λI)=0

這里，I是適當的單位矩陣。左側確定的方程是關于λ的多項式，稱為矩陣A的特征多項式。然后，可以通過求解以下矩陣方程來求出與特征值λi相對應的特征向量：

(A-λiI)x=0

在實踐中，這個過程效率有些低，還有其他策略可以更高效地通過數值計算求得特征值和特征向量。

我們只能計算方陣的特征值和特征向量，對于不是方陣的矩陣，這個定義沒有意義。有一種稱為奇異值的概念可以將特征值和特征向量推廣到非方陣。為了做到這一點，我們必須做出的權衡是，我們必須計算兩個向量u和v以及奇異值σ，然后求解以下方程：

Au=σv

如果A是m×n矩陣，那么u將有n個元素，v將有m個元素。有趣的是，u向量實際上是對稱矩陣A^TA的（正交歸一化的）特征向量，其特征值為σ²。根據這些值，我們可以利用之前的定義方程求出v向量。這將產生所有有趣的組合，但還有其他的向量u和v，使得Au=0和A^Tv=0。

奇異值及對應向量的效用來自奇異值分解（SVD），它將矩陣A寫成以下乘積形式：

A=UΣV^T

這里，U為正交列，V為正交行，Σ為對角矩陣，通常的寫法是數值沿主對角線遞減。我們可以用稍微不同的方式寫出這個公式，如下所示：

這就是說，任何矩陣都可以分解成外積的加權和——假設u和v是n行1列的矩陣，然后將u與向量v的轉置進行矩陣乘法。

一旦完成了這種分解，我們就可以尋找特別小的σ值，這些值對矩陣貢獻很小。如果我們放棄具有小σ值的項，那么我們就可以用更簡單的表示法來有效地近似原始矩陣。這種技術在主成分分析（PCA）中得到了應用——例如，將復雜的高維數據集簡化為對數據整體特征貢獻最大的幾個分量。

在Python中，我們可以使用linalg.svd函數來計算矩陣的SVD。該函數的工作方式與之前描述的eig例程類似，只是它返回的是分解后的三個分量：

該函數返回的數組形狀分別為（2, 2）、（2,）和（4, 4）。顧名思義，U矩陣和VT矩陣是分解中出現的矩陣，而s是包含非零奇異值的一維向量。我們可以通過重構Σ矩陣并計算三個矩陣的乘積來檢查分解是否正確：

請注意，除了矩陣的第一個元素外，矩陣幾乎完全被重構了。左上角元素的值非常接近零——在浮點誤差范圍內可以被視為零。

我們構建矩陣Σ的方法相當不方便。SciPy版本的linalg模塊包含一個特殊例程linalg.diagsvd，用于從一維奇異值數組重構該矩陣。該函數獲取奇異值數組s和原始矩陣的形狀，并構建具有適當形狀的矩陣Σ：

（回想一下，SciPy軟件包是以別名sp導入的。）現在，讓我們改變一下方向，看看如何更有效地描述大多數元素都為零的矩陣，這就是所謂的稀疏矩陣。

1.5.6 稀疏矩陣

前面所討論的線性方程組，在整個數學領域，特別是在數學計算中是非常常見的。在許多應用中，系數矩陣可能會非常龐大，具有成千上萬的行和列，并且很可能來自其他來源而不是簡單地手動輸入。在許多情況下，這也可能是一個稀疏矩陣，其中大多數元素為0。

如果一個矩陣有大量元素為零，則這個矩陣稱為稀疏矩陣。矩陣有多少個元素為零才能稱為稀疏矩陣，這并沒有確切的定義。稀疏矩陣可以通過更高效的表示方式來存儲，例如，只存儲非零元素的索引（i,j）和相應的值ai,j。對于稀疏矩陣，有一整套專門的算法集合，可以在矩陣足夠稀疏的情況下顯著提高性能。

稀疏矩陣出現在許多應用中，并且通常遵循某種模式。特別是，求解偏微分方程（PDE）的幾種技術都涉及求解稀疏矩陣方程（參見第3章），與網絡相關的矩陣通常也是稀疏的。在sparse.csgraph模塊中包含了與網絡（圖）相關的稀疏矩陣的附加例程。我們將在第5章中進一步討論這些內容。

sparse模塊包含幾個不同的類，表示存儲稀疏矩陣的不同方法。存儲稀疏矩陣最基本的方法是存儲三個數組，其中兩個數組包含表示非零元素索引的整數，第三個數組包含相應的元素數據，這是coo_matrix類的格式。此外，還有壓縮稀疏列格式（CSC，使用csc_matrix實現）和壓縮稀疏行格式（CSR，使用csr_matrix實現），分別提供高效的列切片或行切片。在sparse模塊中還有三個額外的稀疏矩陣類，包括dia_matrix類，它可以高效地存儲非零元素沿對角線帶出現的矩陣。

SciPy的sparse模塊包含用于創建和處理稀疏矩陣的例程。我們使用以下導入語句從SciPy導入sparse模塊：

稀疏矩陣可以由滿矩陣（稠密矩陣）或其他類型的數據結構創建。這是使用特定格式（你希望使用這樣的格式存儲稀疏矩陣）的構造函數來完成的。

例如，我們可以使用以下命令將稠密矩陣存儲為CSR格式：

如果手動生成稀疏矩陣，則該矩陣可能遵循某種模式，例如以下三對角矩陣：

這里，非零元素出現在對角線上以及對角線的兩側，并且每行的非零元素遵循相同的模式。要創建這樣的矩陣，我們可以使用sparse中的一個數組創建例程，例如diags，這是一個創建具有對角線模式的矩陣的便捷例程：

這將創建一個如前所述的矩陣T，并將其存儲為CSR格式的稀疏矩陣。第一個參數指定應出現在輸出矩陣中的值，第二個參數是放置數值時相對于對角線的位置。因此，元組中的索引0表示對角線元素，-1表示在行中的對角線左側，+1表示在行中的對角線右側。shape關鍵字參數給出了所生成矩陣的維度，format指定了矩陣的存儲格式。如果沒有使用可選參數來設置格式，則將使用合理的默認值?？梢允褂胻oarray方法將數組T擴展為滿矩陣（稠密矩陣）：

當矩陣很小時（如這里所示），稀疏矩陣求解函數與常規的求解函數之間的性能幾乎沒有差異。

一旦將矩陣以稀疏格式存儲，我們就可以使用linalg中的子模塊sparse求解例程。例如，我們可以使用此模塊中的spsolve例程來求解矩陣方程。如果矩陣不是以稀疏格式提供的，spsolve例程會將矩陣轉換為CSR或CSC格式，這可能會增加計算時間：

sparse.linalg模塊還包含許多可以在NumPy（或SciPy）的linalg模塊中找到的函數，這些函數接受稀疏矩陣而不是完整的NumPy數組，例如eig和inv。

至此，我們對Python及其生態系統中可用的基本數學工具的介紹就結束了。讓我們總結一下我們所學到的內容。