Complete Singular Value Decomposition

[ks147]

Currently sympy has function returning condensed form of singular value decomposition. I was thinking of also implementing the complete SVD algorithm, because for some applications you may require the complete decomposition.
The code I've implemented is not working for some Matrices. If somebody could point out where I am going wrong that would be very helpful

def complete_singular_value_decomposition(A):
    AH = A.H
    m, n = A.shape
    if m >= n:
        
        U, S = (A * AH).diagonalize()  
        
        Singular_vals = S.diagonal()           
        Singular_vals = [sqrt(x) for x in Singular_vals]        
        
        U, _ = U.QRdecomposition()        
        
        V, S = (AH * A).diagonalize()

        Singular_vals = S.diagonal()      
        
        Singular_vals = [sqrt(x) for x in Singular_vals]        
        S = diag(Singular_vals, unpack=True)
        
        V,_ = V.QRdecomposition()       
    
    else:
        
        V, S = (AH * A).diagonalize()

        Singular_vals = S.diagonal()      
        
        Singular_vals = [sqrt(x) for x in Singular_vals]
        
        V,_ = V.QRdecomposition()   
        
        U, S = (A * AH).diagonalize()  
        
        Singular_vals = [sqrt(x) for x in Singular_vals]        
        S = diag(Singular_vals, unpack=True)
        
        U, _ = U.QRdecomposition()

    return U, S, V

In [1]:  C = Matrix([
                    [4, -2, 0], 
                    [-2, 3, 2], 
                    [0, 2, 2]])
          U,S,V = complete_singular_value_decomposition(C)
          assert V * V.T == eye(V.rows)
          assert (V.T * V) == eye(V.cols)
          assert (U.T * U) == eye(U.rows)
          assert(U * U.T) == eye(U.cols)
          assert U * S * V.T == C

Out [1]:

In [2]:  A = Matrix([[1, 2], [2, 1]])
          U, S, V = complete_singular_value_decomposition(A)
          assert V * V.T == eye(V.rows)
          assert (V.T * V) == eye(V.cols)
          assert (U.T * U) == eye(U.rows)
          assert(U * U.T) == eye(U.cols)
          pprint(U * S * V.T)
          U * S * V.T == A

Out [2]:

          ⎡2  1⎤
          ⎢    ⎥
          ⎣1  2⎦
          False

[sylee957]

The difference of your implementation is that you try to find other singular vector in parallel way.
I'm not sure that they can be computed in such parallel way because singular vectors are not unique (You can multiply +1, -1 randomly across each columns to get other set of orthogonal singular vectors)