9

这个很棒的 SO 答案指向了一个很好的稀疏求解器Ax=b，但是我有这样的限制，即x每个元素x都是.>=0<=N

此外，A它很大（大约 2e6x2e6），但<=4每行元素非常稀疏。

有什么想法/建议吗？我正在寻找类似 MATLABlsqlin但具有巨大稀疏矩阵的东西。

我本质上是在尝试解决稀疏矩阵上的大规模有界变量最小二乘问题：

替代文字

编辑： 在CVX中：

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

6 回答 6

3

您的问题类似于非负最小二乘问题 (NNLS)，可以表述为

$$\min_x ||Ax-b||_2^2 \text{ 服从 } x \ge 0$$,

似乎存在许多算法。

实际上，如果除了原始的非负变量 $x$ 之外，您还创建了其他变量 $x'$ 并将它们与线性约束 $x_i+x_i'=N$ 链接，则您的问题或多或少可以转换为 NNLS 问题。这种方法的问题在于，这些额外的线性约束在最小二乘解中可能无法完全满足——然后尝试对它们进行大量加权可能是合适的。

于 2010-06-10T08:32:33.280 回答
3

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

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

于 2010-06-10T09:15:24.210 回答
2

如果您将模型重新表述为：

分钟
受制于：

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

于 2011-10-06T22:23:17.867 回答
1

你的矩阵 A^TA 是半正定的，所以你的问题是凸的；在设置求解器时一定要利用这一点。

大多数首选 QP 求解器都是 Fortran 和/或非免费的；但是，我听说过有关 OOQP 的好消息（http://www.mcs.anl.gov/research/projects/otc/Tools/OOQP/OoqpRequestForm.html），尽管获得副本有点痛苦。

于 2010-06-10T05:04:42.047 回答
1

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 支持稀疏矩阵，因此对您很有用。

于 2010-06-20T16:58:17.637 回答
1

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

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

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

于 2011-10-06T04:48:24.867 回答