QRS разложение вращений - PullRequest
Новая
       9

QRS разложение вращений

1 голос
/ 12 февраля 2020

Я пытаюсь создать функцию, которая вычисляет QR-разложение Givens Rotation, следуя этому псевдокоду.

enter image description here

function [Q,R] = givens(A)
[m,n] = size(A);
indexI = zeros(m,n);
indexJ = zeros(m,n);
C = zeros(m,n);
S = zeros(m,n);
for i = 1:n
    for j = i+1:m
        c = A(i,i)/((A(i,i))^2 + (A(j,i)^2))^0.5;
        s = A(j,i)/((A(i,i))^2 + (A(j,i)^2))^0.5;
        A(i,:) = c*A(i,:) + s*A(j,:);
        A(j,:) = -s*A(i,:) + c*A(j,:);
        indexI(j,i) = i;
        indexJ(j,i) = j;
        C(j,i) = c;
        S(j,i) = s;
    end
end
R = A;
Q = eye(m);
for i = 1:n
    for j= j+1:m
        Q(:,i) = c*Q(:,i) + s*Q(:,j);
        Q(:,j) = -s*Q(:,i) + c*Q(:,j);
    end
end

Однако , матрица R, которую я получаю, не является верхней три angular. Я не могу найти ошибку здесь.

1 Ответ

0 голосов
/ 22 апреля 2020

На слайдах есть некоторая неоднозначность.

Вращение Гивенса фактически выполняет умножение матрицы на две строки за раз.

Предположим, [ri;rj] - это ваши две строки, а Q - соответствующая matirx поворота данных.

Обновление [ri; rj] = Q*[ri; rj], но в вашем коде вы сначала обновляете ri, а затем используете обновленный ri для немедленного обновления rj. 

    A(i,:) = c*A(i,:) + s*A(j,:);
    A(j,:) = -s*A(i,:) + c*A(j,:);

Вот исправить ошибку: 

    tmp = c*A(i,:) + s*A(j,:);
    A(j,:) = -s*A(i,:) + c*A(j,:);
    A(i,:) = tmp;
...