bastardized_svd

def bastard_svd(A, k, p):                                                   # function [U,D,V] = LOCAL_rsft(A,k,p)
    m, n = A.shape                                                          # n          = size(A,2);
    # Create diagonal matrix of random numbers                              #
    dd = np.exp(1j * 2 * np.pi * np.random.random(size=n))                  # dd         = exp(1i*2*pi*rand(1,n));
    Y = A @ np.diag(dd)                                                     # Y          = A*diag(dd);
    # Apply full FFT to rows of AD, extract k + p columns                   #
    Y = np.fft.fft(A, axis=1)  # Axis0 is columns, Axis1 is Rows            # Y          = fft(Y,[],2);
    # Create vector J of length k + p from the set {1, 2, ..., n}           #
    J = np.random.choice(np.arange(n),                                      # [~,ind]    = sort(rand(1,n));
                         size=(k + p),                                      # J          = ind(1:(k+p));
                         replace=False)                                     #
    Y = Y[:, J]                                                             # Y          = Y(:,J);
    # Create Orthonormal matrix Q by orthonomalizing columns of Y           #
    Q, R = slg.qr(Y)                                                        # [Q,~,~]    = qr(Y,0);
    B = Q.T @ A                                                             # B          = Q'*A;
    # SVD decomposition of B, s is vector of singular values.               $
    Uhat, s, Vh = slg.svd(B)                                                # [Uhat,D,V] = svd(B,'econ');
    D = np.zeros((m, n))                                                    #
    D[:len(s), :len(s)] = np.diag(s)                                        #
    V = Vh.T                                                                #
    U = Q @ Uhat[:, :k]                                                     # U          = Q*Uhat(:,1:k);
    D = D[:k, :k]                                                           # D          = D(1:k,1:k);
    V = V[:, :k]                                                            # V          = V(:,1:k);
    # Return decomposition                                                  #
    print('{} {} {}'.format(Uhat.shape, D.shape, V.shape))                  #
    print(np.linalg.norm(A - U @ D @ V.T))                                  #
    return U, D, V                                                          # return