bastard.py

def bastard_svd(A, k, p):
    ########################
    # Stage A              #
    ########################
    m, n = A.shape
    assert(k + p <= n)
    # (1)
    # Create diagonal matrix of random numbers
    D = np.diag(np.random.normal(size=n))
    # Create vector J of length k + p from the set {1, 2, ..., n}
    J = np.random.choice(np.arange(n),
                         size=(k + p),
                         replace=False)
    # (2)
    # Apply full FFT to rows of AD, extract k + p columns
    Z = np.fft.fft(A.dot(D), axis=1)  # Axis0 is columns, Axis1 is Rows
    Y = Z[:, J]
    # (3)
    # Create Orthonormal matrix Q by orthonomalizing columns of Y
    # Q = orth(Y)
    Q = slg.orth(Y)
    ########################
    # Stage B              #
    ########################
    # (4)
    # B = Q* A
    B = Q.T.dot(A)
    # (5)
    # SVD decomposition of B, s is vector of singular values.
    Uhat, s, V = np.linalg.svd(B, full_matrices=False)  # equiv to 'econ'
    #Uhat, s, V = np.linalg.svd(B)  # equiv to 'econ'
    D = np.diag(s)
    print('{} {} {}'.format(Uhat.shape, D.shape, V.shape))
    # (6)
    # U + Q Uhat
    U = Q.dot(Uhat)
    # Return decomposition
    print(np.linalg.norm(A - U.dot(D).dot(V)))
    return U, D, V