Почему в matlab сильно замедляется выполнение алгоритма Винограда-Штрассена?

По учебе надо сравнить классический алгоритм умножения матриц и алгоритм Винограда-Штрассена. Мне непонятно, почему алгоритм Винограда-Штрассена прям сильно замедляется по сравнению с классическим. Я понимаю, что MatLab оптимизированнее выполняет именно классический алгоритм, особенно не с супер супер большими матрицами. И что на использование алгоритма Винограда-Штрассена тратится время на выделение в памяти места для промежуточных матриц и все такое. Но просто замедление при размере матрицы очень большое - приложил на картинке. На картинке классический алгоритм при порядке матрицы 2^10 выполняется пару секунд, а алгоритм Винограда-Штрассена 7 минут.

введите сюда описание изображения

Вот код, мб не так алгоритм понял. Параметр type нужен просто для вызова одной функции с алгоритмами Штрассена и Винограда-Штрассена. Больше вродь нет особенностей

clc;

function AA = split_to_2x2_blocks(matrix)  % разделение матрицы на 4 части
    buf = width(matrix);
    A11 = matrix(1:buf/2, 1:buf/2);
    A12 = matrix(1:buf/2, buf/2+1:buf);
    A21 =  matrix(buf/2+1:buf, 1:buf/2);
    A22 = matrix(buf/2+1:buf, buf/2+1:buf);
    
    AA = {A11, A12
            A21, A22};
end

function CC = strassen_mul_2x2(lb, rb, type) % умножение матрицы после разделения, рекурсивная чушь тут и все такое

    d = strassen_mul(cell2mat(lb(1,1)) + cell2mat(lb(2,2)), cell2mat(rb(1,1)) + cell2mat(rb(2,2)), type);
    d_1 = strassen_mul(cell2mat(lb(1,2)) - cell2mat(lb(2,2)), cell2mat(rb(2,1)) + cell2mat(rb(2,2)), type);
    d_2 = strassen_mul(cell2mat(lb(2,1)) - cell2mat(lb(1,1)), cell2mat(rb(1,1)) + cell2mat(rb(1,2)), type);
    left = strassen_mul(cell2mat(lb(2,2)), cell2mat(rb(2,1)) - cell2mat(rb(1,1)), type);
    right = strassen_mul(cell2mat(lb(1,1)), cell2mat(rb(1,2)) - cell2mat(rb(2,2)), type);
    top = strassen_mul(cell2mat(lb(1,1)) + cell2mat(lb(1,2)), cell2mat(rb(2,2)), type);
    bottom = strassen_mul(cell2mat(lb(2,1)) + cell2mat(lb(2,2)), cell2mat(rb(1,1)), type);
    
    c1 = d + d_1 + left - top;
    c2 = right + top;
    c3 = left + bottom;
    c4 = d + d_2 + right - bottom;
    CC = [c1, c2
            c3, c4];
    
end

function CC = vinograd_strassen_mul_2x2(lb, rb, type) 
    s1 = cell2mat(lb(2,1)) + cell2mat(lb(2,2));
    s2 = s1 - cell2mat(lb(1,1));
    s3 = cell2mat(lb(1,1)) - cell2mat(lb(2,1));
    s4 = cell2mat(lb(1,2)) - s2;
    s5 = cell2mat(rb(1,2)) - cell2mat(rb(1,1));
    s6 = cell2mat(rb(2,2)) - s5;
    s7 = cell2mat(rb(2,2)) - cell2mat(rb(1,2));
    s8 = s6 - cell2mat(rb(2,1));
    p1 = strassen_mul(s2, s6, type);
    p2 = strassen_mul(cell2mat(lb(1,1)), cell2mat(rb(1,1)), type);
    p3 = strassen_mul(cell2mat(lb(1,2)), cell2mat(rb(2,1)), type);
    p4 = strassen_mul(s3, s7, type);
    p5 = strassen_mul(s1, s5, type);
    p6 = strassen_mul(s4, cell2mat(rb(2,2)), type);
    p7 = strassen_mul(cell2mat(lb(2,2)), s8, type);
    t1 = p1 + p2;
    t2 = t1 + p4;
    CC = [p2 + p3, t1 + p5 + p6
            t2 - p7, t2 + p5];
    
end

function c = default_mul(left, right) % классический метод умножения матриц
    c = zeros(length(left));
    
    for ii = 1:length(left)
        for jj = 1:length(right)
            for rr = 1:length(left)

                c(ii, jj) = c(ii, jj) + left(ii, rr)* right(rr, jj);
            end
        end
    end
    
end

function c = strassen_mul(left, right, type ) % типо сюда пихаем матрицу и её будущие части, а потом решаем
    if length(left) == 2
        c = default_mul(left, right);
    else
        %if type == 0
            
        %    c = strassen_mul_2x2(split_to_2x2_blocks(left), split_to_2x2_blocks(right), type);
        %elseif type == 1
            
            c = vinograd_strassen_mul_2x2(split_to_2x2_blocks(left), split_to_2x2_blocks(right), type);
        %else
        %    disp('Не выбран вариант')
        %end
    end
end
% 
% a = [1 2 3 4
%     5 6 7 8
%     9 10 11 12
%     13 14 15 16];
% b = [1 2 3 4
%     5 6 7 8
%     9 10 11 12
%     13 14 15 16];
% a = randi([1, 10], 8);
% b = randi([1, 10], 8);
% count_mas = [0, 0, 0]; % сложение, умножение, умножение матриц
% disp(a)
% disp(b)
% c = strassen_mul(a,b);
% disp(c)
% disp(sum)
% disp(umnoj)
% disp(umnoj_mat)

max_stepen = 6;
all_time = zeros(max_stepen, 2);
for t = 1:max_stepen
    step = 2^t;
    a = randi([1,10], step);
    b = randi([1,10], step);
    tic
    c = default_mul(a, b);
    all_time(t, 1) = toc;

    tic
    c = strassen_mul(a,b, 1);
    all_time(t, 2) = toc;
end
t = 1:max_stepen;
hold on
grid on
title("График времени работы без фоновых процессов")
plot(t, all_time(:, 1),'-gs',  'LineWidth',2, 'MarkerSize',10,'Color', [0 0.4470 0.7410])
plot(t,all_time(:, 2),':gs','LineWidth',2,'MarkerSize',10, 'Color', [0.9290 0.6940 0.1250])
xlabel('Порядок матриц')
ylabel('Время работы')
legend("Классический метод", "Метод Винограда-Штрассена")

Ответы (1 шт):

Автор решения: Сергей Галкин

Вот вариант кода. Я не тестировал его на скорость исполнения.

ma=size(A);
mb=size(B);
k=2^(ceil(log2(max([ma mb]))));
Z=zeros(k);
Ar=Z;
Ar(1:ma(1),1:ma(2))=A;
Br=Z;
Br(1:mb(1),1:mb(2))=B;
A11=Ar(1:k/2,1:k/2);
A21=Ar(k/2+1:end,1:k/2);
A12=Ar(1:k/2,k/2+1:end);
A22=Ar(k/2+1:end,k/2+1:end);
B11=Br(1:k/2,1:k/2);
B21=Br(k/2+1:end,1:k/2);
B12=Br(1:k/2,k/2+1:end);
B22=Br(k/2+1:end,k/2+1:end);
S1=A21+A22;
S2=S1-A11;
S3=A11-A21;
S4=A12-S2;
S5=B12-B11;
S6=B22-S5;
S7=B22-B12;
S8=S6-B21;
P1=S2*S6;
P2=A11*B11;
P3=A12*B21;
P4=S3*S7;
P5=S1*S5;
P6=S4*B22;
P7=A22*S8;
T1=P1+P2;
T2=T1+P4;
C=[P2+P3,T1+P5+P6;T2-P7,T2+P5];
C=C(1:ma(1),1:mb(2));
→ Ссылка