lineales restringida por mínimos cuadrados Escasos solucionador

Question

This great SO answer puntos para un buen solucionador de escasa Ax=b, pero tengo limitaciones en x de tal manera que cada elemento de x es >=0 un <=N.

Además, A es enorme (alrededor 2e6x2e6) pero muy escaso con <=4 elementos por fila.

¿Alguna idea/recomendación? Estoy buscando algo como MATLAB lsqlin pero con enormes matrices dispersas.

estoy esencialmente tratando de resolver el gran escala bounded variable least squares problem en matrices dispersas:

EDIT: En CVX:

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

Answer 1

Está tratando de resolver mínimos cuadrados con limitaciones de caja. Los algoritmos estándar de mínimos cuadrados dispersos incluyen LSQR y, más recientemente, LSMR. Estos solo requieren que aplique productos de matriz vectorial. Para agregar las restricciones, tenga en cuenta que si se encuentra en el interior de la caja (ninguna de las restricciones está "activa"), proceda con el método de punto interior que elija. Para todas las restricciones activas, la siguiente iteración que realice desactivará la restricción o lo obligará a moverse a lo largo del hiperplano de restricción. Con algunas modificaciones (conceptualmente relativamente simples) adecuadas al algoritmo que elija, puede implementar estas restricciones.

Generalmente, sin embargo, puede utilizar cualquier paquete de optimización convexo. Personalmente he resuelto este tipo exacto de problema usando el paquete de Matlab CVX, que usa SDPT3/SeDuMi para un back-end. CVX es simplemente un envoltorio muy conveniente para estos solucionadores de programas semipinos.

Answer 2

Su matriz A^TA es positivo semi-definido, entonces su problema es convexo; asegúrese de aprovecharlo al configurar su solucionador.

La mayoría de los solucionadores de QP se encuentran en Fortran y/o no son gratuitos; sin embargo, he escuchado cosas buenas sobre OOQP (http://www.mcs.anl.gov/research/projects/otc/Tools/OOQP/OoqpRequestForm.html), aunque es un poco doloroso obtener una copia.

Answer 3

Su problema es similar a un mínimos cuadrados no negativos problema (NNLS), que se pueden formular como

$$ \ MIN_X || Ax-b || _2^2 \ text {sujetos a} x \ ge 0 $$,

para el que parece que existen muchos algoritmos.

En realidad, su problema puede convertirse más o menos en un problema NNLS, si, además de sus variables no negativas originales $ x $, crea variables adicionales $ x '$ y las vincula con restricciones lineales $ x_i + x_i' = N $. El problema con este enfoque es que estas restricciones lineales adicionales podrían no satisfacerse exactamente en la solución de mínimos cuadrados; podría ser apropiado entonces intentar ponderarlas con un número grande.

Answer 4

¿Qué tal CVXOPT?Se trabaja con matrices dispersas, y parece que algunos de los solucionadores de programación de cono pueden ser útiles:

http://abel.ee.ucla.edu/cvxopt/userguide/coneprog.html#quadratic-cone-programs

Esta es una simple modificación del código en el documento anterior, para resolver su problema:

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 admite matrices dispersas, por lo que te será útil.

Answer 5

Si tiene Matlab, esto es algo que puede hacer con TFOCS. Esta es la sintaxis que se utiliza para resolver el problema:

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

puede pasar una función como un mango si es demasiado grande para caber en la memoria.

Answer 6

Si a reformular su modelo como:

min
sujeta a:

entonces es un problema de programación cuadrática estándar. Este es un tipo de modelo común que se puede resolver fácilmente con un solucionador comercial como CPLEX o Gurobi. (Descargo de responsabilidad: actualmente trabajo para Gurobi Optimization y anteriormente trabajé para ILOG, que proporcionó CPLEX).