Я попытался реализовать факоризацию 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