bastardized_svd

# Python                                                                    # MATLAB
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(Y, 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, P = slg.qr(Y, pivoting=True, mode='economic')                     # [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, full_matrices=False)                           # [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

With output

(3, 3) (2, 2) (4, 2)
29.0599081337
(15, 15) (5, 5) (140, 5)
2.27895627308
(20, 20) (10, 10) (140, 10)
2.62638440458
(25, 25) (15, 15) (140, 15)
2.69751998152
(30, 30) (20, 20) (140, 20)
2.4971415018
(35, 35) (25, 25) (140, 25)
2.37175522297
(40, 40) (30, 30) (140, 30)
2.75521402538
(45, 45) (35, 35) (140, 35)
2.49978818076
(50, 50) (40, 40) (140, 40)
2.65377354541
(55, 55) (45, 45) (140, 45)
2.46002810262
(60, 60) (50, 50) (140, 50)
2.85140300167
(65, 65) (55, 55) (140, 55)
2.47002074377
(70, 70) (60, 60) (140, 60)
2.74887292659
(75, 75) (65, 65) (140, 65)
2.07012791991
(80, 80) (70, 70) (140, 70)
2.6846967836
(85, 85) (75, 75) (140, 75)
2.43865242545
(90, 90) (80, 80) (140, 80)
2.56051291336