Rzadkie ograniczane liniowych najmniejszych kwadratów Solver

Question

This great SO answer wskazuje na dobrą rozrzedzony solver dla Ax=b, ale mam ograniczenia na x taki sposób, że każdy element w x jest >=0<=N.

Ponadto, A jest ogromny (około 2e6x2e6), ale bardzo rzadki z <=4 elementów w rzędzie.

Wszelkie pomysły/rekomendacje? Szukam czegoś takiego jak MATLAB's lsqlin, ale z ogromnymi rzadkimi matrycami.

jestem zasadniczo próbuje rozwiązać dużą skalę bounded variable least squares problem na macierzach rzadkich:

EDIT: W CVX:

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

Answer 1

Próbujesz rozwiązać najmniejsze kwadraty z ograniczeniami pola. Standardowe algorytmy rzadkich najmniejszych kwadratów obejmują LSQR, a ostatnio LSMR. Wymagają tylko zastosowania produktów matrycowo-wektorowych. Aby dodać ograniczenia, zdaj sobie sprawę, że jeśli znajdujesz się we wnętrzu pudełka (żadne z ograniczeń nie jest "aktywne"), to kontynuuj wybraną metodę punktu wewnętrznego. Dla wszystkich aktywnych więzów, następna iteracja, którą wykonasz, albo dezaktywuje ograniczenie, albo zmusza cię do poruszania się wzdłuż hiperpłaszczyzny więzów. W przypadku niektórych (z założenia względnie prostych) odpowiednich modyfikacji wybranego algorytmu można wdrożyć te ograniczenia.

Generalnie można jednak użyć dowolnego pakietu optymalizacji wypukłej. Osobiście rozwiązałem ten dokładny typ problemu za pomocą pakietu Matlab CVX, który używa SDPT3/SeDuMi dla zaplecza. CVX jest jedynie bardzo wygodnym opakowaniem wokół tych półpełnych programów solverów.

Answer 2

Twój macierzy A^TA jest dodatnia na wpół określony, więc twój problem jest wypukły; pamiętaj, aby skorzystać z tego podczas konfigurowania solwera.

Większość osób, które zdecydowały się na rozwiązanie QP, znajdują się w języku Fortran i/lub non-free; jednak słyszałem dobre rzeczy o OOQP (http://www.mcs.anl.gov/research/projects/otc/Tools/OOQP/OoqpRequestForm.html), choć trochę bólu dostaje kopię.

Answer 3

Problemu jest podobny do nieujemne najmniejszych kwadratów problemu (NNLS), który może być formułowany jako

$$ \ min_x || Ax b || _2^2 \ tekście {zastrzeżeniem} x \ ge 0 $$,

dla których istnieje wiele algorytmów.

W zasadzie problem może zostać mniej lub bardziej przekonwertowany na problem NNLS, jeśli oprócz oryginalnych zmiennych nieujemnych $ x $ tworzysz dodatkowe zmienne $ x '$ i łącz je z ograniczeniami liniowymi $ x_i + x_i' = N $. Problem z tym podejściem polega na tym, że te dodatkowe ograniczenia liniowe mogą nie być spełnione dokładnie w rozwiązaniu najmniejszych kwadratów - może wtedy być właściwe, aby je obciążyć dużą liczbą.

Answer 4

Co z CVXOPT?Współpracuje z nielicznych matryc, i wydaje się, że niektóre z rozwiązują programowania stożek może pomóc:

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

Jest to prosta modyfikacja kodu w dokumencie powyżej, aby rozwiązać swój problem:

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 obsługuje rzadkie macierze, więc może być użyteczne dla ciebie.

Answer 5

Jeśli masz Matlaba, to jest coś, co możesz zrobić z TFOCS. Jest to składnia byłoby użyć do rozwiązania problemu:

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

można przekazać jako uchwyt funkcji, jeśli jest zbyt duży, aby zmieścić się w pamięci.

Answer 6

Jeśli przeformułować swój model jako:

min
zastrzeżeniem:

to jest to norma kwadratowa problemem programowania. Jest to popularny typ modelu, który można łatwo rozwiązać za pomocą komercyjnego rozwiązania, takiego jak CPLEX lub Gurobi. (Zrzeczenie się: obecnie pracuję dla optymalizacji Gurobi i wcześniej pracowałem dla ILOG, który dostarczył CPLEX).