Штрассен: семь умножений вместо восьми

Коротко

Произведение матриц C = A·B определяется так: C[i][j] - это сумма Σ_k A[i][k]·B[k][j], то есть i-я строка первой матрицы скалярно умножается на j-й столбец второй. Элементов n², на каждый n умножений - по определению выходит n³ умножений чисел. Разрежем каждую матрицу на четыре блока n/2 × n/2 - произведение выражается через ВОСЕМЬ произведений блоков, то есть T(n) = 8·T(n/2) + n², и это снова n³. Штрассен в 1969 году нашёл способ обойтись семью: семь хитрых произведений сумм блоков, из которых четыре блока ответа собираются сложениями. Ответ падает до Θ(n^log₂ 7) = Θ(n^2.807). Это был первый алгоритм быстрее n³ за всю историю, и именно он открыл гонку за показателем ω, который сегодня опущен примерно до 2.37.

Если ты знаешь, почему блочное умножение даёт 8·T(n/2), помнишь, что Штрассен экономит одно произведение из восьми, и понимаешь его практические ограничения, листай до «Теперь сам».

Как думаешь, сколько умножений нужно для матриц 2×2

Перемножить две матрицы 2 × 2 по определению - это четыре элемента ответа, каждый сумма двух произведений: восемь умножений и четыре сложения.

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

На что это похоже

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

Механизм: блоки и семь произведений

A = | A₁₁  A₁₂ |      B = | B₁₁  B₁₂ |      каждый блок n/2 × n/2
    | A₂₁  A₂₂ |          | B₂₁  B₂₂ |

наивно:
C₁₁ = A₁₁B₁₁ + A₁₂B₂₁      C₁₂ = A₁₁B₁₂ + A₁₂B₂₂
C₂₁ = A₂₁B₁₁ + A₂₂B₂₁      C₂₂ = A₂₁B₁₂ + A₂₂B₂₂
восемь произведений блоков  →  T(n) = 8·T(n/2) + n²  →  Θ(n^log₂ 8) = Θ(n³)

ШТРАССЕН: семь произведений

M₁ = (A₁₁ + A₂₂)(B₁₁ + B₂₂)        M₅ = (A₁₁ + A₁₂)·B₂₂
M₂ = (A₂₁ + A₂₂)·B₁₁               M₆ = (A₂₁ − A₁₁)(B₁₁ + B₁₂)
M₃ = A₁₁·(B₁₂ − B₂₂)               M₇ = (A₁₂ − A₂₂)(B₂₁ + B₂₂)
M₄ = A₂₂·(B₂₁ − B₁₁)

C₁₁ = M₁ + M₄ − M₅ + M₇            C₁₂ = M₃ + M₅
C₂₁ = M₂ + M₄                      C₂₂ = M₁ − M₂ + M₃ + M₆

T(n) = 7·T(n/2) + n²  →  Θ(n^log₂ 7) = Θ(n^2.807)

Сложения блоков стоят n² и в асимптотике ничего не меняют - решает только число произведений.

Правило

Блочное умножение матриц даёт T(n) = 8·T(n/2) + n², то есть привычный Θ(n³): деление само по себе не помогает. Штрассен собирает четыре блока ответа из СЕМИ произведений сумм блоков, поэтому T(n) = 7·T(n/2) + n² и Θ(n^log₂ 7) = Θ(n^2.807). Плата - восемнадцать сложений блоков вместо четырёх, дополнительная память под промежуточные матрицы и худшая численная устойчивость на вещественных числах. Поэтому Штрассен применяют с порогом (сотни элементов) и редко - но он открыл целое направление: показатель ω сегодня опущен примерно до 2.37.

Разбор: где именно теряется одно умножение

Посчитаем на n = 1024, считая только умножения ЧИСЕЛ.

наивно:      1024³ = 1 073 741 824
Штрассен:    1024^2.807 ≈ 2.8 · 10⁸          (точно: 7^10 = 282 475 249)
отношение:   примерно 3.8 раза

Точное число у Штрассена легко считается по уровням: на каждом уровне число подзадач умножается на семь, а не на восемь, и на n = 2^k элементов выходит ровно 7^k умножений против 8^k. Отношение - (8/7)^k, то есть на k = 10 уже 3.8, а на k = 20 (матрицы миллион на миллион, сугубо теоретически) - 14.4.

Теперь честная оговорка про практику. Матрицы 1024 × 1024 библиотеки перемножают за доли секунды БЕЗ Штрассена - за счёт того, что наивный алгоритм идеально ложится на кеш и векторные инструкции процессора. Штрассен же ломает регулярность доступа к памяти и требует места под промежуточные блоки. Поэтому в реальных библиотеках (BLAS, cuBLAS) он либо не используется вовсе, либо включается на очень больших размерах. Это важный урок про асимптотику: она сравнивает ЧИСЛО операций, а время считается ещё и памятью.

Разбор: почему устойчивость хуже

У наивного умножения каждый элемент ответа - сумма произведений ИСХОДНЫХ чисел. У Штрассена элемент собирается из семи величин, которые сами являются произведениями РАЗНОСТЕЙ и СУММ.

наивно:      C₁₁ = A₁₁B₁₁ + A₁₂B₂₁                     слагаемые того же порядка, что данные
Штрассен:    C₁₁ = M₁ + M₄ − M₅ + M₇                   слагаемые могут быть много больше ответа

Если M₁ и M₅ велики, а их разность мала, то при вычислениях с плавающей точкой значащие цифры сокращаются - это классическая потеря точности при вычитании близких чисел. На целых числах проблемы нет вовсе (вычитание точное), поэтому в компьютерной алгебре и криптографии Штрассен применяют спокойнее, чем в численных расчётах.

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

Что ломается без этого

Первое: ожидать от Штрассена ускорения на матрицах 100 × 100 - на таких размерах наивный алгоритм с хорошей раскладкой по кешу быстрее в разы. Второе: реализовать его без порога и без дополнения матрицы до СТЕПЕНИ двойки: у 6 × 6 блоки выйдут 3 × 3 и дальше не делятся, и нужен либо паддинг, либо отдельная обработка края. Третье: забыть про память: семь промежуточных матриц на каждом уровне рекурсии - это заметный расход, если не переиспользовать буферы. Четвёртое: применять его к вещественным данным, где важна точность, не проверив устойчивость. Пятое: думать, что ω = 2.807 - это и есть сегодняшний рекорд: теоретические алгоритмы опустили показатель до примерно 2.37, но их константы настолько велики, что на практике они бесполезны (их так и называют - «галактические алгоритмы»).

Зачем это в работе

Прямо - почти никогда: матрицы перемножают библиотеками, и там стоит наивный алгоритм, вылизанный под кеш и SIMD. Знать про Штрассена стоит по трём причинам. Первая: он показывает, что «по определению» и «оптимально» - разные вещи, и это работает далеко за пределами матриц. Вторая: он объясняет, почему у задачи бывает показатель ω, про который говорят «неизвестен» - нижняя оценка для умножения матриц до сих пор всего лишь n², и разрыв между 2 и 2.37 открыт полвека. Третья, самая практичная: он даёт готовый шаблон рассуждения «а нельзя ли получить нужные суммы меньшим числом произведений» - тот же приём работает для умножения комплексных чисел: (a+bi)(c+di) наивно требует четырёх вещественных умножений (ac, bd, ad, bc), а хватает трёх - посчитать ac, bd и (a+b)(c+d), тогда мнимая часть равна (a+b)(c+d) − ac − bd; это ровно тождество Карацубы и для сложных схем в цифровой обработке сигналов.

Теперь сам

Матрицы 2 × 2: A = [[1, 2], [3, 4]], B = [[5, 6], [7, 8]]. Посчитай C₁₂ двумя способами - по определению и через M₃ + M₅.

Ответ. По определению: C₁₂ = A₁₁B₁₂ + A₁₂B₂₂ = 1·6 + 2·8 = 6 + 16 = 22.

Через Штрассена (здесь блоки - это просто числа):

M₃ = A₁₁·(B₁₂ − B₂₂) = 1 · (6 − 8) = −2
M₅ = (A₁₁ + A₁₂)·B₂₂ = (1 + 2) · 8  = 24
C₁₂ = M₃ + M₅ = −2 + 24 = 22

Сходится. Обрати внимание на промежуточные величины: M₅ = 24 больше самого ответа 22, а M₃ вообще отрицательное, хотя все исходные числа положительны. Это и есть та самая особенность, из-за которой на вещественных данных может теряться точность: ответ собирается из величин другого масштаба и другого знака.

А ещё заметь, что для матриц 2 × 2 выигрыша нет вовсе: семь умножений и восемнадцать сложений против восьми умножений и четырёх сложений. Выигрыш появляется только когда «числа» в формулах - это БЛОКИ, и одно сэкономленное произведение блоков стоит (n/2)^2.807 операций.

Главное

Блочное умножение матриц даёт T(n) = 8·T(n/2) + n² и всё тот же Θ(n³) - деление не помогает. Штрассен собирает четыре блока ответа из СЕМИ произведений сумм блоков: T(n) = 7·T(n/2) + n², то есть Θ(n^2.807); на n = 2^k это ровно 7^k умножений против 8^k, отношение (8/7)^k. Плата - восемнадцать сложений вместо четырёх, память под промежуточные блоки и худшая численная устойчивость: элементы ответа собираются из величин другого масштаба и знака. На практике библиотеки используют наивный алгоритм, потому что он ложится на кеш и векторные инструкции, а Штрассен включают только на очень больших размерах. Исторически он первым пробил барьер n³ и открыл гонку за показателем ω, опущенным сегодня примерно до 2.37.