numpy教程06---ndarray的進階操作

阿新 • • 發佈：2022-04-06

歡迎關注公眾號【Python開發實戰】, 獲取更多內容!

工具-numpy

numpy是使用Python進行資料科學的基礎庫。numpy以一個強大的N維陣列物件為中心，它還包含有用的線性代數，傅立葉變換和隨機數函式。

線性代數

numpy中二維的ndarray可以在Python中高效地表示矩陣，下面將介紹一些主要的矩陣運算。

匯入numpy

import numpy as np

矩陣轉置

當秩大於等於2時，T屬性相當於呼叫transpose()函式。

m1 = np.arange(10).reshape(2, 5)
m1

輸出：

array([[0, 1, 2, 3, 4],
       [5, 6, 7, 8, 9]])

m1.T

輸出：

array([[0, 5],
       [1, 6],
       [2, 7],
       [3, 8],
       [4, 9]])

m1.transpose()

輸出：

array([[0, 5],
       [1, 6],
       [2, 7],
       [3, 8],
       [4, 9]])

T屬性對秩為0（空）或秩為1的ndarray沒有影響。

m2 = np.arange(5)
m2

輸出：

array([0, 1, 2, 3, 4])

m2.T

輸出：

array([0, 1, 2, 3, 4])

可以將一維的ndarray重塑為單行的二維矩陣，進而得到轉置。

m2r = m2.reshape(1, 5)
m2r

輸出：

array([[0, 1, 2, 3, 4]])

m2r.T

輸出：

array([[0],
       [1],
       [2],
       [3],
       [4]])

矩陣乘法

dot()方法可以計算兩個矩陣的乘法，矩陣乘法需要滿足左邊矩陣的列數必須等於右邊矩陣的行數。

n1 = np.arange(10).reshape(2, 5)
n1

輸出：

array([[0, 1, 2, 3, 4],
       [5, 6, 7, 8, 9]])

n2 = np.arange(15).reshape(5, 3)
n2

輸出：

array([[ 0,  1,  2],
       [ 3,  4,  5],
       [ 6,  7,  8],
       [ 9, 10, 11],
       [12, 13, 14]])

n1.dot(n2)

輸出：

array([[ 90, 100, 110],
       [240, 275, 310]])

注意：n1*n2不是矩陣乘法，而是按元素乘積。

矩陣的逆與偽逆

許多線性代數函式在numpy中都可用。linalg模組中的inv函式可以計算一個方陣的逆。

import numpy.linalg as linalg

m3 = np.array([[1, 2, 3], [5, 7, 11], [21, 29, 31]])
m3

輸出：

array([[ 1,  2,  3],
       [ 5,  7, 11],
       [21, 29, 31]])

linalg.inv(m3)

輸出：

array([[-2.31818182,  0.56818182,  0.02272727],
       [ 1.72727273, -0.72727273,  0.09090909],
       [-0.04545455,  0.29545455, -0.06818182]])

也可以通過pinv函式來計算偽逆。

linalg.pinv(m3)

輸出：

array([[-2.31818182,  0.56818182,  0.02272727],
       [ 1.72727273, -0.72727273,  0.09090909],
       [-0.04545455,  0.29545455, -0.06818182]])

單位矩陣

矩陣與其逆矩陣相乘返回一個單位矩陣，下面的例子會有很小的浮點誤差。

m3.dot(linalg.inv(m3))

輸出：

array([[ 1.00000000e+00, -5.55111512e-17,  0.00000000e+00],
       [-2.98372438e-16,  1.00000000e+00, -5.55111512e-17],
       [ 5.78009862e-15,  1.27675648e-15,  1.00000000e+00]])

可以通過eye函式來建立一個N×N大小的單位矩陣。

np.eye(3)

輸出：

array([[1., 0., 0.],
       [0., 1., 0.],
       [0., 0., 1.]])

QR分解

linalg模組中的qr函式可以計算一個矩陣的QR分解（正交三角分解）。

q, r = linalg.qr(m3)
q

輸出：

array([[-0.04627448,  0.98786672,  0.14824986],
       [-0.23137241,  0.13377362, -0.96362411],
       [-0.97176411, -0.07889213,  0.22237479]])

輸出：

array([[-21.61018278, -29.89331494, -32.80860727],
       [  0.        ,   0.62427688,   1.9894538 ],
       [  0.        ,   0.        ,  -3.26149699]])

q.dot(r)

輸出：

array([[ 1.,  2.,  3.],
       [ 5.,  7., 11.],
       [21., 29., 31.]])

矩陣的行列式

linalg模組中的det函式可以計算矩陣的行列式。

linalg.det(m3)

輸出：

43.99999999999999

特徵值和特徵向量

linalg模組中的eig函式可以計算一個方陣的特徵值和特徵向量。

eigenvalues, eigenvetors = linalg.eig(m3)
eigenvalues   # λ

輸出：

array([42.26600592, -0.35798416, -2.90802176])

eigenvetors # v

輸出：

array([[-0.08381182, -0.76283526, -0.18913107],
       [-0.3075286 ,  0.64133975, -0.6853186 ],
       [-0.94784057, -0.08225377,  0.70325518]])

m3.dot(eigenvetors) - eigenvalues * eigenvetors     # m3.v -λ*v= 0

輸出：

array([[ 9.76996262e-15,  2.22044605e-16, -3.10862447e-15],
       [ 7.10542736e-15,  2.02615702e-15, -1.11022302e-15],
       [ 2.84217094e-14,  5.11049536e-15, -4.88498131e-15]])

奇異值分解

linalg模組中的svd函式可以計算矩陣的奇異值分解。

m4 = np.array([[1, 0, 0, 0, 2], [0, 0, 3, 0, 0], [0, 0, 0, 0, 0], [0, 2, 0, 0, 0]])
m4

輸出：

array([[1, 0, 0, 0, 2],
       [0, 0, 3, 0, 0],
       [0, 0, 0, 0, 0],
       [0, 2, 0, 0, 0]])

U, S_diag, V = linalg.svd(m4)
U

輸出：

array([[ 0.,  1.,  0.,  0.],
       [ 1.,  0.,  0.,  0.],
       [ 0.,  0.,  0., -1.],
       [ 0.,  0.,  1.,  0.]])

S_diag   # Σ對角線上的值

輸出：

array([3.        , 2.23606798, 2.        , 0.        ])

輸出：

array([[-0.        ,  0.        ,  1.        ,  0.        ,  0.        ],
       [ 0.4472136 ,  0.        ,  0.        ,  0.        ,  0.89442719],
       [-0.        ,  1.        ,  0.        ,  0.        ,  0.        ],
       [ 0.        ,  0.        ,  0.        ,  1.        ,  0.        ],
       [-0.89442719,  0.        ,  0.        ,  0.        ,  0.4472136 ]])

svd函式只返回Σ對角線上的值，完整的矩陣可以這樣建立：

S = np.zeros((4, 5))
S[np.diag_indices(4)] = S_diag   # np.diag_indices函式返回索引以訪問（4，4）陣列的主對角線。
S  # Σ

輸出：

array([[3.        , 0.        , 0.        , 0.        , 0.        ],
       [0.        , 2.23606798, 0.        , 0.        , 0.        ],
       [0.        , 0.        , 2.        , 0.        , 0.        ],
       [0.        , 0.        , 0.        , 0.        , 0.        ]])

U.dot(S).dot(V)  # U.Σ.V = m4

輸出：

array([[1., 0., 0., 0., 2.],
       [0., 0., 3., 0., 0.],
       [0., 0., 0., 0., 0.],
       [0., 2., 0., 0., 0.]])

對角線和軌跡

m3

輸出：

array([[ 1,  2,  3],
       [ 5,  7, 11],
       [21, 29, 31]])

np.diag(m3)  # 返回m3對角線上的元素值

輸出：

array([ 1,  7, 31])

np.trace(m3)   # 相當於 np.diag(m3).sum()

輸出：

np.diag(m3).sum()

輸出：

求解線性標量方程組

linalg模組中的solve函式可以求解線性標量方程組, 例如如下方程組

2x + 6y = 6
5x + 3y = -9

coeffs = np.array([[2, 6], [5, 3]])
depvars = np.array([6, -9])
solution = linalg.solve(coeffs, depvars)
solution

輸出：

array([-3.,  2.])

可以驗證一下求解：

coeffs.dot(solution), depvars

輸出：

(array([ 6., -9.]), array([ 6, -9]))

還可以通過另一個方式驗證求解。

np.allclose(coeffs.dot(solution), depvars)   # np.allclose比較兩個ndarray的每一個元素是否都相等

輸出：

True

向量化

如果堅持進行ndarray操作，而不是一次一個地對單個的元素極性操作，那麼程式碼的效率要高的多，這就叫做向量化。這樣，可以從numpy的許多優化中受益。

例如，要根據公式sin(xy/40.5)生成一個768×1024的ndarray，一個不好的選擇就是使用巢狀迴圈進行計算。

import math
data = np.empty((768, 1024))
for y in range(768):
    for x in range(1024):
        data[y, x] = math.sin(x*y/40.5)
data

輸出：

array([[0.        , 0.        , 0.        , ..., 0.        , 0.        ,
        0.        ],
       [0.        , 0.02468885, 0.04936265, ..., 0.07705885, 0.1016508 ,
        0.12618078],
       [0.        , 0.04936265, 0.09860494, ..., 0.15365943, 0.20224852,
        0.25034449],
       ...,
       [0.        , 0.03932283, 0.07858482, ..., 0.6301488 , 0.59912825,
        0.56718092],
       [0.        , 0.06398059, 0.12769901, ..., 0.56844086, 0.51463783,
        0.45872596],
       [0.        , 0.08859936, 0.17650185, ..., 0.50335246, 0.42481591,
        0.34293805]])

上面的方法雖然可行，但是效率非常低，因為迴圈是在純Python中進行的。現在我們把這個方法向量化，首先，使用numpy的meshgrid函式，該函式可以根據座標向量生成座標矩陣。

x_coords = np.arange(1024)
y_coords = np.arange(768)
X, Y = np.meshgrid(x_coords, y_coords)
X

輸出：

array([[   0,    1,    2, ..., 1021, 1022, 1023],
       [   0,    1,    2, ..., 1021, 1022, 1023],
       [   0,    1,    2, ..., 1021, 1022, 1023],
       ...,
       [   0,    1,    2, ..., 1021, 1022, 1023],
       [   0,    1,    2, ..., 1021, 1022, 1023],
       [   0,    1,    2, ..., 1021, 1022, 1023]])

輸出：

array([[  0,   0,   0, ...,   0,   0,   0],
       [  1,   1,   1, ...,   1,   1,   1],
       [  2,   2,   2, ...,   2,   2,   2],
       ...,
       [765, 765, 765, ..., 765, 765, 765],
       [766, 766, 766, ..., 766, 766, 766],
       [767, 767, 767, ..., 767, 767, 767]])

X.shape, Y.shape

輸出：

((768, 1024), (768, 1024))

X和Y都是768×1024的ndarray，X中所有的值對應水平軸的座標， Y中所有的值對應垂直軸的座標。現在可以簡單地使用ndarray運算計算結果。

data = np.sin(X*Y/40.5)
data

輸出：

array([[0.        , 0.        , 0.        , ..., 0.        , 0.        ,
        0.        ],
       [0.        , 0.02468885, 0.04936265, ..., 0.07705885, 0.1016508 ,
        0.12618078],
       [0.        , 0.04936265, 0.09860494, ..., 0.15365943, 0.20224852,
        0.25034449],
       ...,
       [0.        , 0.03932283, 0.07858482, ..., 0.6301488 , 0.59912825,
        0.56718092],
       [0.        , 0.06398059, 0.12769901, ..., 0.56844086, 0.51463783,
        0.45872596],
       [0.        , 0.08859936, 0.17650185, ..., 0.50335246, 0.42481591,
        0.34293805]])

儲存和載入

numpy可以方便地以二進位制或文字格式來儲存和載入ndarray。

二進位制格式 .npy

下面建立一個隨機的ndarray，並儲存。

a = np.random.rand(2, 3)
a

輸出：

array([[0.76953407, 0.79648012, 0.8868019 ],
       [0.88131806, 0.57297333, 0.70059907]])

np.save('my_array', a)

由於檔名不包含副檔名，因此numpy會自動新增副檔名為.npy，下面檢視一下檔案內容：

with open('my_array.npy', 'rb') as f:
    content  = f.read()
content

輸出：

b"\x93NUMPY\x01\x00v\x00{'descr': '<f8', 'fortran_order': False, 'shape': (2, 3), }                                                          \n0L\xbb\xe8\x05\xa0\xe8?\x13s\x00\xdf\xc3|\xe9?\xdd\x05\xf4`\xae`\xec?S\x0e\xad\xee\xc13\xec?\x03\xa6E)\xccU\xe2?\x97\xc6\xdc\xbcNk\xe6?"

想要把該檔案載入到numpy陣列中只需要呼叫load函式。

a_loaded = np.load('my_array.npy')
a_loaded

輸出：

array([[0.76953407, 0.79648012, 0.8868019 ],
       [0.88131806, 0.57297333, 0.70059907]])

文字格式

現在以文字格式來儲存ndarray。

np.savetxt('my_array.csv', a)

現在檢視一下檔案內容：

with open('my_array.csv', 'rt') as f:
    print(f.read())

輸出：

7.695340676822350900e-01 7.964801173689884939e-01 8.868019002549619723e-01
8.813180600777265061e-01 5.729733282179839682e-01 7.005990685194828371e-01

這是以製表符作為分隔符的csv檔案，還可以設定不同的分隔符。

np.savetxt('my_array.csv', a, delimiter=',')

再檢視一下檔案內容

with open('my_array.csv', 'rt') as f:
    print(f.read())

輸出：

7.695340676822350900e-01,7.964801173689884939e-01,8.868019002549619723e-01
8.813180600777265061e-01,5.729733282179839682e-01,7.005990685194828371e-01

載入該檔案，只需要呼叫loadtxt函式。

a_loaded = np.loadtxt('my_array.csv', delimiter=',')
a_loaded

輸出：

array([[0.76953407, 0.79648012, 0.8868019 ],
       [0.88131806, 0.57297333, 0.70059907]])

壓縮格式 .npz

還可以在一個壓縮檔案中儲存多個ndarray。

b = np.arange(24).reshape(2, 3, 4)
b

輸出：

array([[[ 0,  1,  2,  3],
        [ 4,  5,  6,  7],
        [ 8,  9, 10, 11]],

       [[12, 13, 14, 15],
        [16, 17, 18, 19],
        [20, 21, 22, 23]]])

np.savez('my_arrays', my_a=a, my_b=b)

現在看一下檔案內容，檔案的副檔名.npz已經自動新增。

with open('my_arrays.npz', 'rb') as f:
    content = f.read()
repr(content)[:180] + '[...]'

輸出：

'b"PK\\x03\\x04\\x14\\x00\\x00\\x00\\x00\\x00\\x00\\x00!\\x00&\\xa9j\\xe4\\xb0\\x00\\x00\\x00\\xb0\\x00\\x00\\x00\\x08\\x00\\x00\\x00my_a.npy\\x93NUMPY\\x01\\x00v\\x00{\'descr\': \'<f8\', \'fortran_order\': False, \'s[...]'

然後可以這樣載入這個檔案。

my_arrays = np.load('my_arrays.npz')
my_arrays

輸出：

<numpy.lib.npyio.NpzFile at 0x93b8080>

這是一個類似字典的物件，它會懶載入ndarray。

my_arrays.keys()

輸出：

['my_a', 'my_b']

my_arrays['my_a']

輸出：

array([[0.76953407, 0.79648012, 0.8868019 ],
       [0.88131806, 0.57297333, 0.70059907]])

my_arrays['my_b']

輸出：

array([[[ 0,  1,  2,  3],
        [ 4,  5,  6,  7],
        [ 8,  9, 10, 11]],

       [[12, 13, 14, 15],
        [16, 17, 18, 19],
        [20, 21, 22, 23]]])