Python: NumPy and numerical linear algebra
This work by Jephian Lin is licensed under a Creative Commons Attribution 4.0 International License.
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.
np.array( list of lists )
:把 list of lists
中的每個 list
當作矩陣的列。np.array(list).reshape(m,n)
:把 list
重組成 m x n
的矩陣。np.eye(n)
:單位矩陣。np.zeros((m,n))
:全零矩陣。np.diag(list)
:對角矩陣,其對角線上元素由 list
決定。利用 print
來顯示矩陣。
np.array( list of lists )
: construct an array whose rows are the list
in list of lists
.np.array(list).reshape(m,n)
: make list
into an m x n
matrix.np.eye(n)
: identity matrix.np.zeros((m,n))
: zero matrix.np.diag(list)
: diagonal matrix whose entries are determined by list
.Use print
to show the matrix.
A = np.array([[1,2,3],
[4,5,6]])
print(A)
A = np.array([1,2,3,4,5,6]).reshape(2,3)
print(A)
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.
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)
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
是一個矩陣。
A[i,j]
:選取第 ij
項。A[list1, list2]
:選取列在 list1
中行在 list2
中的子矩陣。也可以混合使用﹐如 A[i, list]
或 A[list, j]
。
Let A
be a matrix.
A[i,j]
: the ij
entry of A
.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]
.
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
的格式取代。
a:b
:從 a
到 b
(不包含 b
)。a:
:從 a
到底。:b
:從頭到 b
(不包含 b
)。:
:全部。When selecting a submatrix, the argument list
can be replaced by a:b
.
a:b
: from a
to b
(excluding b
).a:
: from a
to the end.:b
: from the beginning to b
(excluding b
).:
: all.A = np.array([1,2,3,4,5,6]).reshape(2,3)
print(A[:,1:])
可以把選出來的子矩陣設定成給定的矩陣。
We may also assign values to the selected submatrix.
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
為一矩陣。
A.T
:A
的轉置。LA.det
:A
的行列式值。np.trace(A)
:A
的跡。LA.inv
:A
的反矩陣。np.poly
:A
的特徵多項式。Let A
be a matrix.
A.T
: the transpose of A
.LA.det
: the determinant of A
.np.trace(A)
: the trace of A
.LA.inv
: the inverse of A
.np.poly
: the characteristic polynomial of A
.A = np.array([1,2,3,4,5,6]).reshape(2,3)
print(A.T)
A = np.array([1,2,3,4]).reshape(2,2)
print(LA.det(A))
A = np.array([1,2,3,4]).reshape(2,2)
print(np.trace(A))
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))
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
為一矩陣。
eig
或 eigh
:回傳一個列表及一個矩陣,列表為 A
的特徵值,而矩陣的行向量為 A
的特徵向量。eigvals
或 eigvalsh
:回傳 A
的特徵值。其中有 h
的版本是專為對稱矩陣而設計,並有以下特點:
Let A
be a matrix.
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.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:
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)
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])
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))
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)