Golub-Kahan-Lanczos Bidiagonalization Procedure (implementation of http://www.netlib.org/utk/people/JackDongarra/etemplates/node198.html for the case of lower-bidiagonal matrix)

bidiagonalization.py

import numpy as np
import scipy.sparse as scsp

def GKL_bidiagonalization(A):
    A = np.matrix(A)
    m,n = A.shape
    
    p = max(m, n)
    alpha = np.zeros(p)
    beta = np.zeros(p-1)
    U = np.zeros((p, m))  # in fact it's U.H
    V = np.zeros((p, n))  # in fact it's V.H
    # Transposed matrices are used for easier slicing
    U[0, 0] = 1
    
    for i in range(p):
        V[i] = A.H@U[i] - beta[i-1]*V[i-1]
        alpha[i] = np.linalg.norm(V[i])
        V[i] /= alpha[i]

        if i > p - 2: continue
        U[i+1] = A@V[i] - alpha[i]*U[i]
        beta[i] = np.linalg.norm(U[i+1])
        U[i+1] /= beta[i]
    U,V = map(np.matrix, (U, V))
        
    return U.H, (alpha, beta), V.H 

if __name__ == '__main__':
    A = [[1, 2], [2, 3], [3, 5]]
    U, B, V = GKL_bidiagonalization(A)
    B = scsp.diags(B, [0, -1]).toarray()
    print('original\n', np.array(A))
    print('reconstructed\n', U@B@V.H)

Very good algoritm!

But here's a funny thing,
I tested with two different matrices, one with the correct awnser and the other with the "correct awnser*2"

Here's the output from the 'correct awnser*2' tested matrix

original
 [[ 1  2  3]
 [ 4  5  6]
 [ 7  8  9]
 [10 11 12]]
reconstructed
 [[ 2.67551421  3.90617858  5.13684295]
 [ 8.35198497  9.9511132  11.55024143]
 [14.02845573 15.99604782 17.9636399 ]
 [19.70492648 22.04098243 24.37703838]]

This matrix comes from an example of the book "Matrix Computations" from Golub & Van Loan (chapter 5)