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

python/sympy计算施密特正交化向量sympy的符号计算功能很强大,学习矩阵分析,重温了线性代数中施密特正交化的方法,正好可以用sympy解决一些计算问题。施密特正交化,也称Gram-Schmidt正交化过程(Gram–SchmidtOrthogonalizationProcedure).该⽅法以JørgenP.Gram和ErhardSchmidt命名,它更早出现在拉普拉斯和柯西的⽂章中[1],步骤如下:…

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

sympy的符号计算功能很强大,学习矩阵分析,重温了线性代数中施密特正交化的方法,正好可以用sympy解决一些计算问题。施密特正交化,也称 Gram-Schmidt 正交化过程 (Gram–Schmidt Orthogonalization Procedure). 该⽅法以Jørgen P. Gram 和 Erhard Schmidt 命名, 它更早出现在拉普拉斯和柯西的⽂章中[1],步骤如下:

(1) β 1 = α 1 \beta_1=\alpha_1 β1=α1
(2) β j = α j − ( β 1 , α j ) ( β 1 , β 1 ) β 1 − ( β 2 , α j ) ( β 2 , β 2 ) β 2 − ( β j − 1 , α j ) ( β j − 1 , β j − 1 ) β 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} βj=αj(β1,β1)(β1,αj)β1(β2,β2)(β2,αj)β2(βj1,βj1)(βj1,αj)βj1
(3) η j = β j ∣ ∣ β j ∣ ∣ , j = 1 , 2 , . . . , m \eta_j=\frac{\beta_j}{||\beta_j||},j=1,2,…,m ηj=βjβj,j=1,2,...,m

例如求3个无关列向量的正交向量组:
α 1 = ( 1 , 2 , 3 ) T \alpha_1=(1,2,3)^T α1=(1,2,3)T, α 2 = ( 2 , 1 , 3 ) T \alpha_2=(2,1,3)^T α2=(2,1,3)T, α 3 = ( 3 , 2 , 1 ) T \alpha_3=(3,2,1)^T α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 ) = δ i j (\alpha_i,\alpha_j)=\delta_{ij} (αi,αj)=δ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]])

验证:

In [11]: m=b[:,0]

In [12]: n=b[:,1]

In [13]: p=b[:,2]

In [14]: m.dot(n)
Out[14]: -5.551115123125783e-17

In [15]: m.dot(p)
Out[15]: -2.7755575615628914e-17

In [16]: sum([i**2 for i in m])
Out[16]: 0.9999999999999993

关键点

(1)需要正交化的向量写成矩阵的时候应该注意是矩阵每一列代表原来需要正交化的向量;
(2)数学中向量一般都是列向量,书本上为了节省空间通常写作横向量,严谨点的写法应该是横向量后面添加转置符号,但发现目前很多书籍、教材等都没有特意说明,在机器学习里面有很多以向量、矩阵表达的内容,很容易搞得晕头转向;
(3)正交这个概念在数学中很多方面都有涉及,例如

  • 正交向量:向量内积为零, ( x , y ) = 0 (\mathbf{x},\mathbf{y})=0 (x,y)=0
  • 正交矩阵:有 A ∈ R n × n A \in \mathbb{R}^{n \times n} ARn×n,若 A T A = I A^{\rm{T}}A=I ATA=I,即 A − 1 = A T A^{\rm{-1}}=A^{\rm{T}} A1=AT,称 A A A为正交矩阵
  • 正交函数:如果 ( f , g ) = ∫ a b ρ ( x ) f ( x ) g ( x ) d x = 0 (f,g)=\int_a^b\rho(x)f(x)g(x)dx=0 (f,g)=abρ(x)f(x)g(x)dx=0,则称函数 f f f g g g在区间 [ a , b ] [a,b] [a,b]上带权 ρ ( x ) \rho(x) ρ(x)正交
    参考资料
    [1] 黄正华 武汉大学http://aff.whu.edu.cn/huangzh/

今天的文章python/sympy计算施密特正交化向量分享到此就结束了,感谢您的阅读,如果确实帮到您,您可以动动手指转发给其他人。

版权声明:本文内容由互联网用户自发贡献,该文观点仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 举报,一经查实,本站将立刻删除。
如需转载请保留出处:https://bianchenghao.cn/32168.html

(0)
编程小号编程小号

相关推荐

发表回复

您的电子邮箱地址不会被公开。 必填项已用*标注