Для случайной матрицы разве SVD не должен вообще ничего объяснять? Что я делаю неправильно?

Если бы я построил двумерную матрицу, состоящую полностью из случайных данных, я ожидал бы, что компоненты PCA и SVD по существу ничего не объясняют.

Вместо этого кажется, что первый столбец SVD, кажется, объясняет 75% данных. Как это может быть? Что я делаю неправильно?

Вот сюжет:

Вот код 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)

Обновить

Спасибо, Аарон. Исправление, как вы заметили, состояло в том, чтобы добавить масштабирование к матрице так, чтобы числа центрировались вокруг 0 ​​(т.е. среднее значение было 0).

m <- scale(m, scale=FALSE)

Вот исправленное изображение, показывающее для матрицы со случайными данными первый столбец SVD, близкий к 0, как и ожидалось.

Первый компьютер объясняет, что переменные не центрированы вокруг нуля. Масштабирование в первую очередь или центрирование случайных величин вокруг нуля даст ожидаемый результат. Например, любой из них:

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)

Я добавлю более наглядный ответ на ваш вопрос, используя сравнение нулевой модели. Процедура случайным образом перемешивает данные в каждом столбце, чтобы сохранить общую дисперсию, в то время как ковариация между переменными (столбцами) теряется. Это выполняется несколько раз, и полученное распределение сингулярных значений в рандомизированной матрице сравнивается с исходными значениями.

Я использую prcompвместо svdразложения матрицы, но результаты аналогичны:

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

Сравнение нулевой модели выполняется по центрированной матрице ниже:

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)

Ниже приведен блок-график перестановочной матрицы с 95-процентным квантилем каждого сингулярного значения, показанного сплошной линией. Исходные значения PCA mявляются точками. все из которых лежат ниже линии 95% - таким образом, их амплитуда неотличима от случайного шума.

Та же процедура может быть выполнена для нецентрированной версии mс тем же результатом - без значительных единичных значений:

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

Для сравнения, давайте посмотрим на набор данных с неслучайным набором данных: 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)

Здесь 1-е единственное значение является значимым и объясняет более 92% от общей дисперсии:

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