Dziwne korelacje w wynikach SVD danych losowych; czy mają wyjaśnienie matematyczne, czy jest to błąd LAPACK?

21

Obserwuję bardzo dziwne zachowanie w wynikach losowych danych SVD, które mogę odtworzyć zarówno w Matlabie, jak i R. Wygląda to na jakiś problem numeryczny w bibliotece LAPACK; czy to jest

Rysuję próbek z wymiarowego Gaussa z zerową średnią i kowariancją tożsamości: . I zamontować je w dane matrycy . (Opcjonalnie mogę wyśrodkować lub nie, nie ma to wpływu na następujące.) Następnie wykonuję rozkład wartości osobliwych (SVD), aby uzyskać . Weźmy jakieś dwa szczególne elementy , np i , i zapytać, jaka jest korelacja między nimi całej inny czerpie z . Spodziewałbym się tego, gdyby liczbak = 2 X N ( 0 , I ) 1000 × 2 X X X = U S VU U 11 U 22 X N r e pn=1000k=2XN.(0,ja)1000×2)XXX=US.V.UU11U22XN.rmip losowań jest dość duży, wówczas wszystkie takie korelacje powinny wynosić około zera (tj. Korelacje populacji powinny wynosić zero, a korelacje próbek będą małe).

Obserwuję jednak dziwnie silne korelacje (około ) między , , i i tylko między tymi elementami. Wszystkie pozostałe pary elementów mają korelacje wokół zera, zgodnie z oczekiwaniami. Oto jak wygląda macierz korelacji dla „górnych” elementów (pierwsze elementów pierwszej kolumny, a następnie pierwsze elementów drugiej kolumny):U 11 U 12 U 21 U 22±0.2U11U12U21U22U20U101010

Dziwne korelacje SVD

Zwróć uwagę na dziwnie wysokie wartości w lewym górnym rogu każdej ćwiartki.

To komentarz tego @ whubera przyniósł mi ten efekt. @whuber argumentował, że PC1 i PC2 nie są niezależne i przedstawił tę silną korelację jako dowód na to. Mam jednak wrażenie, że przypadkowo odkrył błąd numeryczny w bibliotece LAPACK. Co tu się dzieje?

Oto kod R @ Whubera:

stat <- function(x) {u <- svd(x)$u; c(u[1,1], u[2, 2])};
Sigma <- matrix(c(1,0,0,1), 2);
sim <- t(replicate(1e3, stat(MASS::mvrnorm(10, c(0,0), Sigma))));
cor.test(sim[,1], sim[,2]);

Oto mój kod Matlab:

clear all
rng(7)

n = 1000;     %// Number of variables
k = 2;        %// Number of observations
Nrep = 1000;  %// Number of iterations (draws)

for rep = 1:Nrep
    X = randn(n,k);
    %// X = bsxfun(@minus, X, mean(X));
    [U,S,V] = svd(X,0);

    t(rep,:) = [U(1:10,1)' U(1:10,2)'];
end

figure
imagesc(corr(t), [-.5 .5])
axis square
hold on
plot(xlim, [10.5 10.5], 'k')
plot([10.5 10.5], ylim, 'k')
ameba mówi Przywróć Monikę
źródło
Jeśli użyjesz n = 4 i k = 3, zobaczysz również korelacje.
Aksakal
@Aksakal: tak, rzeczywiście, dzięki. Zedytowałem, aby usunąć deklarowaną różnicę między k = 2 i k = 3.
ameba mówi Przywróć Monikę

Odpowiedzi:

23

To nie jest błąd.

Jak zbadaliśmy (obszernie) w komentarzach, dzieją się dwie rzeczy. Po pierwsze, kolumny U są ograniczone do spełnienia wymagań SVD: każda musi mieć długość jednostki i być prostopadła do wszystkich pozostałych. Przedstawiamy U jako zmiennej losowej utworzonej z losowej macierzy X poprzez konkretnego algorytmu SVD, że w ten sposób zwrócić uwagę, że te k(k+1)/2) funkcjonalnie niezależne ograniczenia tworzenia zależności statystycznych między kolumnami U .

Zależności te można ujawnić w większym lub mniejszym stopniu poprzez badanie korelacji między składowymi U , ale pojawia się drugie zjawisko : rozwiązanie SVD nie jest unikalne. Co najmniej każdą kolumnę U można niezależnie negować, co daje co najmniej 2)k różnych rozwiązań z k kolumnami. Silne korelacje (powyżej 1/2) ) może być wywołane przez zmianę znaku kolumn odpowiednio. (Jednym ze sposobów, aby to zrobić jest podana w moim pierwszym komentarzu do odpowiedzi ameba jest w tym wątku: wymusić wszystko ujaja,ja=1,,k aby mieć ten sam znak, co powoduje, że wszystkie są ujemne lub wszystkie dodatnie z jednakowym prawdopodobieństwem). Z drugiej strony można wyeliminować wszystkie korelacje, wybierając znaki losowo, niezależnie, z jednakowymi prawdopodobieństwami. (Podam przykład poniżej w sekcji „Edycja”).

Ostrożnie, możemy częściowo odróżnić oba te zjawiska podczas czytania macierze rozrzutu składników U . Pewne cechy charakterystyczne - takie jak pojawienie się punktów prawie równomiernie rozmieszczonych w dobrze określonych regionach kołowych - wskazują na brak niezależności. Inne, takie jak wykresy rozrzutu pokazujące wyraźne niezerowe korelacje, oczywiście zależą od wyborów dokonanych w algorytmie - ale takie wybory są możliwe tylko z powodu braku niezależności.

Ostatecznym testem algorytmu rozkładu takiego jak SVD (lub Cholesky, LR, LU itp.) Jest to, czy robi to, co twierdzi. W tej sytuacji wystarczy sprawdzić, czy gdy SVD zwraca potrójną macierz (U,re,V.) , X jest odzyskiwane, aż do przewidywanego błędu zmiennoprzecinkowego, przez iloczyn UreV. ; że kolumny U i V. są ortonormalne; i że re jest ukośne, jego elementy ukośne są nieujemne i są ułożone w porządku malejącym. Zastosowałem takie testy do svdalgorytmu wRi nigdy nie stwierdziłem, że to pomyłka. Chociaż nie jest to żadna gwarancja, że ​​jest całkowicie poprawne, takie doświadczenie - które, jak sądzę, jest podzielane przez wielu ludzi - sugeruje, że każdy błąd wymagałby nadzwyczajnego wkładu, aby się ujawnić.

Poniżej znajduje się bardziej szczegółowa analiza konkretnych kwestii poruszonych w pytaniu.


Korzystając Rz svdprocedury, najpierw możesz sprawdzić, czy wraz ze wzrostem k korelacje między współczynnikami U maleją, ale wciąż są niezerowe. Jeśli po prostu wykonasz większą symulację, przekonasz się, że są one znaczące. (Gdy k=3) powinno wystarczyć 50000 iteracji.) W przeciwieństwie do twierdzenia w pytaniu, korelacje „nie znikają całkowicie”.

Po drugie, lepszym sposobem na zbadanie tego zjawiska jest powrót do podstawowej kwestii niezależności współczynników. Chociaż w większości przypadków korelacje są bliskie zeru, brak niezależności jest wyraźnie widoczny. Jest to najbardziej widoczna przez badania pełnej wielowymiarowego rozkład współczynników U . Charakter rozkładu pojawia się nawet w małych symulacjach, w których niezerowych korelacji nie można (jeszcze) wykryć. Na przykład zbadaj macierz rozrzutu współczynników. Aby było to wykonalne, ustawiłem rozmiar każdego symulowanego zestawu danych na 4 i utrzymałem k=2) , rysując w ten sposób 1000realizacje macierzy U 4×2) , tworząc macierz 1000 × 8 . Oto pełna macierz wykresów rozrzutu, ze zmiennymi wymienionymi według ich pozycji w U :U1000×8U

Postać

Skanowanie pierwszej kolumny ujawnia interesujący brak niezależności między u11 a drugim ujajot : spójrz na przykład, jak górna ćwiartka wykresu rozrzutu z u21 jest na przykład prawie pusta; lub zbadaj eliptyczną chmurę opadającą opisującą relację (u11,u22) i chmurę opadającą w dół dla pary (u21,u12) . Bliższe spojrzenie ujawnia wyraźny brak niezależności wśród prawie wszystkich tych współczynników: bardzo niewiele z nich wygląda na zdalnie niezależnych, chociaż większość z nich wykazuje korelację bliską zeru.

(Uwaga: większość chmur kołowych to rzuty z hipersfery stworzonej przez warunek normalizacji wymuszający jedność kwadratów wszystkich składników każdej kolumny.)

Matryce wykresów rozrzutu o k=3) i k=4 wykazują podobne wzorce: zjawiska te nie są ograniczone do k=2) , ani nie zależą od wielkości każdego symulowanego zestawu danych: po prostu trudniej jest je wygenerować i zbadać.

Wyjaśnienia dotyczące tych wzorców dotyczą algorytmu stosowanego do uzyskania U w rozkładzie liczby pojedynczej, ale wiemy, że takie wzorce nie-niezależności muszą istnieć dzięki bardzo definiującym właściwościom U : ponieważ każda kolejna kolumna jest (geometrycznie) ortogonalna w stosunku do poprzedniego te warunki ortogonalności narzucają zależności funkcjonalne między współczynnikami, co przekłada się na zależności statystyczne między odpowiadającymi zmiennymi losowymi.


Edytować

W odpowiedzi na komentarze warto zwrócić uwagę na to, w jakim stopniu te zjawiska zależności odzwierciedlają podstawowy algorytm (do obliczenia SVD) i jak bardzo są one związane z charakterem procesu.

Te specyficzne wzorce korelacji między współczynnikami zależeć bardzo wiele na arbitralnych wyborów dokonywanych przez algorytm SVD, ponieważ rozwiązanie to nie jest wyjątkowy: kolumny U zawsze może być niezależnie mnożone przez -1 lub 1 . Nie ma wewnętrznego sposobu wyboru znaku. Zatem, gdy dwa algorytmy SVD dokonują różnych (arbitralnych, a może nawet losowych) wyborów znaku, mogą one skutkować różnymi wzorami wykresów rozrzutu wartości (ujajot,ujajot) . Jeśli chcesz to zobaczyć, zamień statfunkcję w poniższym kodzie na

stat <- function(x) {
  i <- sample.int(dim(x)[1]) # Make a random permutation of the rows of x
  u <- svd(x[i, ])$u         # Perform SVD
  as.vector(u[order(i), ])   # Unpermute the rows of u
}

Ten pierwszy losowo porządkuje obserwacje x, wykonuje SVD, a następnie stosuje odwrotną kolejność, uaby dopasować oryginalną sekwencję obserwacji. Ponieważ efektem jest tworzenie mieszanin odbitych i obróconych wersji oryginalnych wykresów rozrzutu, wykresy rozrzutu w matrycy będą wyglądać znacznie bardziej jednolicie. Wszystkie przykładowe korelacje będą bardzo bliskie zeru (z założenia: podstawowe korelacje są dokładnie zerowe). Niemniej jednak brak niezależności będzie nadal oczywisty (w pojawiających się jednolitych okrągłych kształtach, szczególnie między uja,jot i uja,jot ).

Brak danych w niektórych kwadrantach niektórych oryginalnych wykresów rozrzutu (pokazanych na powyższym rysunku) wynika z tego, jak Ralgorytm SVD wybiera znaki dla kolumn.

W konkluzjach nic się nie zmienia. Ponieważ druga kolumna U jest prostopadła do pierwszej, to (uważana za zmienną losową wielowymiarową) jest zależna od pierwszej (również uważana za zmienną losową wielowymiarową). Nie możesz mieć wszystkich składników jednej kolumny niezależnych od wszystkich składników drugiej; wszystko, co możesz zrobić, to spojrzeć na dane w sposób, który przesłania zależności - ale zależność będzie się utrzymywać.


Oto zaktualizowany Rkod do obsługi przypadków k>2) i narysowania części macierzy wykresu rozrzutu.

k <- 2    # Number of variables
p <- 4    # Number of observations
n <- 1e3  # Number of iterations
stat <- function(x) as.vector(svd(x)$u)
Sigma <- diag(1, k, k); Mu <- rep(0, k)
set.seed(17)
sim <- t(replicate(n, stat(MASS::mvrnorm(p, Mu, Sigma))))
colnames(sim) <- as.vector(outer(1:p, 1:k, function(i,j) paste0(i,",",j)))
pairs(sim[, 1:min(11, p*k)], pch=".")
Whuber
źródło
3
Korelacja zachodzi między pierwszymi składowymi kolumn ponieważ tak działa algorytm SVD. To, że rzędy X są Gaussowskie, jest nieistotne: jestem pewien, że zauważyłeś, że współczynniki U nie są Gaussowskie. UXU
whuber
2
Nawiasem mówiąc, właśnie odkryli, że po prostu zastąpienie svd(X,0)przez svds(X)w moim kodu Matlab sprawia, że efekt znika! O ile mi wiadomo, te dwie funkcje używają różnych algorytmów SVD (obie są procedurami LAPACK, ale najwyraźniej różne). Nie wiem, czy R ma funkcję podobną do Matlaba svds, ale zastanawiam się, czy nadal utrzymujesz, że jest to „prawdziwy” efekt, a nie kwestia liczbowa.
ameba mówi Przywróć Monikę
4
Panowie, poczekajcie chwilę. Dlaczego nie mówisz o znaku? Znak wektora własnego jest w zasadzie dowolny. Ale program svd nie przypisuje go losowo, znak zależy od implementacji svd i danych. Jeśli po wypakowaniu Uzdecydujesz losowo, czy każda z kolumn ma pozostać bez zmian, czy też ma zmienić swój znak, to czy korelacje, o których mówisz, nie znikną?
ttnphns
2
@ttnphns To prawda, jak wyjaśniono w mojej edycji. Chociaż powoduje to zniknięcie korelacji, zależności między kolumnami nie znikają. (Dostarczona ulepszona wersja I odpowiada losowej zmianie znaków kolumn.)Ustat
whuber
2
Drobna uwaga (do tego wielkiego wątku!) SVD nie wymaga, aby elementy po przekątnej Sbyły w określonej kolejności; jest to kwestia wygody. Inne procedury gwarantują to (np. MATLAB svds), ale nie jest to ogólny wymóg. @amoeba: Patrząc na svds(co wydaje się wolne od tego problematycznego zachowania) obliczenie opiera się na faktycznym obliczeniu wartości własnych w pierwszej kolejności (więc nie używa standardowych procedur dgesdd/ dgesvdLAPACK - mocno podejrzewam, że używa dsyevr/ dsyevxna początku).
usεr11852 mówi: Przywróć Monic
11

Ta odpowiedź przedstawia replikację wyników @ Whuber w Matlabie, a także bezpośredni dowód, że korelacje są „artefaktem” tego, jak implementacja SVD wybiera znak dla komponentów.

Biorąc pod uwagę długi ciąg potencjalnie mylących komentarzy, chcę podkreślić dla przyszłych czytelników, że w pełni zgadzam się z następującymi kwestiami:

  1. W kontekście tej dyskusji jest z pewnością zmienną losową.U
  2. Kolumny muszą mieć długość 1 . Oznacza to, że elementy wewnątrz każdej kolumny nie są niezależne; ich kwadraty sumują się do jednego. Jednakże, nie oznacza żadnej korelacji pomiędzy U i 1 i U J 1 do ı j i korelacja próbki powinny być małe na dużą liczbę N r e p losowych zwraca.U1Uja1Ujot1jajotN.rmip
  3. Kolumny muszą być ortogonalne. Oznacza to, że elementy z różnych kolumn nie są niezależne; ich iloczyn skalarny wynosi zero. Ponownie, nie oznacza to żadnej korelacji między U i 1 a U j 2 , a korelacja próbki powinna być niewielka.UUja1Ujot2)

Moje pytanie brzmiało: dlaczego widzimy wysokie korelacje nawet dla dużej liczby losowych N r e p = 1000 ?0.2N.rmip=1000

Oto replika przykładu @ Whubera z , k = 2 i N r e p = 1000 w Matlabie:n=4k=2)N.rmip=1000

SVD

Po lewej stronie znajduje się macierz korelacji, po prawej wykresy rozrzutu podobne do @ whuber. Zgodność między naszymi symulacjami wydaje się idealna.

Teraz, kierując się pomysłową sugestią @ttnphns, przypisuję losowe znaki do kolumn , tj. Po tej linii:U

[U,S,V] = svd(X,0);

Dodaję następujące dwa wiersze:

U(:,1) = U(:,1) * sign(randn(1));
U(:,2) = U(:,2) * sign(randn(1));

Oto wynik:

SVD z losowymi znakami

Wszystkie korelacje znikają, dokładnie tak, jak się spodziewałem od samego początku !

11

UU

PS. Gratulacje dla @whuber za przekazanie dziś 100 tys. Reputacji!

ameba mówi Przywróć Monikę
źródło
Jeśli chcesz zobaczyć „wysoki” korelacje, używać tej wersji SVD w miejscu statw moim kodu: stat <- function(x) { u <- svd(x)$u; as.vector(sign(runif(1) - 1/2)*u %*% diag(sign(diag(u)))) }. Wybiera znaki kolumn w taki sposób, aby stworzyć duże korelacje między ( u 11 , u 22 , , uU(u11,u22,,ukk)UU
Właśnie skończyłem edytować tę odpowiedź, @whuber; Wymieniłem svdsprzez svdzastosowaną przez losowy wybór znak dla kolumn . Zgadzam się z tobą: wysokie korelacje były całkowicie spowodowane wyborem znaku, jak postawili hipoteza ttnphns. To zapewnia satysfakcjonującą odpowiedź na moje pytanie. Teraz spróbujmy rozwiązać spory! Zgodziłem się, że elementyUU
R±2/30.2
1
U
1
Intuicyjnie to uczciwe. Gdy tylko pierwsza oś główna zostanie zdefiniowana w przestrzeni, reszta pr. osie mają zmniejszoną swobodę. W przypadku danych 2D drugi (ostatni) komputer jest całkowicie związany, z wyjątkiem znaku. Wolę nazywać to ograniczeniem, a nie zależnością w sensie statystycznym.
ttnphns
0

xy

x2)+y2)=1

doov[x,y]=V.zar[xy]=mi[x2)y2)]-mi[xy]2)

xy

Aksakal
źródło
k(k+1)/2)UUreUUrek(k-1)/2)
U1Unkn>kUnn=1000n=4x2)+y2)=1U
xUyx2)+y2)=1Cov(x,y)=0x=sałata(θ)y=grzech(θ)θ[0,2)π)
UUjajot01/nnn1/n=1
1
UXU11U21ρnρρρ=0
ameba mówi Przywróć Monikę