Czy regresja kalenicy jest bezużyteczna w dużych wymiarach ( )? W jaki sposób OLS może się nie dopasowywać?


50

Rozważ dobry stary problem regresji z predyktorami i wielkością próby . Zazwyczaj mądrość jest taka, że ​​estymator OLS będzie nadrzędny i generalnie będzie lepszy niż estymator regresji grzbietu:Standardowe jest stosowanie weryfikacji krzyżowej w celu znalezienia optymalnego parametru regularyzacji . Tutaj używam 10-krotnego CV. Aktualizacja wyjaśnienia: gdy , przez „estymator OLS” rozumiem „estymator OLS o minimalnej normie” podany przezpβ = ( X X + λ I ) - 1 X Y . λ n < p β OLS = ( X X ) + X Y = X + Y .n

β^=(XX+λI)1Xy.
λn<p
β^OLS=(XX)+Xy=X+y.

Muszę zestawu danych z oraz . Wszystkie predyktory są ustandaryzowane i istnieje kilka takich, które (same) mogą wykonać dobrą robotę w przewidywaniu . Jeśli losowo wybiorę małą ish, powiedzmy , liczbę predyktorów, otrzymam rozsądną krzywą CV: duże wartości dają zero R-kwadrat, małe wartości ujemne R-kwadrat (ponieważ nadmiernego dopasowania), a pomiędzy nimi jest maksimum. Dla krzywa wygląda podobnie. Jednak dla znacznie większej niż ta, np. , nie otrzymuję żadnego maksimum: płaskowyż krzywej, co oznacza, że ​​OLS zn=80p>1000p = 50 < n λ λ p = 100 > n p p = 1000 λ 0 λyp=50<nλλp=100>npp=1000λ0 działa tak dobrze, jak regresja kalenicowa z optymalnym .λ

wprowadź opis zdjęcia tutaj

Jak to możliwe i co mówi o moim zbiorze danych? Czy brakuje mi czegoś oczywistego, czy rzeczywiście jest to sprzeczne z intuicją? Jak może być dowolny jakościowa różnica między i ponieważ obie są większe niż ?p=100p=1000n

W jakich warunkach minimalne rozwiązanie OLS dla nie pasuje?n<p


Aktualizacja: W komentarzach było trochę niedowierzania, więc oto odtwarzalny przykład użycia glmnet. Używam Pythona, ale użytkownicy R z łatwością dostosują kod.

%matplotlib notebook

import numpy as np
import pylab as plt
import seaborn as sns; sns.set()

import glmnet_python    # from https://web.stanford.edu/~hastie/glmnet_python/
from cvglmnet import cvglmnet; from cvglmnetPlot import cvglmnetPlot

# 80x1112 data table; first column is y, rest is X. All variables are standardized
mydata = np.loadtxt('../q328630.txt')   # file is here https://pastebin.com/raw/p1cCCYBR
y = mydata[:,:1]
X = mydata[:,1:]

# select p here (try 1000 and 100)
p = 1000

# randomly selecting p variables out of 1111
np.random.seed(42)
X = X[:, np.random.permutation(X.shape[1])[:p]]

fit = cvglmnet(x = X.copy(), y = y.copy(), alpha = 0, standardize = False, intr = False, 
               lambdau=np.array([.0001, .001, .01, .1, 1, 10, 100, 1000, 10000, 100000]))
cvglmnetPlot(fit)
plt.gcf().set_size_inches(6,3)
plt.tight_layout()

wprowadź opis zdjęcia tutaj wprowadź opis zdjęcia tutaj


2
@DJohnson Bez żartów. Zwykle 10-krotne CV, co oznacza, że ​​każdy zestaw treningowy ma n = 72, a każdy zestaw testowy ma n = 8.
ameba mówi Przywróć Monikę

2
To dalekie od zwykłego CV. Biorąc pod uwagę, że jak można oczekiwać czegoś takiego jak wykrywalny wynik?
Mike Hunter

3
@DJohnson Nie rozumiem, dlaczego mówisz, że jest to dalekie od zwykłych. Oto 10-krotne CV.
ameba mówi Przywróć Monikę

2
@ seanv507 Rozumiem. Proponuję zdefiniować „rozwiązanie z lambda = 0” jako „rozwiązanie o minimalnej normie z lambda = 0”. Wydaje
ameba mówi Przywróć Monikę

3
@amoeba: Dziękuję za to pytanie. Jak dotąd był on bardzo pouczający i interesujący.
usεr11852 mówi: Przywróć Monic

Odpowiedzi:


23

Naturalna regularyzacja zachodzi z powodu obecności wielu małych składników w teoretycznej PCA . Te małe elementy są domyślnie wykorzystywane do dopasowania hałasu przy użyciu małych współczynników. Stosując minimalną normę OLS, dopasowujesz hałas do wielu małych niezależnych komponentów, a to ma efekt regulujący równoważny z regulacją Ridge. Ta regularyzacja jest często zbyt silna i można ją zrekompensować za pomocą „antyregularyzacji” znanej jako grań negatywna . W takim przypadku zobaczysz minimum krzywej MSE pojawiającej się dla ujemnych wartości .λxλ

Przez teoretyczne PCA mam na myśli:

Niech to wielowymiarowy rozkład normalny. Istnieje liniowa izometria taka jak gdzie jest przekątna: składowe są niezależne. jest po prostu uzyskiwane po przekątnej .f u = f ( x ) N ( 0 , D ) D u D ΣxN(0,Σ)fu=f(x)N(0,D)DuDΣ

Teraz można zapisać model (izometria liniowa zachowuje iloczyn skalarny). Jeśli napiszesz , model można zapisać . Ponadtostąd metody dopasowania takie jak Ridge lub minimalna norma OLS są doskonale izomorficzne: estymator jest obrazem estymatora .y=β.x+ϵy=f(β).f(x)+ϵγ=f(β)y=γ.u+ϵβ=γy=γ.u+ϵfy=β.x+ϵ

Teoretyczna PCA przekształca nie-niezależne predyktory w niezależne predyktory. Jest to tylko luźno związane z empiryczną PCA, w której stosuje się empiryczną macierz kowariancji (która bardzo różni się od teoretycznej z małą wielkością próby). Teoretyczna PCA nie jest praktycznie obliczalna, ale jest tu używana tylko do interpretacji modelu w ortogonalnej przestrzeni predykcyjnej.

Zobaczmy, co się stanie, gdy do modelu dołączymy wiele niezależnych predyktorów niezależnych od małej wariancji:

Twierdzenie

Normalizacja grzbietu ze współczynnikiem jest równoważna (gdy ) do:λp

  • dodawanie fałszywych niezależnych predyktorów (wyśrodkowanych i identycznie rozmieszczonych) każdy z wariancjąpλp
  • dopasowanie wzbogaconego modelu do estymatora OLS z minimalną normą
  • zachowując tylko parametry dla prawdziwych predyktorów

(szkic) Dowód

Udowodnimy, że funkcje kosztów są asymptotycznie równe. model na prawdziwe i fałszywe predyktory: . Funkcję kosztu Ridge'a (dla prawdziwych predyktorów) można zapisać:y=βx+βx+ϵ

costλ=β2+1λyXβ2

W przypadku zastosowania minimalnej normy OLS odpowiedź jest idealnie dopasowana: wartość błędu wynosi 0. Funkcja kosztu dotyczy tylko normy parametrów. Można go podzielić na prawdziwe i fałszywe parametry:

costλ,p=β2+inf{β2Xβ=yXβ}

W prawidłowym wyrażeniu minimalne rozwiązanie normy podaje:

β=X+(yXβ)

Teraz używa SVD dla :X

X=UΣV

X+=VΣ+U

Widzimy, że norma zasadniczo zależy od pojedynczych wartości które są odwrotnością pojedynczych wartości . Znormalizowana wersja to . Spojrzałem na literaturę, a pojedyncze wartości dużych losowych macierzy są dobrze znane. Dla wystarczająco dużych wartości i , minimalne wartości i maksimum są w przybliżeniu (patrz twierdzenie 1.1 ):βX+XXp/λXpnsminsmax

smin(p/λX)p(1n/p)
smax(p/λX)p(1+n/p)

Ponieważ w przypadku dużych , dąży do 0, można po prostu stwierdzić, że wszystkie wartości są przybliżone przez pojedyncze . A zatem:pn/pp

β1λyXβ

Wreszcie:

costλ,pβ2+1λyXβ2=costλ

Uwaga : nie ma znaczenia, czy utrzymasz współczynniki fałszywych predyktorów w swoim modelu. Wariancja wprowadzona przez to . W ten sposób zwiększasz MSE tylko o współczynnik który i tak dąży do 1. Jakoś nie musisz traktować fałszywych predyktorów inaczej niż rzeczywistych.βxλpβ21pyXβ2npMSE(β)1+n/p

Wróćmy do danych @ amoeba. Po zastosowaniu teoretycznego PCA do (zakładając, że jest normalny), jest przekształcany za pomocą izometrii liniowej w zmienną której składniki są niezależne i sortowane w malejącym porządku wariancji. Problem jest równoważny przekształconemu problemowi .xxuy=βx+ϵy=γu+ϵ

Teraz wyobraź sobie, jak wygląda wariancja komponentów:

wprowadź opis zdjęcia tutaj

Rozważ wiele ostatnich składników, nazwij sumę ich wariancji . Każdy z nich ma wariancję w przybliżeniu równą i jest niezależny. Odgrywają rolę fałszywych predyktorów w twierdzeniu.pλλ/p

Fakt ten jest wyraźniejszy w modelu @ Jonny: tylko pierwsza składowa teoretycznej PCA jest skorelowana (nie jest proporcjonalna ) i ma ogromną zmienność. Wszystkie pozostałe składniki (proporcjonalne do ) mają stosunkowo bardzo małą wariancję (napisz macierz kowariancji i przekątna, aby to zobaczyć) i odgrywają rolę fałszywych predyktorów. Obliczyłem, że uregulowanie tutaj odpowiada (w przybliżeniu) wcześniejszemu na podczas gdy true . To zdecydowanie przesada. Widać to po tym, że końcowy MSE jest znacznie większy niż idealny MSE. Efekt regularyzacji jest zbyt silny.yx¯xix¯N(0,1p2)γ1γ12=1p

Czasami można poprawić tę naturalną regularyzację przez Ridge. Po pierwsze, czasami potrzebujesz w twierdzeniu naprawdę dużym (1000, 10000 ...), aby poważnie konkurować z Ridge, a skończoność jest jak niedokładność. Ale pokazuje również, że Ridge jest dodatkową regularyzacją w stosunku do naturalnie istniejącej regularyzacji domyślnej, a zatem może mieć tylko bardzo niewielki efekt. Czasami ta naturalna regularyzacja jest już zbyt silna i Ridge może nawet nie być poprawą. Co więcej, lepiej jest zastosować antyregulację: Grzbiet o ujemnym współczynniku. Pokazuje MSE dla modelu @ jonny ( ), używając :ppp=1000λR

wprowadź opis zdjęcia tutaj


2
+1 Bardzo miło, dziękuję za napisanie tego. Myślę, że ważne jest wyjaśnienie, że kiedy mówisz „regularyzacja”, masz na myśli (tj. Grzbiet). Można by mieć nadzieję, że lasso lub elastyczna siatka będą się lepiej zachowywać i rzeczywiście tego używają ludzie w sytuacjach . W takich warunkach nikt nie używa czystego grzbietu, a standardową radą jest stosowanie regularności wymuszających rzadkość; więc zachowanie czystego grzbietu może mieć jedynie akademicki interes. Mimo to wydaje się niesamowite, że odkrywamy to tutaj. Dlaczego to nie jest dobrze znane? L2np
ameba mówi Przywróć Monikę

1
Źle zakodowane proporcjonalne do . Przepraszam, że nie czas na coś właściwego. Moim głównym celem było zachowanie minimalnej normy OLS, aby zobaczyć, że różni się ona od twojego przykładu i że „niektóre niezbyt złe regularyzacje” na 40 pierwszych poziomach były gwałtownie lepsze. λσ2
Benoit Sanchez

3
Wydaje mi się, że zrozumiałem tajemnicę: regularyzacja grzbietu ze współczynnikiem jest równoważna z minimalną normą OLS, dodając fałszywe predyktory, każdy z wariancją (asymptotycznie dla dużego ). W twoich danych i modelu Johny'ego dzieje się to bez robienia niczego dzięki komponentom PCA o najniższej wariancji. Teraz potrzebuję czasu, aby znaleźć sposób, aby wyjaśnić to jasno ...λpλ/pp
Benoit Sanchez,

1
Wyjaśniłem mały punkt: współczynniki fałszywych predyktorów nie zwiększają znacznie błędu (patrz uwaga na końcu dowodu). Jest to ważne, ponieważ w twoich danych / jonach są one nieuchronnie przechowywane.
Benoit Sanchez

3
Próbowałem negatywnej Ridge. Nie wierzę, ale działa !!! (i nie tylko model Jonny'ego ...)
Benoit Sanchez

16

Dziękujemy wszystkim za wspaniałą trwającą dyskusję. Sednem sprawy wydaje się być to, że OLS o minimalnej normie skutecznie wykonuje skurcz podobny do regresji grzbietu. Wydaje się, że dzieje się to za każdym razem, gdy . Jak na ironię, dodawanie predyktorów czystego szumu może być nawet użyte jako bardzo dziwna forma lub regularyzacja.pn


Część I. Demonstracja ze sztucznymi danymi i analitycznym CV

@Jonny (+1) wymyślił naprawdę prosty sztuczny przykład, który nieco tu dostosuję. o wielkości i są wytwarzane tak, że wszystkie zmienne są Gaussa z wariancji jednostkowej, a korelacja pomiędzy każdym predyktor oraz odpowiedź jest . Naprawię .Xn×pyρρ=.2

Użyję CV z pominięciem jednego, ponieważ istnieje błąd analityczny dla błędu kwadratu: jest znany jako PRESS , „przewidywana suma kwadratów”. gdzie są resztkami a jest macierz kapelusza pod względem SVD . Pozwala to replikować wyniki @ Jonny'ego bez użycia i bez przeprowadzania weryfikacji krzyżowej (wykreślam stosunek PRASY do sumy kwadratów ):

PRESS=i(ei1Hii)2,
ei
e=yy^=yHy,
H
H=X(XX+λI)1X=US2S2+λU
X=USVglmnety

wprowadź opis zdjęcia tutaj

To podejście analityczne pozwala obliczyć limit w . Po prostu podłączenie do wzoru PRESS nie działa: gdy i , wszystkie reszty są zerowe, a macierz kapelusza jest macierzą identyczności z jedynymi na przekątnej, co oznacza, że ​​ułamki w PRESS równanie jest niezdefiniowane. Ale jeśli obliczymy limit w , wówczas będzie on odpowiadał minimalnemu standardowi rozwiązania OLS z .λ0λ=0n<pλ=0λ0λ=0

Sztuką jest wykonanie rozszerzenia Taylora macierzy kapelusza, gdy : Tutaj przedstawiłem macierz gramów .λ0

H=U11+λ/S2UU(1λ/S2)U=IλUS2U=IλG1.
G=XX=US2U

Jesteśmy prawie gotowi:Lambda została anulowana, więc tutaj mamy wartość graniczną. Narysowałem go dużą czarną kropką na powyższym rysunku (na panelach, gdzie ) i pasuje idealnie.

PRESS=i(λ[G1y]iλGii1)2=i([G1y]iGii1)2.
p>n

Zaktualizuj 21 lutego. Powyższa formuła jest dokładna, ale możemy uzyskać więcej informacji, dokonując dalszych przybliżeń. Wygląda na to, że ma w przybliżeniu jednakowe wartości na przekątnej, nawet jeśli ma bardzo nierówne wartości (prawdopodobnie dlatego, że całkiem dobrze miesza wszystkie wartości własne). Tak więc dla każdego mamy ten którym nawiasy kwadratowe oznaczają uśrednianie. Korzystając z tego przybliżenia, możemy przepisać:To przybliżenie pokazano na powyższym rysunku za pomocą czerwonych otwartych kół.G1SUiGii1S2

PRESSS2S2Uy2.

Czy to będzie ona większa lub mniejsza niż zależy od wartości singularnych . W tej symulacji jest skorelowane z pierwszym komputerem PC z więc jest duży, a wszystkie inne terminy są małe. (W moich prawdziwych danych jest również dobrze przewidywane przez wiodące komputery PC). Teraz, w przypadku , jeśli kolumny są wystarczająco losowe, wówczas wszystkie liczby osobliwe będą raczej blisko siebie (wiersze w przybliżeniu prostokątny). „Główny” terminy2=Uy2SyXU1yypnXU1yzostaną pomnożone przez współczynnik mniejszy niż 1. Warunki pod koniec zostaną pomnożone przez czynniki większe niż 1, ale niewiele większe. Ogółem norma spada. Natomiast w przypadku pojawią się bardzo małe wartości pojedyncze. Po inwersji staną się dużymi czynnikami, które podniosą ogólną normę.pn

[Ten argument jest bardzo falisty; Mam nadzieję, że można to uściślić.]

W ramach kontroli rozsądku, jeśli zmienię kolejność pojedynczych wartości do S = diag(flipud(diag(S)));tego czasu, przewidywane MSE jest powyżej wszędzie na 2 i 3 panelu.1

figure('Position', [100 100 1000 300])
ps = [10, 100, 1000];

for pnum = 1:length(ps)
    rng(42)
    n = 80;
    p = ps(pnum);
    rho = .2;
    y = randn(n,1);
    X = repmat(y, [1 p])*rho + randn(n,p)*sqrt(1-rho^2);

    lambdas = exp(-10:.1:20);
    press = zeros(size(lambdas));
    [U,S,V] = svd(X, 'econ');
    % S = diag(flipud(diag(S)));   % sanity check

    for i = 1:length(lambdas)
        H = U * diag(diag(S).^2./(diag(S).^2 + lambdas(i))) * U';
        e = y - H*y;
        press(i) = sum((e ./ (1-diag(H))).^2);
    end

    subplot(1, length(ps), pnum)
    plot(log(lambdas), press/sum(y.^2))
    hold on
    title(['p = ' num2str(p)])
    plot(xlim, [1 1], 'k--')

    if p > n
        Ginv = U * diag(diag(S).^-2) * U';
        press0 = sum((Ginv*y ./ diag(Ginv)).^2);
        plot(log(lambdas(1)), press0/sum(y.^2), 'ko', 'MarkerFaceColor', [0,0,0]);

        press0approx = sum((diag(diag(S).^-2/mean(diag(S).^-2)) * U' * y).^2);
        plot(log(lambdas(1)), press0approx/sum(y.^2), 'ro');
    end
end

Część druga. Dodanie czystych predyktorów hałasu jako formy regularyzacji

Dobre argumenty przedstawili @Jonny, @Benoit, @Paul, @Dikran i inni, że zwiększenie liczby predyktorów zmniejszy minimalną normę rozwiązania OLS. Rzeczywiście, gdy , każdy nowy predyktor może jedynie obniżyć normę rozwiązania normy minimalnej. Dodanie predyktorów obniży normę, nieco podobnie do tego, w jaki sposób regresja kalenicy karze tę normę.p>n

Czy można to wykorzystać jako strategię regularyzacji? Rozpoczynamy od oraz , a następnie dodajemy czystych predykcyjnych hałasem próbie regularyzacji. Zrobię LOOCV i porównuję go z LOOCV dla grzbietu (obliczonego jak wyżej). Zauważ, że po uzyskaniu na predyktorach „obcinam” to w ponieważ interesują mnie tylko oryginalne predyktory.n=80p=40qβ^p+qp

wprowadź opis zdjęcia tutaj

TO DZIAŁA!!!

W rzeczywistości nie trzeba „obcinać” wersji beta; nawet jeśli użyję pełnej wersji beta i pełnych predyktorów , mogę uzyskać dobrą wydajność (linia przerywana na prawym wykresie podrzędnym). Myślę, że to naśladuje moje rzeczywiste dane w pytaniu: tylko nieliczne predyktory naprawdę przewidują , większość z nich to czysty szum i służą one jako regularyzacja. W tym systemie dodatkowa regularyzacja grzbietu wcale nie pomaga.p+qy

rng(42)
n = 80;
p = 40;
rho = .2;
y = randn(n,1);
X = repmat(y, [1 p])*rho + randn(n,p)*sqrt(1-rho^2);

lambdas = exp(-10:.1:20);
press = zeros(size(lambdas));
[U,S,V] = svd(X, 'econ');

for i = 1:length(lambdas)
    H = U * diag(diag(S).^2./(diag(S).^2 + lambdas(i))) * U';
    e = y - H*y;
    press(i) = sum((e ./ (1-diag(H))).^2);
end

figure('Position', [100 100 1000 300])
subplot(121)
plot(log(lambdas), press/sum(y.^2))
hold on
xlabel('Ridge penalty (log)')
plot(xlim, [1 1], 'k--')
title('Ridge regression (n=80, p=40)')
ylim([0 2])

ps = [0 20 40 60 80 100 200 300 400 500 1000];
error = zeros(n, length(ps));
error_trunc = zeros(n, length(ps));
for fold = 1:n
    indtrain = setdiff(1:n, fold);
    for pi = 1:length(ps)
        XX = [X randn(n,ps(pi))];
        if size(XX,2) < size(XX,1)
            beta = XX(indtrain,:) \ y(indtrain,:);
        else
            beta = pinv(XX(indtrain,:)) * y(indtrain,:);
        end
        error(fold, pi) = y(fold) - XX(fold,:) * beta;
        error_trunc(fold, pi) = y(fold) - XX(fold,1:size(X,2)) * beta(1:size(X,2));
    end
end

subplot(122)
hold on
plot(ps, sum(error.^2)/sum(y.^2), 'k.--')
plot(ps, sum(error_trunc.^2)/sum(y.^2), '.-')
legend({'Entire beta', 'Truncated beta'}, 'AutoUpdate','off')
legend boxoff
xlabel('Number of extra predictors')
title('Extra pure noise predictors')
plot(xlim, [1 1], 'k--')
ylim([0 2])

@MartijnWeterings W tym eksperymencie zaczynam od n = 80 ip = 40. Gdy całkowita liczba predyktorów (p + q) zbliża się do n = 80, problem staje się źle uwarunkowany, a rozwiązanie OLS drastycznie się przepełnia. Błąd jest ogromny w okolicach q = 40. Gdy tylko p + q> n uruchamia się ograniczenie „minimalnej normy” i błąd zaczyna maleć, ale powrót do miejsca, w którym był z q = 0, zajmuje trochę czasu. Zdarza się to około q = 70, tj. P + q = 130. Następnie błąd zmniejsza się jeszcze bardziej, a ta część wykresu jest podobna do wykresu regresji grzbietu. Czy jest sens?
ameba mówi Przywróć Monikę

@MartijnWeterings W sprawie pierwszego komentarza: jesteśmy na tej samej stronie. W sprawie drugiego komentarza: w moim pytaniu nie obcinam wersji beta, zgadza się. Ale tak naprawdę, jeśli nie obetnę wersji beta w mojej symulacji (użyj y(fold) - XX(fold,:) * betazamiast XX(fold,1:size(X,2)) * beta(1:size(X,2))), to wyniki nie zmieniają się zbytnio. Chyba powinienem dodać to do mojej odpowiedzi. Myślę, że moje oryginalne dane pokazują tego rodzaju zachowanie.
ameba mówi Przywróć Monikę

(1/2): Nadal pracuję nad wszystkimi komentarzami i kodem, aby zrozumieć, ale przychodzi mi do głowy myśl: czy istnieje związek między tym zjawiskiem, które obserwujemy, a związkiem między regresją grzbietu a efektami losowymi?
Ryan Simmons,

(2/2): Według odpowiedzi Randela tutaj ( stats.stackexchange.com/questions/122062/... ) widzimy oszacowanie równoważne między efektami losowymi a regresją grzbietu, gdzie lambda jest równa stosunkowi reszt do wariancji efekt losowy. Tutaj, zgodnie z odpowiedzią Benoit Sanchez, widzimy, że regresja grzbietu jest równoważna dodaniu dowolnej liczby fałszywych niezależnych predyktorów, z których każda ma wariancję równą funkcji lambda i liczbie parametrów. Wydaje mi się, że istnieje związek koncepcyjny.
Ryan Simmons,

@amoeba to był błąd. dodanie wektora skalowanego y do macierzy X reguluje nieco, ale nie to samo, co regresja grzbietu lub wektory szumu. Sprawia jednak, że zastanawiam się, co się stanie, gdy odejmiemy trochę od każdego x, aby każda zmienna była lekko ujemnie skorelowana (lub mniej dodatnio) z wektorem y. Ma to na celu wykonanie pewnej „negatywnej” regularyzacji. Żeby „cofnąć” regularyzację 1000 wektorów (w pewnym momencie może stać się ona zbyt duża, jak widać, gdy szczytowy / optymalny współczynnik regularyzacji jest teraz prawie poza zasięgiem). y
Sextus Empiricus

15

Oto sztuczna sytuacja, w której ma to miejsce. Załóżmy, że każda zmienna predykcyjna jest kopią zmiennej docelowej z dużą ilością zastosowanego szumu gaussowskiego. Najlepszy możliwy model jest średnią wszystkich zmiennych predykcyjnych.

library(glmnet)
set.seed(1846)
noise <- 10
N <- 80
num.vars <- 100
target <- runif(N,-1,1)
training.data <- matrix(nrow = N, ncol = num.vars)
for(i in 1:num.vars){
  training.data[,i] <- target + rnorm(N,0,noise)
}
plot(cv.glmnet(training.data, target, alpha = 0,
               lambda = exp(seq(-10, 10, by = 0.1))))

MSE dla różnych lambda ze 100 predyktorami

100 zmiennych zachowuje się w „normalny” sposób: pewna dodatnia wartość lambda minimalizuje błąd poza próbą.

Ale zwiększ liczbę zmiennych w powyższym kodzie do 1000, a oto nowa ścieżka MSE. (Rozszerzyłem log (Lambda) = -100, aby się przekonać.

MSE dla różnych lambda z 1000 predyktorami

To, co myślę, że się dzieje

Przy dopasowywaniu wielu parametrów o niskiej regularyzacji współczynniki są losowo rozmieszczane wokół ich prawdziwej wartości z dużą zmiennością.

Ponieważ liczba predyktorów staje się bardzo duża, „średni błąd” zmierza w kierunku zera, a lepiej jest pozwolić, aby współczynniki spadły tam, gdzie mogą, i zsumować wszystko, niż przesunąć je w kierunku 0.

Jestem pewien, że ta sytuacja, w której prawdziwe przewidywanie jest średnią wszystkich predyktorów, nie jest jedynym momentem, kiedy to się dzieje, ale nie wiem, jak zacząć tutaj wskazywać największy niezbędny warunek.

EDYTOWAĆ:

Zachowanie „płaskie” dla bardzo niskiej lambda zawsze się zdarza, ponieważ rozwiązanie jest zbieżne z rozwiązaniem OLS o minimalnej normie. Podobnie krzywa będzie płaska dla bardzo wysokiej lambda, ponieważ rozwiązanie zbiega się do 0. Nie będzie minimum, jeśli jedno z tych dwóch rozwiązań jest optymalne.

Dlaczego rozwiązanie OLS o minimalnej normie jest (porównywalnie) dobre w tym przypadku? Myślę, że jest to związane z poniższym zachowaniem, które uważam za bardzo sprzeczne z intuicją, ale refleksja ma wiele sensu.

max.beta.random <- function(num.vars){
  num.vars <- round(num.vars)
  set.seed(1846)
  noise <- 10
  N <- 80
  target <- runif(N,-1,1)
  training.data <- matrix(nrow = N, ncol = num.vars)

  for(i in 1:num.vars){
    training.data[,i] <- rnorm(N,0,noise)
  }
  udv <- svd(training.data)

  U <- udv$u
  S <- diag(udv$d)
  V <- udv$v

  beta.hat <- V %*% solve(S) %*% t(U) %*% target

  max(abs(beta.hat))
}


curve(Vectorize(max.beta.random)(x), from = 10, to = 1000, n = 50,
      xlab = "Number of Predictors", y = "Max Magnitude of Coefficients")

abline(v = 80)

Wykres maksymalnej wielkości współczynników wraz ze wzrostem liczby predyktorów

Przy losowo generowanych predyktorach niezwiązanych z odpowiedzią, gdy p wzrasta, współczynniki stają się większe, ale gdy p jest znacznie większe niż N, kurczą się w kierunku zera. Tak dzieje się również w moim przykładzie. Tak bardzo luźno, nieregularne rozwiązania tych problemów nie wymagają skurczu, ponieważ są już bardzo małe!

Dzieje się tak z trywialnego powodu. może być dokładnie wyrażony jako liniowa kombinacja kolumn . jest wektorem minimalnych norm współczynników. W miarę dodawania kolejnych kolumn norma musi się zmniejszać lub pozostać stała, ponieważ możliwą kombinacją liniową jest utrzymanie poprzednich współczynników na tym samym poziomie i ustawienie nowych współczynników na .yXβ^β^0


1
(+1). Zjawisko to wydaje się zatem mieć miejsce, gdy predyktory są skorelowane. Nie oznacza to formalnie, że krzywa błędu nie ma minimum dla dodatniej , ani że granica przy 0 nie jest duża. Oznacza to po prostu, że krzywa ma tendencję do spłaszczania się i że próg określający, jak mała musi być, aby regularyzacja przestała działać, zmierza w kierunku 0 dla dużego . Tutaj ten próg przekracza limit obliczeniowy, ale odpowiedź Firebuga sugeruje, że zawsze może istnieć. λλp
Benoit Sanchez

1
Dlaczego potrzebujesz glmnetaktualizacji? Jeśli potrzebujesz tylko rozwiązania OLS o minimalnej normie, istnieje bezpośrednia formuła (patrz druga formuła w moim pytaniu), a jeśli obliczymy SVD dla wówczas ta formuła staje się po prostu . Prawdopodobnie jest też funkcja w R, która oblicza to rozwiązanie, ale tak naprawdę nie znam R :)X=USVβ^=VS1Uy
amoeba mówi Przywróć Monikę

2
Myśląc o tym trochę, wcale nie jest to zaskakujące. może być dokładnie wyrażony jako liniowa kombinacja wektorów w . to wektor współczynników o najmniejszej normie. Po dodaniu wektora norma musi się zmniejszyć lub pozostać na tym samym rozmiarze, ponieważ możesz zachować te same współczynniki i ustawić nowe na zero. yXβ^β^
Jonny Lomond

3
Przykład Jonny'ego jest dobry, ponieważ został już przeanalizowany: patrz estymator Jamesa-Steina . Szacując średnią ustalonego wektora o wymiarze 3 lub większym, zawsze możemy poprawić po prostym uśrednieniu poprzez odchylenie do zera, co jest mniej więcej tym, co robi regresja grzbietu. Zastanawiam się, czy może poprawa jest zbyt niewielka, aby można ją było zaobserwować w tym przypadku? θ
Paul

3
Dobrze znany jest fakt, że regresja grzbietu jest równoważna dodaniu dodatkowych „fałszywych” próbek do zestawu danych, przy czym każda próbka ma wartość w jednej funkcji i zera wszędzie indziej, a wszystkie odpowiadające odpowiedzi są równe zero. (Przepisanie funkcji kosztu RR w tej formie jest banalne). Zastanawiam się teraz, czy istnieje sposób na dodanie dodatkowych funkcji (np. Czysty szum?), Które miałyby podobny efekt. Oczywiście dodanie dodatkowych funkcji oznacza, że zwiększy swoją wymiarowość, ale można było spojrzeć na jego wartości tylko w „oryginalnych” predyktorach. @Paulpλβ^
ameba mówi Przywróć Monikę

6

Postanowiłem więc uruchomić zagnieżdżoną weryfikację krzyżową przy użyciu specjalistycznego mlrpakietu w R, aby zobaczyć, co tak naprawdę wynika z podejścia do modelowania.

Kod (uruchomienie zwykłego notebooka zajmuje kilka minut)

library(mlr)
daf = read.csv("https://pastebin.com/raw/p1cCCYBR", sep = " ", header = FALSE)

tsk = list(
  tsk1110 = makeRegrTask(id = "tsk1110", data = daf, target = colnames(daf)[1]),
  tsk500 = makeRegrTask(id = "tsk500", data = daf[, c(1,sample(ncol(daf)-1, 500)+1)], target = colnames(daf)[1]),
  tsk100 = makeRegrTask(id = "tsk100", data = daf[, c(1,sample(ncol(daf)-1, 100)+1)], target = colnames(daf)[1]),
  tsk50 = makeRegrTask(id = "tsk50", data = daf[, c(1,sample(ncol(daf)-1, 50)+1)], target = colnames(daf)[1]),
  tsk10 = makeRegrTask(id = "tsk10", data = daf[, c(1,sample(ncol(daf)-1, 10)+1)], target = colnames(daf)[1])
)

rdesc = makeResampleDesc("CV", iters = 10)
msrs = list(mse, rsq)
configureMlr(on.par.without.desc = "quiet")
bm3 = benchmark(learners = list(
    makeLearner("regr.cvglmnet", alpha = 0, lambda = c(0, exp(seq(-10, 10, length.out = 150))),
    makeLearner("regr.glmnet", alpha = 0, lambda = c(0, exp(seq(-10, 10, length.out = 150))), s = 151)
    ), tasks = tsk, resamplings = rdesc, measures = msrs)

Wyniki

getBMRAggrPerformances(bm3, as.df = TRUE)
#   task.id    learner.id mse.test.mean rsq.test.mean
#1    tsk10 regr.cvglmnet     1.0308055  -0.224534550
#2    tsk10   regr.glmnet     1.3685799  -0.669473387
#3   tsk100 regr.cvglmnet     0.7996823   0.031731316
#4   tsk100   regr.glmnet     1.3092522  -0.656879104
#5  tsk1110 regr.cvglmnet     0.8236786   0.009315037
#6  tsk1110   regr.glmnet     0.6866745   0.117540454
#7    tsk50 regr.cvglmnet     1.0348319  -0.188568886
#8    tsk50   regr.glmnet     2.5468091  -2.423461744
#9   tsk500 regr.cvglmnet     0.7210185   0.173851634
#10  tsk500   regr.glmnet     0.6171841   0.296530437

Robią w zasadzie to samo we wszystkich zadaniach.

A co z optymalnymi lambdami?

sapply(lapply(getBMRModels(bm3, task.ids = "tsk1110")[[1]][[1]], "[[", 2), "[[", "lambda.min")
# [1] 4.539993e-05 4.539993e-05 2.442908e-01 1.398738e+00 4.539993e-05
# [6] 0.000000e+00 4.539993e-05 3.195187e-01 2.793841e-01 4.539993e-05

Zauważ, że lambdy są już przekształcone. Niektóre fold nawet wybrał minimalną lambda .λ=0

Zrobiłem trochę więcej zabawy glmneti nie odkryłem, że tam jest wybierana minimalna lambda. Czek:

EDYTOWAĆ:

Po komentarzach ameby stało się jasne, że ścieżka regularyzacji jest ważnym krokiem w glmnetoszacowaniu, więc kod ją teraz odzwierciedla. W ten sposób większość rozbieżności zniknęła.

cvfit = cv.glmnet(x = x, y = y, alpha = 0, lambda = exp(seq(-10, 10, length.out = 150)))
plot(cvfit)

wprowadź opis zdjęcia tutaj

Wniosek

Zasadniczo naprawdę poprawia dopasowanie ( edytuj: ale nie za bardzo! ).λ>0

Jak to możliwe i co mówi o moim zbiorze danych? Czy brakuje mi czegoś oczywistego, czy rzeczywiście jest to sprzeczne z intuicją?

Prawdopodobnie jesteśmy bliżej prawdziwego rozkładu ustawienia danych na małą wartość większą niż zero. Nie ma w tym jednak nic sprzecznego z intuicją.λ

Edycja: Pamiętaj jednak, że ścieżka regulowania grzbietu korzysta z wcześniejszych oszacowań parametrów, gdy dzwonimy glmnet, ale to jest poza moją wiedzą. Jeśli ustawimy naprawdę niski poziom lambdaizolacji, prawdopodobnie obniży to wydajność.

EDYCJA: Wybór lambda mówi coś więcej o twoich danych. Ponieważ większe lambdy zmniejszają wydajność, oznacza to, że w twoim modelu występują preferencyjne, tj. Większe współczynniki, ponieważ duże lambdas zmniejszają wszystkie współczynniki do zera. Chociaż oznacza, że ​​efektywne stopnie swobody w twoim modelu są mniejsze niż pozorne stopnie swobody, .λ0p

Jak może istnieć jakakolwiek różnica jakościowa między p = 100 ip = 1000, biorąc pod uwagę, że oba są większe niż n?

p=1000 niezmiennie zawiera co najmniej tyle samo informacji, a nawet więcej niż .p=100


Komentarze

Wygląda na to, że otrzymujesz małe minimum dla jakiejś niezerowej lambda (patrzę na twoją figurę), ale krzywa nadal jest naprawdę bardzo płaska na lewo od niej. Tak więc moje główne pytanie pozostaje, dlaczego λ → 0 nie zauważalnie się przenika. Nie widzę tu jeszcze odpowiedzi. Czy spodziewasz się, że będzie to zjawisko ogólne? Czyli dla dowolnych danych z n≪p, lambda = 0 będzie [prawie] tak dobre, jak optymalna lambda? Czy może jest coś specjalnego w tych danych? Jeśli spojrzysz wyżej w komentarzach, zobaczysz, że wiele osób nawet mi nie wierzyło, że to możliwe.

Myślę, że łączysz wydajność sprawdzania poprawności z wydajnością testową i takie porównanie nie jest uzasadnione.

Edycja: zauważ jednak, że kiedy ustawimy lambdana 0 po uruchomieniu całej ścieżki normalizacji, wydajność jako taka nie spada, dlatego ścieżka normalizacji jest kluczem do zrozumienia, co się dzieje!

Nie do końca rozumiem twoją ostatnią linię. Spójrz na wynik cv.glmnet dla p = 100. Będzie miał bardzo inny kształt. Co więc wpływa na ten kształt (asymptota po lewej vs. brak asymptoty), gdy p = 100 lub p = 1000?

Porównajmy ścieżki regularyzacji dla obu:

fit1000 = glmnet(x, y, alpha = 0, lambda = exp(seq(-10,10, length.out = 1001)))
fit100 = glmnet(x[, sample(1000, 100)], y, alpha = 0, lambda = exp(seq(-10,10, length.out = 1001)))
plot(fit1000, "lambda")

wprowadź opis zdjęcia tutaj

x11()
plot(fit100, "lambda")

wprowadź opis zdjęcia tutaj

Staje się jasne, daje większe współczynniki przy wzroście , mimo że ma mniejsze współczynniki dla asymptotycznie grzbietu OLS, po lewej stronie obu wykresów. Zasadniczo więc po lewej stronie wykresu pasuje, a to prawdopodobnie tłumaczy różnicę w zachowaniu między nimi.p=1000λp=100

trudniej jest przeregulować, ponieważ chociaż Ridge zmniejsza współczynniki do zera, nigdy nie osiągają zera. Oznacza to, że moc predykcyjna modelu jest dzielona między wiele innych komponentów, co ułatwia przewidywanie wokół średniej, zamiast być porywanym przez hałas.p=1000


+1 Dziękujemy za wykonanie tych eksperymentów! Wygląda na to, że otrzymujesz małe minimum dla jakiejś niezerowej lambda (patrzę na twoją figurę), ale krzywa nadal jest naprawdę bardzo płaska na lewo od niej. Tak więc moje główne pytanie pozostaje, dlaczego nie zauważalnie się przejmuje. Nie widzę tu jeszcze odpowiedzi. Czy spodziewasz się, że będzie to zjawisko ogólne? Czyli dla dowolnych danych z , lambda = 0 będzie [prawie] tak dobre, jak optymalna lambda? Czy może jest coś specjalnego w tych danych? Jeśli spojrzysz wyżej w komentarzach, zobaczysz, że wiele osób nawet mi nie wierzyło, że to możliwe. n strλ0np
ameba mówi Przywróć Monikę

Nie do końca rozumiem twoją ostatnią linię. Spójrz na cv.glmnetwynik dla p = 100. Będzie miał bardzo inny kształt. Co więc wpływa na ten kształt (asymptota po lewej vs. brak asymptoty), gdy p = 100 lub p = 1000?
ameba mówi Przywróć Monikę

Czy wiesz, czy mlrwybiera lambda.minlub lambda.1se(w cv.glmnetterminologii)?
ameba mówi Przywróć Monikę

@amoeba lambda.min. Jest też regr.cvglmnetuczeń, który prawdopodobnie pozwala wybrać inne reguły.
Firebug

Dzięki. Szczerze mówiąc, nie rozumiem wyników twojego testu porównawczego 1e-100. Np. Dla p = 1100 daje MSE = 1,45. Ale tutaj nie ma strojenia hiperparametrów w wewnętrznej pętli, więc w zasadzie nie trzeba wcale wewnętrznej pętli CV. Oznacza to, że wynik powinien być taki sam, jak w przypadku nie zagnieżdżonego CV przy lambda = 1e-100. Ale widzimy na pierwszej liczbie, że MSE wynosi około 0,7. To nie ma dla mnie sensu.
Ameba mówi Przywróć Monikę

5

W jaki sposób (minimalna norma) OLS może nie pasować?

W skrócie:

Parametry eksperymentalne, które korelują z (nieznanymi) parametrami w prawdziwym modelu, będą bardziej prawdopodobne do oszacowania za pomocą wysokich wartości w procedurze dopasowania minimalnej normy OLS. Wynika to z tego, że będą pasować do „modelu + hałasu”, podczas gdy inne parametry będą pasować tylko do „hałasu” (a zatem będą pasować do większej części modelu o niższej wartości współczynnika i prawdopodobnie będą miały wysoką wartość w minimalnej normie OLS).

Efekt ten zmniejszy ilość przeuczeń w procedurze dopasowania minimalnej normy OLS. Efekt jest bardziej wyraźny, jeśli dostępnych jest więcej parametrów, ponieważ wówczas bardziej prawdopodobne staje się włączenie większej części „prawdziwego modelu” do oszacowania.

Dłuższa część:
(nie jestem pewien, co tu umieścić, ponieważ kwestia nie jest dla mnie całkowicie jasna lub nie wiem z jaką precyzją odpowiedź wymaga odpowiedzi na pytanie)

Poniżej znajduje się przykład, który można łatwo zbudować i pokazuje problem. Efekt nie jest tak dziwny, a przykłady są łatwe do zrobienia.

  • Wziąłem funkcji sin (ponieważ są one prostopadłe) jako zmiennep=200
  • utworzył model losowy z pomiarami. n=50
    • Model jest zbudowany tylko z zmiennych, więc 190 z 200 zmiennych tworzy możliwość generowania nadmiernego dopasowania.tm=10
    • współczynniki modelu są ustalane losowo

W tym przykładzie obserwujemy, że występuje pewne nadmierne dopasowanie, ale współczynniki parametrów należących do prawdziwego modelu mają wyższą wartość. Zatem R ^ 2 może mieć pewną wartość dodatnią.

Poniższy obraz (i kod do jego wygenerowania) pokazują, że nadmierne dopasowanie jest ograniczone. Kropki, które odnoszą się do modelu szacowania 200 parametrów. Czerwone kropki odnoszą się do tych parametrów, które są również obecne w „prawdziwym modelu” i widzimy, że mają wyższą wartość. Zatem istnieje pewien stopień zbliżenia się do modelu rzeczywistego i uzyskania R ^ 2 powyżej 0.

  • Zauważ, że użyłem modelu ze zmiennymi ortogonalnymi (funkcje sinusoidalne). Jeśli parametry są skorelowane, mogą wystąpić w modelu o stosunkowo bardzo wysokim współczynniku i stać się bardziej karane w minimalnej normie OLS.
  • Zauważ, że „zmienne ortogonalne” nie są ortogonalne, gdy weźmiemy pod uwagę dane. Wewnętrzny iloczyn wynosi tylko zero, gdy zintegrujemy całą przestrzeń a nie, gdy mamy tylko kilka próbek . Konsekwencją jest to, że nawet przy zerowym hałasu wystąpi nadmierna montażu (i ^ 2 Wartość R wydaje się zależeć od wielu czynników, oprócz szumu. Oczywiście istnieje zależność i , lecz również istotne jest to, jak wiele zmienne w prawdziwym modelu i ile z nich w dopasowanym modelu).x x n psin(ax)sin(bx)xxnp

przykład ograniczenia nadmiernego dopasowania

library(MASS)

par(mar=c(5.1, 4.1, 9.1, 4.1), xpd=TRUE)

p <- 200       
l <- 24000
n <- 50
tm <- 10

# generate i sinus vectors as possible parameters
t <- c(1:l)
xm <- sapply(c(0:(p-1)), FUN = function(x) sin(x*t/l*2*pi))

# generate random model by selecting only tm parameters
sel <- sample(1:p, tm)
coef <- rnorm(tm, 2, 0.5)

# generate random data xv and yv with n samples
xv <- sample(t, n)
yv <- xm[xv, sel] %*% coef + rnorm(n, 0, 0.1)

# generate model
M <- ginv(t(xm[xv,]) %*% xm[xv,])

Bsol <- M %*% t(xm[xv,]) %*% yv
ysol <- xm[xv,] %*% Bsol

# plotting comparision of model with true model
plot(1:p, Bsol, ylim=c(min(Bsol,coef),max(Bsol,coef)))
points(sel, Bsol[sel], col=1, bg=2, pch=21)
points(sel,coef,pch=3,col=2)

title("comparing overfitted model (circles) with true model (crosses)",line=5)
legend(0,max(coef,Bsol)+0.55,c("all 100 estimated coefficients","the 10 estimated coefficients corresponding to true model","true coefficient values"),pch=c(21,21,3),pt.bg=c(0,2,0),col=c(1,1,2))

Skrócona technika beta w odniesieniu do regresji kalenicy

Przekształciłem kod Pythona z Amoeby w R i połączyłem oba wykresy razem. Dla każdego minimalnego oszacowania normy OLS z dodanymi zmiennymi szumu dopasowuję oszacowanie regresji grzbietu z tą samą (w przybliżeniu) l_2 dla wektora . βl2β

  • Wygląda na to, że model obciętego szumu robi to samo (oblicza tylko trochę wolniej, a może nieco rzadziej).
  • Jednak bez obcięcia efekt jest znacznie mniej silny.
  • Ta zgodność między dodawaniem parametrów a karą kalenicową niekoniecznie jest najsilniejszym mechanizmem za brakiem nadmiernego dopasowania. Widać to zwłaszcza na krzywej 1000p (na zdjęciu pytania) zbliżającej się do prawie 0,3, podczas gdy inne krzywe, z innym p, nie osiągają tego poziomu, bez względu na parametr regresji grzbietu. Dodatkowe parametry, w tym praktycznym przypadku, nie są takie same jak przesunięcie parametru grzbietu (i sądzę, że dzieje się tak, ponieważ dodatkowe parametry stworzą lepszy, bardziej kompletny model).

  • Parametry hałasu zmniejszają z jednej strony normę (podobnie jak regresja kalenicy), ale także wprowadzają dodatkowy hałas. Benoit Sanchez pokazuje, że na granicy, dodając wiele wielu parametrów hałasu przy mniejszym odchyleniu, ostatecznie stanie się taki sam jak regresja kalenicy (rosnąca liczba parametrów hałasu znosi się nawzajem). Ale jednocześnie wymaga znacznie więcej obliczeń (jeśli zwiększymy odchylenie szumu, aby umożliwić użycie mniejszych parametrów i przyspieszenie obliczeń, różnica stanie się większa).

Rho = 0,2 porównanie przyciętego hałasu z regresją grzbietu

Rho = 0,4 porównanie przyciętego hałasu z regresją grzbietu

Rho = 0,2 zwiększając wariancję parametrów hałasu do 2 porównanie przyciętego hałasu z regresją grzbietu

przykład kodu

# prepare the data
set.seed(42)
n = 80
p = 40
rho = .2
y = rnorm(n,0,1)
X = matrix(rep(y,p), ncol = p)*rho + rnorm(n*p,0,1)*(1-rho^2)

# range of variables to add
ps = c(0, 5, 10, 15, 20, 40, 45, 50, 55, 60, 70, 80, 100, 125, 150, 175, 200, 300, 400, 500, 1000)
#ps = c(0, 5, 10, 15, 20, 40, 60, 80, 100, 150, 200, 300) #,500,1000)

# variables to store output (the sse)
error   = matrix(0,nrow=n, ncol=length(ps))
error_t = matrix(0,nrow=n, ncol=length(ps))
error_s = matrix(0,nrow=n, ncol=length(ps))

# adding a progression bar
pb <- txtProgressBar(min = 0, max = n, style = 3)

# training set by leaving out measurement 1, repeat n times 
for (fold in 1:n) {
    indtrain = c(1:n)[-fold]

    # ridge regression
    beta_s <- glmnet(X[indtrain,],y[indtrain],alpha=0,lambda = 10^c(seq(-4,2,by=0.01)))$beta
    # calculate l2-norm to compare with adding variables
    l2_bs <- colSums(beta_s^2)

    for (pi in 1:length(ps)) {
        XX = cbind(X, matrix(rnorm(n*ps[pi],0,1), nrow=80))
        XXt = XX[indtrain,]

        if (p+ps[pi] < n) {
            beta = solve(t(XXt) %*% (XXt)) %*% t(XXt) %*% y[indtrain]
        }
        else {
            beta = ginv(t(XXt) %*% (XXt)) %*% t(XXt) %*% y[indtrain]
        }

        # pickout comparable ridge regression with the same l2 norm      
        l2_b <- sum(beta[1:p]^2)
        beta_shrink <- beta_s[,which.min((l2_b-l2_bs)^2)] 

        # compute errors
        error[fold, pi] = y[fold] - XX[fold,1:p] %*% beta[1:p]
        error_t[fold, pi] = y[fold] - XX[fold,] %*% beta[]
        error_s[fold, pi] = y[fold] - XX[fold,1:p] %*% beta_shrink[]
    }
    setTxtProgressBar(pb, fold) # update progression bar
}

# plotting
plot(ps,colSums(error^2)/sum(y^2) , 
     ylim = c(0,2),
     xlab ="Number of extra predictors",
     ylab ="relative sum of squared error")
lines(ps,colSums(error^2)/sum(y^2))
points(ps,colSums(error_t^2)/sum(y^2),col=2)
lines(ps,colSums(error_t^2)/sum(y^2),col=2)
points(ps,colSums(error_s^2)/sum(y^2),col=4)
lines(ps,colSums(error_s^2)/sum(y^2),col=4)

title('Extra pure noise predictors')

legend(200,2,c("complete model with p + extra predictors",
               "truncated model with p + extra predictors",
               "ridge regression with similar l2-norm",
               "idealized model uniform beta with 1/p/rho"),
       pch=c(1,1,1,NA), col=c(2,1,4,1),lt=c(1,1,1,2))

# idealized model (if we put all beta to 1/rho/p we should theoretically have a reasonable good model)
error_op <- rep(0,n)
for (fold in 1:n) {
  beta = rep(1/rho/p,p)
    error_op[fold] = y[fold] - X[fold,] %*% beta
}
id <- sum(error_op^2)/sum(y^2)
lines(range(ps),rep(id,2),lty=2)

1
(+1) Dzięki. Myślę, że intuicyjny argument na początku twojej odpowiedzi ma sens.
ameba mówi Przywróć Monikę

1

Jeśli znasz operatory liniowe, możesz polubić moją odpowiedź jako najbardziej bezpośrednią drogę do zrozumienia tego zjawiska: dlaczego regresja norm nie zawiedzie wprost? Powodem jest to, że twoim problemem ( ) jest źle postawiony problem odwrotny, a pseudo-odwrotny jest jednym ze sposobów jego rozwiązania. Uregulowanie jest jednak poprawą.np

Ten artykuł jest prawdopodobnie najbardziej zwięzłym i stosownym wyjaśnieniem: Lorenzo Rosasco i in., Uczenie się, regularyzacja i źle odwrócone problemy . Ustawili problem regresji jako uczenie się, patrz równanie 3, gdzie liczba parametrów przekracza liczbę obserwacji: gdzie jest operatorem liniowym w przestrzeni Hilberta, a - zaszumione dane.A g δ

Ax=gδ,
Agδ

Oczywiście jest to źle postawiony odwrotny problem. Możesz więc rozwiązać to za pomocą SVD lub odwrotności Moore-Penrose'a, co faktycznie uczyniłoby rozwiązanie najmniej normalne. Tak więc powinno nie być zaskoczeniem, że najmniej rozwiązanie norma nie zawodzi wprost.

Jeśli jednak posłuchasz artykułu, zobaczysz, że regresja grzbietu byłaby poprawą w stosunku do powyższego. Poprawa jest naprawdę lepszym zachowaniem estymatora, ponieważ rozwiązanie Moore-Penrose'a niekoniecznie jest ograniczone.

AKTUALIZACJA

Uświadomiłem sobie, że nie wyjaśniłem, że źle postawione problemy prowadzą do nadmiernego dopasowania. Oto cytat z pracy Gábor A, Banga JR. Solidne i wydajne oszacowanie parametrów w modelach dynamicznych układów biologicznych . BMC Systems Biology. 2015; 9: 74. doi: 10.1186 / s12918-015-0219-2:

Niewłaściwe uwarunkowanie tych problemów wynika zazwyczaj z (i) modeli o dużej liczbie parametrów (nad parametryzacja), (ii) niedoboru danych eksperymentalnych i (iii) znacznych błędów pomiaru [19, 40]. W rezultacie często uzyskujemy przeregulowanie takich modeli kinetycznych, tj. Modele skalibrowane z rozsądnym dopasowaniem do dostępnych danych, ale słabymi możliwościami generalizacji (niska wartość predykcyjna)

Zatem mój argument można sformułować w następujący sposób:

  • źle postawione problemy prowadzą do przeuczenia
  • Przypadek (n <p) jest wyjątkowo źle postawionym odwrotnym problemem
  • Psudo-inverse Moore-Penrose'a (lub inne narzędzia, takie jak SVD), które określasz w pytaniu jako , rozwiązuje źle postawiony problemX+
  • dlatego dba o przebudowę przynajmniej w pewnym stopniu i nie powinno dziwić, że nie zawiedzie całkowicie, w przeciwieństwie do zwykłego OLS

Ponownie, regularyzacja jest nadal bardziej niezawodnym rozwiązaniem.


1
(+1) Dzięki, ale nie bardzo rozumiem, jak ważny jest ten artykuł. Przyjrzę się temu jutro bardziej szczegółowo. Gdzie dokładnie mówią, że rozwiązanie OLS z minimalną normą nie będzie pasować lub że minimalny wymóg normy można uznać za regularyzację?
ameba mówi Przywróć Monikę

1
Omówmy, kiedy czytasz gazetę. Nie twierdzą, że odwrotność psudo jest regularyzacją. Mówią, że jest to rozwiązanie źle postawionego problemu. Mówię o tym, że nadmierne dopasowanie wynika z niewłaściwej natury problemu, więc zajmując się tym drugim, zajmujesz się tym pierwszym, choć nie tak dobrze, jak z regularyzacją.
Aksakal

1
Myślę, że zagadką nie jest to, że minimalne rozwiązanie norm nie poprawia do pewnego stopnia nadmiernego dopasowania, ale to, że dodanie większej regularności nie poprawi jeszcze bardziej. Także dlaczego minimalne normy są bardziej efektywne, gdy liczba funkcji rośnie. Moją intuicją jest to, że problemy z większą liczbą parametrów wymagają większej regularyzacji (w przeciwnym razie wszystkie rzeczy są takie same), a nie mniej. Jest to naprawdę interesujący problem i może pomóc wyjaśnić, dlaczego np. Nawet nieuregulowane sieci neuronowe nie są nadmiernie dopasowane, jak można się spodziewać.
Dikran Torbacz

1
@Dikran W rzeczywistości inne formy lub regularyzacje mogą nadal poprawiać wydajność: np. Mogę poprawić wydajność (w porównaniu z minimalną normą OLS) z regresją głównych składników lub elastyczną siatką. Po prostu uregulowanie kalenicy staje się bezużyteczne. Analogia do sieci neuronowych to fascynująca myśl, która nie przyszła mi do głowy. Co mi nie myśleć o tym niedawno jednak to, że nic dziwnego, nikt nie rozumie dlaczego skomplikowane głębokie rzeczy uczenia się, jak naprawdę działa partii normalizacji, biorąc pod uwagę, że nawet regresja liniowa grzbiet ze statystyk 101 może być tak zagadkową :-)
ameba mówi dozbrojenie Monica

2
To nie jest główne pytanie, ale myślę, że ta znakomita seria pytań, odpowiedzi i komentarzy została odsunięta na bok od wprowadzenia weryfikacji krzyżowej. W przypadku tej dyskusji o wiele łatwiej byłoby obliczyć predyktor liniowy populacji z , które wykorzystano do symulacji danych, i obliczyć MSE dowolnego estymatora predyktora liniowego. I widziałem przypadek, w którym dla nie mogłem znaleźć optymalnie skorygowanego AIC dla regresji grzbietu ( funkcja pakietu R ). Ale muszę ponownie uruchomić to, używając prawdziwego predyktora liniowego jako złotego standardu. n < < sβn<<prmsols
Frank Harrell,
Korzystając z naszej strony potwierdzasz, że przeczytałeś(-aś) i rozumiesz nasze zasady używania plików cookie i zasady ochrony prywatności.
Licensed under cc by-sa 3.0 with attribution required.