Множення матриць за Штрассеном — від O(n³) до субкубічних алгоритмів
Шкільний алгоритм множення n×n матриць потребує n³ множень. У 1969 Штрассен показав, що матриці 2×2 можна перемножити лише 7 множеннями замість 8, а рекурсія дає складність O(nlog₂7) ≈ O(n2.807) — субкубічну економію, на якій ґрунтуються сучасні бібліотеки числових обчислень.
1. Наївний алгоритм O(n³)
Множення матриць C = A·B для матриць n×n: C[i][j] = Σₖ A[i][k] · B[k][j]. Трициклова реалізація — це O(n³) множень і O(n³) додавань:
function matmul(A, B, n) {
const C = Array.from({length: n}, () => new Float64Array(n));
for (let i = 0; i < n; i++)
for (let k = 0; k < n; k++) // k-цикл усередині задля кешу
for (let j = 0; j < n; j++)
C[i][j] += A[i][k] * B[k][j];
return C;
}
Зверніть увагу на порядок циклів i-k-j (зовнішній-k-внутрішній-j): він
звертається до B рядок за рядком (дружньо до кешу), тоді як
i-j-k звертається до B стовпець за стовпцем (вороже до кешу). Для n
= 1000 наївний алгоритм виконує 10⁹ множень — на що йде
близько 1 секунди на скалярному процесорі.
2. «Розділяй і володарюй»
Розбиваємо кожну матрицю n×n на чотири блочні матриці (n/2)×(n/2):
Розбиття на 4 блоки з рекурсією структурно еквівалентне наївному підходу. Ключова ідея в роботі Штрассена: зменшити кількість рекурсивних підзадач з 8 до 7.
3. Трюк Штрассена із 7 множеннями
Штрассен показав, що чотири добутки C₁₁, C₁₂, C₂₁, C₂₂ можна обчислити лише з 7 проміжних добутків M₁…M₇:
Ці тотожності можна перевірити прямим розкриттям. Вартість — 7 множень і 18 додавань на рівень. Оскільки додавання є Θ(n²) і дешеві порівняно з множеннями (на великих матрицях), сукупний результат загалом швидший.
4. Аналіз складності
5. Кеш-нечутливі алгоритми
Сучасні процесори з ієрархіями кешу L1/L2/L3 штрафують випадковий доступ до пам'яті. Блочне (плиткове) множення матриць покращує локальність даних, працюючи з малими підблоками B×B, що повністю вміщаються в кеш L1 або L2:
Кеш-нечутливий алгоритм досягає оптимальної продуктивності кешу, не знаючи його розміру, шляхом рекурсивного поділу доти, доки підзадача не вміститься в кеш. Формулювання Штрассена у стилі «розділяй і володарюй» природно є кеш-нечутливим.
6. BLAS і практична реалізація
Стандарт Basic Linear Algebra Subprograms (BLAS) визначає інтерфейс для операцій з матрицями. DGEMM (загальне множення матриць подвійної точності) — це базова процедура:
- OpenBLAS, MKL, BLIS: вручну налаштовані бібліотеки BLAS, що використовують SIMD (AVX-512 виконує 16 множень подвійної точності за цикл), багатопотоковість і блокування, оптимізоване під конкретні розміри кешу процесора.
- LAPACK: процедури вищого рівня (LU-розклад, власні значення), побудовані на BLAS. LU-розклад має складність O(n³), у якій домінують виклики DGEMM.
- cuBLAS, rocBLAS: BLAS для GPU, що задіює тисячі ядер CUDA/HIP. Пікова DGEMM на GPU зазвичай досягає 30–100 TFLOPS для великих матриць.
Множення матриць у глибинному навчанні (пакетне, змішаної точності) обробляється cuBLAS-lt та тензорними ядрами — спеціалізованими блоками множення матриць 4×4 у графічних процесорах NVIDIA, що досягають у 2–8× вищої пропускної здатності для операцій FP16/BF16, критичної для навчання трансформерів.
7. Гонка до ω = 2
Показник множення матриць ω визначається як інфімум усіх дійсних чисел α таких, що дві матриці n×n можна перемножити за O(n^α) операцій у полі. Прогрес:
Родина алгоритмів Копперсміт-Віноград (та її нащадки) використовує витончений алгебраїчний апарат — розклад тензорного рангу, лазерний метод, теоретико-групові техніки — породжуючи алгоритми настільки складні, що їх рідко реалізують на практиці. Точка перетину зі Штрассеном астрономічно велика (n ~ 1010 або більше). Але теоретично вони підтверджують, що множення матриць ближче до «лінійної» складності, ніж до «кубічної».