学术, 其他笔记

施密特正交化(附Python代码)

施密特正交化(Schmidt orthogonalization)过程[1]:

Python代码:

"""
This code is supported by the website: https://www.guanjihuan.com
The newest version of this code is on the web page: https://www.guanjihuan.com/archives/10890
"""

import numpy as np


def main():
    A = np.array([[0, 1, 1, -1], [1, 0, -1, 1], [1, -1, 0, 1], [-1, 1, 1, 0]])
    eigenvalue, eigenvector = np.linalg.eig(A)
    print('矩阵:\n', A)
    print('特征值:\n', eigenvalue)
    print('特征向量:\n', eigenvector)

    print('\n判断是否正交:\n', np.dot(eigenvector.transpose(), eigenvector))
    print('判断是否正交:\n', np.dot(eigenvector, eigenvector.transpose()))

    print('对角化验证:')
    print(np.dot(np.dot(eigenvector.transpose(), A), eigenvector))

    # 施密斯正交化
    eigenvector = Schmidt_orthogonalization(eigenvector)

    print('\n施密斯正交化后,特征向量:\n', eigenvector)

    print('施密斯正交化后,判断是否正交:\n', np.dot(eigenvector.transpose(), eigenvector))
    print('施密斯正交化后,判断是否正交:\n', np.dot(eigenvector, eigenvector.transpose()))

    print('施密斯正交化后,对角化验证:')
    print(np.dot(np.dot(eigenvector.transpose(), A), eigenvector))


def Schmidt_orthogonalization(eigenvector):
    num = eigenvector.shape[1]
    for i in range(num):
        for i0 in range(i):
            eigenvector[:, i] = eigenvector[:, i] - eigenvector[:, i0]*np.dot(eigenvector[:, i].transpose().conj(), eigenvector[:, i0])/(np.dot(eigenvector[:, i0].transpose().conj(),eigenvector[:, i0]))
        eigenvector[:, i] = eigenvector[:, i]/np.linalg.norm(eigenvector[:, i])
    return eigenvector


if __name__ == '__main__':
    main()

运行结果:

矩阵:
 [[ 0  1  1 -1]
 [ 1  0 -1  1]
 [ 1 -1  0  1]
 [-1  1  1  0]]
特征值:
 [ 1. -3.  1.  1.]
特征向量:
 [[ 0.8660254   0.5         0.12803516  0.08530895]
 [ 0.28867513 -0.5         0.85020245 -0.16863006]
 [ 0.28867513 -0.5        -0.36108364  0.80962746]
 [-0.28867513  0.5         0.36108364  0.55568845]]

判断是否正交:
 [[ 1.00000000e+00  5.55111512e-17  1.47842273e-01  9.85062898e-02]
 [ 5.55111512e-17  1.00000000e+00  5.55111512e-17 -5.55111512e-17]
 [ 1.47842273e-01  5.55111512e-17  1.00000000e+00 -2.24140369e-01]
 [ 9.85062898e-02 -5.55111512e-17 -2.24140369e-01  1.00000000e+00]]
判断是否正交:
 [[ 1.02367062  0.09447016  0.02283706  0.0936366 ]
 [ 0.09447016  1.08461363 -0.11018839 -0.12004491]
 [ 0.02283706 -0.11018839  1.11921136 -0.0138141 ]
 [ 0.0936366  -0.12004491 -0.0138141   0.77250439]]
对角化验证:
[[ 1.00000000e+00 -1.66533454e-16  1.47842273e-01  9.85062898e-02]
 [-2.22044605e-16 -3.00000000e+00 -2.22044605e-16  0.00000000e+00]
 [ 1.47842273e-01 -1.66533454e-16  1.00000000e+00 -2.24140369e-01]
 [ 9.85062898e-02  1.11022302e-16 -2.24140369e-01  1.00000000e+00]]

施密斯正交化后,特征向量:
 [[ 8.66025404e-01  5.00000000e-01  4.20959579e-17 -3.85081968e-18]
 [ 2.88675135e-01 -5.00000000e-01  8.16496581e-01  2.87496183e-17]
 [ 2.88675135e-01 -5.00000000e-01 -4.08248290e-01  7.07106781e-01]
 [-2.88675135e-01  5.00000000e-01  4.08248290e-01  7.07106781e-01]]
施密斯正交化后,判断是否正交:
 [[1.00000000e+00 0.00000000e+00 4.16333634e-17 0.00000000e+00]
 [0.00000000e+00 1.00000000e+00 2.77555756e-17 0.00000000e+00]
 [4.16333634e-17 2.77555756e-17 1.00000000e+00 0.00000000e+00]
 [0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00]]
施密斯正交化后,判断是否正交:
 [[ 1.00000000e+00  8.32667268e-17  0.00000000e+00  0.00000000e+00]
 [ 8.32667268e-17  1.00000000e+00  0.00000000e+00  2.77555756e-17]
 [ 0.00000000e+00  0.00000000e+00  1.00000000e+00 -8.32667268e-17]
 [ 0.00000000e+00  2.77555756e-17 -8.32667268e-17  1.00000000e+00]]
施密斯正交化后,对角化验证:
[[ 1.00000000e+00 -2.22044605e-16  4.16333634e-17  2.77555756e-17]
 [-1.11022302e-16 -3.00000000e+00 -1.11022302e-16  0.00000000e+00]
 [ 5.55111512e-17 -5.55111512e-17  1.00000000e+00  5.55111512e-17]
 [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  1.00000000e+00]]

另外,这里给的函数可能还不一定稳定,目前还不清楚是什么原因。一个例子为:

"""
This code is supported by the website: https://www.guanjihuan.com
The newest version of this code is on the web page: https://www.guanjihuan.com/archives/10890
"""

import numpy as np

a = [[ 0 , 0  ,  1.5 ,  0.32635182-0.98480775j],
 [0  ,   0  , -0.32635182-0.98480775j, 1.5  ],
 [ 1.5 ,    -0.32635182+0.98480775j ,0, 0 ],
 [ 0.32635182+0.98480775j , 1.5 ,  0, 0 ]]

def Schmidt_orthogonalization(eigenvector):
    num = eigenvector.shape[1]
    for i in range(num):
        for i0 in range(i):
            eigenvector[:, i] = eigenvector[:, i] - eigenvector[:, i0]*np.dot(eigenvector[:, i].transpose().conj(), eigenvector[:, i0])/(np.dot(eigenvector[:, i0].transpose().conj(),eigenvector[:, i0]))
        eigenvector[:, i] = eigenvector[:, i]/np.linalg.norm(eigenvector[:, i])
    return eigenvector

def verify_orthogonality(vectors):
    identity = np.eye(vectors.shape[1])
    product = np.dot(vectors.T.conj(), vectors)
    return np.allclose(product, identity)

# 对 np.linalg.eigh() 的特征向量正交化

E, v = np.linalg.eigh(a)
print(verify_orthogonality(v))

v1 = Schmidt_orthogonalization(v)
print(verify_orthogonality(v1))

from scipy.linalg import orth
v2 = orth(v)
print(verify_orthogonality(v2))

v3, S, Vt = np.linalg.svd(v)
print(verify_orthogonality(v3))

v4, R = np.linalg.qr(v)
print(verify_orthogonality(v4))

print()


# 对 np.linalg.eig() 的特征向量正交化

E, v = np.linalg.eig(a)
print(verify_orthogonality(v))

v1 = Schmidt_orthogonalization(v)
print(verify_orthogonality(v1))

from scipy.linalg import orth
v2 = orth(v)
print(verify_orthogonality(v2))

v3, S, Vt = np.linalg.svd(v)
print(verify_orthogonality(v3))

v4, R = np.linalg.qr(v)
print(verify_orthogonality(v4))

运行结果为:

True
True
True
True
True

False
False
True
True
True

说明:这里的例子 np.linalg.qr() 、np.linalg.svd() 和 scipy.linalg.orth() 都成功完成了正交化,但自定义的 Schmidt_orthogonalization() 函数没有完成正交化,暂时不知道是什么原因导致的。个人比较推荐使用 np.linalg.qr(),有时候 np.linalg.svd() 和 scipy.linalg.orth() 在解决实际的物理问题时可能也会出现问题。

可能的原因是:对于有些例子,矩阵中的特征向量非常接近线性相关,那么导致了某些算法在数值计算上不是很稳定,在具体应用时可以考虑使用多个算法进行计算对比。

参考资料:

[1] 截图自北京科技大学廖福成老师的”高等代数与解析几何“课件

[2] ”高等代数与解析几何“、”线性代数“、”高等数学“等相关教材

3,424 次浏览

【说明:本站主要是个人的一些笔记和代码分享,内容可能会不定期修改。为了使全网显示的始终是最新版本,这里的文章未经同意请勿转载。引用请注明出处:https://www.guanjihuan.com

发表评论

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

Captcha Code