R中的CVX-esque凸优化？答案

【问题标题】：CVX-esque convex optimization in R?R中的CVX-esque凸优化？
【发布时间】：2013-04-26 18:54:07
【问题描述】：

我需要解决（很多时候，对于大量数据，以及一堆其他事情）我认为归结为second order cone program。可以用CVX 简洁地表达如下：

cvx_begin
    variable X(2000);
    expression MX(2000);
    MX = M * X;
    minimize( norm(A * X - b) + gamma * norm(MX, 1) )
  subject to
    X >= 0
    MX((1:500) * 4 - 3) == MX((1:500) * 4 - 2)
    MX((1:500) * 4 - 1) == MX((1:500) * 4)
cvx_end

显示的数据长度和等式约束模式只是来自某些测试数据的任意值，但一般形式将大致相同，有两个客观术语 - 一个最小化错误，另一个鼓励稀疏性 - 以及一个很大的数字对优化变量的变换版本的元素的等式约束（本身被约束为非负）。

这似乎工作得非常好，比我以前的方法要好得多，它使约束变得很糟糕。问题在于，这方面的所有其他事情都发生在 R 中，将它移植到 Matlab 会很麻烦。那么在 R 中这样做是否可行，如果可行，如何实现？

这实际上归结为两个独立的问题：

1) 有没有什么好的 R 资源呢？据我从CRAN task page 得知，SOCP 包选项是CLSCOP 和DWD，其中包括一个 SOCP 求解器作为其分类器的附件。两者都有相似但相当不透明的界面，并且文档和示例有点薄，这让我们：

2) 在这些包使用的约束块格式中表示上述问题的最佳方式是什么？上面的 CVX 语法隐藏了很多繁琐的额外变量等等，我可以看到自己花了周试图做到这一点，所以任何提示或指示可以推动我朝着正确的方向前进会很受欢迎...

【问题讨论】：

重新表述问题（引入松弛变量以去除L^1范数并将L^2范数转化为锥约束）相对容易：将M %*% x的L^1范数替换为@ 987654330@ 并添加约束y >= 0，z >= 0；将A %*% x - b 的L^2 范数替换为t 并添加约束t >= sqrt( t(u) %*% u )、u = A %*% x - b。大多数转换甚至可以是automated，就像CVX一样，但对于这样一个简单的问题，它可能不值得麻烦。
但是，DWD::sqlp 或 CLSOCP::socp 的输入格式未记录：您被告知哪个参数包含约束，但不知道它们是如何编码的......您可以尝试联系作者这些包，以获得有关约束编码的更多信息。您还可以查看Rcsdp 包：它解决了更大类的问题（半确定程序），输入记录在案，但是将您的问题转换为所需的形式并不是那么简单......
@Vincent 谢谢，这很有帮助。（这是一篇相当不错的博文！）DWD::sqlp 看起来是以 SDPT3 为模型的，因此可以通过比较来澄清输入。
我是 CLSOCP 的作者。我的 github 帐户上有一个文档更详细地解释了输入格式：github.com/jcrudy/CLSOCP/blob/master/CLSOCP/inst/doc/manual.pdf

标签： r optimization convex-optimization

【解决方案1】：

您可能会发现 R 包 CVXfromR 很有用。这使您可以将优化问题从 R 传递给 CVX，并将解决方案返回给 R。

【讨论】：

【解决方案2】：

好的，所以对这个问题的简短回答是：在 R 中确实没有非常令人满意的方法来处理这个问题。我最终在 Matlab 中完成了相关部分，在两个系统之间进行了一些尴尬的捏造，并且可能会迁移所有内容最终到Matlab。（我目前的方法早于 user2439686 发布的答案。在实践中，使用 CVXfromR 我的问题同样尴尬，但它看起来确实是一个有用的包，所以我会接受这个答案。）

这方面的R资源非常稀少，但他在cmets中提到的blog post by Vincent Zoonekynd绝对值得一读。

R 包 DWD 中包含的 SOCP 求解器是从 Matlab 求解器SDPT3 移植的（减去 SDP 部分），因此编程接口基本相同。然而，至少在我的测试中，它的运行速度要慢得多，并且几乎会因几千个变量+约束的问题而崩溃，而 SDPT3 在几秒钟内就可以解决它们。（我还没有对此进行完全公平的比较，因为 CVX 对问题进行了一些漂亮的转换以使其更有效，而在 R 中我使用了一个非常幼稚的定义，但仍然如此。）

另一种可能的选择是使用商业的Mosek 求解器，特别是如果您有资格获得学术许可，它有一个 R 接口包Rmosek。我还没有尝试过，但可能会在某个时候试一试。

（顺便说一句，与 CVX 捆绑的另一个求解器 SeDuMi 在同一个问题上完全失败；当 CVX 作者建议尝试多个求解器时，他们并不是在开玩笑。此外，在很大一部分情况下，SDTP3 必须从 Cholesky 切换到 LU 分解，这使得处理速度变慢了几个数量级，与预 LU 步骤相比，目标只有非常微小的改进。我发现降低要求的精度以避免这种情况是值得的，但是 YMMV。）

【讨论】：

【解决方案3】：

有一个新的替代方案：CVXR，它来自同一个人。有website、paper 和github project。

观察到cvxpy（Python）和Convex.jl（Julia），有纪律的凸编程似乎越来越受欢迎，再次得到了同样的人的支持。

【讨论】：