栏目分类:
子分类:
返回
名师互学网用户登录
快速导航关闭
当前搜索
当前分类
子分类
实用工具
热门搜索
名师互学网 > IT > 软件开发 > 后端开发 > Python

Orthogonalization: Gram-Schmidt,Householder,Givens

Python 更新时间: 发布时间: IT归档 最新发布 模块sitemap 名妆网 法律咨询 聚返吧 英语巴士网 伯小乐 网商动力

Orthogonalization: Gram-Schmidt,Householder,Givens

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

Gram-Schmidt 方1
# -*- 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))
转载请注明:文章转载自 www.mshxw.com
本文地址:https://www.mshxw.com/it/269327.html
我们一直用心在做
关于我们 文章归档 网站地图 联系我们

版权所有 (c)2021-2022 MSHXW.COM

ICP备案号:晋ICP备2021003244-6号