Алгоритм Штрассена
Войти  |  Регистрация
Авторизация
» » Алгоритм Штрассена

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



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

В отличие от традиционного алгоритма умножения матриц (по формуле c i , k = ∑ a i , j b j , k {displaystyle c_{i,k}=sum a_{i,j}b_{j,k}} ), работающего за время Θ ( n log 2 ⁡ 8 ) = Θ ( n 3 ) {displaystyle Theta (n^{log _{2}8})=Theta (n^{3})} , алгоритм Штрассена умножает матрицы за время Θ ( n log 2 ⁡ 7 ) = O ( n 2.81 ) {displaystyle Theta (n^{log _{2}7})=O(n^{2.81})} , что даёт выигрыш на больших плотных матрицах начиная, примерно, от 64×64.

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

Алгоритм

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

C = A B A , B , C ∈ R 2 n × 2 n {displaystyle mathbf {C} =mathbf {A} mathbf {B} qquad mathbf {A} ,mathbf {B} ,mathbf {C} in R^{2^{n} imes 2^{n}}}

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

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

A = [ A 1 , 1 A 1 , 2 A 2 , 1 A 2 , 2 ]  ,  B = [ B 1 , 1 B 1 , 2 B 2 , 1 B 2 , 2 ]  ,  C = [ C 1 , 1 C 1 , 2 C 2 , 1 C 2 , 2 ] {displaystyle mathbf {A} ={egin{bmatrix}mathbf {A} _{1,1}&mathbf {A} _{1,2}mathbf {A} _{2,1}&mathbf {A} _{2,2}end{bmatrix}}{mbox{ , }}mathbf {B} ={egin{bmatrix}mathbf {B} _{1,1}&mathbf {B} _{1,2}mathbf {B} _{2,1}&mathbf {B} _{2,2}end{bmatrix}}{mbox{ , }}mathbf {C} ={egin{bmatrix}mathbf {C} _{1,1}&mathbf {C} _{1,2}mathbf {C} _{2,1}&mathbf {C} _{2,2}end{bmatrix}}}

где

A i , j , B i , j , C i , j ∈ R 2 n − 1 × 2 n − 1 {displaystyle mathbf {A} _{i,j},mathbf {B} _{i,j},mathbf {C} _{i,j}in R^{2^{n-1} imes 2^{n-1}}}

тогда

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

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

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

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

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

C 1 , 1 = P 1 + P 4 − P 5 + P 7 {displaystyle mathbf {C} _{1,1}=mathbf {P} _{1}+mathbf {P} _{4}-mathbf {P} _{5}+mathbf {P} _{7}} C 1 , 2 = P 3 + P 5 {displaystyle mathbf {C} _{1,2}=mathbf {P} _{3}+mathbf {P} _{5}} C 2 , 1 = P 2 + P 4 {displaystyle mathbf {C} _{2,1}=mathbf {P} _{2}+mathbf {P} _{4}} C 2 , 2 = P 1 − P 2 + P 3 + P 6 {displaystyle 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 END IF 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) END DO END DO END DO RETURN END SUBROUTINE MUL2 END MODULE STRASSEN_MUL

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

Дальнейшее развитие

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

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

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

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

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

Вычисляются промежуточные матрицы S1S8, P1P7, T1T2

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


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


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

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

C 1 , 1 := P 2 + P 3 {displaystyle mathbf {C} _{1,1}:=mathbf {P} _{2}+mathbf {P} _{3}} C 1 , 2 := T 1 + P 5 + P 6 {displaystyle mathbf {C} _{1,2}:=mathbf {T} _{1}+mathbf {P} _{5}+mathbf {P} _{6}} C 2 , 1 := T 2 − P 7 {displaystyle mathbf {C} _{2,1}:=mathbf {T} _{2}-mathbf {P} _{7}} C 2 , 2 := T 2 + P 5 {displaystyle mathbf {C} _{2,2}:=mathbf {T} _{2}+mathbf {P} _{5}}
Добавлено Admin 8-12-2020, 06:32 Просмотров: 42
Добавить комментарий
Ваше Имя:
Ваш E-Mail:
  • bowtiesmilelaughingblushsmileyrelaxedsmirk
    heart_eyeskissing_heartkissing_closed_eyesflushedrelievedsatisfiedgrin
    winkstuck_out_tongue_winking_eyestuck_out_tongue_closed_eyesgrinningkissingstuck_out_tonguesleeping
    worriedfrowninganguishedopen_mouthgrimacingconfusedhushed
    expressionlessunamusedsweat_smilesweatdisappointed_relievedwearypensive
    disappointedconfoundedfearfulcold_sweatperseverecrysob
    joyastonishedscreamtired_faceangryragetriumph
    sleepyyummasksunglassesdizzy_faceimpsmiling_imp
    neutral_faceno_mouthinnocent