Python实现QR分解

使用Gram-Schmidt正交化,Householder变换,Given旋转三种方法实现了QR分解,但是感觉好像并无卵用,貌似实际生产有更好的改进方法?
矩阵分析或者矩阵论的课都有介绍,网上资料也很多,在此不贴了。

import numpy as np

def SchmitOrth(mat:np.array):
    cols = mat.shape[1]

    Q = np.copy(mat)
    R = np.zeros((cols, cols))

    for col in range(cols):
        for i in range(col):
            k =  np.sum(mat[:, col] * Q[:, i]) / np.sum( np.square(Q[:, i]) )
            Q[:, col] -= k*Q[:, i]
        Q[:, col] /= np.linalg.norm(Q[:, col])

        for i in range(cols):
            R[col, i] = Q[:, col].dot( mat[:, i] )

    return Q, R
def HouseHolder(mat:np.array):
    cols = mat.shape[1]

    Q = np.eye(cols)
    R = np.copy(mat)

    for col in range(cols-1):
        a = np.linalg.norm(R[col:, col])
        e = np.zeros((cols- col))
        e[0] = 1.0
        num = R[col:, col] -a*e
        den = np.linalg.norm(num)

        u = num / den
        H = np.eye((cols))
        H[col:, col:] = np.eye((cols- col))- 2*u.reshape(-1, 1).dot(u.reshape(1, -1))
        R = H.dot(R)

        Q = Q.dot(H)

    return Q, R
def GivenRot(mat:np.array):
    rows, cols = mat.shape

    R = np.copy(mat)
    Q = np.eye(cols)

    for col in range(cols):
        for row in range(col+1, rows):
            if abs(R[row, col]) < 1e-6:
                continue

            f = R[col, col]
            s = R[row, col]
            den = np.sqrt( f*f+ s*s)
            c = f / den
            s = s / den

            T = np.eye(rows)
            T[col, col], T[row, row] = c, c
            T[row, col], T[col, row] = -s, s

            R = T.dot(R)

            Q = T.dot(Q)
    
    return Q.T, R

mat = np.array( [   [1.0, 2.0, 2.0],
                    [2.0, 1.0, 2.0],
                    [1.0, 2.0, 1.0],  ])
q, r = SchmitOrth(mat)
print("Gram-Schmidt Orthogonal")
print("Q: \n", q)
print("R: \n", r)
print("Q x R: \n", q.dot(r))
#
q, r = HouseHolder(mat)
print("HouseHolder")
print("Q: \n", q)
print("R: \n", r)
print("Q x R: \n", q.dot(r))
#
q, r = GivenRot(mat)
print("Givens")
print("Q: \n", q)
print("R: \n", r)
print("Q x R: \n", q.dot(r))      
最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
平台声明:文章内容(如有图片或视频亦包括在内)由作者上传并发布,文章内容仅代表作者本人观点,简书系信息发布平台,仅提供信息存储服务。