Matrix Calculus, часть 2: якобианы и градиенты матричных функций

Это перевод глав 3, 5 и 7 конспекта курса MIT «Matrix Calculus (for Machine Learning and Beyond)». Начало — в части 1, карта курса — в путеводителе. Перевод распространяется под лицензией CC BY-NC-SA 4.0.

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

Эта часть отвечает на вопрос, который возникает у всех, кто пишет обратный проход руками: какой формы градиент по матрице весов и как его вывести, не выписывая сотни частных производных. Ответ главы 5: градиент — это то, на что надо скалярно умножить $dA$, чтобы получить $df$; он всегда той же формы, что и $A$. Именно поэтому в PyTorch W.grad имеет ту же форму, что W (урок 3.3.1). Глава 3 показывает, почему явные якобианы для матриц — плохая идея на практике. Глава 7 — классические формулы для $\det$ и $A^{-1}$. Рисунков в этих главах нет; блоки «Проверьте в NumPy» и примечания переводчика — наше дополнение.

3. Якобианы матричных функций

Для функций, у которых входы и/или выходы — матрицы, мы уже видели в прошлых лекциях, что производную всё равно можно определить как линейный оператор — формулой для $f'$, которая переводит малое изменение входа в соответствующее малое изменение выхода. Однако, когда вы впервые изучали линейную алгебру, большинство линейных операций, вероятно, задавались умножением матриц на векторы, и привыкнуть думать о линейных операциях шире может быть непросто. В этой главе мы обсудим, как всё-таки можно представить $f'$ матрицей Якоби даже для матричных входов и выходов и как самый распространённый способ это сделать использует «векторизацию» матриц и новый вид операции над матрицами — произведение Кронекера. Это даёт ещё один, иногда удобный, способ думать о линейных операторах $f'$; но в то же время важно освоиться и с другими способами записи линейных операторов: иногда явный якобиан скрывает ключевую структуру, а часто ещё и вычислительно неэффективен.

В этой части конспекта мы ссылаемся на Pluto-блокнот курса (notes/2x2Jacobians.jl в репозитории) с вычислительными демонстрациями на Julia, где производная квадрата $A^2$ матриц $2 \times 2$ показана с разных сторон.

3.1. Производные матричных функций: линейные операторы

Как мы уже подчёркивали, производная $f'$ — это линейный оператор, переводящий малое изменение входа в малое изменение выхода. Однако эта идея может принять непривычную форму для функций $f(A)$, переводящих матричные входы $A$ в матричные выходы. Например, мы уже рассматривали такие функции квадратных матриц $m \times m$:

• $f(A) = A^3$, для которой $df = f'(A)[dA] = dA\,A^2 + A\,dA\,A + A^2\,dA$;
• $f(A) = A^{-1}$, для которой $df = f'(A)[dA] = -A^{-1}\,dA\,A^{-1}$.

Пример. Ещё проще функция «квадрат матрицы»

$$f(A) = A^2,$$

для которой по правилу произведения

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

Это можно получить и явно из $df = f(A + dA) - f(A) = (A + dA)^2 - A^2$, отбросив член $(dA)^2$.

Во всех этих примерах $f'(A)$ описывается простой формулой для $f'(A)[dA]$, которая с точностью до первого порядка связывает произвольное изменение $dA$ матрицы $A$ с изменением $df = f(A + dA) - f(A)$ функции $f$. Если дифференциал вас смущает, заметьте, что в эту формулу можно подставить любую матрицу $X$, а не только «бесконечно малое» изменение $dA$. Например, для квадрата матрицы

$$f'(A)[X] = XA + AX$$

для произвольной $X$ (это производная по направлению, см. часть 1). Это выражение линейно по $X$: если масштабировать или складывать входы, выходы соответственно масштабируются или складываются:

$$f'(A)[2X] = 2XA + A\,2X = 2(XA + AX) = 2f'(A)[X],$$ $$\begin{aligned} f'(A)[X + Y] &= (X + Y)A + A(X + Y) = XA + YA + AX + AY = XA + AX + YA + AY \\ &= f'(A)[X] + f'(A)[Y]. \end{aligned}$$

Это совершенно нормальный способ задать линейную операцию! Мы не записываем её здесь в привычной форме $f'(A)[X] = (\text{какая-то матрица?}) \times (X \text{ — вектор?})$, и это нормально! Формулу вроде $XA + AX$ легко записать, легко понять и легко с ней вычислять.

Но иногда всё же хочется думать о $f'$ как об одной матрице Якоби — на самом привычном языке линейной алгебры, — и это возможно! Если вы прошли достаточно абстрактный курс линейной алгебры, то знаете, что любой линейный оператор можно представить матрицей, выбрав базисы во входном и выходном пространствах. Здесь же мы будем гораздо конкретнее, потому что у матриц $A$ есть общепринятый «декартов» базис, который называется «векторизацией», и в этом базисе линейные операторы вроде $AX + XA$ особенно легко записать матрицей, если ввести новый вид произведения матриц, широко применяемый в «многомерной» линейной алгебре.

3.2. Простой пример: квадрат матрицы 2 × 2

Для начала подробнее рассмотрим функцию «квадрат матрицы» $f(A) = A^2$ в простом случае матриц $2 \times 2$: они описываются всего четырьмя числами, так что каждый член производной можно рассмотреть явно.

Пример. Для матрицы $2 \times 2$

$$A = \begin{pmatrix} p & r \\ q & s \end{pmatrix}$$

функция «квадрат матрицы» равна

$$f(A) = A^2 = \begin{pmatrix} p & r \\ q & s \end{pmatrix} \begin{pmatrix} p & r \\ q & s \end{pmatrix} = \begin{pmatrix} p^2 + qr & pr + rs \\ pq + qs & qr + s^2 \end{pmatrix}.$$

Когда функция выписана через элементы матрицы $(p, q, r, s)$, естественно думать о ней как об отображении 4 скалярных входов в 4 скалярных выхода. То есть $f$ можно считать равносильной «векторизованной» функции $\tilde{f}: \mathbb{R}^4 \to \mathbb{R}^4$:

$$\tilde{f}\left(\begin{pmatrix} p \\ q \\ r \\ s \end{pmatrix}\right) = \begin{pmatrix} p^2 + qr \\ pq + qs \\ pr + rs \\ qr + s^2 \end{pmatrix}.$$

Такое превращение матрицы в вектор-столбец называется векторизацией и обычно обозначается операцией «$\operatorname{vec}$»:

$$\operatorname{vec} A = \operatorname{vec} \begin{pmatrix} p & r \\ q & s \end{pmatrix} = \begin{pmatrix} p \\ q \\ r \\ s \end{pmatrix} \begin{matrix} \leftarrow A_{1,1} \\ \leftarrow A_{2,1} \\ \leftarrow A_{1,2} \\ \leftarrow A_{2,2} \end{matrix}, \qquad \operatorname{vec} f(A) = \operatorname{vec} \begin{pmatrix} p^2 + qr & pr + rs \\ pq + qs & qr + s^2 \end{pmatrix} = \begin{pmatrix} p^2 + qr \\ pq + qs \\ pr + rs \\ qr + s^2 \end{pmatrix}.$$

Через $\operatorname{vec}$ наша «векторизованная» функция квадрата $\tilde{f}$ определяется так:

$$\tilde{f}(\operatorname{vec} A) = \operatorname{vec} f(A) = \operatorname{vec}(A^2).$$

Определение. Векторизация $\operatorname{vec} A \in \mathbb{R}^{mn}$ любой матрицы $m \times n$ $A \in \mathbb{R}^{m \times n}$ определяется просто как укладка столбцов $A$ слева направо в один вектор-столбец $\operatorname{vec} A$. То есть если обозначить $n$ столбцов $A$ векторами из $m$ компонент $\vec{a}_1, \vec{a}_2, \ldots \in \mathbb{R}^m$, то

$$\operatorname{vec} \underbrace{\begin{pmatrix} \vec{a}_1 & \vec{a}_2 & \cdots & \vec{a}_n \end{pmatrix}}_{A \in \mathbb{R}^{m \times n}} = \begin{pmatrix} \vec{a}_1 \\ \vec{a}_2 \\ \vdots \\ \vec{a}_n \end{pmatrix} \in \mathbb{R}^{mn}$$

— вектор-столбец из $mn$ компонент, содержащий все элементы $A$.

В компьютере элементы матрицы обычно хранятся подряд в памяти, и это можно считать разновидностью векторизации. Более того, $\operatorname{vec} A$ в точности соответствует так называемому хранению «по столбцам» (column-major), когда подряд хранятся элементы столбца; это формат по умолчанию, например, в Fortran, Matlab и Julia, а почтенное наследие Fortran означает, что хранение по столбцам широко принято в библиотеках линейной алгебры.

Примечание переводчика: NumPy и PyTorch по умолчанию хранят матрицы по строкам (row-major). Поэтому A.flatten() в NumPy — это $\operatorname{vec}(A^T)$, а $\operatorname{vec} A$ получается как A.flatten(order='F') (F — от Fortran).

Задача. Вектор $\operatorname{vec} A$ — это коэффициенты, которые получаются, если разложить матрицу $m \times n$ $A$ по некоторому базису матриц. Что это за базис?

Векторизация превращает непривычное (матричные функции и их производные) в привычное (векторные функции, их якобианы и градиенты). Этим она очень привлекательна — пожалуй, даже слишком: зачем «матричный анализ», если всё можно свести обратно к обычному анализу многих переменных? Но у векторизации есть недостатки. Концептуально она может скрывать математическую структуру (например, $\tilde{f}$ выше совсем не похожа на квадрат матрицы $A^2$), а вычислительно такая потеря структуры иногда ведёт к серьёзной неэффективности (например, к построению огромных якобианов $m^2 \times m^2$, о которых ниже). В целом мы считаем, что основной способ изучать такие матричные функции — видеть у них матричные входы ($A$) и матричные выходы ($A^2$) и так же в общем случае считать их производные линейными операторами на матрицах, а не их векторизованными версиями. Тем не менее с векторизацией полезно быть знакомым, чтобы иметь альтернативную точку зрения.

3.2.1. Якобиан 4 × 4 для квадрата матрицы

Чтобы понять якобианы функций (из матриц в матрицы), начнём с простого вопроса.

Вопрос. Каков размер якобиана функции «квадрат матрицы»?

Если смотреть на функцию через её векторизованный аналог $\tilde{f}$, отображающий $\mathbb{R}^4 \to \mathbb{R}^4$ (векторы из 4 компонент в векторы из 4 компонент), якобиан будет матрицей $4 \times 4$ (из производных каждой компоненты выхода по каждой компоненте входа). Теперь возьмём более общую квадратную матрицу $A$ размера $m \times m$. Чтобы найти якобиан $f(A) = A^2$, можно (символьно) проделать то же самое и получить матрицу $m^2 \times m^2$ (поскольку входов $m^2$ — элементы $A$, — и выходов $m^2$ — элементы $A^2$). Явное вычисление этих $m^4$ частных производных утомительно даже при малых $m$, но с этим справляются средства символьных вычислений, например в Julia или Mathematica. Как видно в блокноте, Julia выдаёт якобиан без труда. Для случая $m = 2$, который мы выписали выше, можно либо продифференцировать $\tilde{f}$ вручную, либо воспользоваться символьными средствами Julia и получить якобиан:

$$\tilde{f}' = \begin{array}{c|cccc} & (1,1) & (2,1) & (1,2) & (2,2) \\ \hline (1,1) & 2p & r & q & 0 \\ (2,1) & q & p + s & 0 & q \\ (1,2) & r & 0 & p + s & r \\ (2,2) & 0 & r & q & 2s \end{array}.$$

Например, первая строка $\tilde{f}'$ состоит из частных производных $p^2 + qr$ (первого выхода) по 4 входам $p$, $q$, $r$ и $s$. Здесь строки подписаны индексами (строка, столбец) $(j_\text{out}, k_\text{out})$ элементов «выходной» матрицы $d(A^2)$, а столбцы — индексами $(j_\text{in}, k_\text{in})$ элементов «входной» матрицы $A$. Хотя мы записали якобиан $\tilde{f}'$ как «двумерную» матрицу, его можно представлять и как «четырёхмерную», индексированную $j_\text{out}, k_\text{out}, j_\text{in}, k_\text{in}$.

Однако взгляд матричного анализа, при котором производная $f'(A)$ — это линейное преобразование матриц (как мы вывели выше),

$$f'(A)[X] = XA + AX,$$

кажется гораздо более информативным, чем явный покомпонентный «векторизованный» якобиан $\tilde{f}'$, и даёт формулу для любой матрицы $m \times m$, не заставляя нас утомительно брать $m^4$ частных производных одну за другой. Если же мы действительно хотим держаться векторизации, нужен способ вернуть часть структуры, которую скрывает утомительное покомпонентное дифференцирование. Ключевой инструмент, связывающий две точки зрения, — операция над матрицами, с которой вы, возможно, не знакомы: произведение Кронекера (обозначается $\otimes$).

3.3. Произведения Кронекера

Линейную операцию вроде $f'(A)[X] = XA + AX$ можно представлять себе как «многомерную матрицу»: обычные «двумерные» матрицы переводят «одномерные» векторы-столбцы в одномерные векторы-столбцы, а чтобы переводить двумерные матрицы в двумерные матрицы, можно представить себе «четырёхмерную» матрицу (её иногда называют тензором). Чтобы превращать двумерные матрицы обратно в одномерные векторы, у нас уже есть векторизация ($\operatorname{vec} A$). Тесно связанный инструмент, превращающий «многомерные» линейные операции над матрицами обратно в «двумерные» матрицы для векторизованных входов и выходов, — произведение Кронекера $A \otimes B$. Во вводных курсах линейной алгебры оно встречается нечасто, но появляется в самых разных приложениях с многомерными данными: в многомерной статистике и анализе данных, в многомерных научных и инженерных задачах.

Определение. Если $A$ — матрица $m \times n$ с элементами $a_{ij}$, а $B$ — матрица $p \times q$, то их произведение Кронекера $A \otimes B$ определяется так:

$$A = \begin{pmatrix} a_{11} & \cdots & a_{1n} \\ \vdots & \ddots & \vdots \\ a_{m1} & \cdots & a_{mn} \end{pmatrix} \Longrightarrow \underbrace{A}_{m \times n} \otimes \underbrace{B}_{p \times q} = \underbrace{\begin{pmatrix} a_{11} B & \cdots & a_{1n} B \\ \vdots & \ddots & \vdots \\ a_{m1} B & \cdots & a_{mn} B \end{pmatrix}}_{mp \times nq},$$

так что $A \otimes B$ — матрица $mp \times nq$, составленная «вклеиванием» копии $B$, умноженной на каждый элемент $A$.

Возьмём, например, матрицы $2 \times 2$

$$A = \begin{pmatrix} p & r \\ q & s \end{pmatrix} \quad \text{и} \quad B = \begin{pmatrix} a & c \\ b & d \end{pmatrix}.$$

Тогда $A \otimes B$ — матрица $4 \times 4$, содержащая все возможные произведения элементов $A$ на элементы $B$. Заметьте, что $A \otimes B \ne B \otimes A$ (хотя они отличаются лишь перестановкой элементов):

$$A \otimes B = \begin{pmatrix} pB & rB \\ qB & sB \end{pmatrix} = \begin{pmatrix} pa & pc & ra & rc \\ pb & pd & rb & rd \\ qa & qc & sa & sc \\ qb & qd & sb & sd \end{pmatrix} \quad \ne \quad B \otimes A = \begin{pmatrix} aA & cA \\ bA & dA \end{pmatrix} = \begin{pmatrix} ap & ar & cp & cr \\ aq & as & cq & cs \\ bp & br & dp & dr \\ bq & bs & dq & ds \end{pmatrix}.$$

(В оригинале одна копия $B$ выделена красным: в $A \otimes B$ это левый верхний блок $pa, pc, pb, pd$, а в $B \otimes A$ те же четыре произведения разбросаны по углам блоков.) Больше примеров произведений Кронекера (в том числе с картинками вместо чисел) — в блокноте.

Выше мы видели, что $f(A) = A^2$ при $A = \begin{pmatrix} p & r \\ q & s \end{pmatrix}$ можно считать равносильной функцией $\tilde{f}(\operatorname{vec} A)$, переводящей векторы из 4 входов в 4 выхода ($\mathbb{R}^4 \to \mathbb{R}^4$), с якобианом $4 \times 4$, который мы (или компьютер) с трудом вычислили как 16 поэлементных частных производных. Оказывается, этот результат получается гораздо изящнее, если лучше понимать произведения Кронекера. Мы увидим, что «векторизованный» якобиан $4 \times 4$ — это просто

$$\tilde{f}' = I_2 \otimes A + A^T \otimes I_2,$$

где $I_2$ — единичная матрица $2 \times 2$. То есть матричный линейный оператор $f'(A)[dA] = dA\,A + A\,dA$ после векторизации равносилен

$$\operatorname{vec} \underbrace{f'(A)[dA]}_{dA\,A + A\,dA} = \underbrace{(I_2 \otimes A + A^T \otimes I_2)}_{\tilde{f}'} \operatorname{vec} dA = \underbrace{\begin{pmatrix} 2p & r & q & 0 \\ q & p + s & 0 & q \\ r & 0 & p + s & r \\ 0 & r & q & 2s \end{pmatrix}}_{\tilde{f}'} \underbrace{\begin{pmatrix} dp \\ dq \\ dr \\ ds \end{pmatrix}}_{\operatorname{vec} dA}.$$

Но чтобы понять, почему это так, сначала нужно разобраться в алгебре произведений Кронекера. Для начала хорошее упражнение — убедиться в нескольких более простых свойствах.

Задача. Из определения произведения Кронекера выведите тождества:

1. $(A \otimes B)^T = A^T \otimes B^T$.
2. $(A \otimes B)(C \otimes D) = (AC) \otimes (BD)$.
3. $(A \otimes B)^{-1} = A^{-1} \otimes B^{-1}$. (Следует из свойства 2.)
4. $A \otimes B$ ортогональна (её транспонированная равна обратной), если $A$ и $B$ ортогональны. (Из свойств 1 и 3.)
5. $\det(A \otimes B) = \det(A)^m \det(B)^n$, где $A \in \mathbb{R}^{n \times n}$ и $B \in \mathbb{R}^{m \times m}$.
6. $\operatorname{tr}(A \otimes B) = (\operatorname{tr} A)(\operatorname{tr} B)$.
7. Если $Au = \lambda u$ и $Bv = \mu v$ — собственные векторы и числа $A$ и $B$, то $\lambda\mu$ — собственное число $A \otimes B$ с собственным вектором $u \otimes v$. (Кроме того, поскольку $u \otimes v = \operatorname{vec} X$, где $X = vu^T$, это можно связать через утверждение ниже с тождеством $BXA^T = Bv(Au)^T = \lambda\mu X$.)

3.3.1. Главное тождество для произведения Кронекера

Чтобы превращать линейные операции вроде $AX + XA$ в произведения Кронекера через векторизацию, нужно главное тождество.

Утверждение. Для матриц $A$, $B$, $C$ согласованных размеров

$$(A \otimes B) \operatorname{vec}(C) = \operatorname{vec}(BCA^T).$$

Значит, $A \otimes B$ можно рассматривать как векторизованный аналог линейной операции $C \mapsto BCA^T$. Хочется ввести параллельное обозначение $(A \otimes B)[C] = BCA^T$ для «невекторизованной» версии этой операции, хотя такое обозначение и не стандартно.

Можно запомнить так: $B$ стоит сразу слева от $C$, а $A$ «обходит кругом» направо и транспонируется.

Откуда берётся это тождество? Его можно разбить на более простые части, рассмотрев сначала случаи, когда $A$ или $B$ — единичная матрица $I$ (подходящего размера). Начнём с $A = I$, так что $BCA^T = BC$. Чему равно $\operatorname{vec}(BC)$? Если обозначить через $\vec{c}_1, \vec{c}_2, \ldots$ столбцы $C$, то, напомним, $BC$ просто умножает $B$ слева на каждый столбец $C$:

$$BC = B \begin{pmatrix} \vec{c}_1 & \vec{c}_2 & \cdots \end{pmatrix} = \begin{pmatrix} B\vec{c}_1 & B\vec{c}_2 & \cdots \end{pmatrix} \Longrightarrow \operatorname{vec}(BC) = \begin{pmatrix} B\vec{c}_1 \\ B\vec{c}_2 \\ \vdots \end{pmatrix}.$$

Как получить этот вектор $\operatorname{vec}(BC)$ как что-то, умноженное на $\operatorname{vec} C$? Сразу видно, что

$$\operatorname{vec}(BC) = \begin{pmatrix} B\vec{c}_1 \\ B\vec{c}_2 \\ \vdots \end{pmatrix} = \underbrace{\begin{pmatrix} B & & \\ & B & \\ & & \ddots \end{pmatrix}}_{I \otimes B} \underbrace{\begin{pmatrix} \vec{c}_1 \\ \vec{c}_2 \\ \vdots \end{pmatrix}}_{\operatorname{vec} C},$$

но эта матрица — в точности произведение Кронекера $I \otimes B$! Итак, мы вывели, что

$$(I \otimes B) \operatorname{vec} C = \operatorname{vec}(BC).$$

А как быть с $A^T$? Это чуть хитрее, но снова упростим задачу до случая $B = I$, когда $BCA^T = CA^T$. Чтобы векторизовать это, нужно посмотреть на столбцы $CA^T$. Что такое первый столбец $CA^T$? Это линейная комбинация столбцов $C$, коэффициенты которой задаёт первый столбец $A^T$ (= первая строка $A$):

$$\text{столбец 1 матрицы } CA^T = \sum_j a_{1j} \vec{c}_j.$$

Аналогично для столбца 2 и т. д., а затем мы «укладываем» эти столбцы друг под другом, чтобы получить $\operatorname{vec}(CA^T)$. Но это в точности формула умножения матрицы $A$ на вектор, если бы «элементами» вектора были столбцы $\vec{c}_j$. Явно:

$$\operatorname{vec}(CA^T) = \begin{pmatrix} \sum_j a_{1j} \vec{c}_j \\ \sum_j a_{2j} \vec{c}_j \\ \vdots \end{pmatrix} = \underbrace{\begin{pmatrix} a_{11} I & a_{12} I & \cdots \\ a_{21} I & a_{22} I & \cdots \\ \vdots & \vdots & \ddots \end{pmatrix}}_{A \otimes I} \underbrace{\begin{pmatrix} \vec{c}_1 \\ \vec{c}_2 \\ \vdots \end{pmatrix}}_{\operatorname{vec} C},$$

и, значит, мы вывели

$$(A \otimes I) \operatorname{vec} C = \operatorname{vec}(CA^T).$$

Полное тождество $(A \otimes B) \operatorname{vec}(C) = \operatorname{vec}(BCA^T)$ получается прямым соединением двух выводов: замените во втором выводе $CA^T$ на $BCA^T$ — тогда $\vec{c}_j$ заменится на $B\vec{c}_j$, а значит, $I$ — на $B$.

3.3.2. Якобиан в записи через произведение Кронекера

Теперь воспользуемся этим утверждением, чтобы выразить якобиан $f(A) = A^2$ через произведение Кронекера. Пусть $dA$ играет роль $C$ в утверждении. Сразу видно, что

$$\operatorname{vec}(A\,dA + dA\,A) = \underbrace{(I \otimes A + A^T \otimes I)}_{\text{якобиан } \tilde{f}'(\operatorname{vec} A)} \operatorname{vec}(dA),$$

где $I$ — единичная матрица того же размера, что $A$. Это можно записать и в нашей «невекторизованной» записи линейного оператора:

$$A\,dA + dA\,A = (I \otimes A + A^T \otimes I)[dA].$$

В примере $2 \times 2$ эти произведения Кронекера вычисляются явно:

$$\underbrace{\begin{pmatrix} 1 & \\ & 1 \end{pmatrix}}_{I} \otimes \underbrace{\begin{pmatrix} p & r \\ q & s \end{pmatrix}}_{A} + \underbrace{\begin{pmatrix} p & q \\ r & s \end{pmatrix}}_{A^T} \otimes \underbrace{\begin{pmatrix} 1 & \\ & 1 \end{pmatrix}}_{I} = \underbrace{\begin{pmatrix} p & r & & \\ q & s & & \\ & & p & r \\ & & q & s \end{pmatrix}}_{I \otimes A} + \underbrace{\begin{pmatrix} p & & q & \\ & p & & q \\ r & & s & \\ & r & & s \end{pmatrix}}_{A^T \otimes I} = \begin{pmatrix} 2p & r & q & 0 \\ q & p + s & 0 & q \\ r & 0 & p + s & r \\ 0 & r & q & 2s \end{pmatrix} = \tilde{f}',$$

что в точности совпадает с якобианом $\tilde{f}'$, который мы с трудом вычислили раньше!

Пример. Для функции «куб матрицы» $A^3$, где $A$ — квадратная матрица $m \times m$, найдите якобиан $m^2 \times m^2$ векторизованной функции $\operatorname{vec}(A^3)$.

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

$$(A^3)'[dA] = dA\,A^2 + A\,dA\,A + A^2\,dA,$$

который векторизуется тройным применением тождества Кронекера:

$$\operatorname{vec}(dA\,A^2 + A\,dA\,A + A^2\,dA) = \underbrace{\left((A^2)^T \otimes I + A^T \otimes A + I \otimes A^2\right)}_{\text{векторизованный якобиан}} \operatorname{vec}(dA).$$

(Примечание переводчика: в оригинале в конце формулы опечатка — $\operatorname{vec}(dX)$ вместо $\operatorname{vec}(dA)$.)

Можно продолжить и найти якобианы $A^4$, $A^5$ и так далее или любой линейной комбинации степеней матрицы. Можно даже представить, как применить похожий процесс к ряду Тейлора любой (аналитической) матричной функции $f(A)$, но это становится неуклюжим. Позже (и в домашних заданиях) мы обсудим более изящные способы дифференцировать другие матричные функции — не как векторизованные якобианы, а как линейные операторы на матрицах.

3.3.3. Вычислительная стоимость произведений Кронекера

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

Пусть $A$, $B$ и $C$ — матрицы $m \times m$. Стоимость умножения двух матриц $m \times m$ (обычными методами) растёт пропорционально $\sim m^3$ — информатики называют это «сложностью» $\Theta(m^3)$. Значит, стоимость линейной операции $C \mapsto BCA^T$ растёт как $\sim m^3$ (два умножения $m \times m$). Если же вычислять тот же ответ через $\operatorname{vec}(BCA^T) = (A \otimes B) \operatorname{vec} C$, придётся:

1. Построить матрицу $A \otimes B$ размера $m^2 \times m^2$. Это требует $m^4$ умножений (все элементы $A$ на все элементы $B$) и $\sim m^4$ памяти. (Сравните с $\sim m^2$ памяти, чтобы хранить $A$ или $B$. При $m = 1000$ это в миллион раз больше памяти — терабайты вместо мегабайт!)
2. Умножить $A \otimes B$ на вектор $\operatorname{vec} C$ из $m^2$ элементов. Умножение матрицы $n \times n$ на вектор требует $\sim n^2$ операций, а здесь $n = m^2$, так что это снова $\sim m^4$ арифметических операций.

Итак, вместо $\sim m^3$ операций и $\sim m^2$ памяти для вычисления $BCA^T$ путь через $(A \otimes B) \operatorname{vec} C$ требует $\sim m^4$ операций и $\sim m^4$ памяти — несравнимо хуже! По сути, это потому, что у $A \otimes B$ много структуры, которую мы не используем (это очень особая матрица $m^2 \times m^2$).

Таких примеров много. Ещё один знаменитый — решение линейных матричных уравнений

$$AX + XB = C$$

относительно неизвестной матрицы $X$ при заданных $A$, $B$, $C$, где все они — матрицы $m \times m$. Это называется уравнением Сильвестра. Это линейные уравнения относительно неизвестной $X$, и их можно превратить в обычную систему из $m^2$ линейных уравнений через произведения Кронекера:

$$\operatorname{vec}(AX + XB) = (I \otimes A + B^T \otimes I) \operatorname{vec} X = \operatorname{vec} C,$$

которую затем можно решить относительно $m^2$ неизвестных $\operatorname{vec} X$ методом Гаусса. Но решение системы $m^2 \times m^2$ методом Гаусса стоит $\sim (m^2)^3 = m^6$. Оказывается, однако, что есть хитрые алгоритмы, решающие $AX + XB = C$ всего за $\sim m^3$ операций (с $\sim m^2$ памяти); при $m = 1000$ это экономит множитель $\sim m^3 = 10^9$ (миллиард) в объёме вычислений.

(Произведения Кронекера могут быть более практичным вычислительным инструментом для разреженных матриц — матриц, которые в основном состоят из нулей, например, с несколькими ненулевыми элементами в строке. Дело в том, что произведение Кронекера двух разреженных матриц тоже разрежено, и это избавляет от огромных затрат памяти, которые требует произведение Кронекера неразреженных, «плотных» матриц. Это удобный способ собирать большие разреженные системы уравнений, например для многомерных уравнений в частных производных.)

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

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

vec = lambda M: M.flatten(order='F')    # vec по столбцам, как в конспекте

# Главное тождество: (A ⊗ B) vec(C) = vec(B C Aᵀ)
A = rng.standard_normal((3, 3)); B = rng.standard_normal((2, 2)); C = rng.standard_normal((2, 3))
print(np.allclose(np.kron(A, B) @ vec(C), vec(B @ C @ A.T)))     # True

# Якобиан 4×4 для A² и проверка малым сдвигом
A = rng.standard_normal((2, 2)); I = np.eye(2)
J = np.kron(I, A) + np.kron(A.T, I)
dA = 1e-7 * rng.standard_normal((2, 2))
print(np.abs(vec((A + dA) @ (A + dA) - A @ A) - J @ vec(dA)).max())   # ~1e-14

5. Производные в общих векторных пространствах

Матричный анализ требует ещё шире обобщить понятия производной и градиента — на функции, у которых входы и/или выходы не просто скаляры или векторы-столбцы. Для этого мы расширим понятие обычного скалярного произведения векторов и обычной евклидовой «длины» вектора до общих скалярных произведений (inner products) и норм на векторных пространствах. Первым примером будут знакомые матрицы, рассмотренные с этой точки зрения.

Вспомним из линейной алгебры: любое множество $V$ можно назвать «векторным пространством», если его элементы можно складывать/вычитать $x \pm y$ и умножать на скаляры $\alpha x$ (с соблюдением некоторых базовых арифметических аксиом, например дистрибутивности). Например, векторное пространство образуют сами матрицы $m \times n$ или даже непрерывные функции $u(x)$ (из $\mathbb{R}$ в $\mathbb{R}$) — главное, что их можно складывать, вычитать, масштабировать и получать элементы того же множества. Оказывается, необычайно полезно распространить дифференцирование на такие пространства — например, на функции, переводящие матрицы в матрицы или функции в числа. Ключевое условие — чтобы у входного и выходного векторных пространств $V$ была норма и, желательно, скалярное произведение.

5.1. Простое скалярное произведение и норма для матриц

Напомним, что для скалярных функций $f(x) \in \mathbb{R}$ с векторными входами $x \in \mathbb{R}^n$ (то есть «векторами-столбцами» из $n$ компонент)

$$df = f(x + dx) - f(x) = f'(x)[dx] \in \mathbb{R}.$$

Значит, $f'(x)$ — линейный оператор, который принимает вектор $dx$ и выдаёт скаляр. Иначе говоря, $f'(x)$ — вектор-строка¹ $(\nabla f)^T$. С этой точки зрения $df$ — скалярное произведение:

$$df = \nabla f \cdot dx.$$

¹ Понятие «вектора-строки» формализуется как «ковектор», «двойственный вектор» или элемент «двойственного (сопряжённого) пространства» — не путать с дуальными числами, которые используются в автоматическом дифференцировании (часть 3).

Это можно обобщить на любое векторное пространство $V$ со скалярным произведением! Для $x \in V$ и скалярной функции $f$ получаем линейный оператор $f'(x)[dx] \in \mathbb{R}$, который называется «линейной формой». Чтобы определить градиент $\nabla f$, нужно скалярное произведение на $V$ — обобщение привычного скалярного произведения векторов на векторные пространства!

Для $x, y \in V$ скалярное произведение $\langle x, y \rangle$ — это отображение, для которого $\langle x, y \rangle \in \mathbb{R}$. Его также часто обозначают $x \cdot y$ или $\langle x \mid y \rangle$. Технически скалярное произведение — это отображение, которое:

1. Симметрично: $\langle x, y \rangle = \langle y, x \rangle$ (или эрмитово-симметрично², $\langle x, y \rangle = \overline{\langle y, x \rangle}$, если бы мы работали с комплексными числами);
2. Линейно: $\langle x, \alpha y + \beta z \rangle = \alpha \langle x, y \rangle + \beta \langle x, z \rangle$;
3. Неотрицательно: $\langle x, x \rangle := \Vert x \Vert^2 \ge 0$, и $= 0$ тогда и только тогда, когда $x = 0$.

² Некоторые авторы для комплексных векторных пространств отличают «dot product» от «inner product»: у первого нет комплексного сопряжения, $x \cdot y = y \cdot x$ (тогда $x \cdot x$ не обязано быть вещественным и не равно $\Vert x \Vert^2$), а второе обязано быть эрмитово-симметричным: $\langle x, y \rangle = \bar{x} \cdot y$. Другой источник путаницы с комплексными пространствами — в некоторых областях математики определяют $\langle x, y \rangle = x \cdot \bar{y}$, то есть сопрягают правый аргумент вместо левого (так что произведение линейно по левому аргументу и сопряжённо-линейно по правому). Ну не рады ли вы, что мы остаёмся с вещественными числами?

Заметьте, что из первых двух свойств следует, что произведение линейно и по левому вектору (или сопряжённо-линейно, если бы мы работали с комплексными числами). Ещё одно полезное следствие этих трёх свойств, которое вывести чуть сложнее, — неравенство Коши — Буняковского — Шварца $|\langle x, y \rangle| \le \Vert x \Vert\, \Vert y \Vert$.

Определение (гильбертово пространство). (Полное) векторное пространство со скалярным произведением называется гильбертовым пространством. (Техническое требование «полноты» по сути означает, что в пространстве можно переходить к пределу; оно важно для строгих доказательств.³)

³ Полнота означает, что у любой последовательности Коши в векторном пространстве — любой последовательности точек, которые сближаются всё сильнее, — есть предел, лежащий в этом же пространстве. Для векторных пространств над вещественными или комплексными скалярами это условие на практике обычно выполняется, но с пространствами функций бывает сложнее: например, пределом последовательности непрерывных функций может оказаться разрывная функция.

Имея гильбертово пространство, можно определить градиент скалярной функции. Пусть $x \in V$ — гильбертову пространству — и $f(x)$ — скаляр; тогда у нас есть линейная форма $f'(x)[dx] \in \mathbb{R}$. При этих предположениях теорема, известная как теорема Рисса о представлении, утверждает, что любая линейная форма (в том числе $f'$) должна быть скалярным произведением с чем-то:

$$f'(x)[dx] = \big\langle \underbrace{\text{(некоторый вектор)}}_{\text{градиент } \nabla f|_x}, dx \big\rangle = df.$$

То есть градиент $\nabla f$ определяется как то, на что надо скалярно умножить $dx$, чтобы получить $df$. Заметьте, что $\nabla f$ всегда имеет «ту же форму», что и $x$.

Первые несколько примеров — обычное гильбертово пространство $V = \mathbb{R}^n$ с разными скалярными произведениями.

Пример. Для $V = \mathbb{R}^n$ из векторов-столбцов длины $n$ есть привычное евклидово скалярное произведение $\langle x, y \rangle = x^T y$. Оно приводит к обычному $\nabla f$.

Пример. На $\mathbb{R}^n$ бывают и другие скалярные произведения. Например,

$$\langle x, y \rangle_W = w_1 x_1 y_1 + w_2 x_2 y_2 + \dots + w_n x_n y_n = x^T \underbrace{\begin{pmatrix} w_1 & & \\ & \ddots & \\ & & w_n \end{pmatrix}}_{W} y$$

с весами $w_1, \dots, w_n > 0$.

Шире, взвешенное скалярное произведение $\langle x, y \rangle_W = x^T W y$ можно определить для любой симметричной положительно определённой матрицы $W$ ($W = W^T$ и $W$ положительно определена — этого достаточно, чтобы произведение было корректным).

Если поменять определение скалярного произведения, поменяется и определение градиента! Например, для $f(x) = x^T A x$ мы ранее нашли $df = x^T (A + A^T)\,dx$. С обычным евклидовым скалярным произведением это дало градиент $\nabla f = (A + A^T)x$. Но если взять взвешенное произведение $x^T W y$, получится другой «градиент» $\nabla^{(W)} f = W^{-1}(A + A^T)x$, такой что $df = \langle \nabla^{(W)} f, dx \rangle_W$.

В этом конспекте мы, если не оговорено иное, используем для $x \in \mathbb{R}^n$ евклидово скалярное произведение и, значит, обычный $\nabla f$. Однако взвешенные скалярные произведения полезны во многих случаях, особенно когда компоненты $x$ имеют разные масштабы или единицы измерения.

Комментарий переводчика: именно эта идея стоит за «предобусловленными» оптимизаторами. Adam, по сути, делит градиент покомпонентно на оценку масштаба — это шаг по градиенту в «взвешенной» геометрии (урок 1.4.2).

Можно рассмотреть и пространство матриц $m \times n$, $V = \mathbb{R}^{m \times n}$. Существует, конечно, изоморфизм векторных пространств $V \ni A \to \operatorname{vec}(A) \in \mathbb{R}^{mn}$. Поэтому в этом пространстве есть аналог привычного евклидова скалярного произведения («фробениусово»), которое удобно переписать через матричные операции с помощью следа.

Определение (скалярное произведение Фробениуса). Скалярное произведение Фробениуса двух матриц $m \times n$ $A$ и $B$:

$$\langle A, B \rangle_F = \sum_{ij} A_{ij} B_{ij} = \operatorname{vec}(A)^T \operatorname{vec}(B) = \operatorname{tr}(A^T B).$$

С ним связана соответствующая норма Фробениуса:

$$\Vert A \Vert_F = \sqrt{\langle A, A \rangle_F} = \sqrt{\operatorname{tr}(A^T A)} = \Vert \operatorname{vec} A \Vert = \sqrt{\sum_{i,j} |A_{ij}|^2}.$$

С её помощью мы теперь можем определить градиент скалярных функций с матричными входами! Это будет наше скалярное произведение матриц по умолчанию (а значит, и наш градиент по матрице по умолчанию) в этом конспекте (иногда индекс $F$ мы опускаем).

Пример. Рассмотрим функцию

$$f(A) = \Vert A \Vert_F = \sqrt{\operatorname{tr}(A^T A)}.$$

Чему равен $df$?

Во-первых, по привычным правилам цепочки и степени для скаляров

$$df = \frac{1}{2\sqrt{\operatorname{tr}(A^T A)}}\, d(\operatorname{tr} A^T A).$$

Далее заметим, что (по линейности следа)

$$d(\operatorname{tr} B) = \operatorname{tr}(B + dB) - \operatorname{tr}(B) = \operatorname{tr}(B) + \operatorname{tr}(dB) - \operatorname{tr}(B) = \operatorname{tr}(dB).$$

Следовательно,

$$\begin{aligned} df &= \frac{1}{2\Vert A \Vert_F} \operatorname{tr}(d(A^T A)) \\ &= \frac{1}{2\Vert A \Vert_F} \operatorname{tr}(dA^T\, A + A^T\, dA) \\ &= \frac{1}{2\Vert A \Vert_F} \left(\operatorname{tr}(dA^T\, A) + \operatorname{tr}(A^T\, dA)\right) \\ &= \frac{1}{\Vert A \Vert_F} \operatorname{tr}(A^T\, dA) = \left\langle \frac{A}{\Vert A \Vert_F}, dA \right\rangle. \end{aligned}$$

Здесь мы воспользовались тем, что $\operatorname{tr} B = \operatorname{tr} B^T$, а на последнем шаге связали $df$ со скалярным произведением Фробениуса. Иными словами,

$$\nabla f = \nabla \Vert A \Vert_F = \frac{A}{\Vert A \Vert_F}.$$

Заметьте, что для векторов-столбцов $x$ получается в точности тот же результат: $\nabla \Vert x \Vert = x / \Vert x \Vert$ (и это на самом деле одно и то же через $x = \operatorname{vec} A$).

Рассмотрим ещё один простой пример.

Пример. Зафиксируем постоянные $x \in \mathbb{R}^m$, $y \in \mathbb{R}^n$ и рассмотрим функцию $f: \mathbb{R}^{m \times n} \to \mathbb{R}$, заданную формулой

$$f(A) = x^T A y.$$

Чему равен $\nabla f$?

Имеем

$$\begin{aligned} df &= x^T\, dA\, y \\ &= \operatorname{tr}(x^T\, dA\, y) \\ &= \operatorname{tr}(y x^T\, dA) \\ &= \big\langle \underbrace{x y^T}_{\nabla f}, dA \big\rangle. \end{aligned}$$

Комментарий переводчика: это ровно градиент линейного слоя. Если выход слоя $z = Wh$, а ошибка (градиент потерь по $z$) — вектор $\delta$, то потери меняются на $\delta^T dW\, h$, и по примеру выше $\nabla_W L = \delta\, h^T$ — внешнее произведение «ошибка × вход», той же формы, что $W$. Именно это считает обратный проход в уроке 3.1.2.

Шире, для любой скалярной функции $f(A)$ из определения скалярного произведения Фробениуса следует:

$$df = f(A + dA) - f(A) = \langle \nabla f, dA \rangle = \sum_{i,j} (\nabla f)_{i,j}\, dA_{i,j},$$

а значит, компоненты градиента — в точности поэлементные производные

$$(\nabla f)_{i,j} = \frac{\partial f}{\partial A_{i,j}},$$

как и в покомпонентном определении градиента из анализа многих переменных! Но для нетривиальных функций от матриц $f(A)$ брать производную по каждому элементу $A$ по отдельности бывает крайне неудобно. С «целостным» определением через скалярное произведение матриц мы вскоре сможем вычислять ещё более сложные градиенты по матрицам, включая $\nabla(\det A)$!

5.2. Производные, нормы и банаховы пространства

Мы весь курс говорим о «норме», но что такое норма технически? Конечно, есть знакомые примеры вроде евклидовой нормы («$\ell^2$») $\Vert x \Vert = \sqrt{\sum_k x_k^2}$ для $x \in \mathbb{R}^n$, но полезно понять, как это понятие обобщается на другие векторные пространства. Оказывается, нормы на самом деле необходимы для определения производной!

Для векторного пространства $V$ норма $\Vert \cdot \Vert$ на $V$ — это отображение $\Vert \cdot \Vert: V \to \mathbb{R}$ с тремя свойствами:

1. Неотрицательность: $\Vert v \Vert \ge 0$ и $\Vert v \Vert = 0 \iff v = 0$;
2. Однородность: $\Vert \alpha v \Vert = |\alpha|\, \Vert v \Vert$ для любого $\alpha \in \mathbb{R}$;
3. Неравенство треугольника: $\Vert u + v \Vert \le \Vert u \Vert + \Vert v \Vert$.

Векторное пространство с нормой называется нормированным векторным пространством. Часто математикам технически нужен чуть более точный вид нормированного пространства с менее очевидным названием — банахово пространство.

Определение (банахово пространство). (Полное) векторное пространство с нормой называется банаховым пространством. (Как и для гильбертовых пространств, «полнота» — техническое требование для некоторых видов строгого анализа; по сути, оно позволяет переходить к пределу.)

Например, по любому скалярному произведению $\langle u, v \rangle$ строится соответствующая норма $\Vert u \Vert = \sqrt{\langle u, u \rangle}$. (Таким образом, любое гильбертово пространство — ещё и банахово.⁴)

⁴ Доказать неравенство треугольника для произвольного скалярного произведения не так просто; для этого используют неравенство Коши — Буняковского — Шварца.

Чтобы определить производные, технически нужно, чтобы банаховыми пространствами были и вход, и выход. Чтобы увидеть это, вспомним наш формализм:

$$f(x + \delta x) - f(x) = \underbrace{f'(x)[\delta x]}_{\text{линейное}} + \underbrace{o(\delta x)}_{\text{меньше}}.$$

Чтобы точно определить, в каком смысле члены $o(\delta x)$ «меньше» или «более высокого порядка», нужны нормы. В частности, обозначение «о-малое» $o(\delta x)$ означает любую функцию, для которой

$$\lim_{\delta x \to 0} \frac{\Vert o(\delta x) \Vert}{\Vert \delta x \Vert} = 0,$$

то есть которая стремится к нулю быстрее, чем линейно по $\delta x$. Для этого нормы нужны и у входа $\delta x$, и у выхода (функции). Такое распространение дифференцирования на произвольные нормированные (банаховы) пространства иногда называют производной Фреше.

7. Производная определителя и обратной матрицы

7.1. Два вывода

Этот раздел конспекта следует Julia-блокноту курса (notes/determinant_and_inverse.jl). Блокнот коротковат, но вычисление важное и полезное.

Теорема. Для квадратной матрицы $A$

$$\nabla(\det A) = \operatorname{cofactor}(A) = (\det A)\, A^{-T} := \operatorname{adj}(A^T) = \operatorname{adj}(A)^T,$$

где $\operatorname{adj}$ — «присоединённая» матрица (adjugate). (Возможно, вы не слышали о присоединённой матрице, но эта формула говорит, что это просто $\operatorname{adj}(A) = \det(A) A^{-1}$, или $\operatorname{cofactor}(A) = \operatorname{adj}(A^T)$.) Кроме того,

$$d(\det A) = \operatorname{tr}(\det(A) A^{-1}\, dA) = \operatorname{tr}(\operatorname{adj}(A)\, dA) = \operatorname{tr}(\operatorname{cofactor}(A)^T dA).$$

Возможно, вы помните, что каждый элемент $(i, j)$ матрицы алгебраических дополнений (cofactor) — это $(-1)^{i+j}$, умноженное на определитель, полученный вычёркиванием из $A$ строки $i$ и столбца $j$. Вот несколько вычислений для $2 \times 2$, чтобы почувствовать эти функции:

$$\begin{aligned} M &= \begin{pmatrix} a & c \\ b & d \end{pmatrix} \\ \implies \operatorname{cofactor}(M) &= \begin{pmatrix} d & -b \\ -c & a \end{pmatrix}, \\ \operatorname{adj}(M) &= \begin{pmatrix} d & -c \\ -b & a \end{pmatrix}, \\ M^{-1} &= \frac{1}{ad - bc} \begin{pmatrix} d & -c \\ -b & a \end{pmatrix}. \end{aligned}$$

Примечание переводчика: в оригинале у матриц cofactor и adj переставлены внедиагональные элементы (а $M^{-1}$ записана через cofactor вместо adj). Для $M$ с первой строкой $(a, c)$ и второй $(b, d)$ дополнение к элементу $c$ (позиция $(1,2)$) равно $-b$, поэтому верно так, как в переводе; проверка: $M \cdot \operatorname{adj}(M) = (ad - bc) I$.

Численно, как сделано в блокноте, можно взять случайную матрицу $A$ размера $n \times n$ (скажем, $9 \times 9$), положить, например, $dA = 0{,}00001\,A$ и увидеть, что

$$\det(A + dA) - \det(A) \approx \operatorname{tr}(\operatorname{adj}(A)\, dA),$$

что численно подтверждает теорему.

Теперь докажем теорему двумя способами. Во-первых, есть прямое доказательство: продифференцировать скаляр по каждому входу, используя разложение определителя по $i$-й строке (разложение Лапласа). Напомним, что

$$\det(A) = A_{i1} C_{i1} + A_{i2} C_{i2} + \dots + A_{in} C_{in}.$$

Значит,

$$\frac{\partial \det A}{\partial A_{ij}} = C_{ij} \implies \nabla(\det A) = C$$

— матрица алгебраических дополнений. (Вычисляя эти частные производные, важно помнить, что дополнение $C_{ij}$ не содержит элементов $A$ из строки $i$ и столбца $j$. Так что, например, $A_{i1}$ явно входит только в первый член и не спрятан ни в одном из $C$ в этом разложении.)

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

$$\det(I + dA) - 1 = \operatorname{tr}(dA),$$

и, значит,

$$\begin{aligned} \det(A + A(A^{-1}\, dA)) - \det(A) &= \det(A)\left(\det(I + A^{-1}\, dA) - 1\right) \\ &= \det(A) \operatorname{tr}(A^{-1}\, dA) = \operatorname{tr}(\det(A) A^{-1}\, dA) \\ &= \operatorname{tr}(\operatorname{adj}(A)\, dA). \end{aligned}$$

Это тоже доказывает теорему.

Почему $\det(I + dA) - 1 = \operatorname{tr}(dA)$ (пояснение переводчика): в разложении определителя $I + dA$ единственное слагаемое без произведений двух малых элементов — произведение диагонали $(1 + dA_{11})(1 + dA_{22}) \cdots$, а оно с точностью до первого порядка равно $1 + dA_{11} + dA_{22} + \dots$

7.2. Применения

7.2.1. Характеристический многочлен

Применим это, чтобы найти производную характеристического многочлена в точке $x$. Пусть $p(x) = \det(xI - A)$ — скалярная функция $x$. Напомним, что через разложение на множители $p(x)$ можно записать через собственные числа $\lambda_i$. Значит, можно спросить: чему равна производная $p(x)$? Средствами анализа первого курса можно просто вычислить

$$\frac{d}{dx} \prod_i (x - \lambda_i) = \sum_i \prod_{j \ne i} (x - \lambda_j) = \prod_i (x - \lambda_i) \left\{ \sum_i (x - \lambda_i)^{-1} \right\},$$

если $x \ne \lambda_i$.

Это вполне хорошее простое доказательство, но с нашей новой техникой есть другое:

$$\begin{aligned} d(\det(xI - A)) &= \det(xI - A) \operatorname{tr}\left((xI - A)^{-1}\, d(xI - A)\right) \\ &= \det(xI - A) \operatorname{tr}\left((xI - A)^{-1}\right) dx. \end{aligned}$$

Здесь мы воспользовались тем, что $d(xI - A) = dx\, I$, когда $A$ постоянна, и $\operatorname{tr}(A\, dx) = \operatorname{tr}(A)\, dx$, поскольку $dx$ — скаляр.

Это снова можно проверить численно, как в блокноте.

7.2.2. Логарифмическая производная

Аналогично по правилу цепочки

$$d(\log(\det A)) = \frac{d(\det A)}{\det A} = \det(A^{-1})\, d(\det A) = \operatorname{tr}(A^{-1}\, dA).$$

Логарифмическая производная часто встречается в прикладной математике. Заметьте: здесь мы воспользовались тем, что $\frac{1}{\det A} = \det(A^{-1})$, поскольку $1 = \det(I) = \det(AA^{-1}) = \det(A)\det(A^{-1})$.

Например, вспомним метод Ньютона для поиска корней $f(x) = 0$ вещественных функций одной переменной $f(x)$: делаем последовательность шагов $x \to x + \delta x$. Ключевая формула метода Ньютона — $\delta x = -f'(x)^{-1} f(x)$, а это то же самое, что $-\frac{1}{(\log f(x))'}$. Поэтому производные логарифма определителя появляются при поиске корней определителей, то есть для $f(x) = \det M(x)$. Когда $M(x) = A - xI$, корни определителя — собственные числа $A$. Для более общих функций $M(x)$ решение $\det M(x) = 0$ поэтому называют нелинейной задачей на собственные значения.

Примечание переводчика: в оригинале в формуле шага Ньютона потерян знак минус.

Комментарий переводчика: $\nabla \log\det A = A^{-T}$ встречается и в машинном обучении — например, в логарифме правдоподобия многомерного нормального распределения (там стоит $\log\det$ матрицы ковариаций) и в нормализующих потоках.

7.3. Якобиан обратной матрицы

Наконец, вычислим производную обратной матрицы (и как линейный оператор, и как явную матрицу Якоби). Есть изящный приём, позволяющий получить эту производную просто из свойства обратной $A^{-1}A = I$. По правилу произведения отсюда следует

$$\begin{aligned} d(A^{-1}A) &= d(I) = 0 = d(A^{-1})A + A^{-1}\, dA \\ &\implies \boxed{d(A^{-1}) = (A^{-1})'[dA] = -A^{-1}\, dA\, A^{-1}}. \end{aligned}$$

Вторая строка задаёт вполне хороший линейный оператор для производной $(A^{-1})'$, но при желании его можно переписать явной матрицей Якоби, используя произведения Кронекера, действующие на «векторизованные» матрицы, как в главе 3:

$$\operatorname{vec}\left(d(A^{-1})\right) = \operatorname{vec}\left(-A^{-1}(dA)A^{-1}\right) = \underbrace{-(A^{-T} \otimes A^{-1})}_{\text{якобиан}} \operatorname{vec}(dA),$$

где $A^{-T}$ обозначает $(A^{-1})^T = (A^T)^{-1}$. Эту формулу можно проверить численно, как в блокноте.

На практике, однако, вы, вероятно, обнаружите, что операторное выражение $-A^{-1}\, dA\, A^{-1}$ полезнее явного якобиана для производных с обратными матрицами. Например, если у вас есть матричная функция $A(t)$ скалярного параметра $t \in \mathbb{R}$, вы сразу получаете $\frac{d(A^{-1})}{dt} = -A^{-1} \frac{dA}{dt} A^{-1}$. Более изощрённое применение — метод сопряжённых уравнений (adjoint method) в главе 6 конспекта.

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

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

A = rng.standard_normal((9, 9))
dA = 1e-6 * rng.standard_normal((9, 9))
inner = lambda X, Y: np.sum(X * Y)          # скалярное произведение Фробениуса = tr(XᵀY)

# ∇ det A = det(A)·A⁻ᵀ
g_det = np.linalg.det(A) * np.linalg.inv(A).T
print(np.linalg.det(A + dA) - np.linalg.det(A), inner(g_det, dA))   # почти равны

# ∇ ‖A‖_F = A / ‖A‖_F
nrm = np.linalg.norm                        # для матрицы по умолчанию — норма Фробениуса
print(nrm(A + dA) - nrm(A), inner(A / nrm(A), dA))

# ∇ (xᵀ A y) = x yᵀ
x, y = rng.standard_normal(9), rng.standard_normal(9)
print(x @ (A + dA) @ y - x @ A @ y, inner(np.outer(x, y), dA))

# d log|det A| = tr(A⁻¹ dA)   (slogdet устойчивее, чем log(det))
print(np.linalg.slogdet(A + dA)[1] - np.linalg.slogdet(A)[1], np.trace(np.linalg.inv(A) @ dA))

Словарик

EN	RU
vectorization, vec	векторизация, vec
column-major / row-major	хранение по столбцам / по строкам
Kronecker product	произведение Кронекера
tensor	тензор
Sylvester equation	уравнение Сильвестра
sparse / dense matrix	разреженная / плотная матрица
inner product	скалярное произведение
Frobenius inner product, Frobenius norm	скалярное произведение Фробениуса, норма Фробениуса
trace	след
Hilbert space, Banach space	гильбертово пространство, банахово пространство
Riesz representation theorem	теорема Рисса о представлении
Cauchy–Schwarz inequality	неравенство Коши — Буняковского — Шварца
triangle inequality, homogeneity	неравенство треугольника, однородность
Fréchet derivative	производная Фреше
cofactor matrix, adjugate	матрица алгебраических дополнений, присоединённая матрица
Laplace (cofactor) expansion	разложение Лапласа (по строке)
characteristic polynomial	характеристический многочлен
logarithmic derivative	логарифмическая производная
nonlinear eigenproblem	нелинейная задача на собственные значения

Дальше: часть 3 — прямой и обратный режим автоматического дифференцирования.