QR-разложение квадратной матрицы

Любая матрица A размера NxN может быть представлена в виде произведения ортонормированной матрицы Q размера NxN (det Q = 1) и верхнетреугольной матрицы R размера NxN.

Такое представление позволяет легко находить обратную матрицу A ^-1, вычислять определитель, решать системы уравнений вида Ax = b. Существуют и другие применения QR-разложения.

QR-разложение матрицы A строится при помощи преобразований вращения. Обходя матрицу, мы применяем к каждому из элементов, расположенных ниже главной диагонали, преобразование вращения, которое обнуляет этот элемент. Если мы обнуляем элемент a_i,j , то осуществляем поворот вокруг i-ой и j-ой координатных осей. За счет поворота элемент a_i,j  обнуляется, а элемент a_j,j  изменяет свое значение. Попутно меняются другие элементы i-ой и j-ой строки. Таким образом мы получаем матрицу R. Это же преобразование мы применяем к единичной матрице (единичной она является только до первого преобразования), чтобы получить матрицу Q.

Немного о реализации алгоритма - по соображениям удобства все преобразования проводятся "на месте", т.е. результирующая матрица R замещает переданную матрицу A, а матрица Q помещается в отдельную переменную.

Если нашли ошибку в алгоритме - сообщите!

Реализации алгоритма на различных языках:

Реализация алгоритма на AlgoPascal:

unit QRDecompositionUnit;

interface QRDecomposition;

implementation

(*******************************************************************
Процедура QR-разложения квадратной матрицы A.

Матрица  А  представляется  в  виде  произведения  ортонормированной
матрицы Q (с определителем, равным 1)  и  верхнетреугольной  матрицы
R: A=QR.

Входные данные:
    A:  array [1..N, 1..N] - разлагаемая матрица
    N:  - размер матрицы

Выходные данные:
    A:  array [1..N, 1..N] - верхнетреугольная матрица R
    Q:  array [1..N, 1..N] - ортонормированная матрица

*******************************************************************)
procedure QRDecomposition(
    var     A:  array of array of Real;
    const   N:  Integer;
    out     Q:  array of array of Real);
var
    Col     :   Integer;
    Row     :   Integer;
    I       :   Integer;
    J       :   Integer;
    X1      :   Real;
    X2      :   Real;
    T       :   Real;
    T1      :   Real;
    T2      :   Real;
    C       :   Real;
    S       :   Real;
    Denomin :   Real;
begin
    //инициализация Q-матрицы
    SetBounds(Q, [1,N], [1,N]);
    for I:=1 to N do
        for J:=1 to N do
            if I=J then
                Q[I,J]:=1
            else
                Q[I,J]:=0;

    //вращение
    for Col:=1 to N do
        for Row:=Col+1 to N do
        begin
            //выбор угла вращения
            X1:=A[Col,Col];
            X2:=A[Row,Col];
            if (X1=0) and (X2=0) then
                Continue;
            if AbsReal(X1)>AbsReal(X2) then
            begin
                T:=X2/X1;
                Denomin:=AbsReal(X1)*Sqrt(1+T*T)
            end
            else
            begin
                T:=X1/X2;
                Denomin:=AbsReal(X2)*Sqrt(1+T*T);
            end;
            C:=X1/Denomin;
            S:=-X2/Denomin;

            //поворот матрицы коэффициентов и обратной матрицы
            for I:=1 to N do
            begin
                T1:=A[Col,I];
                T2:=A[Row,I];
                A[Col,I]:=C*T1-S*T2;
                A[Row,I]:=S*T1+C*T2;
                T1:=Q[I,Col];
                T2:=Q[I,Row];
                Q[I,Col]:= C*T1+S*T2;
                Q[I,Row]:=-S*T1+C*T2;
            end;
        end;
end;

end.

Бочканов Сергей, Быстрицкий Владимир
Copyright © 1999-2004
При поддержке проекта MANUAL.RU