Выпуклая оптимизация в стиле CVX в R?

Мне нужно решить (много раз, для большого количества данных, наряду с множеством других вещей) то, что, по моему мнению, сводится к программа конуса второго порядка. Это можно кратко выразить в 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, пакет SOCP варианты: CLSCOP и DWD, который включает в себя решатель SOCP в качестве дополнения к своему классификатору. Оба имеют схожие, но довольно непрозрачные интерфейсы и немного скудны в документации и примерах, что подводит нас к:

2) Как лучше всего представить указанную выше проблему в формате блока ограничений, используемом этими пакетами? Синтаксис CVX выше скрывает много утомительной возни с дополнительными переменными и т. был бы очень рад ...


r optimization convex-optimization
person walkytalky    schedule 26.04.2013    source источник
comment
Переформулировать проблему (ввести переменные резерва для удаления нормы L ^ 1 и преобразования нормы L ^ 2 в ограничение конуса) относительно легко: заменить норму L ^ 1 для M %*% x на y-z и добавить ограничения y >= 0, z >= 0; замените норму L ^ 2 для A %*% x - b на t и добавьте ограничения t >= sqrt( t(u) %*% u ), u = A %*% x - b. Большинство этих преобразований можно даже автоматизировать, как в случае с CVX. , но для такой простой проблемы, вероятно, не стоит беспокоиться.   -  person Vincent Zoonekynd    schedule 27.04.2013
comment
Однако входной формат DWD::sqlp или CLSOCP::socp недокументирован: вам сообщают, какой аргумент содержит ограничения, но не то, как они закодированы ... Вы можете попытаться связаться с авторами этих пакетов, чтобы получить дополнительную информацию о кодировке ограничения. Вы также можете посмотреть на пакет Rcsdp: он решает более широкий класс проблем (полуопределенные программы), входные данные задокументированы, но преобразование вашей проблемы в желаемую форму не будет таким простым ...   -  person Vincent Zoonekynd    schedule 27.04.2013
comment
@Vincent Спасибо, это полезно. (И это настоящий пост в блоге!) DWD::sqlp выглядит смоделированным на SDPT3, поэтому можно было бы уточнить вводные данные путем сравнения.   -  person walkytalky    schedule 27.04.2013
comment
Я автор CLSOCP. В моей учетной записи github есть документ, в котором более подробно объясняется формат ввода: github.com/jcrudy/CLSOCP/blob/master/CLSOCP/inst/doc/manual.pdf   -  person jcrudy    schedule 13.01.2016

Ответы (3)


arrow_upward
2
arrow_downward

Вы можете найти пакет R CVXfromR полезным. Это позволяет передать задачу оптимизации в CVX из R и вернуть решение в R.

person user2439686    schedule 31.05.2013

arrow_upward
1
arrow_downward

Хорошо, поэтому краткий ответ на этот вопрос: на самом деле нет очень удовлетворительного способа справиться с этим в R. Я закончил тем, что выполнил соответствующие части в Matlab с некоторыми неудобными подделками между двумя системами и, вероятно, в конечном итоге перенесу все в Matlab . (Мой текущий подход предшествует ответу, опубликованному пользователем 2439686. На практике моя проблема была бы столь же неудобной при использовании CVXfromR, но в целом он выглядит как полезный пакет, поэтому я собираюсь принять этот ответ.)

Ресурсов R для этого довольно мало, но сообщение в блоге Винсента Зоонекынд, о котором он упомянул в комментариях, определенно стоит прочитать.

Решатель SOCP, содержащийся в DWD-пакете R, перенесен из решателя Matlab SDPT3 (без частей SDP), поэтому программный интерфейс в основном такой же. Однако, по крайней мере, в моих тестах он работает намного медленнее и в значительной степени решает проблемы с несколькими тысячами варов + ограничений, тогда как SDPT3 решает их за несколько секунд. (Я не проводил полностью честного сравнения по этому поводу, потому что CVX выполняет несколько изящных преобразований по проблеме, чтобы сделать ее более эффективной, в то время как в R я использую довольно наивное определение, но все же.)

Другой возможный вариант, особенно если вы имеете право на академическую лицензию, - это использовать коммерческий решатель Mosek, который имеет пакет интерфейса R Rmosek. Я еще не пробовал это, но в какой-то момент могу попробовать.

(Кстати, другой решатель, связанный с CVX, SeDuMi, полностью не справляется с той же проблемой; авторы CVX не шутят, когда предлагают попробовать несколько решателей. Кроме того, в значительном подмножестве случаев SDTP3 должен переключаться с Cholesky в LU-декомпозицию, которая делает обработку на порядки медленнее, с очень незначительным улучшением цели по сравнению с этапами до LU. Я обнаружил, что для предотвращения этого стоит снизить требуемую точность, но YMMV.)

person walkytalky    schedule 31.05.2013

arrow_upward
1
arrow_downward

Есть новая альтернатива: CVXR, который исходит из те же люди. Существует веб-сайт, paper и проект github.

Дисциплинированное выпуклое программирование, похоже, набирает популярность, наблюдая за cvxpy (Python) и Convex.jl (Джулия), опять же, при поддержке тех же людей.

person sascha    schedule 11.01.2018