Correlações estranhas nos resultados SVD de dados aleatórios; eles têm uma explicação matemática ou é um bug do LAPACK?

21

Observo um comportamento muito estranho no resultado SVD de dados aleatórios, que posso reproduzir tanto no Matlab quanto no R. Parece um problema numérico na biblioteca LAPACK; é isso?

Eu desenho n=1000 amostras do Gaussiano k=2 dimensional com média zero e covariância de identidade: XN(0,I) . I montá-los num 1000×2 matriz de dados X . (Opcionalmente, posso centralizar X ou não, pois isso não influencia o seguinte.) Em seguida, executo a decomposição de valor singular (SVD) para obter X=USV . Vamos pegar alguns dois elementos particulares de U , por exemplo, U11 e U22 , e perguntar o que é a correlação entre eles através de diferentes chama deX . Eu esperaria que, se o númeroNrep de draws for razoavelmente grande, todas essas correlações sejam em torno de zero (ou seja, as correlações da população devem ser zero e as correlações da amostra serão pequenas).

No entanto, observo algumas correlações estranhamente fortes (em torno de ±0.2 ) entre U11 , U12 , U21 e U22 , e apenas entre esses elementos. Todos os outros pares de elementos têm correlações em torno de zero, conforme o esperado. Aqui está como a matriz de correlação dos 20 elementos "superiores" de U parece (primeiros 10 elementos da primeira coluna, depois os primeiros 10 elementos da segunda coluna):

Correlações estranhas SVD

Observe valores estranhamente altos nos cantos superior esquerdo de cada quadrante.

Foi o comentário deste @ whuber que trouxe esse efeito à minha atenção. A @whuber argumentou que o PC1 e o PC2 não são independentes e apresentaram essa forte correlação como uma evidência disso. No entanto, minha impressão é que ele descobriu acidentalmente um bug numérico na biblioteca LAPACK. O que está acontecendo aqui?

Aqui está o código R do @ whuber:

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]);

Aqui está o meu código 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 diz Restabelecer Monica
fonte
Se você usar n = 4 ek = 3, também verá as correlações.
Aksakal
@ Aksakal: sim, de fato, obrigado. Editei para remover a diferença reivindicada entre k = 2 e k = 3.
Ameba diz Reinstate Monica

Respostas:

23

Isso não é um bug.

Como exploramos (extensivamente) nos comentários, há duas coisas acontecendo. A primeira é que as colunas de U são restritas para atender aos requisitos de SVD: cada uma deve ter comprimento unitário e ser ortogonal a todas as outras. Vendo U como uma variável aleatória criada a partir de uma matriz aleatória X através de um algoritmo SVD específico, observamos que esses k(k+1)/2 restrições funcionalmente independentes criar dependências estatísticas entre as colunas deU .

Essas dependências podem ser reveladas em maior ou menor grau, estudando as correlações entre os componentes de U , mas surge um segundo fenômeno : a solução SVD não é única. No mínimo, cada coluna de U pode ser negada independentemente, fornecendo pelo menos 2k soluções distintas com k colunas. Uma forte correlação (superior a 1/2 ) pode ser induzido alterando os sinais das colunas de forma adequada. (Uma maneira de fazer isso é dada no meu primeiro comentário à resposta da Amoeba neste tópico: forço todas asuii,i=1,,k para ter o mesmo sinal, tornando-os todos negativos ou positivos com igual probabilidade.) Por outro lado, todas as correlações podem desaparecer escolhendo os sinais aleatoriamente, independentemente, com probabilidades iguais. (Dou um exemplo abaixo na seção "Editar".)

Com cuidado, podemos parcialmente discernir esses dois fenômenos ao ler matrizes de dispersão dos componentes de U . Certas características - como a aparência de pontos quase uniformemente distribuídos em regiões circulares bem definidas - desmentem uma falta de independência. Outros, como gráficos de dispersão que mostram correlações claras diferentes de zero, obviamente dependem de escolhas feitas no algoritmo - mas tais escolhas são possíveis apenas devido à falta de independência em primeiro lugar.

O teste final de um algoritmo de decomposição como SVD (ou Cholesky, LR, LU, etc.) é se ele faz o que afirma. Nesta circunstância, é suficiente verificar que quando retorna a SVD triplo de matrizes (U,D,V) , que X é recuperado, até erro de ponto flutuante antecipado, pelo produto UDV ; que as colunas de U e de V são ortonormais; e que D é diagonal, seus elementos diagonais não são negativos e estão organizados em ordem decrescente. Eu apliquei esses testes ao svdalgoritmo emRe nunca achei que estivesse errado. Embora isso não seja garantia de que esteja perfeitamente correto, essa experiência - que acredito ser compartilhada por muitas pessoas - sugere que qualquer bug exigiria algum tipo extraordinário de entrada para se manifestar.

O que se segue é uma análise mais detalhada de pontos específicos levantados na questão.


Usando Ro svdprocedimento, primeiro você pode verificar se, à medida que k aumenta, as correlações entre os coeficientes de U tornam mais fracas, mas ainda são diferentes de zero. Se você simplesmente realizasse uma simulação maior, descobriria que são significativas. (Quando k=3 , 50000 iterações devem ser suficientes.) Ao contrário do que afirma a pergunta,as correlações "não desaparecem completamente".

Segundo, uma maneira melhor de estudar esse fenômeno é voltar à questão básica da independência dos coeficientes. Embora as correlações tendam a ser próximas de zero na maioria dos casos, a falta de independência é claramente evidente. Isto é feito mais aparente através do estudo da distribuição multivariada completa dos coeficientes de U . A natureza da distribuição surge mesmo em pequenas simulações nas quais as correlações diferentes de zero ainda não podem ser detectadas. Por exemplo, examine uma matriz de dispersão dos coeficientes. Para tornar isso viável, defino o tamanho de cada conjunto de dados simulado como 4 e mantive k=2 , desenhando1000realizações da matriz 4×2U , criando uma matriz 1000×8 . Aqui está sua matriz completa do gráfico de dispersão, com as variáveis ​​listadas por suas posições em U :

Figura

A varredura da primeira coluna revela uma interessante falta de independência entre u11 e a outra uij : veja como o quadrante superior do gráfico de dispersão com u21 está quase vazio, por exemplo; ou examine a nuvem elíptica inclinada para cima que descreve a relação (u11,u22) e a nuvem inclinada para baixo para a (u21,u12)par 12 ). Um exame mais atento revela uma clara falta de independência entre quase todos esses coeficientes: muito poucos deles parecem remotamente independentes, embora a maioria deles apresente correlação quase zero.

(NB: A maioria das nuvens circulares são projeções de uma hiperesfera criada pela condição de normalização, forçando a soma dos quadrados de todos os componentes de cada coluna a serem unidos.)

Matrizes de gráficos de dispersão com k=3 e k=4 exibem padrões semelhantes: esses fenômenos não estão confinados a k=2 , nem dependem do tamanho de cada conjunto de dados simulado: eles apenas ficam mais difíceis de gerar e examinar.

As explicações para esses padrões vão para o algoritmo usado para obter U na decomposição de valor singular, mas sabemos que esses padrões de não independência devem existir pelas propriedades definidoras deU : uma vez que cada coluna sucessiva é (geometricamente) ortogonal à anterior outras, essas condições de ortogonalidade impõem dependências funcionais entre os coeficientes, o que se traduz em dependências estatísticas entre as variáveis ​​aleatórias correspondentes.


Editar

Em resposta aos comentários, vale a pena observar até que ponto esses fenômenos de dependência refletem o algoritmo subjacente (para calcular um SVD) e o quanto eles são inerentes à natureza do processo.

Os padrões específicos de correlações entre coeficientes dependem muito das escolhas arbitrárias feitas pelo algoritmo SVD, porque a solução não é única: as colunas de U sempre podem ser multiplicadas independentemente por 1 ou 1 . Não há maneira intrínseca de escolher o sinal. Assim, quando dois algoritmos SVD fazem escolhas diferentes (arbitrárias ou talvez aleatórias) de sinal, eles podem resultar em diferentes padrões de gráficos de dispersão dos valores (uij,uij) . Se você quiser ver isso, substitua a statfunção no código abaixo por

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
}

Este primeiro reordena aleatoriamente as observações x, executa SVD e aplica a ordem inversa upara corresponder à sequência de observação original. Como o efeito é formar misturas de versões refletidas e rotacionadas dos gráficos de dispersão originais, os gráficos de dispersão na matriz parecerão muito mais uniformes. Todas as correlações de amostra serão extremamente próximas de zero (por construção: as correlações subjacentes são exatamente zero). No entanto, a falta de independência ainda será óbvia (nas formas circulares uniformes que aparecem, particularmente entre ui,j e ui,j ).

A falta de dados em alguns quadrantes de alguns dos gráficos de dispersão originais (mostrados na figura acima) decorre de como o R algoritmo SVD seleciona sinais para as colunas.

Nada muda sobre as conclusões. Como a segunda coluna de U é ortogonal à primeira, ela (considerada uma variável aleatória multivariada) é dependente da primeira (também considerada como uma variável aleatória multivariada). Você não pode ter todos os componentes de uma coluna independentes de todos os componentes da outra; tudo o que você pode fazer é examinar os dados de maneiras que obscurecem as dependências - mas a dependência persistirá.


Aqui está o Rcódigo atualizado para lidar com os casos k>2 e desenhar uma parte da matriz do gráfico de dispersão.

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
fonte
3
A correlação ocorre entre os primeiros componentes das colunas de porque é assim que o algoritmo SVD funciona. O fato de as linhas de X serem gaussianas é irrelevante: tenho certeza de que você notou que os coeficientes de U não são gaussianos. UXU
whuber
2
A propósito, eu acabei de descobrir que simplesmente substituir svd(X,0)por svds(X)no meu código Matlab faz o efeito desaparecer! Até onde eu sei, essas duas funções usam algoritmos SVD diferentes (ambos são rotinas LAPACK, mas aparentemente diferentes). Não sei se R tem uma função semelhante à do Matlab svds, mas estou imaginando se você ainda vai sustentar que é um efeito "real" e não um problema numérico.
Ameba diz Reinstate Monica
4
Senhores, esperem um pouco. Por que você não está falando do sinal? O sinal de um vetor próprio é basicamente arbitrário. Mas o programa do svd não o atribui aleatoriamente, o sinal depende da implementação do svd e dos dados. Se, após a extração, Uvocê decidir aleatoriamente se cada uma de suas colunas deve permanecer como está ou se deve alterar seu sinal, as correlações de que você está falando desaparecerão?
ttnphns 24/02
2
@ttnphns Isso está correto, como explicado na minha edição. Embora isso faça desaparecer as correlações, as dependências entre as colunas de não desaparecem. (A versão aprimorada de que eu forneci é equivalente a alterar aleatoriamente os sinais das colunas.)Ustat
whuber
2
Um ponto menor (para esta grande discussão!) O SVD não exige que os elementos na diagonal de Sestejam em uma ordem específica; é uma questão de conveniência. Outras rotinas garantem isso (por exemplo, MATLAB svds), mas isso não é um requisito geral. @amoeba: Olhando svds(o que parece livre desse comportamento problemático), a computação é baseada na computação dos valores próprios primeiro (por isso não usa as rotinas padrão dgesdd/ dgesvdLAPACK - eu suspeito fortemente que use dsyevr/ dsyevxa princípio).
usεr11852 diz Reinstate Monic
11

Esta resposta apresenta uma replicação dos resultados do @ whuber no Matlab e também uma demonstração direta de que as correlações são um "artefato" de como a implementação do SVD escolhe assinar para componentes.

Dada a longa cadeia de comentários potencialmente confusos, quero enfatizar para os futuros leitores que concordo plenamente com o seguinte:

  1. No contexto desta discussão, certamente é uma variável aleatória.U
  2. As colunas de devem ter o comprimento 1 . Isso significa que os elementos dentro de cada coluna não são independentes; seus quadrados somam um. No entanto, isso não implica nenhuma correlação entre U i 1 e U j 1 para i j , e a correlação da amostra deve ser pequena para um número grande N r e pU1Ui1Uj1ijNrep de sorteios aleatórios.
  3. As colunas de devem ser ortogonais. Isso significa que os elementos de diferentes colunas não são independentes; seu produto escalar é zero. Novamente, isso não implica nenhuma correlação entre U i 1 e U j 2 , e a correlação da amostra deve ser pequena.UUi1Uj2

A minha pergunta era: por que vemos altas correlações de até mesmo o para grande número de aleatório empates N r e p = 1,000 ?0.2Nrep=1000

Aqui está uma replicação do exemplo de @ whuber com , k = 2 e N r e p = 1000 no Matlab:n=4k=2Nrep=1000

SVD

À esquerda está a matriz de correlação, à direita - gráficos de dispersão semelhantes aos @ whuber's. O acordo entre nossas simulações parece perfeito.

Agora, seguindo uma sugestão engenhosa de @ttnphns, atribuo sinais aleatórios às colunas de , ou seja, após esta linha:você

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

Eu adiciono as duas linhas a seguir:

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

Aqui está o resultado:

SVD com sinais aleatórios

Todas as correlações desaparecem, exatamente como eu esperava desde o início !

11

vocêvocê

PS. Parabéns a @whuber por passar 100k de reputação hoje!

ameba diz Restabelecer Monica
fonte
statstat <- function(x) { u <- svd(x)$u; as.vector(sign(runif(1) - 1/2)*u %*% diag(sign(diag(u)))) }você(você11,você22,...,vocêkk)vocêvocê
whuber
Acabei de editar esta resposta, @whuber; Substituí svdspor svdseguido por uma escolha de sinal aleatório para as colunas devocêvocê
Eu não concordo com isso: a figura mostra com precisão os resultados das configurações padrão do SVD no R. Forneci explicações para esse comportamento, fundamentado no algoritmo e na teoria subjacente. Se algum aspecto disso é "super enganador", é claro que vou me esforçar para corrigir ou esclarecer o problema, mas neste momento não estou ciente de um problema específico a ser corrigido. Você pode estar interessado em observar que a versão modificada fornecida no meu comentário anterior a esta pergunta é capaz de produzir correlações em torno±2/30,2
1
você
1
Intuitivamente, isso é justo. Assim que o primeiro eixo principal for definido no espaço, o restante será pr. os eixos obtêm liberdade reduzida. No caso de dados 2D, o segundo (último) PC está totalmente ligado, exceto o sinal. Prefiro chamá-lo de restrição, não de dependência estatística.
ttnphns
0

xy

x2+y2=1

Cov[x,y]=Vumar[xy]=E[x2y2]-E[xy]2

xy

Aksakal
fonte
k(k+1)/2vocêvocêDvocêvocêDk(k-1)/2
você1Unkn>kUnn=1000n=4x2+y2=1você
xUyx2+y2=1Cov(x,y)=0x=cos(θ)y=sin(θ)θ[0,2π)
vocêvocêEuj0 01/nnn1/n=1
1
vocêXU11U21ρnρρρ=0
Ameba diz Reinstate Monica