Python之Numpy：線性代數/矩陣運算

當你知道工具的用處，理論與工具如何結合的時候，一般會加速我們對二者的學習效率。

零 numpy

那麼，Numpy是什麼？
NumPy(Numerical Python) 是 Python 語言的一個擴展程序庫，支持大量維度的數組與矩陣運算，此外也針對數組運算提供大量的數學函數庫。
NumPy 的前身 Numeric 最先是由 Jim Hugunin 與其它協做者共同開發，2005 年，Travis Oliphant 在 Numeric 中結合了另外一個同性質的程序庫 Numarray 的特點，並加入了其它擴展而開發了 NumPy。NumPy 爲開放源代碼而且由許多協做者共同維護開發。 ---- 摘自 · 菜鳥教程

一 要點

• 假定AX=b，求解未知矩陣X 【線性代數中常遇到的運算問題】
• 矩陣轉置A^(T)
• 矩陣的逆A^(-1)
• 矩陣行列式的值|A|
• 矩陣的秩 rank(A)
• 矩陣的跡 trace(A)
• 其它
  • 單位矩陣
  • 0向量/矩陣
  • ...

二 示例應用

2.1 求解AX=b中的未知參數矩陣X

import numpy as np

# Hypothsis : A*X = b

A = [[2,1,2],
     [3,1,0],
     [1,1,-1]];
b = np.transpose([-3,5,-2])# 轉置
#[or] b = np.transpose(np.array([-3,5,-2]))# 轉置

# 求解未知參數矩陣X
X = np.linalg.solve(A,b) # 方式一：直接使用numpy的solve函數一鍵求解
#A_inv=np.linalg.inv(A) # 方式二：先求逆運算，再點積求值
#X=np.dot(A_inv,b) # a.dot(b) 與 np.dot(a,b) 效果相同；but np.dot(a,b)與np.dot(b,a)效果確定是不一樣的(線性代數/矩陣常識)
print("方程組的解：\n",X);

# [output]
方程組的解：
 [ 4.4 -8.2 -1.8]

2.2 利用最小二乘法擬合函數模型

給出一組數據【5對（Xi,Yi）參數】，用最小二乘法，求形如：f(x)=a+b*x^3的經驗公式。

• 原方程(模型)
  \[ f(x) = a + bx^3 \]
• 其法方程
  \[ A^T A \begin{pmatrix} a \\ b \\ \end{pmatrix} = A^T y \]
  即
  \[ \begin{pmatrix} a \\ b \\ \end{pmatrix} = \left( A^T A \right)^{-1} A^T y \]

import numpy as np
import matplotlib.pyplot as plt  # Python 繪圖工具(業界推薦)

# 數據初始化
A = [
    [1,pow(-3,3)],
    [1,pow(-2,3)],
    [1,pow(-1,3)],
    [1,pow(2,3)],
    [1,pow(4,3)]
];

At =  np.transpose(A); # A的轉置矩陣

y = np.transpose([14.3,8.3,4.7,8.3,22.7]);


# step1:求解
## 令 (a ,b)^T 爲 未知參數X
X = np.dot(np.dot(np.linalg.inv(np.dot(At,A)),At),y)
print(X)
print("a:",X[0])
print("b:",X[1])

# step2:查看擬合效果
x = [-3,-2,-1,2,4];
## 1. 計算擬合數值 fitValue
def fitValue(arg_x):
    a = X[0];
    b = X[1];
    return a + b*pow(arg_x,3);

fitValues = np.zeros([5]); # 建立長爲5的【一維】數組；[1,5]：建立第1行爲5個元素的【二維】數組 
for i in range(0,len(fitValues)):
    fitValues[i] = fitValue(x[i]);
    print(i,":","x(i):",x[i],"fit Y:",fitValues[i]); # just for test
    pass;

## 2. 繪圖可視化
yt = np.transpose(y); # y的轉置
plt.rcParams['figure.dpi'] = 100 #分辨率
plt.scatter(x, yt, marker = '*',color = 'red', s = 10 ,label = 'Actual Dataset') # 真實數據集
plt.scatter(x, fitValues, marker = 'x',color = 'green', s = 10 ,label = 'Fitting Dataset') #[擬合數據集]
plt.legend(loc = 'best')    # 設置 圖例所在的位置 使用推薦位置
plt.show() 


# [output]
[ 10.67505325   0.13679816]
a: 10.6750532504
b: 0.136798159666
0 : x(i): -3 fit Y: 6.98150293942
1 : x(i): -2 fit Y: 9.58066797308
2 : x(i): -1 fit Y: 10.5382550907
3 : x(i): 2 fit Y: 11.7694385277
4 : x(i): 4 fit Y: 19.430135469

2.3 數組建立/初始化

• numpy.linspace(start, stop[, num=50[, endpoint=True[, retstep=False[, dtype=None]]]]])
  • 一維等差數列
  • return 在指定範圍內的均勻間隔的數字（組成的數組），也即返回一個等差數列
  • 參數
    • start - 起始點
    • stop - 結束點
    • num - 元素個數，默認爲50，
    • endpoint - 是否包含stop數值，默認爲True，包含stop值；若爲False，則不包含stop值
    • retstep - 返回值形式，默認爲False，返回等差數列組，若爲True，則返回結果(array([samples, step])),
    • dtype - 返回結果的數據類型，默認無，若無，則參考輸入數據類型。

import numpy as np

a = np.linspace(1,10,5,endpoint= True)
print(a) # [ 1.    3.25  5.5   7.75 10.  ]

b = np.linspace(1,10,5,endpoint= False)
print(b) #[1.  2.8 4.6 6.4 8.2]

c = np.linspace(1,10,5,retstep = False)
print(c) # [ 1.    3.25  5.5   7.75 10.  ]

d = np.linspace(1,10,5,retstep = True)
print(d) # (array([ 1.  ,  3.25,  5.5 ,  7.75, 10.  ]), 2.25)

2.4 線性代數經常使用運算

print("原矩陣A：\n",A);
print("原矩陣b：\n",b);

print("轉置矩陣A^T：\n",np.transpose(A)); # 轉置
print("矩陣的行列式值|A|：\n",np.linalg.det(A)); # 方陣的行列式值：|A|
print("矩陣的跡trace(A)：\n",np.trace(A)); 
print("矩陣的秩rank(A)：\n",np.linalg.matrix_rank(A)); 
print("逆矩陣A^(-1)：\n",np.linalg.inv(A)); #矩陣的逆運算(條件：矩陣A可逆(行列式值不爲0)| 矩陣A爲方陣)

print("*"*30); # 分隔線

print("N階單位矩陣:\n",np.eye(4));
print(np.zeros([5])); # 建立長爲5的【一維】數組；[1,5]：建立第1行爲5個元素的【二維】數組 

# 建立指定的初始化數組
print(np.array([1]))  # 生成 第1行含值爲1的元素的【一維】數組
print(np.array([[56]])) # 生成 第1行含值爲56的元素的【二維】數組
np.full((3,5),3.14) # 建立一個3x5的浮點型數組，數組的值都是3.14

# [output]
原矩陣A：
 [[2, 1, 2], [3, 1, 0], [1, 1, -1]]
原矩陣b：
 [-3  5 -2]
轉置矩陣A^T：
 [[ 2  3  1]
 [ 1  1  1]
 [ 2  0 -1]]
矩陣的行列式值|A|：
 5.0
矩陣的跡trace(A)：
 2
矩陣的秩rank(A)：
 3
逆矩陣A^(-1)：
 [[-0.2  0.6 -0.4]
 [ 0.6 -0.8  1.2]
 [ 0.4 -0.2 -0.2]]
**********************
N階單位矩陣:
 [[ 1.  0.  0.  0.]
 [ 0.  1.  0.  0.]
 [ 0.  0.  1.  0.]
 [ 0.  0.  0.  1.]]
[ 0.  0.  0.  0.  0.]
[1]
[[56]]

三 推薦文獻

• Numpy - 菜鳥教程
• Numpy學習筆記二——初始化數組的10種方法