Dopasowuję regresję liniową regulowaną przez L1 do bardzo dużego zestawu danych (z n >> p.) Zmienne są znane z góry, ale obserwacje pojawiają się w małych porcjach. Chciałbym utrzymać dopasowanie lasso po każdym kawałku.

Mogę oczywiście dopasować cały model po obejrzeniu każdego nowego zestawu obserwacji. Byłoby to jednak dość nieefektywne, biorąc pod uwagę, że jest dużo danych. Ilość nowych danych, które docierają do każdego kroku, jest bardzo mała, a dopasowanie raczej nie zmieni się znacznie między krokami.

Czy mogę coś zrobić, aby zmniejszyć ogólne obciążenie obliczeniowe?

Patrzyłem na algorytm LARS Efrona i wsp., Ale chętnie rozważę jakąkolwiek inną metodę dopasowania, jeśli można ją zastosować do „gorącego startu” w sposób opisany powyżej.

Uwagi:

Szukam głównie algorytmu, ale wskaźniki do istniejących pakietów oprogramowania, które mogą to zrobić, mogą również okazać się wnikliwe.
Oprócz bieżących trajektorii lasso algorytm jest oczywiście mile widziany, aby zachować inny stan.

Bradley Efron, Trevor Hastie, Iain Johnstone i Robert Tibshirani, Regresja najmniejszego kąta , Annals of Statistics (z dyskusją) (2004) 32 (2), 407--499.

regression lasso

— NPE
źródło

Lasso jest dopasowywane za pomocą LARS (iteracyjny proces, który rozpoczyna się przy początkowym oszacowaniu ). Domyślnie ale możesz to zmienić w większości implementacji (i zastąp go optymalnym który już masz). Najbliższa to , tym mniejsza liczba iteracji LARS, którą będziesz musiał wykonać, aby dostać się do . $\beta^0$ $\beta^0=0_p$ $\beta^*_{old}$ $\beta^*_{old}$ $\beta_{new}^*$ $\beta_{new}^*$

EDYTOWAĆ:

Ze względu na komentarze user2763361dodam więcej szczegółów do mojej oryginalnej odpowiedzi.

Z poniższych komentarzy wynika, że użytkownik2763361 sugeruje uzupełnienie mojej oryginalnej odpowiedzi, aby zmienić ją w taką, która może być używana bezpośrednio (z półek), a jednocześnie jest bardzo wydajna.

Aby wykonać pierwszą część, zilustruję rozwiązanie, które proponuję krok po kroku na przykładzie zabawki. Aby zaspokoić drugą część, zrobię to za pomocą najnowszego, wysokiej jakości wewnętrznego narzędzia do rozwiązywania punktów. Wynika to z faktu, że łatwiej jest uzyskać wysokowydajną implementację rozwiązania, które proponuję, używając biblioteki, która może rozwiązać problem lasso poprzez podejście do punktu wewnętrznego, zamiast próbować zhakować algorytm LARS lub simpleks, aby rozpocząć optymalizację od non-non standardowy punkt początkowy (choć to drugie miejsce jest również możliwe).

Zauważ, że czasami twierdzi się (w starszych książkach), że wewnętrzne podejście do rozwiązywania programów liniowych jest wolniejsze niż podejście simpleksowe i może to być prawdą dawno temu, ale ogólnie nie jest prawdą dzisiaj i na pewno nie jest prawdziwe w przypadku problemów na dużą skalę (dlatego większość profesjonalnych bibliotek lubi cplexalgorytm punktu wewnętrznego), a pytanie dotyczy przynajmniej problemów na dużą skalę. Zauważ też, że wewnętrzny solwer punktowy, którego używam, w pełni obsługuje rzadkie macierze, więc nie sądzę, że będzie duża różnica wydajności z LARS (oryginalną motywacją do korzystania z LARS było to, że wiele popularnych solverów w tym czasie nie radziło sobie dobrze z rzadkimi macierzami i są to charakterystyczne cechy problemu LASSO).

(Bardzo) dobra implementacja algorytmu punktu wewnętrznego o otwartym kodzie źródłowym znajduje się ipoptw COIN-ORbibliotece. Innym powodem będę używał ipoptjest to, że posiada interfejs R ipoptr. Znajdziesz bardziej wyczerpujący przewodnik instalacji tutaj , poniżej podaję standardowe polecenia, aby go zainstalować ubuntu.

w bash:

sudo apt-get install gcc g++ gfortran subversion patch wget
svn co https://projects.coin-or.org/svn/Ipopt/stable/3.11 CoinIpopt
cd ~/CoinIpopt
./configure
make 
make install

Następnie, jako root, w Rdo (zakładam, svnże skopiował plik subversion ~/jak to domyślnie):

install.packages("~/CoinIpopt/Ipopt/contrib/RInterface",repos=NULL,type="source")

Stąd podaję mały przykład (głównie z przykładu zabawki podanego przez Jelmera Ypmę jako część jego Ropakowania ipopt):

library('ipoptr')
# Experiment parameters.
lambda <- 1                                # Level of L1 regularization.
n      <- 100                              # Number of training examples.
e      <- 1                                # Std. dev. in noise of outputs.
beta   <- c( 0, 0, 2, -4, 0, 0, -1, 3 )    # "True" regression coefficients.
# Set the random number generator seed.
ranseed <- 7
set.seed( ranseed )
# CREATE DATA SET.
# Generate the input vectors from the standard normal, and generate the
# responses from the regression with some additional noise. The variable 
# "beta" is the set of true regression coefficients.
m     <- length(beta)                           # Number of features.
A     <- matrix( rnorm(n*m), nrow=n, ncol=m )   # The n x m matrix of examples.
noise <- rnorm(n, sd=e)                         # Noise in outputs.
y     <- A %*% beta + noise                     # The outputs.
# DEFINE LASSO FUNCTIONS
# m, lambda, y, A are all defined in the ipoptr_environment
eval_f <- function(x) {
    # separate x in two parts
    w <- x[  1:m ]          # parameters
    u <- x[ (m+1):(2*m) ]

    return( sum( (y - A %*% w)^2 )/2 + lambda*sum(u) )
}
# ------------------------------------------------------------------
eval_grad_f <- function(x) {
    w <- x[ 1:m ]
    return( c( -t(A) %*% (y - A %*% w),  
               rep(lambda,m) ) )
}
# ------------------------------------------------------------------
eval_g <- function(x) {
    # separate x in two parts
    w <- x[  1:m ]          # parameters
    u <- x[ (m+1):(2*m) ]
    return( c( w + u, u - w ) )
}
eval_jac_g <- function(x) {
    # return a vector of 1 and minus 1, since those are the values of the non-zero elements
    return( c( rep( 1, 2*m ), rep( c(-1,1), m ) ) )
}
# ------------------------------------------------------------------
# rename lambda so it doesn't cause confusion with lambda in auxdata
eval_h <- function( x, obj_factor, hessian_lambda ) {
    H <- t(A) %*% A
    H <- unlist( lapply( 1:m, function(i) { H[i,1:i] } ) )
    return( obj_factor * H )
}
eval_h_structure <- c( lapply( 1:m, function(x) { return( c(1:x) ) } ),
                       lapply( 1:m, function(x) { return( c() ) } ) )
# The starting point.
x0 = c( rep(0, m), 
        rep(1, m) )
# The constraint functions are bounded from below by zero.
constraint_lb = rep(   0, 2*m )
constraint_ub = rep( Inf, 2*m )
ipoptr_opts <- list( "jac_d_constant"   = 'yes',
                     "hessian_constant" = 'yes',
                     "mu_strategy"      = 'adaptive',
                     "max_iter"         = 100,
                     "tol"              = 1e-8 )
# Set up the auxiliary data.
auxdata <- new.env()
auxdata$m <- m
    auxdata$A <- A
auxdata$y <- y
    auxdata$lambda <- lambda
# COMPUTE SOLUTION WITH IPOPT.
# Compute the L1-regularized maximum likelihood estimator.
print( ipoptr( x0=x0, 
               eval_f=eval_f, 
               eval_grad_f=eval_grad_f, 
               eval_g=eval_g, 
               eval_jac_g=eval_jac_g,
               eval_jac_g_structure=eval_jac_g_structure,
               constraint_lb=constraint_lb, 
               constraint_ub=constraint_ub,
               eval_h=eval_h,
               eval_h_structure=eval_h_structure,
               opts=ipoptr_opts,
               ipoptr_environment=auxdata ) )

Chodzi mi o to, że jeśli masz nowe dane, po prostu musisz

aktualizuj ( nie zastępuj) macierz ograniczeń i wektor funkcji celu, aby uwzględnić nowe obserwacje.
zmień punkty początkowe punktu wewnętrznego z

x0 = c (rep (0, m), rep (1, m))

$\beta_{new}$ $\beta_{old}$ $\beta_{init}$ x0

| β_{ja n ja t} - β_{n mi w} |_{1} > | β_{n mi w} - β_{o l re} |_{1} (1)

$|\beta_{init}-\beta_{new}|_1>|\beta_{new}-\beta_{old}|_1\quad(1)$

$\beta_{new}$ $\beta_{old}$ $\beta_{init}$ $n$ $p$

Jeśli chodzi o warunki, w jakich obowiązuje nierówność (1), są to:

$\lambda$ $|\beta_{OLS}|_1$ $p$ $n$
gdy nowe obserwacje nie mają wpływu patologicznego, np. gdy są spójne z procesem stochastycznym, który wygenerował istniejące dane.
gdy rozmiar aktualizacji jest niewielki w stosunku do rozmiaru istniejących danych.

— użytkownik603
źródło

p + 1

$p+1$

β

$\beta$

β^{0}

$\beta^0$

0_{p}

$0_p$

@ aix: czy chcesz zaktualizować całą ścieżkę lasso, czy tylko rozwiązanie? (tj. czy kara za rzadkość jest ustalona?).

— user603

λ

$\lambda$

{\hat{β}}^{l za s s o} = \underset{β}{argmin} {\frac{1}{2)} \sum_{ja = 1}^{N.} (y_{ja} - β_{0} - \sum_{jot = 1}^{p} x_{ja jot} β_{jot})^{2)} + λ \sum_{jot = 1}^{p} | β_{jot} |}

$\hat{\beta}^{lasso} = \underset{\beta}{\operatorname{argmin}} \left \{ {1 \over 2} \sum_{i=1}^N(y_i-\beta_0-\sum_{j=1}^p x_{ij} \beta_j)^2 + \lambda \sum_{j=1}^p |\beta_j| \right \}$

β

$\beta$

Czy znasz jakieś biblioteki, które mogą to zrobić bez konieczności edytowania kodu źródłowego?

— user2763361

Aktualizacja dopasowania lasso o nowe obserwacje

EDYTOWAĆ: