稀疏约束线性最小二乘解算器

Question

This great SO answer指向一个良好稀疏求解器用于Ax=b，但是我有制约x使得x每个元素是一个>=0<=N。

另外，A是巨大的（约2e6x2e6），但非常稀疏，每行<=4元素。

任何想法/建议？我正在寻找像MATLAB的lsqlin之类的东西，但是有很大的稀疏矩阵。

我基本上是试图解决大型bounded variable least squares problem稀疏矩阵：

编辑： 在CVX：

cvx_begin 
    variable x(n) 
    minimize(norm(A*x-b)); 
    subject to 
     x <= N; 
     x >= 0; 
cvx_end 

Answer 1

您试图解决方框约束的最小二乘方。标准稀疏最小二乘算法包括LSQR和最近的LSMR。这些只需要你应用矩阵矢量产品。要添加约束条件，请注意，如果您位于方框的内部（约束都不是“激活的”），那么您将继续使用您选择的任何内点方法。对于所有活动约束，您执行的下一次迭代将停用约束，或限制您沿约束超平面移动。通过对您选择的算法进行一些（概念上相对简单的）适当的修改，您可以实现这些约束。

通常情况下，您可以使用任何凸优化软件包。我个人使用Matlab包CVX解决了这个确切类型的问题，它使用SDPT3/SeDuMi作为后端。 CVX仅仅是这些半定义程序解算器的一个非常方便的包装。

Answer 2

你的矩阵A^TA为正半定，所以你的问题是凸的;确保在设置求解器时充分利用这一点。

大多数QP求解器都是Fortran和/或非自由的;但是我听说过OOQP的好东西（http://www.mcs.anl.gov/research/projects/otc/Tools/OOQP/OoqpRequestForm.html），尽管获得副本有点痛苦。

Answer 3

你的问题类似于一个非负最小二乘问题（NNLS），其可以被配制为

$$ \ MIN_X ||斧-B || _2^2 \ {文本受} x \ ge 0 $$，

其中似乎存在许多算法。

实际上，如果除了原始的非负变量$ x $，您还可以创建其他变量$ x'$并将它们与线性约束连接起来，您可以将问题或多或少地转换为NNLS问题$ x_i + x_i' = N $。这种方法的问题在于，这些附加的线性约束可能不能完全满足最小二乘法解决方案 - 那么尝试用大数量来加权它们可能是适当的。

Answer 4

CVXOPT怎么样？它适用于稀疏矩阵，它似乎是一些锥规划求解的可能帮助：

http://abel.ee.ucla.edu/cvxopt/userguide/coneprog.html#quadratic-cone-programs

这是在文档中的代码的简单修改以上，解决你的问题：

from cvxopt import matrix, solvers 
A = matrix([ [ .3, -.4, -.2, -.4, 1.3 ], 
      [ .6, 1.2, -1.7, .3, -.3 ], 
      [-.3, .0, .6, -1.2, -2.0 ] ]) 
b = matrix([ 1.5, .0, -1.2, -.7, .0]) 
N = 2; 
m, n = A.size 
I = matrix(0.0, (n,n)) 
I[::n+1] = 1.0 
G = matrix([-I, I]) 
h = matrix(n*[0.0] + n*[N]) 
print G 
print h 
dims = {'l': n, 'q': [n+1], 's': []} 
x = solvers.coneqp(A.T*A, -A.T*b, G, h)['x']#, dims)['x'] 
print x 

CVXOPT支持稀疏矩阵，所以它对你很有用。

Answer 5

如果你有Matlab，这是你可以用TFOCS做的事情。这是您用来解决问题的语法：

x = tfocs(smooth_quad, { A, -b }, proj_box(0, N)); 

如果它太大而无法放入内存，则可以将A作为函数句柄传递。

Answer 6

如果重新制定模型为：

分钟
受：

那么它是一个标准的二次规划问题。这是一种常见的模型，可以使用CPLEX或Gurobi等商业求解器轻松解决。 （免责声明：我目前为Gurobi Optimization工作，之前曾为提供CPLEX的ILOG工作）。