Sparso vincolata lineare dei minimi quadrati risolutore

Question

Questa grande risposta SO punta a un buon risolutore sparse per Ax=b, ma ho vincoli x tale che ogni elemento in x è >=0 un <=N.

Inoltre, A è enorme (circa 2e6x2e6) ma molto scarsa con elementi <=4 per riga.

Tutte le idee / suggerimenti? Sto cercando qualcosa di simile a MATLAB di lsqlin ma con enormi matrici sparse.

sto essenzialmente cercando di risolvere la grande scala variabile delimitata almeno piazze problema su matrici sparse:

Modifica In CVX :

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

Answer 1

Si sta cercando di risolvere i minimi quadrati con vincoli box. Standard minimi sparse algoritmi quadrati includono LSQR e, più recentemente, LSMR. Questi richiedono solo di applicare i prodotti matrice-vettore. Per aggiungere i vincoli, si rende conto che se si è all'interno della scatola (nessuno dei vincoli sono "attivi"), poi si procede con qualunque interni punto metodo scelto. Per tutti i vincoli attivi, la prossima iterazione si esegue sarà o disattivare il vincolo, oppure limitare a muoversi lungo all'iperpiano vincolo. Con alcuni (concettualmente relativamente semplici) opportune modifiche all'algoritmo che si sceglie, è possibile implementare questi vincoli.

In generale tuttavia, è possibile utilizzare qualsiasi pacchetto di ottimizzazione convessa. Ho personalmente risolto questo tipo esatto di problema utilizzando il Matlab pacchetto CVX, che utilizza SDPT3 / SeDuMi per un back-end. CVX è solo un molto conveniente involucro intorno a questi risolutori programma semidefinita.

Answer 2

Il tuo problema è simile ad un non negativi minimi quadrati problema (NNLS), che può essere formulato come

$$ \ min_x || Ax-b || _2 ^ 2 \ text {} oggetto di x \ ge 0 $$,

per i quali non sembra esistere molti algoritmi.

In realtà, è problema può essere più o meno trasformato in un problema di NNLS, se, oltre alle variabili non negative originali $ x $ si crea ulteriori variabili $ x '$ e il loro collegamento con vincoli lineari $ x_i + x_i' = N $. Il problema di questo approccio è che questi vincoli lineari aggiuntivi non potrebbero essere soddisfatte esattamente nella soluzione dei minimi quadrati -. Potrebbe essere opportuno quindi per cercare di peso loro con un gran numero

Answer 3

Se si riformula il proprio modello come:

min
soggetto a:

, allora è un problema di programmazione quadratica standard. Questo è un tipo comune di modello che può essere facilmente risolto con un risolutore commerciale come CPLEX o Gurobi. . (Disclaimer: io attualmente lavoro per Gurobi Ottimizzazione e in precedenza ha lavorato per ILOG, che ha fornito CPLEX)

Answer 4

Il tuo matrice A ^ T A è positiva semi-definita, in modo che il problema è convesso; Assicurati di approfittare di questo prima di impostare il risolutore.

La maggior parte di go-to risolutori QP sono in Fortran e / o non-free; tuttavia Ho sentito parlare bene di OOQP ( http: //www.mcs.anl.gov/research/projects/otc/Tools/OOQP/OoqpRequestForm.html ), anche se è un po 'un dolore che ottiene una copia.

Answer 5

Come su CVXOPT? Funziona con matrici sparse, e sembra che alcuni dei risolutori di programmazione cono possono aiutare:

http: //abel.ee.ucla .edu / cvxopt / userguide / coneprog.html #-cono-programmi quadratica

Questa è una semplice modifica del codice nel documento di cui sopra, per risolvere il 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

supporto CVXOPT matrici sparse, in modo da essere utile per voi.

Answer 6

Se si dispone di Matlab, questo è qualcosa che si può fare con TFOCS . Questa è la sintassi si usa per risolvere il problema:

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

È possibile passare una come impugnatura funzione se è troppo grande per entrare nella memoria.