Разложение LU с частичным поворотом Matlab

Я пытаюсь реализовать свою собственную декомпозицию LU с частичным поворотом. Мой код ниже и вроде работает нормально, но для некоторых матриц дает разные результаты при сравнении со встроенной функцией [L, U, P] = lu(A) в матлабе

Может ли кто-нибудь заметить, где это неправильно?

function [L, U, P] = lu_decomposition_pivot(A)
    n = size(A,1);
    Ak = A;
    L = zeros(n);
    U = zeros(n);
    P = eye(n);
    for k = 1:n-1
        for i = k+1:n
            [~,r] = max(abs(Ak(:,k)));

            Ak([k r],:) = Ak([r k],:);
            P([k r],:) = P([r k],:);

            L(i,k) = Ak(i,k) / Ak(k,k);
            for j = k+1:n
                U(k,j-1) = Ak(k,j-1);
                Ak(i,j) = Ak(i,j) - L(i,k)*Ak(k,j);
            end
        end
    end
    L(1:n+1:end) = 1;
    U(:,end) = Ak(:,end);
return

Вот две матрицы, с которыми я тестировал. Первый правильный, а во втором некоторые элементы перевернуты.

A = [1 2 0; 2 4 8; 3 -1 2];

A = [0.8443 0.1707 0.3111;
     0.1948 0.2277 0.9234;
     0.2259 0.4357 0.4302];

ОБНОВЛЕНИЕ

Я проверил свой код и исправил некоторые ошибки, но с частичным поворотом чего-то не хватает. В первом столбце последние две строки всегда инвертируются (по сравнению с результатом lu() в Matlab)

function [L, U, P] = lu_decomposition_pivot(A)
    n = size(A,1);
    Ak = A;
    L = eye(n);
    U = zeros(n);
    P = eye(n);
    for k = 1:n-1
        [~,r] = max(abs(Ak(k:end,k)));
        r = n-(n-k+1)+r;
        Ak([k r],:) = Ak([r k],:);
        P([k r],:) = P([r k],:);
        for i = k+1:n
            L(i,k) = Ak(i,k) / Ak(k,k);
            for j = 1:n
                U(k,j) = Ak(k,j);
                Ak(i,j) = Ak(i,j) - L(i,k)*Ak(k,j);
            end
        end
    end
    U(:,end) = Ak(:,end);
return

matlab matrix decomposition
person BRabbit27    schedule 08.03.2013    source источник

Ответы (3)


arrow_upward
3
arrow_downward

Я забыл, что если в матрице P была замена, мне пришлось поменять местами и матрицу L. Так что просто добавьте следующую строку после замены P, и все будет работать отлично.

L([k r],:) = L([r k],:);
person BRabbit27    schedule 09.03.2013
comment
Можете ли вы указать, где именно вы поместили эту строку? Я проверил ваш код, и есть проблема с матрицей L - person CTZStef; 10.01.2014
comment
Должно быть сразу после P([k r],:) = P([r k],:); - person BRabbit27; 10.01.2014

arrow_upward
3
arrow_downward

Обе функции некорректны. Вот правильный.

function [L, U, P] = LU_pivot(A)
    [m, n] = size(A); L=eye(n); P=eye(n); U=A;
    for k=1:m-1
        pivot=max(abs(U(k:m,k)))
        for j=k:m
            if(abs(U(j,k))==pivot)
                ind=j
                break;
            end
        end
        U([k,ind],k:m)=U([ind,k],k:m)
        L([k,ind],1:k-1)=L([ind,k],1:k-1)
        P([k,ind],:)=P([ind,k],:)
        for j=k+1:m
            L(j,k)=U(j,k)/U(k,k)
            U(j,k:m)=U(j,k:m)-L(j,k)*U(k,k:m)
        end
        pause;
    end
end
person panch    schedule 06.08.2014

arrow_upward
1
arrow_downward

Мой ответ здесь:

function [L, U, P] = LU_pivot(A)
[n, n1] = size(A); L=eye(n); P=eye(n); U=A;
for j = 1:n
  [pivot m] = max(abs(U(j:n, j)));     
  m = m+j-1;
  if m ~= j
    U([m,j],:) =  U([j,m], :);   % interchange rows m and j in U
    P([m,j],:) =  P([j,m], :);   % interchange rows m and j in P
    if j >= 2;    % very_important_point
      L([m,j],1:j-1) =  L([j,m], 1:j-1);   % interchange rows m and j in columns 1:j-1 of L
    end;
  end
  for i = j+1:n      
    L(i, j) = U(i, j) / U(j, j);
    U(i, :) =  U(i, :) - L(i, j)*U(j, :);
  end
end
person Dr Karatay    schedule 04.03.2014