Многомерное нормальное распределение

Совместное распределение $d$ величин — это в общем случае функция $d$ аргументов, которую не вместит никакая таблица. Одно семейство устроено иначе. Оно полностью описывается вектором средних и ковариационной матрицей, не выходит за пределы семейства при линейных отображениях, маргинализации и переходе к условным распределениям, и всё это — несколько строк линейной алгебры. Это семейство — многомерное нормальное распределение (multivariate Gaussian) — служит моделью шума в линейной регрессии, априорным распределением гауссовских процессов (Gaussian processes), состоянием фильтра Калмана (Kalman filter) и латентным пространством большинства генеративных моделей. В этой заметке мы строим его из одномерного нормального распределения и ковариации из предыдущей заметки.

Обозначения. $X \in \mathbb R^d$ — случайный вектор с вектором средних $\mu = E[X]$ и ковариационной матрицей $\Sigma$, $\Sigma_{ij} = \operatorname{Cov}(X_i, X_j)$, симметричной и положительно определённой. Пишем $X \sim N(\mu, \Sigma)$; $Z \sim N(0, I)$ — стандартный нормальный вектор. $\lvert \Sigma \rvert$ — определитель, $A^\top$ — транспонированная матрица. Заметка использует немного линейной алгебры: произведения матриц, обратные матрицы и определители, квадратичную форму $x^\top A x$ и собственные векторы симметричной матрицы.

Плотность#

Вектор $X \sim N(\mu, \Sigma)$ имеет плотность

$$f(x) = \frac{1}{(2\pi)^{d/2} \lvert \Sigma \rvert^{1/2}} \exp\Big(-\frac12 (x - \mu)^\top \Sigma^{-1} (x - \mu)\Big) .$$

При $d = 1$ это нормальная плотность из заметки о непрерывных распределениях с $\Sigma = \sigma^2$. Квадратичная форма в показателе экспоненты, $\Delta^2(x) = (x - \mu)^\top \Sigma^{-1} (x - \mu)$, — это квадрат расстояния Махаланобиса (Mahalanobis distance): это расстояние от среднего, измеренное в стандартных отклонениях вдоль каждого направления. Плотность постоянна на множествах, где постоянно $\Delta$, а эти множества — эллипсоиды.

Если $\Sigma$ диагональна, квадратичная форма — это сумма $\sum_i (x_i - \mu_i)^2 / \sigma_i^2$, экспонента суммы — произведение, и плотность распадается в произведение нормальных плотностей, по одной на координату. Стандартный нормальный вектор $Z \sim N(0, I)$ — это просто $d$ независимых координат с распределением $N(0, 1)$.

Ковариация как геометрия#

У симметричной положительно определённой матрицы есть ортонормированный базис из собственных векторов $q_1, \dots, q_d$ с положительными собственными значениями $\lambda_1, \dots, \lambda_d$: $\Sigma = Q \Lambda Q^\top$. В координатах этого базиса ковариация диагональна, поэтому плотность в них распадается в произведение, с дисперсией $\lambda_i$ вдоль $q_i$. Линии уровня — эллипсы с осями вдоль собственных векторов и полуосями $k\sqrt{\lambda_i}$ на расстоянии Махаланобиса $k$. Для

$$\Sigma = \begin{pmatrix} 2 & 1.2 \\ 1.2 & 1 \end{pmatrix}, \qquad \lambda = 2.8,\ 0.2, \qquad \rho = \frac{1.2}{\sqrt{2 \cdot 1}} = 0.85 ,$$

длинная ось направлена вдоль $(0.83, 0.55)$:

Рисунок

Ковариационная матрица как эллипс

На длинную ось приходится $2.8 / 3 = 93\%$ полной дисперсии $\operatorname{tr} \Sigma = 3$, и проекция на неё — первая главная компонента метода главных компонент (principal component analysis, PCA). Эллипсы исправляют и одну распространённую интуицию. В одномерном случае интервал в одно стандартное отклонение вокруг среднего содержит 68% вероятности; в двумерном эллипс расстояния Махаланобиса 1 содержит лишь $1 - e^{-1/2} = 39.3\%$, а эллипс расстояния 2 — $86.5\%$, потому что $\Delta^2$ — сумма $d$ квадратов стандартных нормальных величин и имеет распределение хи-квадрат (chi-squared) с $d$ степенями свободы.

Аффинные отображения#

Линейные отображения переводят гауссовские векторы в гауссовские. Если $X \sim N(\mu, \Sigma)$, то для любой матрицы $B$ и вектора $c$

$$BX + c \sim N(B\mu + c,\ B \Sigma B^\top) .$$

Среднее следует из линейности математического ожидания, ковариация — из билинейности ковариации; то, что результат снова гауссовский, для обратимой квадратной матрицы $B$ следует из формулы замены переменных из заметки о функциях от случайных величин, записанной в $d$ измерениях, а для любой $B$ — из характеристической функции; если матрица $B \Sigma B^\top$ вырождена, результат — вырожденное гауссовское распределение, сосредоточенное на подпространстве и не имеющее плотности. Три следствия используются постоянно:

• Сэмплирование. $X = \mu + LZ$ с $LL^\top = \Sigma$ (например, с множителем Холецкого, Cholesky factor) имеет распределение $N(\mu, \Sigma)$; заметка о гауссовских распределениях курса «Генеративные модели» разбирает это на PyTorch.
• Отбеливание (whitening). Обратно, $L^{-1}(X - \mu) \sim N(0, I)$: любое нормальное распределение в подходящих координатах становится стандартным.
• Линейные комбинации. $a^\top X \sim N(a^\top \mu, a^\top \Sigma a)$. Для примера выше (с нулевым средним) $X_1 + X_2 \sim N(0, 2 + 1 + 2 \cdot 1.2) = N(0, 5.4)$.

Маргинальные и условные распределения#

Разобьём $X$ на два блока $X_1$, $X_2$ и соответственно разобьём среднее и ковариацию:

$$\mu = \begin{pmatrix} \mu_1 \\ \mu_2 \end{pmatrix}, \qquad \Sigma = \begin{pmatrix} \Sigma_{11} & \Sigma_{12} \\ \Sigma_{21} & \Sigma_{22} \end{pmatrix} .$$

Маргинальное распределение $X_2$ просто считывается: $X_2 \sim N(\mu_2, \Sigma_{22})$ — частный случай аффинного отображения, которое выбирает координаты. Условное распределение $X_2$ при $X_1 = x_1$ тоже гауссовское, с

$$E[X_2 \mid X_1 = x_1] = \mu_2 + \Sigma_{21} \Sigma_{11}^{-1} (x_1 - \mu_1), \qquad \operatorname{Cov}(X_2 \mid X_1 = x_1) = \Sigma_{22} - \Sigma_{21} \Sigma_{11}^{-1} \Sigma_{12} .$$

Для двух скалярных координат блоки — числа, и формулы принимают вид $E[X_2 \mid X_1 = x_1] = \mu_2 + (\Sigma_{21} / \Sigma_{11})(x_1 - \mu_1)$ и $\operatorname{Var}(X_2 \mid X_1 = x_1) = \Sigma_{22} - \Sigma_{21}^2 / \Sigma_{11}$. В нашем примере $\Sigma_{21}/\Sigma_{11} = 1.2/2 = 0.6$, поэтому при $X_1 = x_1$ вторая координата имеет распределение $N(0.6\, x_1,\ 1 - 1.2^2/2) = N(0.6\, x_1, 0.28)$:

Рисунок

Условное распределение гауссовского вектора

Два свойства этих формул важны далеко за пределами примера. Условное среднее — линейная функция наблюдённого значения, поэтому для совместно гауссовских (jointly Gaussian) данных линейная регрессия — не приближение, а точный наилучший прогноз; это тема заметки об условном математическом ожидании. А условная ковариация — дополнение Шура (Schur complement) $\Sigma_{22} - \Sigma_{21}\Sigma_{11}^{-1}\Sigma_{12}$ — не зависит от наблюдённого значения: наблюдение $X_1$ уменьшает неопределённость относительно $X_2$ на одну и ту же величину, какое бы значение ни приняла $X_1$. Регрессия на основе гауссовских процессов и фильтр Калмана — это те же две формулы, применённые к большим ковариационным матрицам.

Для гауссовских векторов некоррелированность означает независимость#

Для совместно гауссовского вектора нулевая ковариация между двумя блоками делает $\Sigma$ блочно-диагональной, и плотность распадается в произведение: из некоррелированности следует независимость. Это свойство совместного распределения, а не маргинальных. Пусть $X \sim N(0, 1)$ и $Y = SX$, где $S = \pm 1$ — независимый случайный знак, оба значения которого равновероятны. Тогда $Y \sim N(0, 1)$, и $\operatorname{Cov}(X, Y) = E[S]\, E[X^2] = 0$, но $\lvert Y \rvert = \lvert X \rvert$: величины сильно зависимы. Пара не является совместно гауссовской: любая линейная комбинация координат гауссовского вектора нормальна, а $X + Y$ равна 0 с вероятностью $1/2$, чего не бывает ни с одной нормальной величиной.

Высокая размерность: оболочка#

Дополнительный материал: почему сэмплы многомерного нормального распределения избегают центра. Остальная часть заметки на него не опирается.

Плотность $N(0, I)$ максимальна в начале координат, и всё же в высокой размерности сэмплы практически никогда не оказываются рядом с ним. Квадрат длины $\lVert Z \rVert^2 = Z_1^2 + \dots + Z_d^2$ имеет среднее $d$ и дисперсию $2d$, поэтому длина сосредоточена около $\sqrt d$ с разбросом, который остаётся около $0.7$:

Рисунок

Длина стандартного нормального вектора в d измерениях

В размерности 100 длина равна $10 \pm 0.7$, а длина меньше 5 имеет вероятность около $10^{-15}$. Плотность около центра высока, но объёма у центра почти нет; объём, а с ним и вероятность, сосредоточен в тонкой оболочке. Отсюда два практических следствия. Линейная интерполяция между двумя латентными сэмплами генеративной модели проходит через область, которую модель никогда не видела при обучении, — поэтому вместо неё часто используют сферическую интерполяцию (spherical interpolation). А в высокой размерности расстояния между независимыми гауссовскими точками почти одинаковы, около $\sqrt{2d}$, поэтому ближайшие соседи становятся менее информативными.

Что изменилось#

	Совместное распределение общего вида	Многомерное нормальное распределение
Параметры	функция $d$ аргументов	$\mu$ и $\Sigma$: $d + d(d+1)/2$ чисел
Линии уровня	какие угодно	эллипсоиды, оси — собственные векторы $\Sigma$
Линейные отображения, маргинальные распределения	новые распределения, которые нужно вычислять	$N(B\mu + c, B\Sigma B^\top)$; блоки считываются
Условные распределения	перенормированные срезы	гауссовские, среднее линейно, ковариация — дополнение Шура
Некоррелированность	недостаточна для независимости	равносильна независимости

Ограничения#

• Лёгкие хвосты и одна мода. Реальные данные высокой размерности редко бывают гауссовскими: у изображений, эмбеддингов текстов и финансовых доходностей бывают тяжёлые хвосты, кластеры и искривлённая структура; именно поэтому существуют смеси, нормализующие потоки и диффузионные модели.
• Ковариационная матрица обходится дорого. В ней $d(d+1)/2$ свободных элементов, 500 500 при $d = 1\,000$; чтобы её оценить, нужно гораздо больше сэмплов, чем размерностей, а обращение стоит $O(d^3)$. На практике используют диагональную или низкоранговую структуру либо разреженную обратную ковариационную матрицу (sparse precision).
• Нормальность маргинальных распределений ничего не говорит о совместном. Каждая координата может быть нормальной, а вектор — нет, как показывает пример со знаком; формулы для условных распределений и для независимости применимы, только если вектор совместно гауссовский.

Что нужно запомнить#

• $N(\mu, \Sigma)$ имеет плотность $\propto \exp(-\frac12 (x - \mu)^\top \Sigma^{-1} (x - \mu))$; линии уровня — эллипсоиды постоянного расстояния Махаланобиса.
• Собственные векторы $\Sigma$ задают оси, квадратные корни её собственных значений — полуоси; в двумерном случае эллипс расстояния Махаланобиса 1 содержит 39%, а не 68%.
• $BX + c \sim N(B\mu + c, B\Sigma B^\top)$; $X = \mu + LZ$ с $LL^\top = \Sigma$ позволяет сэмплировать любое нормальное распределение.
• Маргинальные распределения задаются блоками $\mu$ и $\Sigma$; условные — гауссовские, со средним, линейным по наблюдению, и ковариацией $\Sigma_{22} - \Sigma_{21}\Sigma_{11}^{-1}\Sigma_{12}$.
• Для совместно гауссовских векторов некоррелированность означает независимость.
• Стандартный нормальный вектор в $d$ измерениях имеет длину около $\sqrt d \pm 0.7$.

Код#

Пример с рисунков, проверенный сэмплированием на NumPy:

import numpy as np
from scipy import stats

S = np.array([[2.0, 1.2], [1.2, 1.0]])
w, V = np.linalg.eigh(S)
print(w.round(3), np.abs(V[:, 1]).round(3))                         # [0.2 2.8] [0.832 0.555]: the long axis
print(round(1.2 / np.sqrt(2), 2), [round(float(stats.chi2(2).cdf(k * k)), 3) for k in (1, 2)])   # 0.85 [0.393, 0.865]

rng = np.random.default_rng(0)
L = np.linalg.cholesky(S)                                            # sampling: X = L Z
X = rng.normal(size=(1_000_000, 2)) @ L.T
print(np.cov(X.T).round(2).tolist())                                 # [[2.0, 1.2], [1.2, 1.0]]
d2 = np.einsum('ij,jk,ik->i', X, np.linalg.inv(S), X)                # squared Mahalanobis distances
print(round(float((d2 <= 1).mean()), 3), round(float((d2 <= 4).mean()), 3))   # 0.393 0.865
print(round(float(X.sum(axis=1).var()), 2))                          # 5.39, against 2 + 1 + 2 * 1.2 = 5.4

# conditioning on X1 = 2: the samples with X1 close to 2
near = X[np.abs(X[:, 0] - 2) < 0.02, 1]
print(round(1.2 / 2 * 2, 2), round(1 - 1.2**2 / 2, 2), round(float(near.mean()), 2), round(float(near.var()), 2))   # 1.2 0.28 1.2 0.28

# normal marginals, zero correlation, and still dependent: Y = S X with a random sign S
x = rng.normal(size=1_000_000)
y = rng.choice([-1, 1], size=x.size) * x
print(round(float(np.corrcoef(x, y)[0, 1]), 3), round(float((x + y == 0).mean()), 3))   # -0.001 0.501

# the shell: the length of a standard normal vector in 100 dimensions
r = np.linalg.norm(rng.normal(size=(100_000, 100)), axis=1)
print(round(float(r.mean()), 2), round(float(r.std()), 2), f'{stats.chi(100).cdf(5):.1e}')   # 9.97 0.71 1.1e-15

Дальше#

Для гауссовских векторов наилучший прогноз $X_2$ по $X_1$ оказался линейным, а оставшаяся дисперсия — известной. Вопрос имеет смысл для любой пары случайных величин: какая функция от $X$ лучше всего предсказывает $Y$ и какая ошибка при этом остаётся? Следующая заметка отвечает на него с помощью условного математического ожидания $E[Y \mid X]$ — функции регрессии, которую приближает модель, обученная с квадратичной функцией потерь; другие функции потерь нацелены на другие характеристики условного распределения, например на медиану или квантиль.

Литература#

• C. M. Bishop. Pattern Recognition and Machine Learning. Springer, 2006. Раздел 2.3 (гауссовское распределение, условные и маргинальные гауссовские распределения).
• M. P. Deisenroth, A. A. Faisal, C. S. Ong. Mathematics for Machine Learning. Cambridge University Press, 2020. Раздел 6.5 (гауссовское распределение).
• J. K. Blitzstein, J. Hwang. Introduction to Probability, 2nd ed. CRC Press, 2019. Раздел 7.5 (многомерное нормальное распределение).
• T. White. Sampling Generative Networks. arXiv:1609.04468, 2016. Сферическая интерполяция в латентных пространствах.