<style>
.markdown-body table{
display: unset;
}
</style>
# 使用 Python 及 C++ 計算矩陣乘法、反矩陣
> 作者：王一哲
> 日期：2024年7月7日
<br />
## 前言
矩陣乘法及反矩陣的運算會出現在高中數學教材之中，同時也可以運用在物理課程之中，例如克希荷夫定律求分支電流。我之前都是直接用 NumPy 處理，為了更了解運算的過程，我再花了一些時間研究不使用 NumPy 的寫法。
<br /><br />
## 矩陣乘法
### 運算原則
矩陣橫的部分為**列 (row)**、直的部分為**欄 (column)**，如果某個矩陣共有 $m$ 列、$n$ 欄，則這個矩陣的**維度 (dimension)** 為 $m \times n$。假設 $A$ 矩陣的維度為 $m \times n$，$B$ 矩陣的維度為 $n \times p$，則兩個矩陣相乘 $AB$ 的維度為 $m \times p$。如果兩個矩陣要相乘，第一個矩陣的欄與第二個矩陣的列數量必須相等。需要特別注意，矩陣乘法沒有交換性。
假設矩陣
$$
A =
\begin{bmatrix}
a_{0, 0} & a_{0, 1} & a_{0, 2} & \dots \\
a_{1, 0} & a_{1, 1} & a_{1, 2} & \dots \\
\vdots & \vdots & \vdots & \ddots
\end{bmatrix} =
\begin{bmatrix}
A_1 & A_2 & \dots
\end{bmatrix}
$$
$$
B =
\begin{bmatrix}
b_{0, 0} & b_{0, 1} & b_{0, 2} & \dots \\
b_{1, 0} & b_{1, 1} & b_{1, 2} & \dots \\
\vdots & \vdots & \vdots & \ddots
\end{bmatrix} =
\begin{bmatrix}
B_1 \\ B_2 \\ \vdots
\end{bmatrix}
$$
相乘的結果
$$
AB =
\begin{bmatrix}
a_{0, 0} b_{0, 0} + a_{0, 1} b_{1, 0} + \dots & a_{0, 0} b_{0, 1} + a_{0, 1} b_{1, 1} + \dots & \dots \\
a_{1, 0} b_{0, 0} + a_{1, 1} b_{1, 0} + \dots & a_{1, 0} b_{1, 1} + a_{0, 1} b_{1, 1} + \dots & \dots \\
\vdots & \vdots & \ddots
\end{bmatrix} = A_1 B_1 + A_2 B_2 + \dots
$$
<br /><br />
### 使用 Python 的 NumPy 套件計算矩陣乘法
先用 numpy.asarray 將二維串列轉成 ndarray，再用內建的矩陣乘法即可，如果這兩個矩陣維度不合、無法相乘，執行時會回傳 ValueError。
```python=
import numpy as np
A = np.asarray([[0, 1, 2, 3],
[4, 5, 6, 7]])
B = np.asarray([[0, 1, 2],
[3, 4, 5],
[-5, -4, -3],
[-2, -1, 0]])
print(A @ B)
print(np.dot(A, B))
print(A.dot(B))
print(np.matmul(A, B))
```
<br />
以上四種寫法效果相同，輸出都是
```python
[[-13 -7 -1]
[-29 -7 15]]
```
<br /><br />
### 使用 Python 但不使用 NumPy 套件計算矩陣乘法
我習慣上會先建立儲存相乘結果的串列 C，並將所有元素都預設為0。由於計算時會依序處理矩陣 C 的 $c_{0, 0}, c_{0, 1}, c_{0, 2}, \dots$，再處理 $c_{1, 0}, c_{1, 1}, c_{1, 2}, \dots$，因此第一層 for 迴圈的變數 i 依序為 $0, 1, 2, \dots, m-1$，第二層 for 迴圈的變數 j 依序為 $0, 1, 2, \dots, p-1$，第三層 for 迴圈的變數 k 依序為 $0, 1, 2, \dots, n-1$。
```python=
def matmul(A, B):
m, n, p = len(A), len(A[0]), len(B[0])
C = [[0]*p for _ in range(m)]
for i in range(m):
for j in range(p):
for k in range(n):
C[i][j] += A[i][k]*B[k][j]
return C
if __name__ == "__main__":
A = [[0, 1, 2, 3],
[4, 5, 6, 7]]
B = [[0, 1, 2],
[3, 4, 5],
[-5, -4, -3],
[-2, -1, 0]]
C = matmul(A, B)
for r in C: print(*r)
```
<br />
輸出為
```python
-13 -7 -1
-29 -7 15
```
<br /><br />
### 使用 C++ 計算矩陣乘法
程式運作的邏輯與上一個寫法一樣，以下的程式碼是用 vector 作為容器，也可以使用 array，但是 array 沒有 size 屬性，寫法會略有不同。
```cpp=
#include <iostream>
#include <vector>
using namespace std;
void matmul(vector<vector<int>> &A, vector<vector<int>> &B, vector<vector<int>> &C) {
size_t m = A.size(), n = A[0].size(), p = B[0].size();
for(size_t i=0; i<m; i++)
for(size_t j=0; j<p; j++)
for(size_t k=0; k<n; k++)
C[i][j] += A[i][k]*B[k][j];
}
int main() {
vector<vector<int>> A = {{0, 1, 2, 3},
{4, 5, 6, 7}};
vector<vector<int>> B = {{0, 1, 2},
{3, 4, 5},
{-5, -4, -3},
{-2, -1, 0}};
vector<vector<int>> C (A.size(), vector<int> (B[0].size(), 0));
matmul(A, B, C);
for(size_t i=0; i<C.size(); i++)
for(size_t j=0; j<C[0].size(); j++)
cout << C[i][j] << " \n"[j == C[0].size()-1];
return 0;
}
```
<br />
輸出為
```cpp
-13 -7 -1
-29 -7 15
```
<br /><br />
### 執行時間
我在 Ubuntu 指令介面中使用 time 測試這幾個寫法的執行時間，很意外地發現用 Python list 的寫法執行速度居然有點快，C++ 的部分是測試編輯後的可執行檔，速度明顯快很多。
<div style="text-align:center">
<table>
| <center>程式碼</center> | real (s) | user (s) | sys (s) |
| ----- | -------- | -------- | ------- |
| A @ B | 0.196 | 0.548 | 1.804 |
| np.dot(A, B) | 0.122 | 0.471 | 1.164 |
| A.dot(B) | 0.150 | 0.417 | 1.216 |
| np.matmul(A, B) | 0.156 | 0.477 | 1.239 |
| Python list | 0.035 | 0.023 | 0.012 |
| C++ | 0.002 | 0.002 | 0.002 |
</table>
</div>
<br /><br />
## 反矩陣
### 定義
如果一個$n$階方陣（$n \times n$ 矩陣）$A$，與另一個$n$階方陣$B$，兩者相乘後為$n$階單位矩陣
$$
AB = BA = I_n
$$
則矩陣$B$為$A$的反矩陣，符號為$A^{-1}$。在高中數學課本中，通常會用高斯消去法求反矩陣，先將矩陣$A$與單位矩陣$I$並列，再透過整列同加或同減的運算，將左側的$A$化為單位矩陣，此時右側就是$A^{-1}$。詳細的計算過程可以參考[英文版維基百科的 Invertible matrix 條目](https://en.wikipedia.org/wiki/Invertible_matrix#Gaussian_elimination)。以下的程式碼是以這兩個矩陣為例。
$$
A =
\begin{bmatrix}
2 & 2 & 0 \\
2 & 0 & 2 \\
0 & 2 & 0
\end{bmatrix}~~~~~A^{-1} =
\begin{bmatrix}
\frac{1}{2} & 0 & -\frac{1}{2} \\
0 & 0 & \frac{1}{2} \\
-\frac{1}{2} & \frac{1}{2} & \frac{1}{2}
\end{bmatrix}
$$
<br /><br />
### 使用 Python 的 NumPy 套件計算反矩陣
直接用 numpy.linalg.inv 即可。
```python=
import numpy as np
A = np.asarray([[2, 2, 0],
[2, 0, 2],
[0, 2, 0]])
print(np.linalg.inv(A))
```
<br />
輸出為
```python
[[ 0.5 0. -0.5]
[-0. -0. 0.5]
[-0.5 0.5 0.5]]
```
<br/><br />
### 使用 Python 但不使用 NumPy 套件計算反矩陣
以下的程式碼使用高斯消去法求反矩陣，還有其它更有效率的作法，有興趣的同學可以參考[這個網頁](https://www.scicoding.com/how-to-calculate-matrix-inverse-in-python/)。
```python=
import copy
def inv(A):
# 取 A 的列、欄數量，如果不是方陣，回傳 ValueError 並中止程式
n, m = len(A), len(A[0])
if n != m: raise ValueError
# 建立右側的單位矩陣 B
B = [[0]*n for _ in range(n)]
for i in range(n): B[i][i] = 1.0
# 使用 copy.deepcopy 複製 A
test = copy.deepcopy(A)
# 消除 test 左下角的數值
for i in range(n): # 依序處理第 0 ~ (n-1) 列
coef = test[i][i] # 以對角線上的元素作為係數
for j in range(i+1, n): # 處理第 (i+1) ~ (n-1) 列
mul = -test[j][i] / coef # 計算要將 test[i] 同乘的倍數 mul
for k in range(n): # 處理第 0 ~ (n-1) 欄
test[j][k] += mul*test[i][k] # 化簡 test 的左下角
B[j][k] += mul*B[i][k] # 同時處理 B
# 消除 test 右下角的數值
for i in range(n-1, -1, -1): # 依序處理第 (n-1) ~ 0 列
coef = test[i][i] # 以對角線上的元素作為係數
for j in range(i-1, -1, -1): # 處理第 (i-1) ~ (0) 列
mul = -test[j][i] / coef # 計算要將 test[i] 同乘的倍數 mul
for k in range(n-1, -1, -1): # 處理第 (n-1) ~ 0 欄
test[j][k] += mul*test[i][k] # 化簡 test 的右下角
B[j][k] += mul*B[i][k] # 同時處理 B
# 將 test 對角線上的元素都變為1
for i in range(n): # 依序處理第 0 ~ (n-1) 列
coef = test[i][i] # 以對角線上的元素作為係數
test[i][i] /= coef # 處理對角線上的元素
for j in range(n): # 依序處理 B 的第 0 ~ (n-1) 欄
B[i][j] /= coef
# 回傳反矩陣 B
return B
if __name__ == "__main__":
A = [[2, 2, 0],
[2, 0, 2],
[0, 2, 0]]
print(inv(A))
```
<br />
輸出為
```python
[[0.5, 0.0, -0.5], [-0.0, -0.0, 0.5], [-0.5, 0.5, 0.5]]
```
<br /><br />
### 使用 C++ 計算矩陣乘法
一樣使用高斯消去法求反矩陣。
```cpp=
#include <iostream>
#include <vector>
using namespace std;
void inv(vector<vector<float>> &A, vector<vector<float>> &B) {
int n = (int)A.size(); // A 的列、欄數量
for(int i=0; i<n; i++) B[i][i] = 1.0; // 將 B 的對角線元素都設定為 1.0
vector<vector<float>> test (A.begin(), A.end()); // 複製 A 的資料到 test
/* 消除 test 左下角的數值 */
for(int i=0; i<n; i++) { // 依序處理第 0 ~ (n-1) 列
float coef = test[i][i]; // 以對角線上的元素作為係數
for(int j=i+1; j<n; j++) { // 處理第 (i+1) ~ (n-1) 列
float mul = -test[j][i] / coef; // 計算要將 test[i] 同乘的倍數 mul
for(int k=0; k<n; k++) { // 處理第 0 ~ (n-1) 欄
test[j][k] += mul*test[i][k]; // 化簡 test 的左下角
B[j][k] += mul*B[i][k]; // 同時處理 B
}
}
}
/* 消除 test 右上角的數值 */
for(int i=n-1; i>=0; i--) { // 依序處理第 (n-1) ~ 0 列
float coef = test[i][i]; // 以對角線上的元素作為係數
for(int j=i-1; j>=0; j--) { // 處理第 (i-1) ~ (0) 列
float mul = -test[j][i] / coef; // 計算要將 A[i] 同乘的倍數 mul
for(int k=n-1; k>=0; k--) { // 處理第 (n-1) ~ 0 欄
test[j][k] += mul*test[i][k]; // 化簡 test 的右下角
B[j][k] += mul*B[i][k]; // 同時處理 B
}
}
}
/* 將 test 對角線上的元素都變為1 */
for(int i=0; i<n; i++) { // 依序處理第 0 ~ (n-1) 列
float coef = test[i][i]; // 以對角線上的元素作為係數
test[i][i] /= coef; // 處理對角線上的元素
for(int j=0; j<n; j++) B[i][j] /= coef; // 依序處理 B 的第 0 ~ (n-1) 欄
}
}
int main() {
vector<vector<float>> A = {{2.0, 2.0, 0.0},
{2.0, 0.0, 2.0},
{0.0, 2.0, 0.0}};
vector<vector<float>> B (A.size(), vector<float> (A[0].size(), 0.0));
inv(A, B);
for(int i=0; i<(int)B.size(); i++)
for(int j=0; j<(int)B[0].size(); j++)
cout << B[i][j] << " \n"[j == (int)B[0].size()-1];
return 0;
}
```
<br />
輸出為
```cpp=
0.5 0 -0.5
-0 -0 0.5
-0.5 0.5 0.5
```
<br /><br />
### 執行時間
我在 Ubuntu 指令介面中使用 time 測試這幾個寫法的執行時間，很意外地發現用 Python list 的寫法執行速度居然有點快，C++ 的部分是測試編輯後的可執行檔，速度明顯快很多。
<div style="text-align:center">
<table>
| <center>程式碼</center> | real (s) | user (s) | sys (s) |
| ----- | -------- | -------- | ------- |
| np.linalg.inv(A) | 0.165 | 0.579 | 1.276 |
| Python list | 0.038 | 0.034 | 0.004 |
| C++ | 0.003 | 0.003 | 0.000 |
</table>
</div>
<br /><br />
## 結語
由於矩陣運算可以運用在解線性的聯立方程式，例如克希荷夫定律求分支電流，對物理、電機科系的學生很有用。實務上直接用 NumPy 會比較方便，不使用 NumPy 的寫法，是為了更了解運算的過程，有興趣的同學可以參考看看。
<br /><br />
## 參考資料
1. [WikiPedia Matrix multiplication](https://en.wikipedia.org/wiki/Matrix_multiplication)
2. [NumPy for MATLAB users](https://numpy.org/devdocs/user/numpy-for-matlab-users.html)
3. [numpy.dot](https://numpy.org/doc/stable/reference/generated/numpy.dot.html)
4. [numpy.matmul](https://numpy.org/doc/stable/reference/generated/numpy.matmul.html)
5. [WikiPeida Invertible matrix](https://en.wikipedia.org/wiki/Invertible_matrix#Gaussian_elimination)
6. [numpy.linalg.inv](https://numpy.org/doc/stable/reference/generated/numpy.linalg.inv.html)
7. [How to Calculate Matrix Inverse in Python](https://www.scicoding.com/how-to-calculate-matrix-inverse-in-python/)
---
###### tags:`Python`、`C++`、`Math`
Sign in with Wallet
Connect another wallet