python - 用于 python 广义奇异值分解的 GSVD

Question

MATLAB 有一个 gsvd 函数来执行广义 SVD。自 2013 年以来，我认为关于 github pages 的讨论很多。关于把它放在 scipy 中，有些页面有我可以使用的代码，例如 here这对于像我这样的新手来说 super 复杂(让它运行起来)。

我还找到了 LJWilliams github 页面 implementation .这是不好的，因为在转移到 python 3 时有很多错误。尝试更正简单的错误，例如 assert 和 print。它很快就会变得复杂。

谁能帮我提供 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)

File "/home/eghx/agent18/master_thesis/AMfe/amfe/gsvd.py", line 260, in gsvd U, V, Z, C, S = csd(Q[0:m,:],Q[m:m+n,:])

File "/home/eghx/agent18/master_thesis/AMfe/amfe/gsvd.py", line 107, in csd Q,R = scipy.linalg.qr(S[q:n,m:p])

File "/home/eghx/anaconda3/lib/python3.5/site-packages/scipy/linalg/decomp_qr.py", line 141, in qr overwrite_a=overwrite_a)

File "/home/eghx/anaconda3/lib/python3.5/site-packages/scipy/linalg/decomp_qr.py", line 19, in safecall ret = f(*args, **kwargs)

ValueError: failed to create intent(cache|hide)|optional array-- must have defined dimensions but got (0,)

Answer 1

如果您想使用 github 上的 LJWillams 实现，则有几个错误。但是，要完全理解该技术，我可能会建议您自己动手实现它。我查看了 Octave(相当于 MATLAB 的免费软件)的作用及其 "code is a wrapper to the corresponding Lapack dggsvd and zggsvd routines." ，恕我直言，这是 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) (注意 * 是 element-wise multiplication在 numpy 中，不是矩阵乘法。)