# Python: NumPy 與數值線性代數

Python: NumPy and numerical linear algebra

This work by Jephian Lin is licensed under a [Creative Commons Attribution 4.0 International License](http://creativecommons.org/licenses/by/4.0/).

import numpy as np
from numpy import linalg as LA

## 符號運算與數值運算

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

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

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

<!-- eng start -->
**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.
<!-- eng end -->

## 建構矩陣

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` 來顯示矩陣。

<!-- eng start -->
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.
<!-- eng end -->

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)`。

<!-- eng start -->
Different from Sage, most of the operations in NumPy is entrywise. We may use `A.dot(B)` for the usual matrix multiplication.
<!-- eng end -->

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` 是一個矩陣。

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

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

<!-- eng start -->
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]` .
<!-- eng end -->

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. `:`：全部。

<!-- eng start -->
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.
<!-- eng end -->

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

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

<!-- eng start -->
We may also assign values to the selected submatrix.
<!-- eng end -->

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` 的特徵多項式。

<!-- eng start -->
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` .
<!-- eng end -->

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` 為一矩陣。

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

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

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

<!-- eng start -->
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.
<!-- eng end -->

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)