Разреженный линейный решатель наименьших квадратов - PullRequest
Новая
       9

Разреженный линейный решатель наименьших квадратов

8 голосов
/ 10 июня 2010

Этот отличный ответ SO указывает на хороший разреженный решатель для Ax=b, но у меня есть ограничения на x, такие что каждый элемент в x равен >=0 и <=N ,

Кроме того, A - это огромный (около 2e6x2e6), но очень редкий элемент <=4 в строке.

Есть идеи / рекомендации? Я ищу что-то вроде MATLAB lsqlin, но с огромными разреженными матрицами.

Я, по сути, пытаюсь решить проблему наименьших квадратов в больших масштабах на разреженных матрицах:

alt text

EDIT: В CVX : 

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

Ответы [ 6 ]

3 голосов
/ 10 июня 2010

Вы пытаетесь решить наименьших квадратов с ограничениями коробки. Стандартные алгоритмы разреженных наименьших квадратов включают LSQR и, в последнее время, LSMR. Это требует только от вас применения матрично-векторных произведений. Чтобы добавить ограничения, имейте в виду, что если вы находитесь внутри поля (ни одно из ограничений не является «активным»), то вы переходите к любому выбранному вами методу внутренней точки. Для всех активных ограничений следующая итерация, которую вы выполняете, будет либо деактивировать ограничение, либо заставит вас двигаться по гиперплоскости ограничения. С некоторыми (концептуально относительно простыми) подходящими модификациями выбранного вами алгоритма вы можете реализовать эти ограничения.

Как правило, вы можете использовать любой выпуклый пакет оптимизации. Я лично решил эту проблему именно с помощью пакета CVX Matlab, который использует SDPT3 / SeDuMi для бэкэнда. CVX - просто очень удобная оболочка для этих полуопределенных программных решателей.

3 голосов
/ 10 июня 2010

Ваша проблема похожа на неотрицательную задачу наименьших квадратов (NNLS), которую можно сформулировать как

$$ \ min_x || Ax-b || _2 ^ 2 \text {подчиненный} x \ ge 0 $$,

, для которого, кажется, существует много алгоритмов.

На самом деле, ваша проблема может быть более или менее преобразована в проблему NNLS, если,В дополнение к вашим исходным неотрицательным переменным $ x $ вы создаете дополнительные переменные $ x '$ и связываете их с линейными ограничениями $ x_i + x_i' = N $.Проблема с этим подходом состоит в том, что эти дополнительные линейные ограничения могут не удовлетворяться точно в решении наименьших квадратов - тогда может быть целесообразно попытаться взвесить их с большим числом.

1 голос
/ 07 октября 2011

Если вы переформулируете свою модель как:

мин image
при условии:
image
image

, тогда это стандартная квадратичнаяпроблема программирования.Это распространенный тип модели, который может быть легко решен с помощью коммерческого решателя, такого как CPLEX или Gurobi.(Отказ от ответственности: в настоящее время я работаю в Gurobi Optimization и ранее работал в ILOG, которая предоставила CPLEX).

1 голос
/ 06 октября 2011

Если у вас есть Matlab, это то, что вы можете сделать с TFOCS .Синтаксис, который вы использовали бы для решения проблемы: 

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

Вы можете передать A в качестве дескриптора функции, если он слишком велик, чтобы поместиться в память.

1 голос
/ 10 июня 2010

Ваша матрица A ^ TA положительно-полуопределена, поэтому ваша проблема выпуклая;убедитесь, что воспользовались этим при настройке вашего солвера.

Большинство решателей QP идут в Фортран и / или не являются бесплатными;однако я слышал хорошие отзывы об OOQP (http://www.mcs.anl.gov/research/projects/otc/Tools/OOQP/OoqpRequestForm.html),, хотя получение копии немного затруднительно.

0 голосов
/ 20 июня 2010

Как насчет 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 поддерживает разреженные матрицы, так что это будет полезно для вас.

...