Python: NumPy 與數值線性代數¶

Python: NumPy and numerical linear algebra

Creative Commons License
This work by Jephian Lin is licensed under a Creative Commons Attribution 4.0 International License.

In [ ]:
import numpy as np
from numpy import linalg as LA

符號運算與數值運算¶

符號運算 通常比較慢，但帶來精確的結果；
數值運算 通常比較快，但必須處理計算誤差。

在線性代數裡，
符號運算擅長的包含：最簡階梯型、零解空間、秩與零維度、喬丹標準型等等；
數值運算擅長的包含：矩陣指數、特徵值、特徵向量、奇異值分解、QR 分解等等。

Sage 為建立在 Python 上的一套代數系統，
主要目的在於處理符號運算；
而 NumPy 為純 Python 環境中的一個套件，
主要目的在於處理大型陣列
（把矩陣看成兩個維度的話，NumPy 可以處理任意維度的陣列）
的數值運算。

Symbolic computation is slower but leads to an exact solution; while
numerical computation is faster but leads to a solution with numerical errors.

In linear algebra, symbolic computation is good for finding the reduced echelon form, the kernel, the rank and the nullity, and the Jordan canonical form, etc.; while numerical computation is good for finding the matrix exponential, the eigenvalues, the eigenvectors, the singular value decomposition, and the QR decomposition, etc.

Sage is a algebra system built on Python, which is good at symbolic computation; while NumPy is a package in Python, which is good at numerical computation on huge arrays. (A matrix can be viewed as an array of two dimension, while NumPy can deal with arrays of arbitrary dimensions.

建構矩陣¶

  1. np.array( list of lists )：把 list of lists 中的每個 list 當作矩陣的列。
  2. np.array(list).reshape(m,n)：把 list 重組成 m x n 的矩陣。
  3. np.eye(n)：單位矩陣。
  4. np.zeros((m,n))：全零矩陣。
  5. np.diag(list)：對角矩陣，其對角線上元素由 list 決定。

利用 print 來顯示矩陣。

  1. np.array( list of lists ) : construct an array whose rows are the list in list of lists.
  2. np.array(list).reshape(m,n) : make list into an m x n matrix.
  3. np.eye(n) : identity matrix.
  4. np.zeros((m,n)) : zero matrix.
  5. np.diag(list) : diagonal matrix whose entries are determined by list .

Use print to show the matrix.

In [ ]:
A = np.array([[1,2,3], 
              [4,5,6]])
print(A)
In [ ]:
A = np.array([1,2,3,4,5,6]).reshape(2,3)
print(A)
In [ ]:
A = np.eye(3)
# A = np.zeros((3,3))
# A = np.diag([1,2,3])
print(A)

陣列之間的運算¶

Operations between matrices

有別於 Sage，在 NumPy 中的所有常見運算都是逐項處理。
若要處理矩陣相乘，可以使用 A.dot(B)。

Different from Sage, most of the operations in NumPy is entrywise. We may use A.dot(B) for the usual matrix multiplication.

In [ ]:
A = np.array([1,2,3,4,5,6]).reshape(2,3)
B = np.array([6,5,4,3,2,1]).reshape(2,3)

C = A + B
# C = A - B
# C = A * B
# C = A / B
# C = A ** B
print(C)
In [ ]:
A = np.array([1,2,3,4,5,6]).reshape(2,3)
B = np.array([6,5,4,3,2,1]).reshape(3,2)

C = A.dot(B)
print(C)

從矩陣中選取各項或子矩陣¶

Selecting an entry or a submatrix from a matrix

若 A 是一個矩陣。

  1. A[i,j]：選取第 ij 項。
  2. A[list1, list2]：選取列在 list1 中行在 list2 中的子矩陣。

也可以混合使用﹐如 A[i, list] 或 A[list, j]。

Let A be a matrix.

  1. A[i,j] : the ij entry of A .
  2. A[list1, list2] : the submatrix of $A$ induced on rows in list1 and columns in list2 .

One may also mix the two usages, such as A[i, list] or A[list, j] .

In [ ]:
A = np.array([1,2,3,4,5,6]).reshape(2,3)
print(A)
print(A[0,1])
print(A[[0,1],[1,2]])

選取子矩陣中 list 的可以用 a:b 的格式取代。

  1. a:b：從 a 到 b（不包含 b）。
  2. a:：從 a 到底。
  3. :b：從頭到 b（不包含 b）。
  4. :：全部。

When selecting a submatrix, the argument list can be replaced by a:b .

  1. a:b : from a to b (excluding b ).
  2. a: : from a to the end.
  3. :b : from the beginning to b (excluding b ).
  4. : : all.
In [ ]:
A = np.array([1,2,3,4,5,6]).reshape(2,3)
print(A[:,1:])

可以把選出來的子矩陣設定成給定的矩陣。

We may also assign values to the selected submatrix.

In [ ]:
A = np.zeros((2,3) )
print(A)
A[0,0] = 100
print(A)
A[:,1:] = np.eye(2)
print(A)

線性代數上的運算¶

Operations in linear algebra

若 A 為一矩陣。

  1. A.T：A 的轉置。
  2. LA.det：A 的行列式值。
  3. np.trace(A)：A 的跡。
  4. LA.inv：A 的反矩陣。
  5. np.poly：A 的特徵多項式。

Let A be a matrix.

  1. A.T : the transpose of A .
  2. LA.det : the determinant of A .
  3. np.trace(A) : the trace of A .
  4. LA.inv : the inverse of A .
  5. np.poly : the characteristic polynomial of A .
In [ ]:
A = np.array([1,2,3,4,5,6]).reshape(2,3)
print(A.T)
In [ ]:
A = np.array([1,2,3,4]).reshape(2,2)
print(LA.det(A))
In [ ]:
A = np.array([1,2,3,4]).reshape(2,2)
print(np.trace(A))
In [ ]:
A = np.array([1,2,3,4]).reshape(2,2)
Ainv = LA.inv(A)
print("A^{-1} =")
print(Ainv)
print("A A^{-1} =")
print(A.dot(Ainv))
In [ ]:
A = np.array([1,2,3,4]).reshape(2,2)
p = np.poly(A)
print("characteristic polynomial has coefficients")
print(p)
### Cayley--Hamilton Theorem
print("p_A(A) =")
print(p[0] * A.dot(A) + p[1] * A + p[2] * np.eye(2))

特徵值、特徵向量、對角化¶

Eigenvalues, eigenvectors, and diagonalization

若 A 為一矩陣。

  1. eig 或 eigh：回傳一個列表及一個矩陣，列表為 A 的特徵值，而矩陣的行向量為 A 的特徵向量。
  2. eigvals 或 eigvalsh：回傳 A 的特徵值。

其中有 h 的版本是專為對稱矩陣而設計，並有以下特點：

  • 回傳的特徵值都是實數，且由小到大排列。
  • 回傳的特徵向量會互向垂直，且長度是一。
  • 只適用對稱矩陣。

Let A be a matrix.

  1. eig or eigh : returns a list and a matrix, where the list consists of the eigenvalues of A , while the matrix contains the eigenvectors of A as the columns.
  2. eigvals or eigvalsh : only the eigenvalues of A .

Here the version with an extra h is customized for Hermitian matrices and has the following features:

  • The eigenvalues are all real, and they are sorted from small to large in the output.
  • The eigenvectors are orthonormal.
  • It only works for Hermitian matrices.
In [ ]:
A = np.array([1,2,3,4]).reshape(2,2)
vals, vecs = LA.eig(A)
print("eigenvalues =")
print(vals)
print("eigenvectors are the columns of")
print(vecs)
In [ ]:
A = np.array([1,2,3,4]).reshape(2,2)
vals, vecs = LA.eig(A)
### A v = lambda v
k = 0
print("A v")
print(A.dot(vecs[:,k]))
print("lambda v")
print(vals[k] * vecs[:,k])
In [ ]:
A = np.array([1,2,3,4]).reshape(2,2)
vals, vecs = LA.eig(A)
### D = Qinv A Q
print("D =")
print(np.diag(vals))
print("Q^{-1} A Q =")
print(LA.inv(vecs).dot(A).dot(vecs))
In [ ]:
A = np.array([1,2,3,4]).reshape(2,2)
vals, vecs = LA.eig(A)
### e^A = Q e^D Qinv
eD = np.diag(np.exp(vals))
print("e^D =")
print(eD)
eA = vecs.dot(eD).dot(LA.inv(vecs))
print("e^A =")
print(eA)