Para uma matriz aleatória, um SVD não deveria explicar nada? O que estou fazendo de errado?

Se eu construir uma matriz 2D composta inteiramente de dados aleatórios, esperaria que os componentes PCA e SVD não explicassem nada.

Em vez disso, parece que a primeira coluna SVD parece explicar 75% dos dados. Como isso pode ser possível? O que estou fazendo de errado?

Aqui está o enredo:

insira a descrição da imagem aqui

Aqui está o código R:

set.seed(1)
rm(list=ls())
m <- matrix(runif(10000,min=0,max=25), nrow=100,ncol=100)
svd1 <- svd(m, LINPACK=T)
par(mfrow=c(1,4))
image(t(m)[,nrow(m):1])
plot(svd1$d,cex.lab=2, xlab="SVD Column",ylab="Singluar Value",pch=19)

percentVarianceExplained = svd1$d^2/sum(svd1$d^2) * 100
plot(percentVarianceExplained,ylim=c(0,100),cex.lab=2, xlab="SVD Column",ylab="Percent of variance explained",pch=19)

cumulativeVarianceExplained = cumsum(svd1$d^2/sum(svd1$d^2)) * 100
plot(cumulativeVarianceExplained,ylim=c(0,100),cex.lab=2, xlab="SVD column",ylab="Cumulative percent of variance explained",pch=19)

Atualizar

Obrigado @Aaron. A correção, como você observou, foi adicionar escala à matriz para que os números fiquem centralizados em torno de 0 (ou seja, a média é 0).

m <- scale(m, scale=FALSE)

Aqui está a imagem corrigida, mostrando para uma matriz com dados aleatórios, a primeira coluna SVD está perto de 0, como esperado.

Imagem corrigida

r pca svd Contango
fonte

[0, 1]^{100}

$[0,1]^{100}$

R^{100}

$\mathbb{R}^{100}$

R^{n}

$\mathbb{R}^n$

n

$n$

1 / 3

$1/3$

n / 3 - (n - 1) / 12

$n/3-(n-1)/12$

1 / 12

$1/12$

(n / 3 - (n - 1) / 12) / (n / 3) = 3 / 4 + 1 / (4 n)

$(n/3-(n-1)/12)/(n/3)=3/4+1/(4n)$

n = 100

$n=100$

75.25

$75.25$

Respostas:

O primeiro PC está explicando que as variáveis não estão centralizadas em torno de zero. Escalar primeiro ou centralizar suas variáveis aleatórias em torno de zero terá o resultado esperado. Por exemplo, um destes:

m <- matrix(runif(10000,min=0,max=25), nrow=100,ncol=100)
m <- scale(m, scale=FALSE)

m <- matrix(runif(10000,min=-25,max=25), nrow=100,ncol=100)

Aaron deixou Stack Overflow
fonte

Você levanta um bom argumento, mas acho que isso conta apenas uma parte da história. Na verdade, eu acho que o OP tentaria cada uma delas e ainda ficaria surpreso com o resultado. O fato é que, como os valores singulares são ordenados na saída, eles não aparecerão (e de fato não serão) distribuídos uniformemente, como seria ingenuamente esperado a partir de dados "aleatórios". A distribuição Marchenko-Pastur rege seu comportamento neste caso.

cardinal

@ Aaron Obrigado, você estava absolutamente certo. Adicionei um gráfico da saída corrigida acima para mostrar como é bonito o resultado.

Contango

@ cardinal Obrigado pelo seu comentário, você está absolutamente certo (veja o gráfico produzido pelo código corrigido acima). Acredito que os valores SVD se tornariam menos uniformemente distribuídos à medida que a matriz diminuísse, pois uma matriz menor teria mais chances de ter padrões que não seriam esmagados pela lei de grandes números.

Contango

Adicionarei uma resposta mais visual à sua pergunta, usando uma comparação de modelo nulo. O procedimento aleatoriamente embaralha os dados em cada coluna para preservar a variação geral enquanto a covariância entre variáveis (colunas) é perdida. Isso é realizado várias vezes e a distribuição resultante de valores singulares na matriz aleatória é comparada aos valores originais.

Eu uso em prcompvez de svdpara a decomposição da matriz, mas os resultados são semelhantes:

set.seed(1)
m <- matrix(runif(10000,min=0,max=25), nrow=100,ncol=100)

S <- svd(scale(m, center = TRUE, scale=FALSE))
P <- prcomp(m, center = TRUE, scale=FALSE)
plot(S$d, P$sdev) # linearly related

A comparação do modelo nulo é realizada na matriz centralizada abaixo:

library(sinkr) # https://github.com/marchtaylor/sinkr

# centred data
Pnull <- prcompNull(m, center = TRUE, scale=FALSE, nperm = 100)
Pnull$n.sig
boxplot(Pnull$Lambda[,1:20], ylim=range(Pnull$Lambda[,1:20], Pnull$Lambda.orig[1:20]), outline=FALSE, col=8, border="grey50", log="y", main=paste("m (center=FALSE); n sig. =", Pnull$n.sig))
lines(apply(Pnull$Lambda, 2, FUN=quantile, probs=0.95))
points(Pnull$Lambda.orig[1:20], pch=16)

A seguir, é apresentado um gráfico de caixa da matriz permutada com o quantil de 95% de cada valor singular mostrado como a linha sólida. Os valores originais de PCA de msão os pontos. todos os quais estão abaixo da linha de 95% - Portanto, sua amplitude é indistinguível de ruído aleatório.

O mesmo procedimento pode ser feito na versão não centralizada mcom o mesmo resultado - sem valores singulares significativos:

# centred data
Pnull <- prcompNull(m, center = FALSE, scale=FALSE, nperm = 100)
Pnull$n.sig
boxplot(Pnull$Lambda[,1:20], ylim=range(Pnull$Lambda[,1:20], Pnull$Lambda.orig[1:20]), outline=FALSE, col=8, border="grey50", log="y", main=paste("m (center=TRUE); n sig. =", Pnull$n.sig))
lines(apply(Pnull$Lambda, 2, FUN=quantile, probs=0.95))
points(Pnull$Lambda.orig[1:20], pch=16)

Para comparação, vejamos um conjunto de dados com um conjunto de dados não aleatório: iris

# iris dataset example
m <- iris[,1:4]
Pnull <- prcompNull(m, center = TRUE, scale=FALSE, nperm = 100)
Pnull$n.sig
boxplot(Pnull$Lambda, ylim=range(Pnull$Lambda, Pnull$Lambda.orig), outline=FALSE, col=8, border="grey50", log="y", main=paste("m (center=FALSE); n sig. =", Pnull$n.sig))
lines(apply(Pnull$Lambda, 2, FUN=quantile, probs=0.95))
points(Pnull$Lambda.orig[1:20], pch=16)

Aqui, o 1º valor singular é significativo e explica mais de 92% da variação total:

P <- prcomp(m, center = TRUE)
P$sdev^2 / sum(P$sdev^2)
# [1] 0.924618723 0.053066483 0.017102610 0.005212184

Marc na caixa
fonte

+1. O exemplo do conjunto de dados Iris é interessante, porque olhando para os dois primeiros PCs (como, por exemplo, em seu próprio post aqui stats.stackexchange.com/a/88092 ), é bem claro que o segundo transmite algum sinal. O teste de permutação (também conhecido como shuffling) indica que apenas o 1º é "significativo". Está claro que o embaralhamento tende a subestimar o número de PCs: a grande variação do primeiro PC real se espalhará pelos PCs embaralhados e elevará todos eles a partir do segundo. Pode-se conceber testes mais sensíveis que explicam isso, mas isso raramente é feito.

Ameba diz Reinstate Monica

@amoeba - Excelente comentário. Há algum tempo que me pergunto sobre o efeito "espalhar". Suponho que um teste de validação cruzada possa ser um dos mais sensíveis a que você se refere (por exemplo, sua resposta aqui )? Seria ótimo se você pudesse fornecer um exemplo / referência.

Marc na caixa

Normalmente, prefiro usar a validação cruzada (com base no erro de reconstrução, conforme minha resposta aqui ), mas na verdade não tenho certeza se não está sofrendo o mesmo tipo de insensibilidade ou não. Pode fazer sentido experimentá-lo no conjunto de dados Iris. Com relação às abordagens baseadas em embaralhamento, não conheço nenhuma referência para explicar essa "disseminação", apenas conheço algumas pessoas que estavam trabalhando nela recentemente. Eu acho que eles queriam escrever em breve. A idéia é introduzir alguns fatores de redução de escala para as variações de PCs mais embaralhados.

Ameba diz Reinstate Monica

@amoeba - Obrigado por esse link. Isso explica muito para mim. Achei especialmente interessante ver que a validação cruzada no PCA utiliza métodos que podem operar em conjuntos de dados com valores ausentes. Fiz algumas tentativas nessa abordagem e (como você declara) a abordagem de shuffling de modelo nulo de fato tende a subestimar o número de PCs significativos. No entanto, para o conjunto de dados da íris, eu sempre retorno um único PC para erro de reconstrução. Interessante, dado o que você mencionou sobre o enredo. Pode ser que, se estivéssemos avaliando o erro com base nas previsões das espécies, os resultados pudessem ser diferentes.

Marc na caixa

Por curiosidade, tentei nos dados do Iris. Na verdade, recebo dois PCs significativos com o método de validação cruzada. Eu atualizei meu post vinculado, por favor, veja lá.

Ameba diz Reinstate Monica