使用QR_Basic, QR_Heisenberg和经典Jacobi三种方法实现了SVD分解，迭代使用常数退出，没用误差退出是因为写QR分解时懒的写，到SVD更不想写了囧。
具体资料网上很多，对着代码瞅瞅就完事了。

import numpy as np

def QR_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 QR_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
#---------------------------------------
def QR_Basic(mat:np.array):
    eig = np.copy(mat)

    rows, cols = mat.shape
    eigV = np.eye(cols)
    for _ in range(50):
        q, r = QR_HouseHolder(eig)
        eig = r.dot(q)
        eigV = eigV.dot(q)

    return np.diag(eig), eigV
#---------------------------------------
def QR_Heisenberg(mat:np.array):

    rows, cols = mat.shape
    Hsbg = np.copy(mat)

    eigv = np.eye(cols)
    for row in range(1, rows-1):
        array = mat[row:, row-1]

        #householder
        a = np.linalg.norm(array)

        e = np.zeros((1, rows- row))[0]
        e[0] = 1.0
        num = array-a*e
        u = num/ np.linalg.norm(num)

        H = np.eye(rows- row)- 2*u.reshape(-1, 1).dot(u.reshape(1, -1))

        #Hsbg
        T = np.eye(rows)
        T[row:, row:] = H

        Hsbg = T.T.dot(Hsbg).dot(T)
        eigv = eigv.dot(T.T)

    eigv = eigv.T
    # Given
    for _ in range(50):
        q, r = QR_GivenRot(Hsbg)
        Hsbg = r.dot(q)
        eigv = eigv.dot(q)
    return np.diag(Hsbg), eigv
    
#-----------------------------
def Jacobi_Basic(mat:np.array):
    rows, cols = mat.shape


    eig = np.copy(mat)
    eigV = np.eye(rows, cols)
    for _ in range(8):
        maxRow, maxCol = 0, 0
        maxValue = -1
        for row in range(rows):
            for col in range(cols):
                if row == col:
                    continue
                if abs(eig[row, col]) > maxValue:
                    maxValue = abs(eig[row, col])
                    maxRow = row
                    maxCol = col
        a = 0
        if abs(eig[maxRow, maxRow] - eig[maxCol, maxCol]) < 1e-5:
            a = 0.25*np.pi 
        else:
            a = 0.5* np.arctan(  (2*eig[maxRow, maxCol]) /(eig[maxRow, maxRow]- eig[maxCol, maxCol]) )

        P = np.eye(rows)
        P[maxRow, maxRow] = np.cos(a)
        P[maxCol, maxCol] = np.cos(a)
        P[maxRow, maxCol] = -np.sin(a)
        P[maxCol, maxRow] = np.sin(a)
        eig = P.T.dot(eig).dot(P)
        eigV = eigV.dot(P)
    return np.diag(eig), eigV

#----------------
def SortEig(eig, eigV):
    e = np.copy(eig)
    v = np.copy(eigV)
    indices = np.argsort(e)
    e= np.sort(e)
    for i, j in enumerate(indices):
        v[:, i] = eigV[:, j]

    return e, v

def SVD_Solver(mat, Func):
    aTa = mat.T.dot(mat)
    aaT = mat.T.dot(mat)

    eig, eigv = Func(aTa)
    eig, eigv = SortEig(eig, eigv)

    rows, cols = mat.shape
    U = eigv
    S = np.zeros_like(mat)
    for i in range(rows):
        S[i, i] = eig[i]
    
    eig, eigv = Func(aaT)
    eig, eigv = SortEig(eig, eigv)
    V = eigv

    return U, np.sqrt(S), V

def SVD(mat:np.array, Method:int):
    #Method: 0: QR_Basic, 1:QR_Heisenberg, 2:Jacobi_Basic
    if Method == 0:
        U, S, V = SVD_Solver(mat, QR_Basic)
    elif Method == 1:
        U, S, V = SVD_Solver(mat, QR_Heisenberg)
    else:
        U, S, V = SVD_Solver(mat, Jacobi_Basic)
    return U, S, V


mat = np.array( [   [3.0, 2.0, 2.0],
                    [2.0, 5.0, 1.0],
                    [2.0, 1.0, 4.0],  ])

u, s, v=SVD(mat, 0)
print(u.dot(s).dot(v.T))
u, s, v=SVD(mat, 1)
print(u.dot(s).dot(v.T))
u, s, v=SVD(mat, 2)
print(u.dot(s).dot(v.T))