Странные корреляции в результатах SVD случайных данных; у них есть математическое объяснение или это ошибка LAPACK?

21

Я наблюдаю очень странное поведение в результате SVD случайных данных, которое я могу воспроизвести как в Matlab, так и в R. Это похоже на некоторую числовую проблему в библиотеке LAPACK; это?

Я рисую выборок из мерного гауссиана с нулевым средним и единичной ковариацией: . Я собрать их в данных матрица . (При желании я могу центрировать или нет, это не влияет на следующее.) Затем я выполняю разложение по сингулярным значениям (SVD), чтобы получить . Давайте возьмем некоторые два отдельных элементов , например и , и спросить , что корреляция между ними по разным розыгрышам . Я ожидаю, что если числоNзнак равно1000Кзнак равно2Икс~N(0,я)1000×2ИксИксИксзнак равноUSВUU11U22ИксNреп розыгрышей достаточно велик, тогда все такие корреляции должны быть около нуля (т. Е. Корреляции популяций должны быть нулевыми, а выборочные корреляции будут небольшими).

Однако я наблюдаю некоторые странно сильные корреляции (около ) между , , и и только между этими элементами. Все остальные пары элементов имеют корреляции около нуля, как и ожидалось. Вот как выглядит матрица корреляции для «верхних» элементов (первые элементов первого столбца, затем первые элементов второго столбца):±0.2U11U12U21U2220U1010

SVD странные корреляции

Обратите внимание на странно высокие значения в верхнем левом углу каждого квадранта.

Именно этот комментарий @ whuber привлек мое внимание к этому эффекту. @whuber утверждал, что ПК1 и ПК2 не являются независимыми, и представил эту сильную корреляцию в качестве доказательства этого. Однако у меня сложилось впечатление, что он случайно обнаружил числовую ошибку в библиотеке LAPACK. Что здесь происходит?

Вот код @ whuber's:

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

Вот мой код 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')
амеба говорит восстановить монику
источник
Если вы используете n = 4 и k = 3, вы также увидите корреляции.
Аксакал
@Aksakal: да, действительно, спасибо. Я отредактировал, чтобы убрать заявленную разницу между k = 2 и k = 3.
амеба говорит восстановить Монику

Ответы:

23

Это не ошибка.

Как мы подробно рассмотрели в комментариях, происходит две вещи. Во-первых, столбцы ограничены для соответствия требованиям SVD: каждый должен иметь длину блока и быть ортогональным по отношению ко всем остальным. Просмотр как случайная величина , созданная из случайной матрицы через конкретный алгоритм SVD, мы тем самым отметить , что это функционально независимым ограничения создают статистические зависимости между столбцами .UUИксК(К+1)/2U

Эти зависимости могут быть выявлены в большей или меньшей степени путем изучения корреляций между компонентами , но возникает второе явление : решение SVD не является уникальным. Как минимум, каждый столбец можно независимо отрицать, давая как минимум различных решений с столбцами. Сильные корреляции (превышающие ) могут быть вызваны соответствующим изменением знаков столбцов. (Один из способов сделать это дан в моем первом комментарии к ответу Амебы в этой теме: я принудительно заставляю всеUU 2 к к +1 / 2 у я я , я = 1 , ... , KU2КК1/2Uяя,язнак равно1,...,Киметь один и тот же знак, что делает их всех отрицательными или положительными с одинаковой вероятностью.) С другой стороны, все корреляции можно сделать равными нулю, выбрав знаки случайным образом, независимо с равными вероятностями. (Я приведу пример ниже в разделе «Редактировать».)

С осторожностью, мы можем частично различить оба эти явления при чтении диаграммы рассеивания матрицы компонентов . Определенные характеристики - такие как появление точек, почти равномерно распределенных в четко определенных круглых областях, - свидетельствуют об отсутствии независимости. Другие, такие как диаграммы рассеяния, показывающие четкие ненулевые корреляции, очевидно, зависят от выбора, сделанного в алгоритме - но такие выборы возможны только из-за отсутствия независимости в первую очередь.U

Окончательный тест алгоритма декомпозиции, такого как SVD (или Cholesky, LR, LU и т. Д.), Заключается в том, выполняет ли он то, что заявляет. В этом случае достаточно проверить, что когда SVD возвращает тройку матриц , то восстанавливается с точностью до ожидаемой ошибки с плавающей запятой произведением ; что столбцы и ортонормированы; и что диагонально, его диагональные элементы неотрицательны и расположены в порядке убывания. Я применил такие тесты к алгоритму в(U,D,В)ИксUDВ'UВDsvdRи никогда не находил это ошибочным. Хотя это и не является гарантией, что это совершенно правильно, такой опыт, который, как я полагаю, разделяют очень многие люди, предполагает, что любая ошибка потребует какого-то необычного вклада для проявления.

Далее следует более подробный анализ конкретных вопросов, затронутых в этом вопросе.


Используя R«s svdпроцедуры, сначала вы можете проверить , что в качестве увеличивается, корреляция между коэффициентами ослабевать, но они по - прежнему отличны от нуля. Если бы вы просто выполнили большую симуляцию, вы бы обнаружили, что они значительны. (Когда , должно быть достаточно 50000 итераций.) Вопреки утверждению в вопросе, корреляции не «полностью исчезают».КUk = 3Кзнак равно3

Во-вторых, лучший способ изучить это явление - вернуться к основному вопросу независимости коэффициентов. Хотя корреляции, как правило, близки к нулю в большинстве случаев, отсутствие независимости ясно видно. Это стало наиболее очевидным при изучении полного многомерного распределения коэффициентов . Характер распределения проявляется даже в небольших моделях, в которых ненулевые корреляции не могут (пока) быть обнаружены. Например, рассмотрим матрицу рассеяния для коэффициентов. Чтобы сделать это практически возможным, я установил размер каждого моделируемого набора данных равным и сохранил , тем самым нарисовав реализацийU4Кзнак равно210004×2 матрицы , создав матрицу . Вот его полная матрица рассеяния с переменными, перечисленными по их позициям в :U1000×8U

фигура

Сканирование вниз по первому столбцу показывает интересное отсутствие независимости между и другим : посмотрите, например, как верхний квадрант диаграммы рассеяния с почти свободен; или исследуйте эллиптическое наклонное вверх облако, описывающее соотношение и наклонное облако вниз для пары . При внимательном рассмотрении обнаруживается явное отсутствие независимости почти всех этих коэффициентов: очень немногие из них выглядят дистанционно независимыми, хотя большинство из них демонстрируют почти нулевую корреляцию.U11UяJU21(U11,U22)(U21,U12)

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

Матрицы диаграммы рассеяния с и демонстрируют сходные закономерности: эти явления не ограничиваются и не зависят от размера каждого моделируемого набора данных: их просто сложнее генерировать и исследовать.Кзнак равно3Кзнак равно4Кзнак равно2

Объяснения для этих шаблонов идут к алгоритму, используемому для получения в разложении по сингулярному значению, но мы знаем, что такие шаблоны несамостоятельности должны существовать с помощью определяющих свойств : поскольку каждый последующий столбец (геометрически) ортогонален предыдущему Кроме того, эти условия ортогональности налагают функциональные зависимости между коэффициентами, которые тем самым переводят в статистические зависимости среди соответствующих случайных величин.UUU


редактировать

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

В конкретных моделях корреляций между коэффициентами зависят много от произвольных выборов , сделанных с помощью алгоритма SVD, потому что решение не является уникальным: столбцы всегда могут быть независимы умноженные на или . Не существует внутреннего способа выбрать знак. Таким образом, когда два алгоритма SVD делают разные (произвольные или, возможно, даже случайные) варианты выбора знака, они могут привести к различным шаблонам диаграмм рассеяния значений . Если вы хотите увидеть это, замените функцию в коде ниже наU-11(UяJ,Uя'J')stat

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
}

Это сначала случайным образом xпереупорядочивает наблюдения , выполняет SVD, затем применяет обратный порядок, uчтобы соответствовать исходной последовательности наблюдения. Поскольку эффект состоит в формировании смеси отраженных и повернутых версий исходных диаграмм рассеяния, диаграммы рассеяния в матрице будут выглядеть намного более однородными. Все выборочные корреляции будут очень близки к нулю (по построению: лежащие в основе корреляции точно равны нулю). Тем не менее, отсутствие независимости все еще будет очевидным (в однородных круглых формах, которые появляются, особенно между и ).Uя,JUя,J'

Отсутствие данных в некоторых квадрантах некоторых из исходных диаграмм рассеяния (показано на рисунке выше) возникает из-за того, как Rалгоритм SVD выбирает знаки для столбцов.

В выводах ничего не меняется. Поскольку второй столбец ортогонален первому, он (рассматриваемый как многомерная случайная величина) зависит от первого (также рассматриваемого как многомерная случайная величина). Вы не можете иметь все компоненты одного столбца независимыми от всех компонентов другого; все, что вы можете сделать, это посмотреть на данные таким образом, чтобы скрыть зависимости, но зависимость сохранится.U


Здесь обновлен Rкод для обработки случаев и рисования части матрицы рассеяния.К>2

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
источник
3
Корреляция происходит между первыми компонентами столбцов потому что так работает алгоритм SVD. То , что ряды гауссово несущественно: Я уверен , вы заметили , что коэффициенты являются не гауссовы. X UUИксU
whuber
2
Кстати, я только что обнаружил, что простая замена svd(X,0)на svds(X)в моем коде Matlab эффект исчезает! Насколько я знаю, эти две функции используют разные алгоритмы SVD (обе являются подпрограммами LAPACK, но, очевидно, разными). Я не знаю, имеет ли R функцию, аналогичную функции Matlab svds, но мне интересно, если вы все еще собираетесь утверждать, что это «реальный» эффект, а не числовая проблема.
говорит амеба, восстанови Монику
4
Господа, подождите минутку. Почему ты не говоришь о знаке? Знак собственного вектора в основном произвольный. Но программа svd не назначает его случайным образом, знак зависит от реализации svd и данных. Если после извлечения Uвы случайным образом решите, должен ли каждый из его столбцов оставаться в том виде, в каком он есть, или будет менять свой знак, не исчезнут ли корреляции, о которых вы говорите?
ttnphns
2
@ttnphns Это правильно, как объяснено в моем редактировании. Хотя это делает корреляции исчезающими, зависимости между столбцами тем самым не исчезают . (Расширенная версия, которую я предоставил, эквивалентна случайному изменению знаков столбцов.)Ustat
whuber
2
Незначительный момент (для этой великой темы!) SVD не требует, чтобы элементы в диагонали Sбыли в определенном порядке; это вопрос удобства. Другие процедуры гарантируют это (например, MATLAB svds), но это не является общим требованием. @amoeba: Смотря svds(что кажется свободным от этого проблемного поведения), вычисления основаны на фактическом вычислении собственных значений вначале (поэтому он не использует стандартные процедуры dgesdd/ dgesvdLAPACK - я сильно подозреваю, что он использует dsyevr/ dsyevxсначала).
usεr11852 говорит восстановить Monic
11

Этот ответ представляет репликацию результатов @ whuber в Matlab, а также прямую демонстрацию того, что корреляции являются «артефактом» того, как реализация SVD выбирает знак для компонентов.

Учитывая длинную цепочку потенциально запутанных комментариев, я хочу подчеркнуть для будущих читателей, что я полностью согласен со следующим:

  1. В контексте этого обсуждения определенно является случайной величиной.U
  2. Столбцы должны иметь длину . Это означает, что элементы внутри каждого столбца не являются независимыми; их квадраты равны единице. Однако это не подразумевает какой-либо корреляции между и для , и выборочная корреляция должна быть крошечной для большого числа случайных отрисовок. 1 U i 1 U j 1U1Uя1UJ1N r e pяJNреп
  3. Столбцы должны быть ортогональны. Это означает, что элементы из разных столбцов не являются независимыми; их скалярное произведение равно нулю. Опять же, это не подразумевает никакой корреляции между и , и выборочная корреляция должна быть крошечной.UU U j 2Uя1UJ2

Мой вопрос был: почему мы видим высокие корреляции даже для большого числа случайных дро ?~0.2Nрепзнак равно1000

Вот репликация примера @ whuber с , и в Matlab:k = 2Nзнак равно4Кзнак равно2Nрепзнак равно1000

SVD

Слева - корреляционная матрица, справа - диаграммы рассеяния, подобные @ whuber's. Соглашение между нашими симуляциями кажется идеальным.

Теперь, следуя гениальному предложению @ttnphns, я назначаю случайные знаки столбцам , то есть после этой строки:U

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

Я добавляю следующие две строки:

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

Вот результат:

СВД со случайными признаками

Все корреляции исчезают, как я и ожидал с самого начала !

Как говорит @whuber, отсутствие независимости можно увидеть в идеальной круглой форме некоторых диаграмм рассеяния (поскольку длина каждого столбца должна равняться , сумма квадратов любых двух элементов не может превышать ). Но корреляции исчезают.111

Подводя итог всей проблеме, мы видим, что появляются сильные корреляции, потому что LAPACK выбирает знаки для столбцов определенным образом, который, кажется, зависит от первых двух точек данных. UЭто, конечно, не ошибка, потому что декомпозиция правильная. Но LAPACK по сути создает эти «артефактные» корреляции, используя свободу присваивать знаки. Эти корреляции не отражают зависимость элементов от ; вместо этого они отражают свободу в решении SVD и конкретное соглашение LAPACK, чтобы исправить это.UU

PS. Поздравляем @whuber за то, что сегодня он получил 100 тыс. Репутации!

амеба говорит восстановить монику
источник
statstat <- function(x) { u <- svd(x)$u; as.vector(sign(runif(1) - 1/2)*u %*% diag(sign(diag(u)))) }U(U11,U22,...,UКК)UU
svdssvdUU
R±2/30.2
1
U
1
Интуитивно, это честно. Как только первая главная ось определена в пространстве, остальные пр. топоры получают уменьшенную свободу. В случае двумерных данных второй (последний) ПК привязан целиком, кроме знака. Я бы скорее назвал это ограничением, а не зависимостью в статистическом смысле.
ttnphns
0

ИксY

Икс2+Y2знак равно1

Соv[Икс,Y]знак равноВaр[ИксY]знак равноЕ[Икс2Y2]-Е[ИксY]2

ИксY

Аксакал
источник
К(К+1)/2UUDUUDК(К-1)/2
U1UNКN>КUNNзнак равно1000Nзнак равно4Икс2+Y2знак равно1U
ИксUYИкс2+Y2знак равно1Cov(Икс,Y)знак равно0Иксзнак равносоз(θ)Yзнак равногрех(θ)θ[0,2π)
UUяJ01/NNN1/Nзнак равно1
1
UИксU11U21ρNρρρзнак равно0
говорит амеба, восстанови Монику