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 V ⊤ U U 11 U 22 X N r e p 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 22U10
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')
źródło
Odpowiedzi:
To nie jest błąd.
Jak zbadaliśmy (obszernie) w komentarzach, dzieją się dwie rzeczy. Po pierwsze, kolumnyU 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ładowymiU , 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 uja ja, i = 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ówU . 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, D , V) , X jest odzyskiwane, aż do przewidywanego błędu zmiennoprzecinkowego, przez iloczyn UD V.′ ; ż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
svd
algorytmu wR
i 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ąck 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”.
R
zsvd
procedury, najpierw możesz sprawdzić, czy wraz ze wzrostemPo 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ówU . 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 1000 realizacje 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 :U 1000 × 8 U
Skanowanie pierwszej kolumny ujawnia interesujący brak niezależności międzyu11 a drugim uI j : 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 ok = 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 uzyskaniaU 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: kolumnyU 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 ( uI j, uja′jot′) . Jeśli chcesz to zobaczyć, zamień
stat
funkcję w poniższym kodzie naTen pierwszy losowo porządkuje obserwacjeuja , j i uja , j′ ).
x
, wykonuje SVD, a następnie stosuje odwrotną kolejność,u
aby 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ędzyBrak danych w niektórych kwadrantach niektórych oryginalnych wykresów rozrzutu (pokazanych na powyższym rysunku) wynika z tego, jak
R
algorytm SVD wybiera znaki dla kolumn.W konkluzjach nic się nie zmienia. Ponieważ druga kolumnaU 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 zaktualizowanyk > 2 i narysowania części macierzy wykresu rozrzutu.
R
kod do obsługi przypadkówźródło
svd(X,0)
przezsvds(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 Matlabasvds
, ale zastanawiam się, czy nadal utrzymujesz, że jest to „prawdziwy” efekt, a nie kwestia liczbowa.U
zdecydujesz 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ą?stat
S
były w określonej kolejności; jest to kwestia wygody. Inne procedury gwarantują to (np. MATLABsvds
), ale nie jest to ogólny wymóg. @amoeba: Patrząc nasvds
(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 procedurdgesdd
/dgesvd
LAPACK - mocno podejrzewam, że używadsyevr
/dsyevx
na początku).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:
Moje pytanie brzmiało: dlaczego widzimy wysokie korelacje nawet dla dużej liczby losowych N r e p = 1000 ?∼ 0,2 N.r e p= 1000
Oto replika przykładu @ Whubera z , k = 2 i N r e p = 1000 w Matlabie:n = 4 k = 2 N.r e p= 1000
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
Dodaję następujące dwa wiersze:
Oto wynik:
Wszystkie korelacje znikają, dokładnie tak, jak się spodziewałem od samego początku !
PS. Gratulacje dla @whuber za przekazanie dziś 100 tys. Reputacji!
źródło
stat
w 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 , … , usvds
przezsvd
zastosowaną 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 elementyR
źródło