В эпоху глубокого обучения (DL) все больше и больше людей сталкивались и использовали (сознательно или нет) случайные матрицы. В большинстве случаев это использование ограничивается инициализацией весов сетей, что может быть выполнено одной строкой кода в вашей любимой структуре DL. Тем не менее, случайные матрицы имеют богатую математическую теорию с далеко идущими приложениями в физике, теории сетей, машинном обучении, финансах и т. Д. Этот увлекательный диапазон приложений также означает, что каждая область часто разрабатывала свою специальную терминологию для описания одной и той же математической концепции, часто с запутанные последствия. Цель этой статьи (статей) - представить теорию случайных матриц (RMT) с точки зрения физика, пытаясь связать различные инструменты и терминологию, особенно в литературе по математике и физике. Хотя вам действительно нужно знать основы исчисления, я постараюсь представить новые концепции и дать ссылки на случай, если вы захотите углубиться.

Этот пост дополнен численным анализом, который вы можете найти в этом репозитории GitHub в формате блокнота Jupyter. Некоторое содержимое записной книжки также представлено в последнем разделе.

Немного предыстории

Исторически случайные матрицы были введены в статистику Джоном Вишартом (1928) и впервые появились в физике, когда Юджин Вигнер [1] использовал их для моделирования ядер тяжелых атомов (1955). Возможно, вы знаете, что на микроскопическом уровне, где преобладают квантово-механические эффекты, энергетический спектр атома не является непрерывным, а имеет дискретные уровни энергии. Оценка этих уровней для сложных атомов на современных компьютерах требует больших численных затрат, но в то время, когда Вигнер занимался проблемой, это было практически невозможно. Вместо этого он предположил, что эти спектры могут быть сродни собственным значениям случайной матрицы, свойства которой должны зависеть только от ее класса симметрии (подробнее об этом позже). Это пример универсальности: сложные системы, которые различаются в каком-то микроскопическом масштабе, имеют тенденцию иметь схожие свойства в «макроскопическом». При переходе от «микро к макро» многие детали вымываются из описания системы, а оставшиеся релевантные степени свободы, как правило, являются общими для более широкого класса систем, обладающих некоторой определенной симметрией. Эта идея классов симметрии в RMT была развита Фрименом Дайсоном [2]. Существует три классических ансамбля RM: ортогональный, унитарный и симплектический; Дайсон называет их тройственным путем. Если все это звучит немного неясно, не волнуйтесь, мы вернемся к этому позже.

Обновление нескольких ключевых понятий линейной алгебры

Прежде чем углубиться в RMT, позвольте мне резюмировать некоторые основные концепции стандартных матриц. Цель этого раздела - определить терминологию, познакомить с основными инструментами, с которыми мы будем работать, и освежить вашу память на случай, если вы забыли некоторые из этих концепций. Поступая таким образом, я не буду делать никаких попыток дать подробное объяснение, но предполагаю некоторые базовые знания в области исчисления и линейной алгебры. Рассмотрим матрицу M размером N × N; типичная установка проблемы:

где x и b - два вектора: первый - это переменная, для которой мы хотим найти, а второй - «исходный» термин. Вы можете рассматривать эту формулу как компактный способ представить систему связанных линейных уравнений. Типичная задача линейной алгебры - найти собственные значения и собственные векторы этой матрицы путем нахождения корней характеристического уравнения:

где λ - собственные значения, а I - единичная матрица (в уравнениях я буду использовать шляпу для обозначения матрицы). Отсюда мы можем найти матрицу U, которая диагонализирует M, то есть M = U 𝚲 U ^ (- 1), где 𝚲 = diag (λ1, λ2,…, λN). Обратите внимание, что M может иметь прямоугольную форму, т.е. быть матрицей M × N, и в этом случае необходимо вычислить так называемые сингулярные значения; хотя эти прямоугольные матрицы не очень популярны в физике, они широко используются в машинном обучении (входные данные, веса и т. д.) или финансах, но подробнее об этом в следующем посте. N собственных значений M образуют его «спектр», содержащий много важной информации о системе. Например, как упоминалось во введении, они могут описывать уровни энергии квантово-механической системы.

Часто бывает удобно оценить спектр M с помощью метода резольвенты; Преимущество этого метода заключается в том, что он довольно легко обобщается на бесконечномерные пространства и может использоваться для решения неоднородных, линейных или нелинейных дифференциальных уравнений. В этом последнем контексте резольвента известна как функция Грина, и это один из центральных объектов для оценки при выполнении вычислений в квантовой теории поля, см., Например, Ref. [4] для общего введения. Возвращаясь к уравнению [1], мы хотим решить x = M ^ (- 1) b. Этот обратный оператор представлен резольвентой:

Идея состоит в том, что у резольвенты есть полюсы (то есть особенности), расположенные в реальных собственных значениях M. Затем мы можем использовать набор инструментов комплексного анализа, определив z = λ + i ε так, чтобы сингулярности были удалены от действительной оси на бесконечно малую величину ε. Эта процедура обычно известна как регуляризация. (Нормализованный) след резольвенты известен в математической литературе как преобразование Стилтьеса для M, и это очень знакомый объект для физиков:

где мы определили распределение собственных значений ρ (λ), также известное как плотность состояний в квантовой механике. Здесь δ (.) - дельта-функция Дирака, см. [5]. Наша последняя задача - оценить распределение собственных значений резольвенты и, следовательно, матрицы M. Последний кусок информации, который нам нужен, - это то, как связать распределение собственных значений со следом резольвенты, то есть с помощью формулы Племеля, см. [4]:

В случайном случае записи M будут выбраны из распределения, поэтому указанный выше объект необходимо усреднить по всем возможным реализациям M. Как только среднее значение Tr [G (z)] будет оценено, мы продолжим аналитически результат, т. Е. Используем z → λ + i ε, возьмем мнимую часть и предел ε → 0 ^ (+) в конце ( знак + означает, что мы приближаемся к 0 от положительной вещественной оси). В следующем разделе мы пересмотрим несколько ключевых концепций теории вероятностей, которые будут полезны для выполнения нашей задачи.

Напоминание о нескольких ключевых концепциях вероятности

Этот раздел предназначен для освежения некоторых концепций вероятности, которые будут полезны для понимания следующего раздела. Я предполагаю, что вы знакомы с концепцией случайной величины x. В зависимости от области определения x может принимать ряд значений в зависимости от лежащего в основе распределения вероятностей P (x). Скажем, x обозначает положение объекта в пространстве, тогда P (x) будет вероятностью того, что объект окажется в позиции x. Матрица случайных чисел X является расширением этой концепции; вы можете думать об этом как о наборе случайных величин. У нас может быть два объекта, расположенных в позициях x1 и x2; в общем, две позиции могут не быть независимыми, поэтому теперь нам нужно совместное распределение вероятностей для описания нашей проблемы. Если два объекта независимы, тогда:

это одно из фундаментальных свойств, которым удовлетворяют вероятности двух независимых случайных процессов. Подумайте о наиболее распространенных примерах - бросании кости. Вероятность того, что, скажем, число 3 выпадет в одном броске, равна 1/6; вероятность того, что мы получим две подряд тройки, равна 1/6 × 1/6 и так далее. Это неверно, если два события коррелированы, то есть если они могут влиять друг на друга. Таким образом, положения двух объектов могут не быть независимыми, но как только мы наблюдаем один, мы можем спросить, какова вероятность того, что второй объект находится в определенной позиции с учетом того факта, что мы наблюдали его в другом; математически это обозначается как:

где | выражает условное утверждение. Чтобы увидеть, как это связано со случайными матрицами, позвольте мне показать явное представление двух вышеупомянутых уравнений, предполагая, что случайная матрица X является гауссовой, но имея в виду, что эти два вышеупомянутых соотношения действительны в целом. . Чтобы быть явным, X представляет собой матрицу 2 × 2 с унитарной дисперсией:

Обратите внимание, что я использую Tr как сокращение для суммы по всем индексам; это не обычное определение следа (сумма по диагональным элементам), но оно часто используется в физической литературе. Из уравнения (8) видно, что даже если отдельные элементы RM являются i.i.d. (независимо, индивидуально распределенный), сам RM включает в себя определенную степень корреляции из-за наличия недиагональных членов. Если последние равны нулю, мы явно имеем:

С другой стороны, мы можем переписать матричное распределение в терминах его элементов как условную вероятность:

Вторые объекты, которые мне нужно кратко представить, - это производящие функции момента и кумулянта. Первый определяется как:

где x - случайная величина, P (x) - ее распределение вероятностей, а ξ - параметр, используемый для генерации всех моментов. В физике этот параметр имеет разные имена, в зависимости от того, в какой области вы работаете, но обычно идет под именем исходный термин. Формально все моменты функции можно получить, многократно дифференцируя по ξ:

Однако часто бывает полезнее работать с кумулянтами, поскольку они часто упрощают вычисления; Кумулянтная производящая функция определяется как:

Как это связано с RM или статистической механикой? Центральным объектом статистической механики является статистическая сумма распределения Больцмана (см. Ниже). Чтобы получить средние физические величины, мы вводим источник ξ и продифференцируем по нему. Однако, поскольку нам нужна правильно нормализованная величина, мы берем среднее значение логарифма статистической суммы, то есть кумулянтной производящей функции, определенной выше. Вы можете обобщить эту концепцию на матрицы или даже на функционалы, как это делается в квантовой теории поля… но это выходит за рамки данной статьи :)

Собираем все вместе

В этом разделе я дам «грязную» оценку плотности собственных значений. Хотя эта оценка дает правильный результат, она скрывает несколько тонких моментов, о которых я сейчас лишь кратко упомяну. Я буду придерживаться подхода статистической механики, в основном сосредоточенного на расчетах. Кроме того, пока мы сосредоточимся только на простейшем из ансамблей RMT, гауссовском ортогональном ансамбле или GOE для простоты. Все матрицы, принадлежащие этому ансамблю, имеют элементы, которые являются гауссовскими переменными; другое условие - это класс симметрии, при котором эти матрицы инвариантны, инвариантность относительно вращения. Пусть будет O ортогональная матрица, тогда O ^ T = O ^ {- 1}, где верхний индекс «T» означает транспонирование операция. Эта матрица воздействует на элементы матрицы M и вращает их. Например, в двух и трех измерениях эти матрицы имеют известный вид:

где суффикс «z» в трехмерной матрице означает, что вращения выполняются вокруг оси z.

Что означает, что M инвариантен относительно вращения? Повернутая матрица M ' = O M O ^ T и M имеет те же собственные значения, но повернутые собственные векторы w = O v. Итак, спектр M ’ и M одинаков. Первый шаг в вычислении заключается в манипулировании следом резольвенты следующим образом:

Глядя на это выражение, вы можете подумать, что мы сделали шаг назад, а не вперед! Однако идея состоит в том, чтобы переписать след резольвенты как более привычный объект в статистической механике. Рассмотрим следующий m-мерный матричный интеграл Гаусса:

идея состоит в том, чтобы использовать эту идентичность «в обратном порядке», чтобы переписать:

где Z теперь имеет форму статистической суммы статистической механики для конкретной реализации M и меры 𝒟 𝜑 = ∏ d𝜑 / (2 𝜋) ^ (N / 2 ), где произведение превышает i = (1,2…, N). До сих пор мы только что выполнили некоторые формальные манипуляции, чтобы преобразовать наше исходное выражение в форму, которую можно обработать с помощью инструментов статистической механики. Мы также можем идентифицировать F_N (z) = - 1 / N log Z (z) как свободную энергию системы (с обратной температурой β = 1 / T); обратите внимание, что с точки зрения вероятности это (средняя) кумулянтная производящая функция, определенная в уравнении (14). На этом этапе нам нужно усреднить след резольвенты по различным реализациям M, используя инвариантное относительно вращения гауссово распределение вероятностей:

где первый член в правой части - нормировочная константа, и мы использовали единичное стандартное отклонение, масштабированное по N как σ = 1 / √N и нулевое среднее. Это подводит нас к первой (и только для этого поста) тонкости, о которой я упоминал ранее: когда мы берем математическое ожидание следа резольвенты в уравнении (17), нам нужно вычислить ‹log Z_M› _M, где ‹. ›_M обозначает ожидание относительно M. Проблема в следующем: как мы можем оценить математическое ожидание логарифма функции и каково его значение? Подробный ответ на этот вопрос будет дан в следующем посте, а пока позвольте мне изложить здесь несколько фактов без доказательств:

• ‹log Z_M› _M и log ‹Z_M› _M известны, соответственно, как средние значения после закалки и отжига.

  • Для классических ансамблей случайных матриц оказывается, что ‹log Z_M› _M = log ‹Z_M› _M; Представленные здесь результаты действительно дают правильный ответ в рамках более простой схемы расчета. Однако неясно, почему эти два типа средних должны быть равными (или разными) на данный момент.

Далее мы оценим отожженное среднее статистической суммы. Чтобы получить все предварительные факторы, удобно ли записывать все в компонентах и ​​отделить диагональ от недиагональной части. При этом мы используем тот факт, что для симметричной матрицы X мы можем разделить сумму по элементам следующим образом:

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

Обратите внимание на вновь сгенерированный негауссовский член в экспоненте, пропорциональной 𝜑⁴; В физической теории этот термин, по-видимому, описывает взаимодействия какого-то типа, например взаимодействие плотности-плотности между «спинами» на сайте i и j. В общем, этот термин описывает, как колебания 𝜑 (дисперсии) на сайте i соотносятся с колебаниями на сайте j. Чтобы отделить этот член, мы вводим дополнительный параметр, который также играет роль параметра порядка в статистической механике. Вы можете думать об этом как об изменении переменной в интеграле или как об обобщенном преобразовании Хаббарда-Стратоновича в теории конденсированного состояния. Идея состоит в том, чтобы использовать дельта-функцию для обеспечения скалярного «ограничения»:

Прежде чем переходить к расчету, давайте разберемся, что означает этот параметр. Конечно, он переопределяет квадратичную переменную как линейную. Однако его значение состоит в том, чтобы измерить колебания на сайте i; ясно, что это не простая замена переменной! Параметры порядка являются центральными объектами статистической механики и квантовой теории поля, поскольку они используются для изучения переходов между различными фазами системы.

Это определенно выглядит лучше, чем наше первоначальное выражение! Во второй строке уравнения (23) мы использовали фурье-представление дельты Дирака и таким образом ввели скалярный множитель Лагранжа ξ [6]. Это выглядит как изящный трюк, но это невероятно полезный и далеко идущий способ наложения ограничений в статистической механике и квантовой теории поля, где в последней дельта становится функционалом, который можно использовать, например регуляризовать калибровочные теории [7]! Прежде чем продолжить, удобно переопределить множитель Лагранжа следующим образом: η = 2 i ξ / N и выполнить интеграл Гаусса по:

Если вам интересно, член log η происходит от повторного возведения в степень результата гауссова интеграла по 𝜑. Нас интересует решение этих двух интегралов в пределе больших N, т.е. мы ищем непрерывное приближение к плотности собственных значений. В этом пределе вычисление двух интегралов становится особенно простым, поскольку мы можем использовать метод наискорейшего спуска [8]. Если вы знакомы с концепцией теории среднего поля в физике, то самый крутой спуск дает вам именно такой результат (пока просто имейте это в виду, я вернусь к этому в другом посте). Идея состоит в том, чтобы вычислить интеграл в стационарных точках (η, q):

Это дает квадратное уравнение для q и, следовательно, две стационарные точки:

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

Обратите внимание, что множитель N в определении следа резольвенты сокращается с нашим решением, и мы получаем конечный и независимый от N (помните, что это непрерывное приближение) результат. Чтобы получить плотность собственных значений, нам нужно сначала выполнить аналитическое продолжение z → λ + i ε, взять мнимую часть, а затем предел ε → 0 ^ (+). Обратите внимание, что порядок этих операций важен! Сложная часть - это мнимая часть, входящая в квадратный корень; квадратный корень имеет не изолированные полюса, а ветвь в комплексной плоскости. Изящный и простой способ решить эту проблему без использования расширенного комплексного анализа приведен в работе. [3] посредством следующего тождества:

где я определил:

Используя этот результат в определении плотности собственных значений и взяв предел ε → 0 ^ (+), мы получаем:

При ближайшем рассмотрении видно, что если | λ | ›2 спектр равен нулю (это действительная функция!), А конечен при | λ |‹ 2. Поскольку спектр определяется как положительная величина, знак решения должен быть выбран как + для положительного λ и - в противном случае. В итоге мы получаем знаменитый закон полукруга Вигнера:

Так много пота для такой маленькой формулы :) В следующем разделе я построю приведенное выше выражение в сравнении с явной числовой оценкой.

Численные эксперименты

В заключение мы сравниваем результаты с прямой численной оценкой. Подробности можно найти в этом репозитории G ithub, здесь я пройду только основные шаги.

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

Одна из проблем, с которыми вы сталкиваетесь при моделировании, - это конечный размер системы. Помните, что уравнение (31) было получено в пределе больших N. Но насколько велико? Насколько большим должна быть наша симуляция, чтобы соответствовать аналитической формуле? Чтобы ответить на этот вопрос, давайте напишем функцию для выполнения всей оценки:

Псевдокод выглядит следующим образом:

  1. Сгенерируйте экземпляр X матрицы размером N × N путем выборки из нормального распределения с нулевым средним и дисперсией σ = 1 / √N.
  2. Создайте симметричную версию матрицы, используя: Xs = (X + X ^ T) / √2.
  3. Получите собственные значения с помощью numpy.eigvalsh(Xs), обратите внимание, что это численно оптимизированный метод получения собственных значений симметричных матриц.
  4. Повторите n_samples раз, чтобы попробовать различные реализации X.

Для небольшой матрицы 10 × 10 вы найдете что-то такое:

Здесь вы можете увидеть две вещи: первая заключается в том, что хвосты численного распределения превышают область [-2,2] объемного полукруга. Во-вторых, бины не слишком хорошо согласуются с аналитическим выражением. Мы можем повторить оценку для N = 500:

Вы видите, что результаты определенно улучшаются. Мы можем оценить ошибку между аналитическими и численными значениями как функцию размера матрицы. Здесь я вычисляю среднеквадратичную ошибку (RMSE):

В приведенном выше уравнении первая плотность собственных значений - это плотность, вычисленная численно, а вторая - это выражение, которое мы вычислили аналитически, а S - это количество выборок. На графике ниже я отделяю хвосты от основной ошибки:

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

Заключительное слово

В этой статье я представил простое пошаговое введение в теорию случайных матриц и, в частности, в гауссовский ортогональный ансамбль. Представленный мною расчет не является стандартным, который вы найдете в книгах, потому что некоторые из шагов, которые я выполнил, например среднее значение отжига трудно оправдать таким образом. В следующей статье я представлю более стандартный вывод, основанный на так называемой «аналогии с кулоновским газом». Другой стандартный способ оценки закаленного среднего - использование реплик [9], введение и обоснование которых требует отдельного места.

Надеюсь, вам понравился этот первый взгляд на случайные матрицы, и оставляйте комментарии и вопросы в комментариях.

использованная литература

[1] Э. Вигнер, Характеристические векторы матриц с краями бесконечной размерности, Annals of Mathematics. 62 (3): 548–564 (1955).

[2] Ф. Дайсон, Тройной путь. Алгебраическая структура групп и ансамблей симметрии в квантовой механике, Журнал математической физики. 3 (6) (1962).

[3] Г. Ливан, М. Новаес, П. Виво, Введение в теорию и практику случайных матриц, arXiv: 1712.07903v1 (2017).

[4] М. Стоун, П. Голдбарт, Математика для физики, Cambridge University Press (2010).

[5] Если вы раньше сталкивались с функцией Дирака, то знаете, что это действительно распределение; в противном случае предлагаю прочитать этот средний пост.

[6] Если вы не знакомы с использованием множителей Лагранжа для оценки оптимумов системы ограничений, посмотрите этот средний пост.

[7] В контексте калибровочных теорий это известно как метод Фадеева-Попова, а поля множителей лагранжа - как призраки.

[8] Если вы не знакомы с этим методом, загляните сюда; для интегралов, определенных в реальном пространстве, этот метод известен как метод Лапласа, в то время как он имеет разные названия для интегралов в комплексном пространстве, например седловая точка, стационарная фаза и т. д.

[9] М. Мезар, Г. Паризи и М. А. Вирасоро, Теория спинового стекла и не только, World Scientific, Сингапур (1987).

Дополнительные ссылки

Стандартная книга по теории случайных матриц - книга Мехты:

  • М.Л. Мехта, Случайные матрицы (Academic Press, 1967)

Мне было трудно следовать этой книге для начинающих, но она содержит много продвинутого материала. Более удобное введение можно найти в [1]. Еще одно хорошее введение с интересным анализом геометрии, лежащей в основе RMT, можно найти в:

  • Ян В. Федоров, Введение в теорию случайных матриц: гауссовский унитарный ансамбль и не только, arXiv: 0412017v2

Наконец, для математически склонного читателя я рекомендую конспекты лекций Теренса Тао: