秩亏NNLS

计算科学 线性代数 最小二乘 二次规划
2021-12-15 20:02:03

我想找到秩不足的最小二乘问题的最小范数解,受正性约束,例如 其中A是大的、稀疏的、矩形的和秩不足的。

minx x2s.t.Ax=b, x0
A

稳健地解决这个问题的最有效方法是什么?有推荐的代码吗?如果A是满秩的,我可以使用 NNLS;如果没有不等式约束,我可以使用稀疏 QR 分解。上面的优化问题显然是一个 QP,所以在最坏的情况下,我假设内点方法应该可以很好地工作,尽管在过去我很难让例如 IpOpt 与秩不足的等式约束矩阵一起工作。

编辑:典型A大小为 50,000 × 70,000,有 1,000,000 个非零×我更喜欢 C 代码,但 FORTRAN、python 等也很好。出于性能原因,我想尽可能避免使用 Matlab。

3个回答

我们使用Split-Bregman 算法解决了几乎相同的问题(有关实现细节,请参阅this )。

本质上,您将问题转换为其中移除所有负值添加了残差然后重复解决这个问题,每次迭代都会改变minxx2+λAxb+γuxuxbbub

如果您愿意,我可以为您编写整个算法,但简单地查看参考中的伪代码可能更容易。

其他人已经为这个问题提供了两个最有可能的答案,但我将添加一些比较和一种方法来帮助在这两种方法之间做出决定。我建议要么

  1. 一种用于二次规划的原始对偶内点法。对于这种方法,您需要使用其他人开发的库。很大程度上取决于您是否想要/需要开源软件,或者是否可以免费研究(但不是开源)或商业而不是免费的研究。要考虑的一些代码包括 CPLEX、Gurobi、LOQO 和 Clp。

要么

  1. 一阶方法(Augmented Lagrangian、ADMM、split Bregman 等)这些可以自己实现,而无需使用打包的库。

一阶方法的最大优点是它们可以解决非常大的问题,这些问题太大而无法通过原始对偶内点方法解决,但这些方法不会产生非常精确的解决方案并且收敛速度非常慢。

您的问题的大小应该可以合理地使用原始双内点方法解决(假设您有一台具有足够内存和 CPU 的计算机。)如果是这样,您可能会对解决方案的质量感到满意由内点法产生。

您还可以考虑使用更通用的 LP/SOCP/SDP 求解器,例如 Mosek、SeDuMi 或 SDPT3,但由于您的问题实际上是 QP,因此将您的问题重新表述为 SOCP 没有多大意义。

这是一个快速测试,可以帮助您确定原始对偶内点法是否可用于此问题。计算(使用例如 MATLAB)然后计算的 Cholesky 分解(使用合适的稀疏 Cholesky 例程 - MATLAB 内置了这个。)如果 Cholesky 分解可以在您的机器上合理计算,那么原始对偶内点代码应该有足够的内存来解决 QP。如果太大太密集以至于无法实现,那么使用原始对偶方法是不可能的,您需要查看一阶方法。 H=AAT+IHH

可能能够形成 A 的(稀疏)然后使用矩阵和系统形成等效的满秩线性系统。QRARRx=QTb

除此之外,您可以尝试其他 QP 求解器,而不仅仅是 IPOPT。CPLEX 和 Gurobi 具有免费的学术许可证并包含 QP 求解器。CVX 之类的东西,或者 SeDuMi 和 SDPT3 之类的 SDP 求解器也可以工作。

像 LKlevin 提出的分裂增强拉格朗日方法对于特定问题非常有效(例如,Boyd 有一篇很好的论文使用乘法器的交替方向方法进行回归),并且它们可以比内点求解器更快利用问题结构。1