Matrix Calculus, часть 1: производная как линейный оператор

Это перевод первых двух глав конспекта курса MIT «Matrix Calculus (for Machine Learning and Beyond)» — матричное дифференциальное исчисление для машинного обучения и не только. Карта всего курса — в путеводителе. Перевод распространяется под лицензией CC BY-NC-SA 4.0 (как материалы курса на MIT OpenCourseWare).

Как читать вместе с курсом ChinaAI

Сначала пройдите уроки 1.2.1 «Производные», 1.2.2 «Правило цепочки» и 1.2.3 «Градиенты». Эта глава показывает то же самое «сверху»: производная — не число и не таблица частных производных, а линейное отображение малых изменений входа в малые изменения выхода. С этой точки зрения формулы для матриц выводятся в две строки. Раздел о стоимости умножения матриц объясняет, почему обучение нейросетей использует именно обратное распространение (урок 3.1.2). Рисунки оригинала описаны словами; блоки «Проверьте в NumPy» и примечания переводчика — наше дополнение (в курсе используется язык Julia).

1. Обзор и мотивация

Для начала: где матричный анализ стоит в каталоге курсов MIT? Там есть 18.01 (анализ функций одной переменной) и 18.02 (векторный анализ), которые студенты MIT обязаны пройти. Но кажется, что эта последовательность обрывается произвольно:

$$\text{скаляр} \to \text{вектор} \to \text{матрица} \to \text{массивы большей размерности?}$$

Ведь именно так эта последовательность устроена во многих языках программирования, включая Julia! Почему анализ должен останавливаться на векторах?

За последнее десятилетие линейная алгебра стала всё важнее во многих областях: машинном обучении, статистике, инженерии и т. д. Она постепенно заняла гораздо бо́льшую часть современного инструментария — сейчас линейная алгебра нужна всем. Поэтому естественно хотеть заниматься анализом на этих массивах высших порядков, и обобщение будет не простым и не очевидным (например, $\frac{d}{dA} A^2 \ne 2A$ для нескалярных матриц $A$).

Шире говоря, дифференцирование и анализ чувствительности (sensitivity analysis) гораздо глубже, чем можно подумать по простым правилам первого-второго семестра. Дифференцирование функций, у которых входы и/или выходы лежат в более сложных векторных пространствах (матрицы, функции и т. д.), — одна часть предмета. Другая — эффективное вычисление производных функций, заданных очень сложными вычислениями, от нейросетей до огромных инженерных симуляций. Это приводит к «сопряжённому» (adjoint) или «обратному» (reverse-mode) дифференцированию, оно же «обратное распространение ошибки» (backpropagation). Автоматическое дифференцирование (AD) программ компиляторами — ещё одна неожиданная тема: компьютер делает нечто совсем иное, чем человек, который сначала выписывает явную символьную формулу, а потом проводит через неё правило цепочки. Это лишь несколько примеров; главное — дифференцирование сложнее, чем вы, возможно, думаете, и эти сложности становятся важны для всё большего числа приложений.

1.1. Приложения (кратко)

Машинное обучение. С ним связано много модных слов: оптимизация параметров, стохастический градиентный спуск, автоматическое дифференцирование, обратное распространение. Во всём этом видна лишь часть того, как матричный анализ применяется в машинном обучении.

Физическое моделирование. Большие физические симуляции, например в инженерном проектировании, всё чаще описываются огромным числом параметров, и производные результатов симуляции по этим параметрам нужны, чтобы оценить чувствительность к неопределённостям и применить крупномасштабную оптимизацию. Форму крыла самолёта можно задать тысячами параметров; если вычислить производную силы сопротивления (из большой симуляции потока) по этим параметрам, можно оптимизировать форму крыла. Крайний случай — «топологическая оптимизация», где материал в «каждой точке» пространства — степень свободы.

Анализ данных и многомерная статистика. В многомерной статистике модели часто формулируются через матричные входы и выходы. Например, «простая» линейная многомерная модель $Y(X) = XB + U$, где $B$ — неизвестная матрица коэффициентов, а $U$ — неизвестная матрица случайного шума. Регрессия сводится к минимизации некоторой функции ошибки $U(B) = Y - XB$, например матричной нормы $\Vert U \Vert_F^2 = \operatorname{tr} U^T U$, определителя $\det U^T U$ или чего-то сложнее. Для этого нужно дифференцировать такие функции по $B$.

Автоматическое дифференцирование. Обычные курсы анализа построены на символьных вычислениях: студенты по сути учатся делать то, что умеют Mathematica или Wolfram Alpha. Сегодняшнее автоматическое дифференцирование — ни то, ни другое (ни символьные формулы, ни конечные разности). Оно ближе к теории компиляторов, чем к математике. Но математика, лежащая в его основе, интересна, и мы её изучим!

Даже приближённое дифференцирование на компьютере сложнее, чем кажется. Для функции одной переменной производная определяется как предел разности $[f(x + \delta x) - f(x)]/\delta x$ при $\delta x \to 0$. Грубая аппроксимация «конечной разностью» просто берёт эту формулу при малом $\delta x$, но это поднимает много интересных вопросов: баланс ошибок усечения и округления, аппроксимации высших порядков, численная экстраполяция.

1.2. Первые производные

Производная функции одной переменной сама является функцией одной переменной — грубо говоря, это линеаризация функции. То есть она имеет вид $(f(x) - f(x_0)) \approx f'(x_0)(x - x_0)$. В этом смысле со скалярными функциями скаляров (функциями, которые принимают одно число и выдают одно число) «всё просто».

Для этой линеаризации встречаются разные обозначения:

• $\delta y \approx f'(x) \delta x$,
• $dy = f'(x)\, dx$,
• $(y - y_0) \approx f'(x_0)(x - x_0)$,
• $df = f'(x)\, dx$.

В этом курсе мы предпочитаем последнее. Можно думать о $dx$ и $dy$ как о «совсем маленьких числах». В математике их называют бесконечно малыми (infinitesimals) и строго определяют через пределы. Заметьте: здесь мы не хотим делить на $dx$. Со скалярами это совершенно нормально, но когда дойдём до векторов и матриц, делить можно будет не всегда!

С числами таких производных легко поиграть. Возьмём $f(x) = x^2$ и точку $(x_0, f(x_0)) = (3, 9)$. Рядом с $(3, 9)$ получаем такие значения:

$$\begin{aligned} f(3.\mathbf{0001}) &= 9.\mathbf{0006}0001 \\ f(3.\mathbf{00001}) &= 9.\mathbf{00006}00001 \\ f(3.\mathbf{000001}) &= 9.\mathbf{000006}000001 \\ f(3.\mathbf{0000001}) &= 9.\mathbf{0000006}0000001. \end{aligned}$$

Здесь выделенные цифры слева — это $\Delta x$, а выделенные цифры справа — $\Delta y$. Заметьте, что $\Delta y = 6 \Delta x$. Следовательно,

$$f(3 + \Delta x) = 9 + \Delta y = 9 + 6\Delta x \implies f(3 + \Delta x) - f(3) = 6\Delta x \approx f'(3)\Delta x.$$

Значит, линеаризация $x^2$ в точке $x = 3$ — это функция $f(x) - f(3) \approx 6(x - 3)$.

---

Теперь мы покидаем мир скалярного анализа и входим в мир векторного и матричного! Профессор Эдельман предлагает думать о матрицах целостно — не просто как о таблице чисел.

Понятие линеаризации функции концептуально переносится на производные функций, которые принимают или выдают больше одного числа. Конечно, у производной тогда будет другая «форма», чем у одного числа. Вот таблица формы первой производной. Входы функции — по строкам, выходы — по столбцам.

вход ↓, выход →	скаляр	вектор	матрица
скаляр	скаляр	вектор (например, скорость)	матрица
вектор	градиент = вектор (столбец)	матрица (матрица Якоби)	массив высшего порядка
матрица	матрица	массив высшего порядка	массив высшего порядка

В этом курсе вы научитесь подробно разбираться с каждым случаем! Цель таблицы — посеять мысль, что дифференциал — это линеаризация. Рассмотрим пример.

Пример. Пусть $f(x) = x^T x$, где $x$ — матрица $2 \times 1$, и выход, значит, — матрица $1 \times 1$. Убедитесь, что $2x_0^T dx$ действительно дифференциал $f$ в точке $x_0 = \begin{pmatrix} 3 & 4 \end{pmatrix}^T$.

Сначала вычислим $f(x_0)$:

$$f(x_0) = x_0^T x_0 = 3^2 + 4^2 = 25.$$

Пусть $dx = [0{,}001;\ 0{,}002]$. Тогда

$$f(x_0 + dx) = (3{,}001)^2 + (4{,}002)^2 = 25{,}\mathbf{022}005.$$

Заметим, что $2x_0^T dx = 2 \begin{pmatrix} 3 & 4 \end{pmatrix} dx = 0{,}022$. Следовательно,

$$f(x_0 + dx) - f(x_0) \approx 2x_0^T dx = 0{,}022.$$

Как мы сейчас увидим, $2x_0^T dx$ взялось не из ниоткуда!

1.3. Введение: правило произведения для матриц и векторов

Для матриц правило произведения на самом деле тоже есть! Подробно мы обсудим его в следующих главах, а здесь — небольшая проба.

Теорема (правило произведения для дифференциалов). Пусть $A$, $B$ — две матрицы. Тогда для $AB$ выполняется

$$d(AB) = (dA)B + A(dB).$$

Под дифференциалом матрицы $A$ мы понимаем малое (ничем не ограниченное) изменение матрицы $A$. Позже на допустимые возмущения можно будет наложить ограничения.

Заметьте, однако, что (по нашей таблице) производная матрицы — матрица! Так что в общем случае множители переставлять нельзя.

Если $x$ — вектор, то по правилу произведения

$$d(x^T x) = (dx^T)x + x^T(dx).$$

Но это скалярное произведение, а скалярные произведения коммутируют (поскольку $\sum a_i \cdot b_i = \sum b_i \cdot a_i$), так что

$$d(x^T x) = (2x)^T dx.$$

Замечание. Для векторов как матриц правило произведения работает так, что транспонирование «едет вместе» с множителем. См. пример ниже.

Пример. По правилу произведения

1. $d(u^T v) = (du)^T v + u^T(dv) = v^T du + u^T dv$, поскольку скалярные произведения коммутируют;
2. $d(uv^T) = (du)v^T + u(dv)^T$.

Как доказывать такие утверждения, показано в главе 2.

2. Производные как линейные операторы

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

2.1. Снова анализ функций одной переменной

Рисунок 1 оригинала: суть производной — линеаризация: предсказать малое изменение $\delta f$ выхода $f(x)$ по малому изменению $\delta x$ входа $x$ с точностью до первого порядка по $\delta x$. На рисунке — график кривой $f(x)$, точка $x$ и сдвинутая точка $x + \delta x$; касательная в точке $x$ предсказывает изменение $f'(x)\,\delta x$, а настоящее изменение $\delta f$ отличается от него на малую величину более высокого порядка.

В первом семестре анализа (вроде 18.01 в MIT) производную $f'(x)$ вводят как наклон касательной в точке $(x, f(x))$, что можно рассматривать и как линейное приближение $f$ около $x$. Как показано на рисунке 1, это равносильно тому, чтобы предсказать изменение $\delta f$ «выхода» $f(x)$ по малому изменению $\delta x$ «входа» с точностью до первого (линейного) порядка по $\delta x$:

$$\delta f = f(x + \delta x) - f(x) = f'(x)\, \delta x + \underbrace{(\text{члены высших порядков})}_{o(\delta x)}.$$

Точнее записать члены высших порядков можно асимптотическим обозначением «о-малое» $o(\delta x)$: так обозначают любую функцию, модуль которой убывает намного быстрее, чем $|\delta x|$, при $\delta x \to 0$, так что при достаточно малых $\delta x$ она пренебрежимо мала по сравнению с линейным членом $f'(x)\,\delta x$. (Варианты этого обозначения часто используются в информатике; формальное определение мы опускаем.¹) Примеры таких членов высших порядков: $(\delta x)^2$, $(\delta x)^3$, $(\delta x)^{1{,}001}$ и $(\delta x)/\log(\delta x)$.

¹ Коротко: функция $g(\delta x)$ есть $o(\delta x)$, если $\lim_{\delta x \to 0} \frac{\Vert g(\delta x) \Vert}{\Vert \delta x \Vert} = 0$. Мы вернёмся к этому в главе 5.

Замечание. Здесь $\delta x$ — не бесконечно малая, а просто маленькое число. Обратите внимание: наш символ «$\delta$» (греческая строчная «дельта») — не то же самое, что символ «$\partial$», которым обычно обозначают частные производные.

Это понятие производной может напомнить вам первые два члена ряда Тейлора $f(x + \delta x) = f(x) + f'(x)\,\delta x + \cdots$ (хотя на самом деле оно гораздо базовее ряда Тейлора!), и это обозначение хорошо обобщается на высшие размерности и другие векторные пространства. В дифференциальных обозначениях та же идея выражается так:

$$df = f(x + dx) - f(x) = f'(x)\, dx.$$

В этой записи мы неявно отбрасываем член $o(\delta x)$, который исчезает в пределе, когда $\delta x$ становится бесконечно малым.

Это и будет нашим обобщённым определением производной. В такой формулировке мы избегаем деления на $dx$, потому что скоро позволим $x$ (а значит, и $dx$) быть чем-то кроме числа — а если $dx$ — вектор, разделить на него мы просто не сможем!

2.2. Линейные операторы

С точки зрения линейной алгебры производной функции $f$ мы считаем линейный оператор $f'(x)$, такой что

$$df = f(x + dx) - f(x) = f'(x)[dx].$$

Как и выше, дифференциал $dx$ обозначает произвольное малое изменение $x$, и мы неявно отбрасываем все члены $o(dx)$, то есть убывающие быстрее линейных при $dx \to 0$. Часто мы будем опускать квадратные скобки и писать просто $f'(x)\,dx$ вместо $f'(x)[dx]$, но понимать это нужно как линейный оператор $f'(x)$, действующий на $dx$. Не пишите $dx\,f'(x)$ — в общем случае это бессмыслица!

Такое определение позволит распространить дифференцирование на произвольные векторные пространства входов $x$ и выходов $f(x)$. (Технически нам понадобятся векторные пространства с нормой $\Vert x \Vert$ — «банаховы пространства», — чтобы точно определить отбрасываемые члены $o(\delta x)$. Мы вернёмся к банаховым пространствам позже.)

Напоминание (векторное пространство). Нестрого: векторное пространство (над $\mathbb{R}$) — это множество элементов, для которых определены сложение и вычитание, а также умножение на вещественные числа. Например, перемножать произвольные векторы из $\mathbb{R}^n$ бессмысленно, но складывать их и масштабировать на постоянный множитель — безусловно можно.

Примеры векторных пространств:

• $\mathbb{R}^n$, как описано выше. Шире — $\mathbb{R}^{n \times m}$, пространство вещественных матриц $n \times m$. Снова заметьте: если $n \ne m$, умножение между элементами не определено.
• $C^0(\mathbb{R}^n)$ — множество непрерывных функций на $\mathbb{R}^n$ с поточечным сложением.

Напоминание (линейный оператор). Линейный оператор — это отображение $L$ вектора $v$ из векторного пространства $V$ в вектор $L[v]$ (иногда пишут просто $Lv$) в каком-то другом векторном пространстве. Точнее, $L$ линеен, если

$$L[v_1 + v_2] = Lv_1 + Lv_2 \quad \text{и} \quad L[\alpha v] = \alpha L[v]$$

для скаляров $\alpha \in \mathbb{R}$.

Замечание. В этом курсе $f'$ — это отображение, которое принимает $x$ и выдаёт линейный оператор $f'(x)$ (производную $f$ в точке $x$). В свою очередь, $f'(x)$ — линейное отображение, которое принимает направление $v$ и выдаёт вектор $f'(x)[v]$ (позже мы истолкуем его как производную по направлению, см. раздел 2.2.1). Когда направление $v$ — бесконечно малое $dx$, выход $f'(x)[dx] = df$ — это дифференциал $f$ (соответствующее бесконечно малое изменение $f$).

Обозначения производных. В ходу несколько обозначений производных и несколько родственных понятий: производная, дифференцирование, дифференциалы. Таблица ниже их сводит; обозначения, принятые в курсе, в оригинале взяты в рамку — здесь они выделены жирным.

Название	Обозначения	Пояснение
производная	$\boldsymbol{f'}$, а также $\frac{df}{dx}$, $Df$, $f_x$, $\partial_x f$, …	линейный оператор $f'(x)$, который переводит малое изменение входа $dx$ в малое изменение выхода $df = f'(x)[dx]$. В анализе одной переменной этот оператор задаётся одним числом — «наклоном»: если $f(x) = \sin x$, то $f'(x) = \cos x$ — число, на которое умножаем $dx$, чтобы получить $dy = \cos(x)\, dx$. В анализе многих переменных оператор $f'(x)$ задаётся матрицей — якобианом $J$, так что $df = f'(x)[dx] = J\,dx$. Но мы увидим, что выражать $f'$ матрицей не всегда удобно, даже когда это возможно
дифференцирование	$\boldsymbol{{}'}$, $\frac{d}{dx}$, $D$, …	линейный оператор, который переводит функцию $f$ в её производную $f'$
разность	$\boldsymbol{\delta x}$ и $\boldsymbol{\delta f} = f(x + \delta x) - f(x)$	малое (но не бесконечно малое) изменение входа $x$ и выхода $f$ соответственно (зависит от $x$ и $\delta x$): элемент векторного пространства, а не линейный оператор
дифференциал	$\boldsymbol{dx}$ и $\boldsymbol{df} = f(x + dx) - f(x)$	сколь угодно малое («бесконечно малое» — члены высших порядков отброшены) изменение входа $x$ и выхода $f$: элемент векторного пространства, а не линейный оператор2
градиент	$\boldsymbol{\nabla f}$	вектор, скалярное произведение которого $df = \langle \nabla f, dx \rangle$ с малым изменением входа $dx$ даёт малое изменение выхода $df$. «Транспонированная производная»: $\nabla f = (f')^T$ (см. раздел 2.3)
частная производная	$\boldsymbol{\frac{\partial f}{\partial x}}$, $f_x$, $\partial_x f$	линейный оператор, который переводит малое изменение $dx$ одного аргумента функции многих аргументов в соответствующее изменение выхода; например, для $f(x, y)$ имеем $df = \frac{\partial f}{\partial x}[dx] + \frac{\partial f}{\partial y}[dy]$

² Нестрого можно считать, что бесконечно малые $dx$ живут в том же пространстве, что и $x$ (малое изменение вектора — всё равно вектор). Формально «пространство бесконечно малых» определяют по-разному, например как кокасательное пространство в дифференциальной геометрии, но подробнее мы на этом останавливаться не будем.

Примеры линейных операторов:

• Умножение на скаляр $\alpha$, то есть $Lv = \alpha v$. Также умножение векторов-столбцов $v$ на матрицы $A$, то есть $Lv = Av$.
• Некоторые функции, вроде $f(x) = x^2$, очевидно нелинейны. А как насчёт $f(x) = x + 1$? На графике она выглядит линейной, но это не линейная операция, потому что $f(2x) = 2x + 1 \ne 2f(x)$. Такие функции — линейные плюс ненулевая константа — называются аффинными.
• Есть много других линейных операций, которые неудобно или непросто записать как произведение матрицы на вектор. Например, если $A$ — матрица $3 \times 3$, то $L[A] = AB + CA$ — линейный оператор при заданных матрицах $3 \times 3$ $B$ и $C$. Транспонирование $f(x) = x^T$ вектора-столбца $x$ линейно, но не задаётся никакой матрицей, умноженной на $x$. А если рассматривать векторные пространства функций, то операции анализа — дифференцирование и интегрирование — тоже линейные операторы!

2.2.1. Производные по направлению

Есть равносильный способ истолковать производную как линейный оператор, который вы, возможно, видели в анализе многих переменных: как производную по направлению.

Если у нас есть функция $f(x)$ от произвольных векторов $x$, то производная в точке $x$ по направлению (вектору) $v$ определяется так:

$$\left. \frac{\partial}{\partial \alpha} f(x + \alpha v) \right|_{\alpha = 0} = \lim_{\delta\alpha \to 0} \frac{f(x + \delta\alpha\, v) - f(x)}{\delta\alpha},$$

где $\alpha$ — скаляр. Это возвращает производные из произвольных векторных пространств обратно в анализ одной переменной. Величина измеряет скорость изменения $f$ из точки $x$ в направлении $v$. Но оказывается, что она очень просто связана с нашим линейным оператором $f'(x)$, потому что (отбрасывая члены высших порядков в пределе $\delta\alpha \to 0$):

$$f(x + \underbrace{d\alpha\, v}_{dx}) - f(x) = f'(x)[dx] = d\alpha\, f'(x)[v],$$

где на последнем шаге мы вынесли скаляр $d\alpha$, потому что $f'(x)$ — линейный оператор. Сравнивая с формулой выше, сразу получаем, что производная по направлению равна

$$\boxed{\left. \frac{\partial}{\partial \alpha} f(x + \alpha v) \right|_{\alpha = 0} = f'(x)[v]}.$$

Это в точности наш прежний $f'(x)$! (Это можно увидеть и как частный случай правила цепочки из раздела 2.5.) Один урок отсюда: на вход $f'(x)[v]$ вполне можно подавать произвольный не бесконечно малый вектор $v$ — результат будет не $df$, а просто производная по направлению.

2.3. Снова анализ многих переменных, часть 1: скалярные функции

Рисунок 2 оригинала: для вещественной $f(x)$ градиент $\nabla f$ определён так, что указывает направление «в гору» в точке $x$, перпендикулярное линиям уровня $f$. Он не обязательно смотрит прямо на ближайший локальный максимум (если линии уровня не окружности), но «идти в гору» — отправная точка многих алгоритмов численной оптимизации при поиске максимума. На рисунке — линии уровня функции на плоскости и стрелка градиента, перпендикулярная линии уровня.

Пусть $f$ — скалярная функция, которая принимает векторы-«столбцы» $x \in \mathbb{R}^n$ и выдаёт скаляр (из $\mathbb{R}$). Тогда

$$df = f(x + dx) - f(x) = f'(x)[dx] = \text{скаляр}.$$

Поскольку $dx$ — вектор-столбец (в произвольном направлении, произвольное малое изменение $x$), линейный оператор $f'(x)$, который выдаёт скаляр $df$, должен быть вектором-строкой («матрицей из одной строки», или, формальнее, тем, что называют ковектором, «двойственным» вектором или линейной формой)! Эту строку мы называем транспонированным градиентом $(\nabla f)^T$, так что $df$ — это скалярное произведение $dx$ на градиент:

$$df = \nabla f \cdot dx = \underbrace{(\nabla f)^T}_{f'(x)} dx, \quad \text{где} \quad dx = \begin{pmatrix} dx_1 \\ dx_2 \\ \vdots \\ dx_n \end{pmatrix}.$$

Некоторые авторы считают градиент вектором-строкой (отождествляя его с $f'$ или якобианом), но считать его «вектором-столбцом» (транспонированной $f'$), как в этом курсе, — распространённый и полезный выбор. Как вектор-столбец градиент можно рассматривать как направление «в гору» (наискорейшего подъёма) в пространстве $x$, как на рисунке 2. Кроме того, такой градиент проще обобщить на скалярные функции в других векторных пространствах. Как бы то ни было, в этом курсе мы всегда определяем $\nabla f$ так, чтобы он имел ту же «форму», что и $x$, — тогда $df$ есть скалярное произведение $dx$ на градиент.

Это полностью согласуется с тем, что вы, возможно, помните из анализа многих переменных, где градиент — вектор из компонент

$$\nabla f = \begin{pmatrix} \frac{\partial f}{\partial x_1} \\ \frac{\partial f}{\partial x_2} \\ \vdots \\ \frac{\partial f}{\partial x_n} \end{pmatrix},$$

или, равносильно,

$$df = f(x + dx) - f(x) = \nabla f \cdot dx = \frac{\partial f}{\partial x_1} dx_1 + \frac{\partial f}{\partial x_2} dx_2 + \cdots + \frac{\partial f}{\partial x_n} dx_n.$$

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

Посмотрим на примере, как вычислять такой дифференциал.

Пример. Пусть $f(x) = x^T A x$, где $x \in \mathbb{R}^n$, а $A$ — квадратная матрица $n \times n$, так что $f(x) \in \mathbb{R}$. Найдите $df$, $f'(x)$ и $\nabla f$.

Сделаем это прямо по определению:

$$\begin{aligned} df &= f(x + dx) - f(x) \\ &= (x + dx)^T A (x + dx) - x^T A x \\ &= \cancel{x^T A x} + dx^T A x + x^T A\, dx + \underbrace{dx^T A\, dx}_{\text{высший порядок}} - \cancel{x^T A x} \\ &= \underbrace{x^T (A + A^T)}_{f'(x) = (\nabla f)^T} dx \implies \nabla f = (A + A^T)x. \end{aligned}$$

Здесь мы отбросили члены с более чем одним множителем $dx$, поскольку они асимптотически пренебрежимы. Ещё один приём — объединить $dx^T A x$ и $x^T A\, dx$, заметив, что это скаляры, а значит, равны своему транспонированию: $dx^T A x = (dx^T A x)^T = x^T A^T dx$. Итак, мы нашли $f'(x) = x^T (A + A^T) = (\nabla f)^T$, или, равносильно, $\nabla f = [x^T (A + A^T)]^T = (A + A^T)x$.

Разумеется, тот же градиент можно найти и покомпонентно, как вас, вероятно, учили в анализе многих переменных. Сначала пришлось бы явно записать $f(x)$ через компоненты $x$: $f(x) = x^T A x = \sum_{i,j} x_i A_{i,j} x_j$. Затем — найти $\partial f / \partial x_k$ для каждого $k$, не забыв, что $x$ входит в сумму дважды. Но такой подход неуклюж, подвержен ошибкам, трудоёмок и быстро становится ещё хуже, когда мы переходим к более сложным функциям. Мы считаем, что гораздо лучше привыкнуть обращаться с векторами и матрицами целиком, а не как с простыми наборами чисел.

2.4. Снова анализ многих переменных, часть 2: векторные функции

Теперь пусть $f$ — векторная функция: принимает векторы $x \in \mathbb{R}^n$ и выдаёт векторы $f(x) \in \mathbb{R}^m$. Тогда $df$ — вектор-столбец из $m$ компонент, $dx$ — вектор-столбец из $n$ компонент, и мы должны получить линейный оператор $f'(x)$, для которого

$$\underbrace{df}_{m \text{ компонент}} = \underbrace{f'(x)}_{m \times n} \underbrace{dx}_{n \text{ компонент}},$$

так что $f'(x)$ должна быть матрицей $m \times n$ — она называется матрицей Якоби (якобианом) $f$!

Матрица Якоби $J$ представляет линейный оператор, который переводит $dx$ в $df$:

$$df = J\, dx.$$

Элементы матрицы $J$ — это $J_{ij} = \frac{\partial f_i}{\partial x_j}$ ($i$-я строка, $j$-й столбец $J$).

Пусть, например, $f: \mathbb{R}^2 \to \mathbb{R}^2$. Разберёмся, как вычислить дифференциал $f$:

$$df = \begin{pmatrix} \frac{\partial f_1}{\partial x_1} & \frac{\partial f_1}{\partial x_2} \\ \frac{\partial f_2}{\partial x_1} & \frac{\partial f_2}{\partial x_2} \end{pmatrix} \begin{pmatrix} dx_1 \\ dx_2 \end{pmatrix} = \begin{pmatrix} \frac{\partial f_1}{\partial x_1} dx_1 + \frac{\partial f_1}{\partial x_2} dx_2 \\ \frac{\partial f_2}{\partial x_1} dx_1 + \frac{\partial f_2}{\partial x_2} dx_2 \end{pmatrix}.$$

Пример. Рассмотрим функцию $f(x) = Ax$, где $A$ — постоянная матрица $m \times n$. Применяя дистрибутивность умножения матрицы на вектор, получаем

$$df = f(x + dx) - f(x) = A(x + dx) - Ax = \cancel{Ax} + A\,dx - \cancel{Ax} = A\,dx = f'(x)\,dx.$$

Следовательно, $f'(x) = A$.

Заметьте: линейный оператор $A$ — сам себе матрица Якоби!

Рассмотрим теперь некоторые правила дифференцирования.

• Правило суммы. Для $f(x) = g(x) + h(x)$ получаем

  $$df = dg + dh \implies f'(x)\,dx = g'(x)\,dx + h'(x)\,dx.$$

  Значит, $f' = g' + h'$, как и следовало ожидать. Это линейный оператор $f'(x)[v] = g'(x)[v] + h'(x)[v]$; заметьте, что линейные операторы (вроде $g'$ и $h'$) можно складывать так же, как матрицы! В этом смысле линейные операторы сами образуют векторное пространство.

• Правило произведения. Пусть $f(x) = g(x)h(x)$. Тогда

  $$\begin{aligned} df &= f(x + dx) - f(x) \\ &= g(x + dx)h(x + dx) - g(x)h(x) \\ &= (g(x) + \underbrace{g'(x)\,dx}_{dg})(h(x) + \underbrace{h'(x)\,dx}_{dh}) - g(x)h(x) \\ &= gh + dg\,h + g\,dh + \underbrace{dg\,dh}_{\to 0} - gh \\ &= dg\,h + g\,dh, \end{aligned}$$

  где член $dg\,dh$ — высшего порядка и потому отброшен в записи через бесконечно малые. Заметьте, как обычно, что $dg$ и $h$ могут не коммутировать, поскольку они больше не обязательно скаляры!

Несколько коротких примеров того, как удобно применять правило произведения.

Пример. Пусть $f(x) = Ax$ (отображает $\mathbb{R}^n \to \mathbb{R}^m$), где $A$ — постоянная матрица $m \times n$. Тогда

$$df = d(Ax) = \underbrace{dA}_{=0}\, x + A\,dx = A\,dx \implies f'(x) = A.$$

Здесь $dA = 0$, потому что $A$ не меняется, когда мы меняем $x$.

Пример. Пусть $f(x) = x^T A x$ (отображает $\mathbb{R}^n \to \mathbb{R}$). Тогда

$$df = dx^T (Ax) + x^T d(Ax) = \underbrace{dx^T A x}_{=\, x^T A^T dx} + x^T A\,dx = x^T (A + A^T)\,dx = (\nabla f)^T dx,$$

и, значит, $f'(x) = x^T (A + A^T)$. (В частом случае симметричной $A$ это упрощается до $f'(x) = 2x^T A$.) Заметьте, что здесь мы применили предыдущий пример, вычисляя $d(Ax) = A\,dx$. Кроме того, $f$ — скалярная функция, так что можно получить и градиент $\nabla f = (A + A^T)x$, как раньше (для симметричной $A$ это упрощается до $2Ax$).

Пример (поэлементные произведения). Для $x, y \in \mathbb{R}^m$ определим

$$x \mathbin{.*} y = \begin{pmatrix} x_1 y_1 \\ \vdots \\ x_m y_m \end{pmatrix} = \underbrace{\begin{pmatrix} x_1 & & & \\ & x_2 & & \\ & & \ddots & \\ & & & x_m \end{pmatrix}}_{\operatorname{diag}(x)} \begin{pmatrix} y_1 \\ \vdots \\ y_m \end{pmatrix}$$

— поэлементное произведение векторов (его также называют произведением Адамара); для удобства мы определили $\operatorname{diag}(x)$ как диагональную матрицу $m \times m$ с $x$ на диагонали. Далее для $A \in \mathbb{R}^{m \times n}$ определим $f: \mathbb{R}^n \to \mathbb{R}^m$ как

$$f(x) = A(x \mathbin{.*} x).$$

В качестве упражнения проверьте:

• (a) $x \mathbin{.*} y = y \mathbin{.*} x$;
• (b) $A(x \mathbin{.*} y) = A \operatorname{diag}(x)\, y$;
• (c) $d(x \mathbin{.*} y) = (dx) \mathbin{.*} y + x \mathbin{.*} (dy)$. Так что если взять $y$ постоянным и определить $g(x) = y \mathbin{.*} x$, её матрица Якоби — $\operatorname{diag}(y)$;
• (d) $df = A(2x \mathbin{.*} dx) = 2A \operatorname{diag}(x)\, dx = f'(x)[dx]$, так что матрица Якоби — $J = 2A \operatorname{diag}(x)$;
• (e) производная $f$ в точке $x$ по направлению $v$ (раздел 2.2.1) равна просто $f'(x)[v] = 2A(x \mathbin{.*} v)$. Можно численно проверить для произвольных $A$, $x$, $v$, что $f(x + 10^{-8} v) - f(x) \approx 10^{-8}\,(2A(x \mathbin{.*} v))$.

Примечание переводчика: .* — это поэлементное умножение в Julia и MATLAB; в NumPy и PyTorch это просто x * y.

2.5. Правило цепочки

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

• Правило цепочки. Пусть $f(x) = g(h(x))$. Тогда

  $$\begin{aligned} df = f(x + dx) - f(x) &= g(h(x + dx)) - g(h(x)) \\ &= g'(h(x))[h(x + dx) - h(x)] \\ &= g'(h(x))[h'(x)[dx]] \\ &= g'(h(x))\, h'(x)[dx], \end{aligned}$$

  где $g'(h(x))h'(x)$ — композиция $g'$ и $h'$ как матриц.

  Иными словами, $f'(x) = g'(h(x))h'(x)$: якобиан (линейный оператор) $f'$ — это просто произведение (композиция) якобианов $g'h'$. Порядок важен, потому что линейные операторы в общем случае не коммутируют: слева направо = от выходов к входам.

Посмотрим внимательнее на формы этих матриц Якоби в примере, где каждая функция переводит вектор-столбец в вектор-столбец.

Пример. Пусть $x \in \mathbb{R}^n$, $h(x) \in \mathbb{R}^p$ и $g(h(x)) \in \mathbb{R}^m$. Пусть $f(x) = g(h(x))$ отображает $\mathbb{R}^n$ в $\mathbb{R}^m$. Тогда правило цепочки утверждает, что

$$f'(x) = g'(h(x))\,h'(x),$$

и это осмысленно: $g'$ — матрица $m \times p$, $h'$ — матрица $p \times n$, так что произведение — матрица $m \times n$, как и $f'$! Но заметьте, что это не то же самое, что $h'(x)g'(h(x))$: матрицы $p \times n$ и $m \times p$ (при $n \ne m$) перемножить нельзя, а даже если $n = m$, ответ будет неверным, потому что они, скорее всего, не коммутируют.

Важен не только порядок умножения — практически важна и ассоциативность умножения матриц. Рассмотрим функцию

$$f(x) = a(b(c(x))),$$

где $c: \mathbb{R}^n \to \mathbb{R}^p$, $b: \mathbb{R}^p \to \mathbb{R}^q$ и $a: \mathbb{R}^q \to \mathbb{R}^m$. По правилу цепочки

$$f'(x) = a'(b(c(x)))\, b'(c(x))\, c'(x).$$

По ассоциативности это то же самое, что

$$f' = (a'b')c' = a'(b'c')$$

(аргументы функций опущены для краткости). Слева — умножение слева направо, справа — умножение справа налево.

Ну и что? Оказывается, ассоциативность глубоко важна. Настолько, что у двух порядков есть названия: умножение слева направо в области автоматического дифференцирования (AD) называют «обратным режимом» (reverse mode), а умножение справа налево — «прямым режимом» (forward mode). Обратный режим дифференцирования в некоторых контекстах известен также как «метод сопряжённых уравнений» (adjoint method) или «обратное распространение ошибки» (backpropagation); подробнее мы разберём его позже. Почему это важно? Подумаем о вычислительной стоимости умножения матриц.

2.5.1. Стоимость умножения матриц

Рисунок 3 оригинала: умножение матриц ассоциативно — $(AB)C = A(BC)$ для любых $A$, $B$, $C$, — но умножать слева направо может быть гораздо выгоднее, чем справа налево, если в самой левой матрице всего одна (или несколько) строк. Соответственно, порядок применения правила цепочки резко меняет объём вычислений. Порядок «слева направо» называют «обратным режимом» или «обратным распространением»; он лучше всего подходит, когда выходов намного меньше, чем входов. На рисунке — строка, умноженная на две большие квадратные матрицы: если сначала умножить строку на первую матрицу, всё время получаются строки; если сначала перемножить квадратные матрицы, получается дорогая большая матрица.

Умножая матрицу $m \times q$ на матрицу $q \times p$, обычно вычисляют $mp$ скалярных произведений длины $q$ (или делают какую-то равносильную перестановку этих операций). Скалярное произведение длины $q$ требует $q$ умножений и $q - 1$ сложений чисел. Всего это примерно $2mpq$ скалярных операций. В информатике записали бы, что это «$\Theta(mpq)$»: объём вычислений асимптотически пропорционален $mpq$ при больших $m, p, q$.

Так почему порядок в правиле цепочки важен? Рассмотрим два примера.

Пример. Пусть у вас много входов $n \gg 1$ и только один выход $m = 1$, а промежуточных значений тоже много, то есть $q = p = n$. Тогда обратный режим (слева направо) обойдётся в $\Theta(n^2)$ скалярных операций, а прямой (справа налево) — в $\Theta(n^3)$! Это огромная разница в стоимости, схематично показанная на рисунке 3.

Наоборот, пусть у вас много выходов $m \gg 1$ и только один вход $n = 1$, а промежуточных значений $q = p = m$. Тогда обратный режим будет стоить $\Theta(m^3)$ операций, а прямой — всего $\Theta(m^2)$!

Мораль. Если входов много, а выходов мало (обычный случай в машинном обучении и оптимизации), вычисляйте правило цепочки слева направо (обратный режим). Если выходов много, а входов мало, — справа налево (прямой режим). Мы вернёмся к этому в главе 8.

Числовой пример переводчика: пусть $n = 1000$. Слева направо: строка $1 \times 1000$ на матрицу $1000 \times 1000$ — около $2 \cdot 10^6$ операций, и ещё раз столько же: всего $\approx 4 \cdot 10^6$. Справа налево: сначала две матрицы $1000 \times 1000$ — около $2 \cdot 10^9$ операций, то есть в 500 раз больше. У нейросети с миллиардами параметров и одной функцией потерь разница ещё драматичнее — поэтому градиент считают обратным проходом.

flowchart LR
X["x: n входов"] -->|"c′: p×n"| H1["c(x)"] -->|"b′: q×p"| H2["b(c(x))"] -->|"a′: 1×q"| Y["f(x): 1 выход"]

Якобиан f′ = a′·b′·c′. Обратный режим начинает с a′ (одна строка) и идёт к входу: на каждом шаге строка × матрица. Прямой режим начинает с c′ и перемножает большие матрицы.

2.6. За пределами анализа многих переменных

Теперь вычислим производные, которые выходят за рамки первого курса анализа: входы и выходы лежат в более общих векторных пространствах. Например:

Пример. Пусть $A$ — матрица $n \times n$. Возможны такие матричнозначные функции:

• $f(A) = A^3$;
• $f(A) = A^{-1}$, если $A$ обратима;
• или $U$ — матрица, которая получается после применения метода Гаусса к $A$!

Возможны и скалярные выходы, например:

• $f(A) = \det A$;
• $f(A) = \operatorname{tr} A$ (след);
• или $f(A) = \sigma_1(A)$ — наибольшее сингулярное число $A$.

Сосредоточимся в этой лекции на двух более простых примерах.

Пример. Пусть $f(A) = A^3$, где $A$ — квадратная матрица. Найдите $df$.

Применяем правило произведения шаг за шагом:

$$df = dA\,A^2 + A\,dA\,A + A^2\,dA = f'(A)[dA].$$

Заметьте: это не равно $3A^2\,dA$ (если только $dA$ и $A$ не коммутируют, а в общем случае это не так, ведь $dA$ — произвольное малое изменение $A$). Правая часть — линейный оператор $f'(A)$, действующий на $dA$, но её не так просто истолковать как одну «матрицу Якоби», умноженную на $dA$!

Пример. Пусть $f(A) = A^{-1}$, где $A$ — квадратная обратимая матрица. Найдите $df = d(A^{-1})$.

Здесь мы применим небольшую хитрость. Заметим, что $AA^{-1} = I$ — единичная матрица. Значит, дифференциал можно найти по правилу произведения (учитывая, что $dI = 0$: изменение $A$ не меняет $I$):

$$d(AA^{-1}) = dA\,A^{-1} + A\,d(A^{-1}) = d(I) = 0 \implies d(A^{-1}) = -A^{-1}\,dA\,A^{-1}.$$

Как это читать: для числа $a$ это знакомое $d(1/a) = -da/a^2$. Для матриц «$a^2$ в знаменателе» расщепляется на $A^{-1}$ слева и $A^{-1}$ справа, и переставить их нельзя.

Проверьте в NumPy (дополнение переводчика)

import numpy as np
rng = np.random.default_rng(0)

# 1) Пример из раздела 1.2: f(x) = xᵀx, x0 = (3, 4), dx = (0.001, 0.002)
x0 = np.array([3., 4.]); dx = np.array([0.001, 0.002])
print((x0 + dx) @ (x0 + dx) - x0 @ x0)   # 0.022005 — настоящее изменение
print(2 * x0 @ dx)                       # 0.022    — линейное предсказание 2x0ᵀdx

# 2) Градиент xᵀAx равен (A + Aᵀ)x
A = rng.standard_normal((3, 3)); x = rng.standard_normal(3)
dx = 1e-6 * rng.standard_normal(3)
f = lambda x: x @ A @ x
print(f(x + dx) - f(x), ((A + A.T) @ x) @ dx)    # два почти одинаковых числа

# 3) d(A³) = dA·A² + A·dA·A + A²·dA, а «3A²dA» — неверно
dA = 1e-6 * rng.standard_normal((3, 3))
true = np.linalg.matrix_power(A + dA, 3) - np.linalg.matrix_power(A, 3)
good = dA @ A @ A + A @ dA @ A + A @ A @ dA
bad = 3 * A @ A @ dA
print(np.abs(true - good).max())   # ~1e-12: ошибка второго порядка малости
print(np.abs(true - bad).max())    # ~1e-6: ошибка того же порядка, что сам сдвиг

# 4) d(A⁻¹) = -A⁻¹ dA A⁻¹
Ai = np.linalg.inv(A)
print(np.abs(np.linalg.inv(A + dA) - Ai - (-Ai @ dA @ Ai)).max())   # ~1e-11

Словарик

EN	RU
linearization	линеаризация
linear operator, linear map	линейный оператор, линейное отображение
affine function	аффинная функция
differential, infinitesimal	дифференциал, бесконечно малая
little-o notation $o(\delta x)$	«о-малое»
directional derivative	производная по направлению
gradient, steepest ascent	градиент, наискорейший подъём
covector, dual vector, linear form	ковектор, двойственный вектор, линейная форма
Jacobian matrix	матрица Якоби, якобиан
chain rule	правило цепочки (дифференцирование сложной функции)
forward mode / reverse mode	прямой режим / обратный режим
backpropagation, adjoint method	обратное распространение ошибки, метод сопряжённых уравнений
elementwise (Hadamard) product	поэлементное произведение (произведение Адамара)
sensitivity analysis	анализ чувствительности
Banach space	банахово пространство

Дальше: часть 2 — градиенты матричных функций.