1. 程式人生 > >python/sympy計算施密特正交化向量

python/sympy計算施密特正交化向量

sympy的符號計算功能很強大,學習矩陣分析,重溫了線性代數中施密特正交化的方法,正好可以用sympy解決一些計算問題。施密特正交化,也稱 Gram-Schmidt 正交化過程 (Gram–Schmidt Orthogonalization Procedure). 該⽅法以Jørgen P. Gram 和 Erhard Schmidt 命名, 它更早出現在拉普拉斯和柯西的⽂章中[1],步驟如下:

(1) β1=α1\beta_1=\alpha_1
(2) βj=αj(β1,αj)(β1,β1)β1(β2,αj)(β2,β2)β2(βj1,αj)(βj1,βj1)βj

1\beta_j=\alpha_j-\frac{(\beta_1,\alpha_j)}{(\beta_1,\beta_1)}\beta_1-\frac{(\beta_2,\alpha_j)}{(\beta_2,\beta_2)}\beta_2-\frac{(\beta_{j-1},\alpha_j)}{(\beta_{j-1},\beta_{j-1})}\beta_{j-1}
(3)ηj=βjβj,j=1,2,...,
m\eta_j=\frac{\beta_j}{||\beta_j||},j=1,2,...,m

例如求3個無關列向量的正交向量組:
α1=(1,2,3)T\alpha_1=(1,2,3)^T,α2=(2,1,3)T\alpha_2=(2,1,3)^T,α3=(3,2,1)T\alpha_3=(3,2,1)^T
手動計算可以按照前面的步驟求解,下面介紹sympy相關函式,很簡單,如下面程式碼所示:

from sympy import *

L = [Matrix([1,
2,3]), Matrix([2,1,3]), Matrix([3,2,1])] o1=GramSchmidt(L) o1 Out[4]: [Matrix([ [1], [2], [3]]), Matrix([ [15/14], [ -6/7], [ 3/14]]), Matrix([ [ 4/3], [ 4/3], [-4/3]])] o2=GramSchmidt(L,True) # 標準化 o2 Out[6]: [Matrix([ [ sqrt(14)/14], [ sqrt(14)/7], [3*sqrt(14)/14]]), Matrix([ [ 5*sqrt(42)/42], [-2*sqrt(42)/21], [ sqrt(42)/42]]), Matrix([ [ sqrt(3)/3], [ sqrt(3)/3], [-sqrt(3)/3]])]

函式GramSchmidt(vlist, orthonormal=False),將引數orthonormal設為True 計算結果便是標準正交基,記作:
(αi,αj)=δij(\alpha_i,\alpha_j)=\delta_{ij}

注意,scipy.linalg包也提供了函式orth()來計算標準正交基,但是根據文件介紹其使用的方法是SVD奇異值分解,所以求解結果和sympy的不一樣

from scipy.linalg import *
import numpy as np

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

a
Out[7]: 
array([[1, 2, 3],
       [2, 1, 3],
       [3, 2, 1]])

a=a.T

a
Out[9]: 
array([[1, 2, 3],
       [2, 1, 2],
       [3, 3, 1]])

orthogonal_procrustes?

orth(a)
Out[11]: 
array([[-0.56525513,  0.68901653,  0.45358886],
       [-0.47238331,  0.18041382, -0.86273105],
       [-0.67626966, -0.70193096,  0.22350007]])

參考資料
[1] 黃正華 武漢大學http://aff.whu.edu.cn/huangzh/