如果你想在github上从LJWillams实现工作，有几个bug。然而，为了完全理解这项技术，我可能会建议您自己尝试一下实现它。我查了Octave（MATLAB免费软件等价物）做什么和他们的"code is a wrapper to the corresponding Lapack dggsvd and zggsvd routines."，这是Scipy应该做的IMHO。
我会发布我发现的bug，但我不打算发布完整的工作顺序，因为我不确定这对版权有何影响，因为它是从受版权保护的MATLAB实现翻译而来的。

注意事项：我不是广义SVD的Maven，并且仅从调试的Angular 来处理这个问题，而不是底层算法是否正确。我已经在你的原始随机数组和Python文件中已经存在的测试用例上工作了。

Bug

设置k

在第63行附近，设置k的条件和对numpy.argparse的误解（特别是与MATLAB的find相比）似乎在某些情况下将k设置为错误。把代码改成

if q == 1:
    k = 0
elif m < p:
    k = n;
else:
    k = max([0,sum((np.diag(C) <= 1/np.sqrt(2)))])

第79行
我认为S[1，1]应该是S[0，0]（Python 0索引数组）

第83行及以后

这里的神经质矩阵切片似乎是错的。我通过将第83-95行改为：

UT, ST, VT = scipy.linalg.svd(slice_matrix(S,i,j))
    ST = add_zeros(ST,np.zeros([n-k,r-k]))

    if k > 0: 
        print('Zeroing elements of S in row indices > r, to be replaced by ST')
        S[0:k,k:r] = 0
    S[k:n,k:r] = ST
    C[:,j] = np.dot(C[:,j],VT)
    V[:,i] = np.dot(V[:,i],UT)
    Z[:,j] = np.dot(Z[:,j],VT)
    i = np.arange(k,q)
    Q,R = scipy.linalg.qr(C[k:q,k:r])

    C[i,j] = np.diag(diagf(R))
    U[:,k:q] = np.dot(U[:,k:q],Q)

在diagp()中

有两个使用X*Y的矩阵乘法应该是np.dot(X,Y)（注意*是numpy中的element-wise multiplication，而不是矩阵乘法）。

广义SVD是一个很重要的问题，我很惊讶它没有直接在numpy或scipy中实现！实际上，正则SVD和广义特征值分解都是GSVD的特例！您可能需要检查scipy.linalg.eig中可行的Generalized Eigen Vaue解决方案。我对LAPACK的理解是，它是用Fortran 77编写的。有一个C++ implementation of LAPACK，但他们确实说这不是一个完整的实现，我不确定他们是否实现了广义SVD。在python中实现LAPACK和Numerical Recipies的人将成为英雄：-）