# Python: NumPy 與數值線性代數

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 可以處理任意維度的陣列）
的數值運算。

## 建構矩陣

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

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)

## 陣列之間的運算

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

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)

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

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

也可以混合使用﹐如 `A[i, list]` 或 `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` 的格式取代。
1. `a:b`：從 `a` 到 `b`（不包含 `b`）。
2. `a:`：從 `a` 到底。
3. `:b`：從頭到 `b`（不包含 `b`）。
4. `:`：全部。

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

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

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

## 線性代數上的運算

若 `A` 為一矩陣。

1. `A.T`：`A` 的轉置。
2. `LA.det`：`A` 的行列式值。
3. `np.trace(A)`：`A` 的跡。
4. `LA.inv`：`A` 的反矩陣。
5. `np.poly`：`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))

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

若 `A` 為一矩陣。

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

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

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

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)