Я наблюдаю очень странное поведение в результате SVD случайных данных, которое я могу воспроизвести как в Matlab, так и в R. Это похоже на некоторую числовую проблему в библиотеке LAPACK; это?
Я рисую выборок из мерного гауссиана с нулевым средним и единичной ковариацией: . Я собрать их в данных матрица . (При желании я могу центрировать или нет, это не влияет на следующее.) Затем я выполняю разложение по сингулярным значениям (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')
источник
Ответы:
Это не ошибка.
Как мы подробно рассмотрели в комментариях, происходит две вещи. Во-первых, столбцы ограничены для соответствия требованиям SVD: каждый должен иметь длину блока и быть ортогональным по отношению ко всем остальным. Просмотр как случайная величина , созданная из случайной матрицы через конкретный алгоритм SVD, мы тем самым отметить , что это функционально независимым ограничения создают статистические зависимости между столбцами .U U X k(k+1)/2 U
Эти зависимости могут быть выявлены в большей или меньшей степени путем изучения корреляций между компонентами , но возникает второе явление : решение SVD не является уникальным. Как минимум, каждый столбец можно независимо отрицать, давая как минимум различных решений с столбцами. Сильные корреляции (превышающие ) могут быть вызваны соответствующим изменением знаков столбцов. (Один из способов сделать это дан в моем первом комментарии к ответу Амебы в этой теме: я принудительно заставляю всеU U 2 к к +1 / 2 у я я , я = 1 , ... , KU 2k k 1/2 uii,i=1,…,k иметь один и тот же знак, что делает их всех отрицательными или положительными с одинаковой вероятностью.) С другой стороны, все корреляции можно сделать равными нулю, выбрав знаки случайным образом, независимо с равными вероятностями. (Я приведу пример ниже в разделе «Редактировать».)
С осторожностью, мы можем частично различить оба эти явления при чтении диаграммы рассеивания матрицы компонентов . Определенные характеристики - такие как появление точек, почти равномерно распределенных в четко определенных круглых областях, - свидетельствуют об отсутствии независимости. Другие, такие как диаграммы рассеяния, показывающие четкие ненулевые корреляции, очевидно, зависят от выбора, сделанного в алгоритме - но такие выборы возможны только из-за отсутствия независимости в первую очередь.U
Окончательный тест алгоритма декомпозиции, такого как SVD (или Cholesky, LR, LU и т. Д.), Заключается в том, выполняет ли он то, что заявляет. В этом случае достаточно проверить, что когда SVD возвращает тройку матриц , то восстанавливается с точностью до ожидаемой ошибки с плавающей запятой произведением ; что столбцы и ортонормированы; и что диагонально, его диагональные элементы неотрицательны и расположены в порядке убывания. Я применил такие тесты к алгоритму в(U,D,V) X UDV′ U V D
svd
R
и никогда не находил это ошибочным. Хотя это и не является гарантией, что это совершенно правильно, такой опыт, который, как я полагаю, разделяют очень многие люди, предполагает, что любая ошибка потребует какого-то необычного вклада для проявления.Далее следует более подробный анализ конкретных вопросов, затронутых в этом вопросе.
Используяk U k = 3к = 3
R
«ssvd
процедуры, сначала вы можете проверить , что в качестве увеличивается, корреляция между коэффициентами ослабевать, но они по - прежнему отличны от нуля. Если бы вы просто выполнили большую симуляцию, вы бы обнаружили, что они значительны. (Когда , должно быть достаточно 50000 итераций.) Вопреки утверждению в вопросе, корреляции не «полностью исчезают».Во-вторых, лучший способ изучить это явление - вернуться к основному вопросу независимости коэффициентов. Хотя корреляции, как правило, близки к нулю в большинстве случаев, отсутствие независимости ясно видно. Это стало наиболее очевидным при изучении полного многомерного распределения коэффициентов . Характер распределения проявляется даже в небольших моделях, в которых ненулевые корреляции не могут (пока) быть обнаружены. Например, рассмотрим матрицу рассеяния для коэффициентов. Чтобы сделать это практически возможным, я установил размер каждого моделируемого набора данных равным и сохранил , тем самым нарисовав реализацийU 4 к = 2 1000 4 × 2 матрицы , создав матрицу . Вот его полная матрица рассеяния с переменными, перечисленными по их позициям в :U 1000 × 8 U
Сканирование вниз по первому столбцу показывает интересное отсутствие независимости между и другим : посмотрите, например, как верхний квадрант диаграммы рассеяния с почти свободен; или исследуйте эллиптическое наклонное вверх облако, описывающее соотношение и наклонное облако вниз для пары . При внимательном рассмотрении обнаруживается явное отсутствие независимости почти всех этих коэффициентов: очень немногие из них выглядят дистанционно независимыми, хотя большинство из них демонстрируют почти нулевую корреляцию.U11 Uя ж U21 ( ты11, у22) ( ты21, у12)
(Примечание: большинство круговых облаков являются проекциями из гиперсферы, созданной условием нормализации, в результате чего сумма квадратов всех компонентов каждого столбца равна единице.)
Матрицы диаграммы рассеяния с и демонстрируют сходные закономерности: эти явления не ограничиваются и не зависят от размера каждого моделируемого набора данных: их просто сложнее генерировать и исследовать.к = 3 к = 4 к = 2
Объяснения для этих шаблонов идут к алгоритму, используемому для получения в разложении по сингулярному значению, но мы знаем, что такие шаблоны несамостоятельности должны существовать с помощью определяющих свойств : поскольку каждый последующий столбец (геометрически) ортогонален предыдущему Кроме того, эти условия ортогональности налагают функциональные зависимости между коэффициентами, которые тем самым переводят в статистические зависимости среди соответствующих случайных величин.U UU
редактировать
В ответ на комментарии, возможно, стоит отметить степень, в которой эти явления зависимости отражают лежащий в основе алгоритм (для вычисления SVD) и насколько они присущи природе процесса.
В конкретных моделях корреляций между коэффициентами зависят много от произвольных выборов , сделанных с помощью алгоритма SVD, потому что решение не является уникальным: столбцы всегда могут быть независимы умноженные на или . Не существует внутреннего способа выбрать знак. Таким образом, когда два алгоритма SVD делают разные (произвольные или, возможно, даже случайные) варианты выбора знака, они могут привести к различным шаблонам диаграмм рассеяния значений . Если вы хотите увидеть это, замените функцию в коде ниже наU - 1 1 ( тыя ж, уя'J')
stat
Это сначала случайным образомUя , дж Uя , дж'
x
переупорядочивает наблюдения , выполняет SVD, затем применяет обратный порядок,u
чтобы соответствовать исходной последовательности наблюдения. Поскольку эффект состоит в формировании смеси отраженных и повернутых версий исходных диаграмм рассеяния, диаграммы рассеяния в матрице будут выглядеть намного более однородными. Все выборочные корреляции будут очень близки к нулю (по построению: лежащие в основе корреляции точно равны нулю). Тем не менее, отсутствие независимости все еще будет очевидным (в однородных круглых формах, которые появляются, особенно между и ).Отсутствие данных в некоторых квадрантах некоторых из исходных диаграмм рассеяния (показано на рисунке выше) возникает из-за того, как
R
алгоритм SVD выбирает знаки для столбцов.В выводах ничего не меняется. Поскольку второй столбец ортогонален первому, он (рассматриваемый как многомерная случайная величина) зависит от первого (также рассматриваемого как многомерная случайная величина). Вы не можете иметь все компоненты одного столбца независимыми от всех компонентов другого; все, что вы можете сделать, это посмотреть на данные таким образом, чтобы скрыть зависимости, но зависимость сохранится.U
Здесь обновленk > 2
R
код для обработки случаев и рисования части матрицы рассеяния.источник
svd(X,0)
наsvds(X)
в моем коде Matlab эффект исчезает! Насколько я знаю, эти две функции используют разные алгоритмы SVD (обе являются подпрограммами LAPACK, но, очевидно, разными). Я не знаю, имеет ли R функцию, аналогичную функции Matlabsvds
, но мне интересно, если вы все еще собираетесь утверждать, что это «реальный» эффект, а не числовая проблема.U
вы случайным образом решите, должен ли каждый из его столбцов оставаться в том виде, в каком он есть, или будет менять свой знак, не исчезнут ли корреляции, о которых вы говорите?stat
S
были в определенном порядке; это вопрос удобства. Другие процедуры гарантируют это (например, MATLABsvds
), но это не является общим требованием. @amoeba: Смотряsvds
(что кажется свободным от этого проблемного поведения), вычисления основаны на фактическом вычислении собственных значений вначале (поэтому он не использует стандартные процедурыdgesdd
/dgesvd
LAPACK - я сильно подозреваю, что он используетdsyevr
/dsyevx
сначала).Этот ответ представляет репликацию результатов @ whuber в Matlab, а также прямую демонстрацию того, что корреляции являются «артефактом» того, как реализация SVD выбирает знак для компонентов.
Учитывая длинную цепочку потенциально запутанных комментариев, я хочу подчеркнуть для будущих читателей, что я полностью согласен со следующим:
Мой вопрос был: почему мы видим высокие корреляции даже для большого числа случайных дро ?∼ 0,2 Nт е р= 1000
Вот репликация примера @ whuber с , и в Matlab:k = 2n = 4 к = 2 Nт е р= 1000
Слева - корреляционная матрица, справа - диаграммы рассеяния, подобные @ whuber's. Соглашение между нашими симуляциями кажется идеальным.
Теперь, следуя гениальному предложению @ttnphns, я назначаю случайные знаки столбцам , то есть после этой строки:U
Я добавляю следующие две строки:
Вот результат:
Все корреляции исчезают, как я и ожидал с самого начала !
Как говорит @whuber, отсутствие независимости можно увидеть в идеальной круглой форме некоторых диаграмм рассеяния (поскольку длина каждого столбца должна равняться , сумма квадратов любых двух элементов не может превышать ). Но корреляции исчезают.11 1
Подводя итог всей проблеме, мы видим, что появляются сильные корреляции, потому что LAPACK выбирает знаки для столбцов определенным образом, который, кажется, зависит от первых двух точек данных.U Это, конечно, не ошибка, потому что декомпозиция правильная. Но LAPACK по сути создает эти «артефактные» корреляции, используя свободу присваивать знаки. Эти корреляции не отражают зависимость элементов от ; вместо этого они отражают свободу в решении SVD и конкретное соглашение LAPACK, чтобы исправить это.UU
PS. Поздравляем @whuber за то, что сегодня он получил 100 тыс. Репутации!
источник
stat
stat <- function(x) { u <- svd(x)$u; as.vector(sign(runif(1) - 1/2)*u %*% diag(sign(diag(u)))) }
svds
svd
R
источник