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 amostras do Gaussiano dimensional com média zero e covariância de identidade: . I montá-los num matriz de dados . (Opcionalmente, posso centralizar ou não, pois isso não influencia o seguinte.) Em seguida, executo a decomposição de valor singular (SVD) para obter . Vamos pegar alguns dois elementos particulares de , por exemplo, e , e perguntar o que é a correlação entre eles através de diferentes chama de . Eu esperaria que, se o número 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 ) entre , , e , 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 elementos "superiores" de parece (primeiros elementos da primeira coluna, depois os primeiros elementos da segunda coluna):
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')
fonte
Respostas:
Isso não é um bug.
Como exploramos (extensivamente) nos comentários, há duas coisas acontecendo. A primeira é que as colunas deU 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 deU , 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 asvocêeu eu, 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 devocê . 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 vocêD V′ ; que as colunas de você 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
svd
algoritmo emR
e 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.
Usandok aumenta, as correlações entre os coeficientes de você 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".
R
osvd
procedimento, primeiro você pode verificar se, à medida queSegundo, 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 devocê . 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 , desenhando1000 realizações da matriz 4 × 2 você , 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 você :
A varredura da primeira coluna revela uma interessante falta de independência entrevocê11 e a outra vocêeu j : veja como o quadrante superior do gráfico de dispersão com você21 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 comk = 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 obtervocê na decomposição de valor singular, mas sabemos que esses padrões de não independência devem existir pelas propriedades definidoras devocê : 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 devocê 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 ( ueu j, uEu′j′) . Se você quiser ver isso, substitua a
stat
função no código abaixo porEste primeiro reordena aleatoriamente as observaçõesvocêi , j e vocêi , j′ ).
x
, executa SVD e aplica a ordem inversau
para 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 entreA 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 devocê é 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á ok > 2 e desenhar uma parte da matriz do gráfico de dispersão.
R
código atualizado para lidar com os casosfonte
svd(X,0)
porsvds(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 Matlabsvds
, mas estou imaginando se você ainda vai sustentar que é um efeito "real" e não um problema numérico.U
você 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?stat
S
estejam em uma ordem específica; é uma questão de conveniência. Outras rotinas garantem isso (por exemplo, MATLABsvds
), mas isso não é um requisito geral. @amoeba: Olhandosvds
(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ãodgesdd
/dgesvd
LAPACK - eu suspeito fortemente que usedsyevr
/dsyevx
a princípio).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:
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,2 Nr e p= 1000
Aqui está uma replicação do exemplo de @ whuber com , k = 2 e N r e p = 1000 no Matlab:n = 4 k = 2 Nr e p= 1000
À 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ê
Eu adiciono as duas linhas a seguir:
Aqui está o resultado:
Todas as correlações desaparecem, exatamente como eu esperava desde o início !
PS. Parabéns a @whuber por passar 100k de reputação hoje!
fonte
stat
stat <- function(x) { u <- svd(x)$u; as.vector(sign(runif(1) - 1/2)*u %*% diag(sign(diag(u)))) }
svds
porsvd
seguido por uma escolha de sinal aleatório para as colunas deR
. 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 tornofonte