我如何矢量化矩阵/输入，以便scipy.optimize.minimize可以使用它？

Posted 2021-04-26

我在为scipy.optimize.minimize转换无约束问题时遇到了一个问题。我想运行方法L-BFGS。

基本问题看起来像这样：

min || A - XY ||

S.T. X，Y

A是给定的矩阵和X $ in R ^ {nxl} $和Y $ in R ^ {lxm}由于scipy只接受Vector输入，我试图将XY解释为更大的变量：Z =（X ，Y）我已经将X和Y的列放在彼此之下。

首先，我试图编写它将转换我的输入向量的函数。对于基本示例，它工作正常（可能是因为矩阵是密集的？idk）

这是我的代码：

import numpy as np
from scipy.optimize import minimize
R=np.array(np.arange(12)).reshape(3, 4)
Z0 = np.array(np.random.random(14))
#X=3x2 = 6
#Y=2x4 = 8
def whatineed (Z):
    return np.linalg.norm(R - np.dot(Z[:6].reshape(3,2),Z[6:].reshape(2,4)))

A = minimize(fun=whatineed, x0=Z0, method='L-BFGS-B', options={'disp': 1})

#print A

上面只是一个（看似？）工作的虚拟代码。它给了我/结果：

x: array([ 1.55308851, -0.50000733,  1.89812395,  1.44382572,  2.24315938, 3.38765876,  0.62668062,  1.23575295,  1.8448253 ,  2.45389762, 1.94655245,  1.83844053,  1.73032859,  1.62221667])

如果我用一个大的它运行它根本不起作用。

RUNNING THE L-BFGS-B CODE

           * * *

Machine precision = 2.220D-16
 N =       377400     M =           10
 This problem is unconstrained.

At X0         0 variables are exactly at the bounds

没有进一步。 R实际上是一个或多或少的稀疏矩阵。但我真的不知道从哪里开始。这是我的功能码吗？这是R的稀疏性吗？都？有什么工作？

更新：解决方案适用于非常小的维度。如果我变得更大，则会发生此错误：

ABNORMAL_TERMINATION_IN_LNSRCH                              

 Line search cannot locate an adequate point after 20 function
  and gradient evaluations.  Previous x, f and g restored.
 Possible causes: 1 error in function or gradient evaluation;
                  2 rounding error dominate computation.

 Cauchy                time 0.000E+00 seconds.
 Subspace minimization time 0.000E+00 seconds.
 Line search           time 0.000E+00 seconds.

 Total User time 0.000E+00 seconds.

正如您在用户时间可以看到的那样，问题非常小并且停止工作。手写（L-BFGS）运行导致完全没有步骤/下降。

答案

如果我们试图解决

min_{B is rank <= k} ‖ A - B ‖_2

众所周知，solution是A的等级k截断SVD。

相反，我们正试图解决

min_{X,Y} ‖ A - XY ‖_2

其中X的形状为n×k，Y的形状为k×n（我使用k因为它比l更容易阅读）

声明：这些是同等的问题。要看到这一点，我们需要证明：

1. XY是秩<= k（具有上述形状）。
2. 任何秩k矩阵都可以写成乘积XY（具有上述形状）。

证明：

1. 第一个事实是Y是秩<= k并且XY的零空间包含在Y的零空间中
2. 第二个是通过写入秩<= k矩阵B = UDV *的SVD并且观察到（UD）具有形状n×k并且V具有形状k×n，其中除了来自分解的前k个奇异值之外我们已经丢弃了所有因为它们保证为零。

实现为了解决OP所述的问题，我们只需要计算A的SVD并将其截断为k级。您可以使用np.linalg.svd或sp.sparse.linalg.svds，具体取决于您的矩阵是否稀疏。对于numpy版本，排名k svd可以计算为：

m,n = 10,20
A = np.random.randn(m,n)

k = 6
u,s,vt = np.linalg.svd(A)

X = u[:,:k]*s[:k]
Y = vt[:k]

print(X.shape, Y.shape)
print(np.linalg.norm(A-X@Y,2))

sp.sparse.linalg.svds的语法几乎相同，除非您能够提前指定所需的等级。