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

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 $n=1000$$k=2$$X\sim \mathcal N (0, \mathbf I)$$1000 \times 2$$\mathbf X$$\mathbf X$$\mathbf X=\mathbf{USV}^\top$$\mathbf U$$U_{11}$$U_{22}$$\mathbf X$$N_\mathrm{rep}$ 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 $\pm0.2$$U_{11}$$U_{12}$$U_{21}$$U_{22}$$20$$\mathbf U$$10$$10$

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')

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 $u_{ii},i=1,\ldots, 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, D, V)$ , $X$ jest odzyskiwane, aż do przewidywanego błędu zmiennoprzecinkowego, przez iloczyn $UDV^\prime$ ; że kolumny $U$ i $V$ są ortonormalne; i że $D$ 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 $1000$realizacje macierzy U $4\times 2$ , tworząc macierz 1000 × 8 . Oto pełna macierz wykresów rozrzutu, ze zmiennymi wymienionymi według ich pozycji w U :$U$$1000\times 8$$U$

Skanowanie pierwszej kolumny ujawnia interesujący brak niezależności między $u_{11}$ a drugim $u_{ij}$ : spójrz na przykład, jak górna ćwiartka wykresu rozrzutu z $u_{21}$ jest na przykład prawie pusta; lub zbadaj eliptyczną chmurę opadającą opisującą relację $(u_{11}, u_{22})$ i chmurę opadającą w dół dla pary $(u_{21}, u_{12})$ . 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 $(u_{ij}, u_{i^\prime j^\prime})$ . 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 $u_{i,j}$ i $u_{i,j^\prime}$ ).

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\gt 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=".")

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ą.$\mathbf 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.$\mathbf U$$1$$U_{i1}$$U_{j1}$$i\ne j$$N_\mathrm{rep}$
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.$\mathbf U$$U_{i1}$$U_{j2}$

Moje pytanie brzmiało: dlaczego widzimy wysokie korelacje nawet dla dużej liczby losowych N r e p = 1000 ?$\sim 0.2$$N_\mathrm{rep}=1000$

Oto replika przykładu @ Whubera z , k = 2 i N r e p = 1000 w Matlabie:$n=4$$k=2$$N_\mathrm{rep}=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$

[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:

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

$1$$1$

$\mathbf U$$\mathbf U$

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

$x$$y$

$x$$y$