Para una matriz aleatoria, ¿no debería una SVD explicar nada? ¿Qué estoy haciendo mal?

Si construyo una matriz 2D compuesta completamente de datos aleatorios, esperaría que los componentes PCA y SVD esencialmente no explicaran nada.

En cambio, parece que la primera columna SVD parece explicar el 75% de los datos. ¿Cómo puede ser esto posible? ¿Qué estoy haciendo mal?

Aquí está la trama:

Aquí está el 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)

Actualizar

Gracias @ Aaron. La solución, como notó, era agregar una escala a la matriz para que los números se centren alrededor de 0 (es decir, la media es 0).

m <- scale(m, scale=FALSE)

Aquí está la imagen corregida, que muestra para una matriz con datos aleatorios, la primera columna SVD está cerca de 0, como se esperaba.

La primera PC explica que las variables no están centradas alrededor de cero. Escalar primero o centrar sus variables aleatorias alrededor de cero tendrá el resultado que espera. Por ejemplo, cualquiera de estos:

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)

Agregaré una respuesta más visual a su pregunta, mediante el uso de una comparación de modelo nulo. El procedimiento baraja aleatoriamente los datos en cada columna para preservar la varianza general mientras se pierde la covarianza entre variables (columnas). Esto se realiza varias veces y la distribución resultante de valores singulares en la matriz aleatoria se compara con los valores originales.

Yo uso en prcomplugar de svdpara la descomposición de la matriz, pero los resultados son similares:

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

La comparación del modelo nulo se realiza en la matriz centrada a continuación:

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)

El siguiente es un diagrama de caja de la matriz permutada con el cuantil del 95% de cada valor singular que se muestra como la línea continua. Los valores originales de PCA de mson los puntos. todo lo cual se encuentra debajo de la línea del 95%. Por lo tanto, su amplitud es indistinguible del ruido aleatorio.

Se puede realizar el mismo procedimiento en la versión no centrada de mcon el mismo resultado. Sin 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)

A modo de comparación, veamos un conjunto de datos con un conjunto de datos no aleatorio: 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)

Aquí, el primer valor singular es significativo y explica más del 92% de la varianza total:

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