qr_givens.m

基于Givens transformation的QR分解。其使用需要调用Givens transformation程序。

function [r,q] = qr_givens(x)

% purpose: 用Givens平面旋转变换计算QR分解，x = qr，size(x)=[m,n]
% 目前只适用于 m ≥ n 的情况，且size(q) = size(x)
% Data: 31/3/2008
% Version 1.0
% Copyright (C) 2008, Yang Chun
% 
% 输入 - x 为列满秩矩阵或标量
% 输出 - q 为x的正交基
%      - r 为正交基的坐标，为上三角矩阵
% 
% r = qr_givens(x) 返回正交基坐标
% [r,q] = qr_givens(x) 返回正交基与正交基坐标

[m,n] = size(x);
if m > n
    iterend = n;
else
    iterend = n - 1;
end
for k1 = 1:iterend
    for k2 = m:-1:k1+1
        [c,s,p] = givens([x(k1,k1);x(k2,k1)]);
        x([k1,k2],k1:n) = [c,s;-s,c]*x([k1,k2],k1:n);
        x(k2,k1) = p;
    end
end

if nargout == 1 
    r = triu(x(1:n,:));
else
    r = triu(x(1:n,:));
    q = eye(m,n);
    for k1 = iterend:-1:1
        for k2 = k1+1:m
            if x(k2,k1) == 1
                c = 0; s = 1;
            elseif abs(x(k2,k1)) < 1
                s = 2*x(k2,k1);
                c = (1-s^2)^0.5;
            else
                c = 2/x(k2,k1);
                s = (1-c^2)^0.5;
            end
            q([k1,k2],k1:n) = [c,-s;s,c]*q([k1,k2],k1:n);
        end
    end
end