GSVD for python 广义奇异值分解

Question

MATLAB 有一个 gsvd 函数来执行广义 SVD。自 2013 年以来，我认为github 页面上已经有很多关于将其放入 scipy 的讨论，并且有些页面有我可以使用的代码，例如这里，这对于像我这样的新手来说非常复杂（要让它运行）。

我还找到了 LJWilliams github 页面及其实现。这是没有好处的，因为当转移到 python 3 时有很多错误。尝试纠正简单的错误，例如断言和打印。它很快就会变得复杂。

有人可以帮助我使用 python 的 gsvd 代码或向我展示如何使用在线代码吗？

另外，一旦打印和断言语句被更正，这就是我通过 LJWilliams 实现得到的结果。代码看起来很复杂，我不确定花时间在上面是最好的事情！另外，有些人在同一个 github 页面上报告了问题，我不确定这些问题是否已修复或已连接。

n = 10
m = 6
p = 6

A = np.random.rand(m,n)
B = np.random.rand(p,n)
gsvd(A,B)

文件“/home/eghx/agent18/master_thesis/AMfe/amfe/gsvd.py”，第 260 行，在 gsvd U、V、Z、C、S = csd(Q[0:m,:],Q[m: m+n，：]）

文件“/home/eghx/agent18/master_thesis/AMfe/amfe/gsvd.py”，第 107 行，csd Q,R = scipy.linalg.qr(S[q:n,m:p])

文件“/home/eghx/anaconda3/lib/python3.5/site-packages/scipy/linalg/decomp_qr.py”，第 141 行，qr overwrite_a=overwrite_a）

文件“/home/eghx/anaconda3/lib/python3.5/site-packages/scipy/linalg/decomp_qr.py”，第 19 行，在 safecall ret = f(*args, **kwargs)

ValueError：无法创建意图（缓存|隐藏）|可选数组 - 必须已定义尺寸但得到（0，）

Answer 1

如果您想使用 github 上的 LJWillams 实现进行工作，则存在一些错误。然而，为了充分理解该技术，我可能建议您尝试自己实现它。我查了一下 Octave（MATLAB 等价的免费软件）的功能，他们的“代码是相应 Lapack dggsvd 和 zggsvd 例程的包装器”。，恕我直言，这就是 scipy 应该做的事情。

我将发布我发现的错误，但我不会发布完整工作顺序的代码，因为考虑到翻译它的受版权保护的 MATLAB 实现，我不确定这在版权方面如何。

警告：我不是广义 SVD 方面的专家，并且只是从调试的角度来解决这个问题，而不是从底层算法是否正确的角度来解决这个问题。我已经在你的原始随机数组上进行了这项工作，并且测试用例已经存在于 Python 文件中。

虫子

环境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 行开始

这里的 numpy 矩阵切片似乎是错误的。我通过将第 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，而不是矩阵乘法。）