Алгоритм Штрассена

Материал из Википедии — свободной энциклопедии
Перейти к: навигация, поиск

Алгоритм Штрассена предназначен для быстрого умножения матриц. Он был разработан Фолькером Штрассеном в 1969 году и является обобщением метода умножения Карацубы на матрицы.

В отличие от традиционного алгоритма умножения матриц (по формуле cik = Σaijbjk), работающего за время Θ(n³) = Θ(nlog2 8), алгоритм Штрассена умножает матрицы за время Θ(nlog2 7) = Θ(n2.81), что даёт выигрыш на больших плотных матрицах начиная, примерно, от 64×64.

Несмотря на то, что алгоритм Штрассена является асимптотически не самым быстрым из существующих алгоритмов быстрого умножения матриц, он проще программируется и эффективнее при умножении матриц относительно малого размера, поэтому именно он чаще используется на практике.

Алгоритм[править | править вики-текст]

Пусть A, B — две квадратные матрицы над кольцом R. Матрица C получается по формуле:

\mathbf{C} = \mathbf{A} \mathbf{B} \qquad \mathbf{A},\mathbf{B},\mathbf{C} \in R^{2^n \times 2^n}

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

Разделим матрицы A, B и C на равные по размеру блочные матрицы

 
\mathbf{A} =
\begin{bmatrix}
\mathbf{A}_{1,1} & \mathbf{A}_{1,2} \\
\mathbf{A}_{2,1} & \mathbf{A}_{2,2}
\end{bmatrix}
\mbox { , }
\mathbf{B} =
\begin{bmatrix}
\mathbf{B}_{1,1} & \mathbf{B}_{1,2} \\
\mathbf{B}_{2,1} & \mathbf{B}_{2,2}
\end{bmatrix}
\mbox { , }
\mathbf{C} =
\begin{bmatrix}
\mathbf{C}_{1,1} & \mathbf{C}_{1,2} \\
\mathbf{C}_{2,1} & \mathbf{C}_{2,2}
\end{bmatrix}

где

\mathbf{A}_{i,j}, \mathbf{B}_{i,j}, \mathbf{C}_{i,j} \in R^{2^{n-1} \times 2^{n-1}}

тогда

\mathbf{C}_{1,1} = \mathbf{A}_{1,1} \mathbf{B}_{1,1} + \mathbf{A}_{1,2} \mathbf{B}_{2,1}
\mathbf{C}_{1,2} = \mathbf{A}_{1,1} \mathbf{B}_{1,2} + \mathbf{A}_{1,2} \mathbf{B}_{2,2}
\mathbf{C}_{2,1} = \mathbf{A}_{2,1} \mathbf{B}_{1,1} + \mathbf{A}_{2,2} \mathbf{B}_{2,1}
\mathbf{C}_{2,2} = \mathbf{A}_{2,1} \mathbf{B}_{1,2} + \mathbf{A}_{2,2} \mathbf{B}_{2,2}

Однако с помощью этой процедуры нам не удалось уменьшить количество умножений. Как и в обычном методе, нам требуется 8 умножений.

Теперь определим новые элементы

\mathbf{P}_{1} := (\mathbf{A}_{1,1} + \mathbf{A}_{2,2}) (\mathbf{B}_{1,1} + \mathbf{B}_{2,2})
\mathbf{P}_{2} := (\mathbf{A}_{2,1} + \mathbf{A}_{2,2}) \mathbf{B}_{1,1}
\mathbf{P}_{3} := \mathbf{A}_{1,1} (\mathbf{B}_{1,2} - \mathbf{B}_{2,2})
\mathbf{P}_{4} := \mathbf{A}_{2,2} (\mathbf{B}_{2,1} - \mathbf{B}_{1,1})
\mathbf{P}_{5} := (\mathbf{A}_{1,1} + \mathbf{A}_{1,2}) \mathbf{B}_{2,2}
\mathbf{P}_{6} := (\mathbf{A}_{2,1} - \mathbf{A}_{1,1}) (\mathbf{B}_{1,1} + \mathbf{B}_{1,2})
\mathbf{P}_{7} := (\mathbf{A}_{1,2} - \mathbf{A}_{2,2}) (\mathbf{B}_{2,1} + \mathbf{B}_{2,2})

которые затем используются для выражения Ci, j. Таким образом, нам нужно всего 7 умножений на каждом этапе рекурсии. Элементы матрицы C выражаются из Pk по формулам

\mathbf{C}_{1,1} = \mathbf{P}_{1} + \mathbf{P}_{4} - \mathbf{P}_{5} + \mathbf{P}_{7}
\mathbf{C}_{1,2} = \mathbf{P}_{3} + \mathbf{P}_{5}
\mathbf{C}_{2,1} = \mathbf{P}_{2} + \mathbf{P}_{4}
\mathbf{C}_{2,2} = \mathbf{P}_{1} - \mathbf{P}_{2} + \mathbf{P}_{3} + \mathbf{P}_{6}

Рекурсивный процесс продолжается n раз, до тех пор пока размер матриц Ci, j не станет достаточно малым, далее используют обычный метод умножения матриц. Это делают из-за того, что алгоритм Штрассена теряет эффективность по сравнению с обычным на малых матрицах в силу большего числа сложений. Оптимальный размер матриц для перехода к обычному умножению зависит от характеристик процессора и языка программирования и на практике лежит в пределах от 32 до 128.

Пример реализации на Фортране[править | править вики-текст]

Приведён пример реализации алгоритма Штрассена на Фортране. Предполагается, что все матрицы квадратные, их размер является степенью 2.

MODULE STRASSEN_MUL
 
CONTAINS
! X = A * B
! V - dimension of matrices
RECURSIVE SUBROUTINE MUL(A, B, V, C)
 
INTEGER, INTENT(IN) :: V
DOUBLE PRECISION, INTENT(IN) :: A( : , : ), B( : , : )
INTEGER :: H ! H = V/2 (see below)
DOUBLE PRECISION, INTENT(OUT) :: C(V, V)
DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: P1, P2, P3, P4, P5, P6, P7
DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: A11, A12, A21, A22, B11, B12, B21, B22
 
IF (V <= 64) THEN ! if dimension equals 64 use MUL2
CALL MUL2 (A, B, V, C)
RETURN
ENDIF
 
H = V/2
 
ALLOCATE (P1(H,H), P2(H,H), P3(H,H), P4(H,H), P5(H,H), P6(H,H), P7(H,H))
ALLOCATE (A11(H,H), A12(H,H), A21(H,H), A22(H,H), B11(H,H), B12(H,H), B21(H,H), B22(H,H))
 
A11 (1:H, 1:H) = A (1:H, 1:H)
A12 (1:H, 1:H) = A (1:H, H+1:V)
A21 (1:H, 1:H) = A (H+1:V, 1:H)
A22 (1:H, 1:H) = A (H+1:V, H+1:V)
 
B11 (1:H, 1:H) = B (1:H, 1:H)
B12 (1:H, 1:H) = B (1:H, H+1:V)
B21 (1:H, 1:H) = B (H+1:V, 1:H)
B22 (1:H, 1:H) = B (H+1:V, H+1:V)
 
!$OMP PARALLEL
CALL MUL(A11 + A22, B11 + B22, H, P1) ! P1 = (A11 + A22) * (B11 + B22)
CALL MUL(A21 + A22, B11, H, P2) ! etc. ...
CALL MUL(A11, B12 - B22, H, P3)
CALL MUL(A22, B21 - B11, H, P4)
CALL MUL(A11 + A12, B22, H, P5)
CALL MUL(A21 - A11, B11 + B12, H, P6)
CALL MUL(A12 - A22, B21 + B22, H, P7)
!$OMP END PARALLEL
 
DEALLOCATE (B11, B12, B21, B22)
DEALLOCATE (A11, A12, A21, A22)
 
C (1:H, 1:H) = P1 + P4 + P7 - P5
C (1:H, H+1:V) = P3 + P5
C (H+1:V, 1:H) = P2 + P4
C (H+1:V, H+1:V) = P1 - P2 + P3 + P6
 
DEALLOCATE (P1, P2, P3, P4, P5, P6, P7)
 
RETURN
END SUBROUTINE MUL
 
! X = A * B using standard method
SUBROUTINE MUL2 (A, B, V, X)
IMPLICIT NONE
INTEGER, INTENT(IN) :: V
DOUBLE PRECISION, INTENT(IN) :: A( : , : ), B( : , : )
DOUBLE PRECISION, INTENT(OUT), DIMENSION (:,:) :: X
INTEGER :: I, J, K
DO I = 1, V
DO J = 1, V
X (I, J) = 0
DO K = 1, V
X (I, J) = X (I, J) + A (I, K) * B (K, J)
ENDDO
ENDDO
ENDDO
RETURN
END SUBROUTINE MUL2
 
END MODULE STRASSEN_MUL

Вычисления промежуточных матриц P1 — P7 можно проводить параллельно при использовании таких библиотек как OpenMP или MPI, что позволяет значительно повысить скорость работы алгоритма на многопроцессорных машинах.

Дальнейшее развитие[править | править вики-текст]

Штрассен был первым, кто показал возможность умножения матриц более эффективным способом, чем стандартный. После публикации его работы в 1969 начались активные поиски более быстрого алгоритма. Самым асимптотически быстрым алгоритмом на сегодняшний день является алгоритм Копперсмита — Винограда, позволяющий перемножать матрицы за {\rm O}(n^{2.376}) операций[1], предложенный в 1987 году и усовершенствованный в 2011 году до уровня {\rm O}(n^{2.3727})[1]. Этот алгоритм не представляет практического интереса в силу астрономически большой константы в оценке арифметической сложности. Вопрос о предельной в асимптотике скорости перемножения матриц не решен. Существует гипотеза Штрассена о том, что для достаточно больших n существует алгоритм перемножения двух матриц размера n \times n за {\rm O}(n^{2+\varepsilon}) операций, где  \varepsilon сколь угодно малое наперед заданное положительное число. Эта гипотеза имеет сугубо теоретический интерес, так как размер матриц для которых  \varepsilon действительно мало по всей видимости очень велик.

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

Алгоритм Винограда-Штрассена[править | править вики-текст]

Существует модификация алгоритма Штрассена, для которой требуется 7 умножений и 15 сложений (вместо 18 для обычного алгоритма Штрассена).

Матрицы A, B и C делятся на блочные подматрицы как показано выше.

Вычисляются промежуточные матрицы S1 — S8, P1 — P7, T1 — T2

\mathbf{S}_{1} := (\mathbf{A}_{2,1} + \mathbf{A}_{2,2})
\mathbf{S}_{2} := (\mathbf{S}_{1} - \mathbf{A}_{1,1})
\mathbf{S}_{3} := (\mathbf{A}_{1,1} - \mathbf{A}_{2,1})
\mathbf{S}_{4} := (\mathbf{A}_{1,2} - \mathbf{S}_{2})
\mathbf{S}_{5} := (\mathbf{B}_{1,2} - \mathbf{B}_{1,1})
\mathbf{S}_{6} := (\mathbf{B}_{2,2} - \mathbf{S}_{5})
\mathbf{S}_{7} := (\mathbf{B}_{2,2} - \mathbf{B}_{1,2})
\mathbf{S}_{8} := (\mathbf{S}_{6} - \mathbf{B}_{2,1})


\mathbf{P}_{1} := \mathbf{S}_{2} \mathbf{S}_{6}
\mathbf{P}_{2} := \mathbf{A}_{1,1} \mathbf{B}_{1,1}
\mathbf{P}_{3} := \mathbf{A}_{1,2} \mathbf{B}_{2,1}
\mathbf{P}_{4} := \mathbf{S}_{3} \mathbf{S}_{7}
\mathbf{P}_{5} := \mathbf{S}_{1} \mathbf{S}_{5}
\mathbf{P}_{6} := \mathbf{S}_{4} \mathbf{B}_{2,2}
\mathbf{P}_{7} := \mathbf{A}_{2,2} \mathbf{S}_{8}


\mathbf{T}_{1} := \mathbf{P}_{1} + \mathbf{P}_{2}
\mathbf{T}_{2} := \mathbf{T}_{1} + \mathbf{P}_{4}

Элементы матрицы C вычисляются по формулам

\mathbf{C}_{1,1} := \mathbf{P}_{2} + \mathbf{P}_{3}
\mathbf{C}_{1,2} := \mathbf{T}_{1} + \mathbf{P}_{5} + \mathbf{P}_{6}
\mathbf{C}_{2,1} := \mathbf{T}_{2} - \mathbf{P}_{7}
\mathbf{C}_{2,2} := \mathbf{T}_{2} + \mathbf{P}_{5}

Примечания[править | править вики-текст]

Литература[править | править вики-текст]

  • Volker Strassen: Gaussian Elimination is not Optimal. In: Numerische Mathemetik, Bd. 13 (1969), S. 354–356, ISSN 00298-599X.
  • Ананий В. Левитин Глава 4. Метод декомпозиции: Умножение больших целых чисел и алгоритм умножения матриц Штрассена // Алгоритмы: введение в разработку и анализ = Introduction to The Design and Analysis of Aigorithms. — М.: «Вильямс», 2006. — С. 189-195. — ISBN 0-201-74395-7.
  • Кормен, Томас Х., Лейзерсон, Чарльз И., Ривест, Рональд Л., Штайн, Клифорд Глава 28. Работа с матрицами // Алгоритмы: построение и анализ = Introduction to Algorithms. — 2-e издание. — М.: «Вильямс», 2005. — С. 833 - 839. — ISBN 5-8459-0857-4.