https://wenku.baidu.com/view/b3baedc949649b6648d7475d.html https://wenku.baidu.com/view/8730e159804d2b160b4ec01a.html https://stackoverflow.com/questions/47349233/modified-gram-schmidt-in-python-for-complex-vectors https://stackoverflow.com/questions/53489237/how-can-you-implement-householder-based-qr-decomposition-in-python https://stackoverflow.com/questions/64432175/python-using-givens-rotation-for-qr-decomposition https://en.wikipedia.org/wiki/Householder_transformation https://en.wikipedia.org/wiki/QR_decomposition#Using_Givens_rotations
Gram–Schmidt process, which uses projection
Householder transformation, which uses reflection
Givens rotation
Symmetric orthogonalization, which uses the Singular value decomposition
# -*- coding: utf-8 -*-
"""
@time: 2021-09-23 上午 11:30
"""
import numpy as np
A = np.array([[1,1,1],[1,1,0],[1,0,0]])
B = np.zeros([A.shape[0], A.shape[1]])
for i in range(A.shape[1]):
bi = A[:,i]
# 当前向量与标准正交向量内积
# 求差后 当前向量转为正交向量
for j in range(i):
rij = np.vdot(bi, B[:,j])
bi = bi - rij*B[:,j]
# 标准化
rij = np.linalg.norm(bi)
B[:,i] = bi/rij
print(B)
'''
[[ 0.57735027 0.40824829 0.70710678]
[ 0.57735027 0.40824829 -0.70710678]
[ 0.57735027 -0.81649658 0. ]]
'''
方2
import numpy as np
# 正交化
A = np.array([[1,1,1],[1,1,0],[1,0,0]])
Q, R = np.linalg.qr(A)
B = Q.copy()
# 标准化
for i in range(Q.shape[1]):
B[:,i] = Q[:,i]/np.linalg.norm(Q[:,i])
print(B)
'''
[[-5.77350269e-01 -4.08248290e-01 -7.07106781e-01]
[-5.77350269e-01 -4.08248290e-01 7.07106781e-01]
[-5.77350269e-01 8.16496581e-01 5.55111512e-17]]
'''
Householder变换
Householder transformations are widely used in numerical linear algebra, for example, to annihilate the entries below the main diagonal of a matrix; to perform QR decompositions and in the first step of the QR algorithm.;They are also widely used for transforming to a Hessenberg form.
QR分解,将矩阵A分解为一个正交矩阵Q与一个上三角矩阵R,通过施密特,Householder变换,Givens旋转可实现。
与施密特正交化每次生产一个正交向量不同,householder正交化最终一次生成所有正交向量
import numpy as np
from typing import Union
def householder_vectorized(a):
# copysign将第一个元素的符号赋给其他元素
v = a / (a[0] + np.copysign(np.linalg.norm(a), a[0]))
v[0] = 1
tau = 2 / (v.T @ v)
return v,tau
def qr_decomposition(A: np.ndarray) -> Union[np.ndarray, np.ndarray]:
m,n = A.shape
R = A.copy()
Q = np.identity(m)
for j in range(0, n):
# Apply Householder transformation.
v, tau = householder_vectorized(R[j:, j, np.newaxis])
H = np.identity(m)
H[j:, j:] -= tau * (v @ v.T)
R = H @ R
Q = H @ Q
return Q[:n].T, np.triu(R[:n])
m = 5
n = 4
A = np.array([[1,1,1],[1,1,0],[1,0,0]])
q, r = np.linalg.qr(A)
Q, R = qr_decomposition(A)
with np.printoptions(linewidth=9999, precision=20, suppress=True):
print("**** Q from qr_decomposition")
print(Q)
print("**** Q from np.linalg.qr")
print(q)
print()
print("**** R from qr_decomposition")
print(R)
print("**** R from np.linalg.qr")
print(r)
Givens旋转
import numpy as np
def QRrot(A):
m,n = A.shape
indexI = np.zeros([m,n])
indexJ = np.zeros([m,n])
C = np.zeros([m,n])
S = np.zeros([m,n])
for i in range(1,n):
for j in range(i+1, m):
c = A[i,i]/((A[i,i])**2 + (A[j,i]**2))**0.5
s = A[j,i]/((A[i,i])**2 + (A[j,i]**2))**0.5
A[i,:] = c*A[i,:] + s*A[j,:]
A[j,:] = -s*A[i,:] + c*A[j,:]
indexI[j,i] = i
indexJ[j,i] = j
C[j,i] = c
S[j,i] = s
R = A.copy()
Q = np.eye(m)
Q[:,i] = c*Q[:,i] + s*Q[:,j]
Q[:,j] = -s*Q[:,i] + c*Q[:,j]
return Q,R
# 原矩阵
B = np.random.uniform(3,10,[3,3])
print(B)
# 分解后的矩阵
Q, R = QRrot(B)
print(Q, R)
# 恢复出原矩阵
print(np.dot(Q, R))



