Как повысить эффективность кода умножения матриц в MATLAB при выполнении факторизации LU? - PullRequest
0 голосов
/ 06 октября 2018

Я попытался реализовать факоризацию LU на страницах 3-6 этого документа:

http://www4.ncsu.edu/~kksivara/ma505/handouts/lu-pivot.pdf

Пример выполнения приведенного ниже кода для простой матрицы 4x4:

A=[2 1 1 0; 4 3 3 1; 8 7 9 5 ; 6 7 9 8];
intch=[];flag=0;
>> [Bp]=Gauss2(A,4,intch,flag);
>> B=lu(A);
>> Bp

Bp =

    8.0000    7.0000    9.0000    5.0000
    0.7500    1.7500    2.2500    4.2500
    0.5000   -0.2857   -0.8571   -0.2857
    0.2500   -0.4286    0.3333    0.6667

>> B

B =

    8.0000    7.0000    9.0000    5.0000
    0.7500    1.7500    2.2500    4.2500
    0.5000   -0.2857   -0.8571   -0.2857
    0.2500   -0.4286    0.3333    0.6667

Как вы можете видеть, матрицы Bp и B одинаковы, что доказывает, что мой алгоритм факторизации LU работает правильно для матрицы 4x4 A. Я изменил A на 60x60 и 256x256, оба из которых такжеработать правильно (т.е. Bp == B).Однако сейчас я запускаю этот код в наборе данных с размером 4000x4000.Моя программа не заканчивается (даже после 10 часов работы в течение ночи).В MATLAB, когда программа работает, я могу нажать кнопку «Пауза».Это переводит программу в режим отладки и сообщает мне точную строку кода, на которую программа тратит много времени.Оказывается, в приведенной ниже функции «Gauss2» программа тратит много времени на два вложенных цикла for в конце функции (см. Ниже, где я определил «BOTTLENECK» в комментариях).

Как мне улучшить скорость умножения матриц здесь?Если вы идете в самом низу стр.5 из http://www4.ncsu.edu/~kksivara/ma505/handouts/lu-pivot.pdf вы увидите уравнения 21.6 и 21.7, которые говорят:

L_k'=P(m-1)...P(k+1)L_k*P(k+1)...*P(m+1)

и

L=inverse(L(m-1)' .... L(2)'L(1)')

Это именно то, что две вложенные для циклов ниже пытаются сделатьделать.Есть ли более умный способ вычислить L_k '?

function [Bp]= Gauss2(A, n, intch, flag)

    intch=zeros(n,1);
    flag=0;
    a=zeros(n,n);

    for k=1:n-1
        [amax,amax_index] = max(abs(A(k:n,k)));
        amax_index=amax_index+k-1;
        if (amax == 0)
            flag=1; % matrix is singular
            intch(k)=0;
        else
            intch(k)=amax_index;
            if (amax_index ~= k)
                % pivoting
                % interchange rows k and m
                A([k,amax_index],:) = A([amax_index,k],:);
            end

            % calculate multipliers store, and store in "a"
            for i=k+1:n
                a(i,k)=-A(i,k)/A(k,k);
            end

            % perform elimination for column "k"
            for i=k+1:n
                A(i,:)=A(i,:)+a(i,k)*A(k,:);

            end 
        end   
    end
    if (A(n,n)==0)
        flag=1;
        intch(n)=0;
    else
        intch(n)=n;
    end

    U=A; % this is the uppper matrix 

    % next, need to determine "L_k'" from algorithm 
    % for all k
    prevX=eye(n);
    for k=1:n-1
        X=eye(n);
        for j=(n-1):-1:(k+1)
            X=X*P(j,intch(j),n); % <---- BOTTLENECK HERE
        end
        X=X*L(a,k,n);
        for j=(k+1):(n-1)
            X=X*P(j,intch(j),n); % < --- BOTTLENECK HERE
        end

        prevX=X*prevX;

    end


    Lower=inv(prevX);

    Bp=U+Lower-eye(n);

end

Вот мой код для Pm и Lm ниже.Эти функции используются в коде выше для генерации матриц перестановок P_1, P_2, ..., P_k и L_1, L_2, .... L_k.Мне нужно вычислить их, чтобы выяснить матрицы L_k ', которые затем могут быть использованы для определения окончательной L нижней матрицы.

function [Y]=L(a,k,n)
    Y=eye(n);
    Y((k+1):n,k)=a(k+1:n,k);
end


function [X]=P(i,j,n)
    X=eye(n);
    X([i,j],:) = X([j,i],:);
end
...