2013-04-26 6 views
6

Мне нужно решить (много раз, для большого количества данных, рядом с кучей других вещей), что я думаю сводится к second order cone program. Это может быть кратко выражено в CVX что-то вроде этого:CVX-esque выпуклая оптимизация в R?

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 выше скрывает много утомительных ошибок с дополнительными переменными и т. Д., И я могу просто увидеть, как я трачу недель, пытаясь понять это правильно, поэтому любые советы или указатели, чтобы подтолкнуть меня в правильном направлении, были бы очень желанными. ..

+1

Переформулировки проблемы (вводя слабину переменных , чтобы удалить L^1 норму и преобразовать L 2 нормы^в ограничение конуса) относительно легко: заменить L^1 норму 'M% *% x' с 'yz' и добавьте ограничения' y> = 0', 'z> = 0'; замените норму L^2 'A% *% x - b' с' t' и добавьте ограничения 't> = sqrt (t (u)% *% u)', 'u = A % *% x - b'. Большинство из этих преобразований может быть даже [автоматизировано] (http://zoonek.free.fr/blosxom/R/2012-06-01_Optimization.html), как с CVX, , но для простой такой проблемы, это, вероятно, не стоит того. –

+1

Однако формат ввода из 'DWD :: sqlp' или' CLSOCP :: socp' недокументирован: вам сообщается, какой аргумент содержит ограничения, но не как они закодированы ... Вы можете попытаться связаться с авторам этих пакетов, , чтобы получить дополнительную информацию о кодировании ограничений. Вы также можете посмотреть пакет 'Rcsdp': он обращается к большему классу проблем (полуопределенные программы), Входы задокументированы, , но преобразование вашей проблемы в нужную форму будет не так просто ... –

+0

@ Vincent Спасибо, это полезно. (И это довольно сообщение в блоге!) 'DWD :: sqlp' выглядит, чтобы быть смоделированным на SDPT3, поэтому было бы возможно прояснить входы путем сравнения. – walkytalky

ответ

2

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

1

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

R ресурсов для этого довольно тонкий на земле, но blog post by Vincent Zoonekynd, о котором он упомянул в комментариях, определенно стоит прочитать.

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

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

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

1

Существует новая альтернатива: CVXR, которая исходит от тех же людей. Существует website, a paper и github project.

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

 Смежные вопросы

  • Нет связанных вопросов^_^