Givens旋转QR分解

问题描述 投票:0回答:1

我想按照这个伪代码,创建一个计算Givens Rotation QR分解的函数。

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矩阵,并不是上三角。我似乎找不到这里的错误。

numerical-methods
1个回答
0
投票

幻灯片上有一些模糊不清的地方。

吉文斯旋转实际上是每次对两行进行矩阵乘法。

假设 [ri;rj] 是你的两行和 Q 是对应的givens旋转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;
© www.soinside.com 2019 - 2024. All rights reserved.