Урок 37. Умножение матриц и цена вычисления

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

Три вложенных цикла

Чтобы перемножить $A\in\mathbb R^{m\times k}$ и $B\in\mathbb R^{k\times n}$ , каждый из $mn$ элементов результата собирают как скалярное произведение строки на столбец:

C_{ij}=\sum_{\ell=1}^{k}A_{i\ell}B_{\ell j}.

Это $mkn$ умножений и примерно столько же сложений. Для квадратных матриц $n\times n$ получается порядка $n^3$ операций.

Восьмикратный рост при удвоении стороны — вот что делает большие матрицы дорогими. Но у той же формулы есть скрытая степень свободы, которая меняет цену, не трогая ответ.

Порядок скобок меняет цену

Умножение матриц ассоциативно: $(AB)C=A(BC)$ , ответ один. А вот стоимость двух порядков может отличаться в разы. Возьмём цепочку $A\,(40\times100)$ , $B\,(100\times5)$ , $C\,(5\times60)$ .

(AB)C:\quad 40\cdot100\cdot5+40\cdot5\cdot60=32\,000,

A(BC):\quad 100\cdot5\cdot60+40\cdot100\cdot60=270\,000.

Разница почти в девять раз — на ровном месте, просто из-за расстановки скобок.

УпражнениеВыберите порядок

Даны матрицы $A\,(2\times50)$ , $B\,(50\times2)$ , $C\,(2\times100)$ . Какой порядок умножения дешевле — $(AB)C$ или $A(BC)$ — и во сколько раз?

Показать ответ

$(AB)C$ : сначала $A B$ стоит $2\cdot50\cdot2=200$ , затем результат $2\times2$ на $C$ стоит $2\cdot2\cdot100=400$ ; всего $600$ . $A(BC)$ : сначала $BC$ стоит $50\cdot2\cdot100=10\,000$ , затем $A$ на результат $50\times100$ стоит $2\cdot50\cdot100=10\,000$ ; всего $20\,000$ . Первый порядок дешевле примерно в $33$ раза: выгодно сперва свернуть матрицы в маленькую $2\times2$ .

Слева три матрицы цепочки A 40 на 100, B 100 на 5, C 5 на 60 разного размера. Справа два столбца: порядок скобок A B в скобках умножить на C стоит 32 тысячи умножений, а A умножить на B C в скобках — 270 тысяч; разница почти в восемь раз — Рис. 37.1. Порядок умножений меняет цену в разы
Одна цепочка, один ответ, но $(AB)C$ и $A(BC)$ стоят разного числа умножений. Выгоднее сначала свернуть то, что резко уменьшает размер.

Замерь куб сам

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

Почему замер важнее формулы

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

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

Квадратные матрицы увеличили с $n=256$ до $n=512$ . Во сколько раз выросло число скалярных умножений и во сколько — число элементов результата?

Показать ответ

Умножений становится в $(512/256)^3=2^3=8$ раз больше: кубический рост. Элементов результата — только в $(512/256)^2=4$ раза (их $n^2$ ). Работа на один элемент результата растёт линейно по $n$ (сумма длины $n$ ), поэтому арифметика и объём результата растут в разном темпе — отсюда и смысл арифметической интенсивности.

Лог-лог график измеренного времени numpy.matmul от размера n. Синие точки — реальные измерения — ложатся ниже красной пунктирной линии закона n в кубе, потому что оптимизированная библиотека распараллеливает работу; общий наклон около 2,5 — Рис. 37.2. Кубический закон, измеренный на самом деле
Реальные измерения на этой машине. Оптимизированная библиотека даже обгоняет чистый $n^3$ : её точки уходят ниже эталонной линии, потому что она распараллеливает работу и лучше использует память. Наивная же реализация даёт наклон ровно три.

Блоки и кэш

Спасение от лишних чтений — блочное умножение. Матрицы режут на плитки $b\times b$ , которые целиком помещаются в быструю память. Загруженная плитка используется много раз, прежде чем её вытеснят:

C_{IJ}\mathrel{+}=A_{IK}B_{KJ}.

Арифметика не меняется, а обмен с медленной памятью резко падает.

ИдеяПлитка живёт в кэше

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

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

ПрактикаПорядок циклов

Один и тот же результат, но порядок ijk читает по столбцам (прыжки по памяти), а ikj — по строкам (подряд). Второй быстрее, хотя число операций то же.

Семь произведений Штрассена

Можно ли перебить кубический закон? Штрассен в 1969 году показал, что да. Обычное блочное умножение матриц $2\times2$ требует восьми произведений блоков; Штрассен обошёлся семью — ценой большего числа сложений. Рекурсивно это даёт

T(n)=7\,T(n/2)+O(n^2)\ \Longrightarrow\ T(n)=O\!\left(n^{\log_2 7}\right)\approx O(n^{2{,}807}).

Рисунок шире экрана — проведите по немуОткрыть целиком ↗

Рис. 37.3. Семь произведений вместо восьми — но не бесплатно

Показатель степени у Штрассена меньше, поэтому на больших $n$ его линия уходит ниже. Но добавочные сложения дают перевес лишь после перелома (здесь около $n\approx128$ ); на малых размерах обычный метод выигрывает.

Поэтому на практике берут порог: выше него — рекурсия Штрассена, ниже — оптимизированное обычное умножение. Слишком ранняя рекурсия проигрывает из-за временных матриц и лишних сложений.

УпражнениеПочему Штрассен проигрывает на малом n

Показатель степени у Штрассена ( $2{,}81$ ) меньше, чем у обычного метода ( $3$ ). Почему же на матрице $n=64$ он часто оказывается медленнее?

Показать ответ

Меньший показатель побеждает лишь асимптотически, на больших $n$ . На малых размерах решает постоянный множитель: каждый уровень рекурсии добавляет $18$ сложений и вычитаний блоков и заводит временные матрицы, а сама рекурсия тратит память и мешает кэшу. На $n=64$ этот накладной расход перевешивает экономию одного умножения из восьми. Именно поэтому вводят порог и ниже него считают обычным методом.

Точность быстрого умножения

За скорость Штрассен платит и точностью: вычитания близких чисел усиливают относительную ошибку округления. Поэтому быстрый алгоритм нельзя принимать на веру — результат сверяют с эталоном повышенной точности по относительной норме:

\frac{\lVert C_{\text{тест}}-C_{\text{эталон}}\rVert_F}{\lVert C_{\text{эталон}}\rVert_F}.

Один случайный пример с равномерными числами не проверяет устойчивость: трудные тесты — это матрицы с числами очень разных масштабов и почти взаимной компенсацией.

Свёртка как умножение матриц

Умножение матриц — не абстракция: в него сворачивается сама свёртка из урока 28. Приём im2col превращает каждое локальное окно картинки в строку матрицы, а ядра — в столбцы; одно матричное умножение сразу выдаёт все позиции и каналы. Платят за это временной памятью: пиксели окон повторяются.

МостСвёртка становится GEMM

im2col разворачивает окна свёртки в строки матрицы, ядра — в столбцы. Тогда свёртка — это одно матричное умножение (GEMM), которое ускоритель считает очень быстро, платя памятью за повтор пикселей.

Для слоя $32\times32\times16$ с ядром $3\times3$ матрица im2col имеет $1024$ строки и $144$ столбца на изображение — заметно больше исходных $16\,384$ значений. Именно поэтому архитектуру CNN оценивают по двум числам сразу: параметрам и операциям, ведь один вес свёртки применяется во многих позициях.

УпражнениеОперации свёрточного слоя

Вход $32\times32\times16$ , ядро $3\times3$ , паддинг сохраняет размер, выход $32$ канала. Оцените число умножений на одно изображение через GEMM после im2col.

Показать ответ

Матрица im2col имеет по одной строке на каждую из $32\cdot32=1024$ выходных позиций и по $3\cdot3\cdot16=144$ столбца (окно на все входные каналы). Матрица ядер — $144\times32$ . GEMM делает $1024\cdot144\cdot32\approx4{,}7$ млн умножений. Это и есть цена одного слоя на одну картинку; на пакет из $128$ — умножаем на $128$ .

::::

Русская линия здесь ведёт к тензорным методам. Академик Е. Е. Тыртышников и его школа развивают малоранговые и тензорные разложения (в том числе тензорный поезд), которые позволяют работать с матрицами и данными огромной размерности, недоступными прямому кубическому умножению. Поиск быстрых произведений и тензорные разложения — две стороны одной идеи: заменить лобовой счёт учётом скрытой структуры.

Матрицы в обучении

При обучении сети мало одного $AB$ : обратный проход из урока 25 требует ещё двух умножений на транспонированные матрицы, $\bar A=\bar C B^\top$ и $\bar B=A^\top\bar C$ . Поэтому один полносвязный слой запускает несколько GEMM за шаг, и ускорение только прямого прохода не ускоряет эпоху целиком. Цена матриц — это цена обучения.

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

Что забрать с собой

Умножение матриц стоит порядка $n^3$ операций, и на нём держится время обучения. Но цена прячется в двух местах: в арифметике и в движении данных по памяти, поэтому её мерят, а не только считают. Ассоциативность даёт бесплатную экономию — правильный порядок скобок меняет цену в разы. Штрассен перебивает кубический показатель, но лишь после порога и ценой точности. А свёртка через im2col сводится к тому же матричному умножению, связывая эту главу с CNN и обратным распространением.

Задачи

4баллов

Классический алгоритм вычисляет $C=AB$ скалярными произведениями. Для $A$ размера $10\times20$ и $B$ размера $20\times30$ найдите число умножений и сложений, считая сумму длины $k$ за $k-1$ сложений. Повторите для размеров $20\times40$ на $40\times60$ и для $10\times40$ на $40\times30$ . Для каждого случая приведите отношение к исходному числу операций и объясните, почему удвоение всех трёх размеров и удвоение только внутреннего размера дают разные множители.

5баллов

Матрицы $A,B,C$ размера $4\times4$ хранятся построчно (row-major) в float64; строка кэша равна $64$ байтам и вмещает восемь соседних чисел. Считайте только чтения $B$ : кэш пуст, полностью ассоциативен, вмещает ровно одну строку кэша, LRU, без предвыборки. Для порядков for i: for j: for k: C[i,j]+=A[i,k]*B[k,j] и for i: for k: for j: C[i,j]+=A[i,k]*B[k,j] выпишите при $i=0$ последовательность линейных индексов $B$ , номера строк кэша и hit/miss. Найдите общее число промахов за все четыре $i$ и объясните различие через расположение в памяти и повторное использование.

6баллов

Для каждого $n\in\{64,128,256,512\}$ один раз с seed 17 сгенерируйте float64-матрицы $A,B\sim U[-1,1]^{n\times n}$ . Реализуйте порядки ijk, ikj и тройное блочное умножение без matmul, dot, einsum или @; для блоков $b\in\{8,16,32,64,128\}$ пропускайте $b>n$ . Эталон — numpy.matmul(A,B); критерий корректности $\lVert C-C_{\text{ref}}\rVert_F/\lVert C_{\text{ref}}\rVert_F<10^{-12}$ . Зафиксируйте один поток, три неучитываемых прогрева и $15$ запусков; измеряйте только умножение. Приведите медиану и интервал $10$ – $90\%$ , ускорение относительно ijk и лучший $b$ для каждого $n$ ; объясните, почему лучший $b$ зависит от машины.

6баллов

Реализуйте float64-Штрассен для $n\in\{64,128,256,512,1024\}$ (дополняя нулями до степени двойки) с переходом на numpy.matmul при пороге $c\in\{16,32,64,128,256\}$ . Для каждого $n$ и seed 17 постройте два теста: $A,B\sim N(0,1)$ и «трудный» $A=D_1G_1D_2$ , $B=D_2^{-1}G_2D_3$ с диагоналями показателей $10^{\mathrm{linspace}(-6,6,n)}$ . Для всех $n$ возьмите $C_{\text{ref}}$ повышенной точности и измерьте относительную ошибку Штрассена и обычного numpy.matmul. Время измеряйте после трёх прогревов по медиане семи запусков при одном потоке BLAS. Найдите точку перелома — минимальный $n$ , где лучший порог быстрее numpy.matmul, — и укажите модель процессора, версии библиотек и число потоков.

Три вложенных цикла

Порядок скобок меняет цену

Замерь куб сам

Измерь цену умножения матриц

Почему замер важнее формулы

Блоки и кэш

Точность быстрого умножения

Свёртка как умножение матриц

Матрицы в обучении

Что забрать с собой

Задачи