Matlab: QR разложение комплексных матриц

[WinterLord1]

Добрый день. Занимаюсь QR разложением корреляционных матриц (симметрична, с единичками по диагонали), использую matlab без применения стандартной функции $[Q,R] = qr(A)$. Во время многочисленных попыток написания кода столкнулся с тем, что для входной матрицы вещественных чисел QR преобразование считается правильно (результат совпадает с матлабовской функцией), а при комплексных входных значениях получается неверный результат.
Подскажите пожалуйста, есть ли какая-то хитрость для работы с комплексными значениями, которую я упускаю?
Ниже привожу использованный мной способ для QR разложения:

Используется синтаксис Matlab M

function [ Q, R ] = houseTriangular( A )

[r, ~] = size(A);
Q = eye(r);
k = 1;
for i = 1 : r - 1
    v = zeros(1,r);
    v(k) = A(k,k) + sign(A(k,k)) * Sum(A(:,k), k);
    for j = k + 1 : r
        v(j) = A(j,k);
    end
    H = eye(r) - 2 * ( (v' * v) / (v * v'));
    Q = Q * H;
    A = H * A;
    k = k + 1;
end
R = A;
end

% Функция для расчета суммы и затем взятия квадратного корня 
% элементов столбца матрицы А, по формуле Хаусхолдера, где 
% Входные параметры:  а - вектор-столбец матрицы А
%                     k - номер эл-та столбца, с которого начинается сумма
% Выходные параметры: S - значение квадратного корня из суммы 
%                   
function [S] = Sum(a, k)

S = 0;
l = length(a);
v = zeros(1, l);
for i = k : l
    v(i) = a(i);
end

S = norm(v);

end

[Lia]

Тема перемещена из форума «Помогите решить / разобраться (М)» в форум «Околонаучный софт»

[WinterLord1]

Ок, оказалось не все так просто. Метод отражений в комплексном простанстве принимает несколько иной вид $Q = I -(1+\omega)vv^*$, $\omega = \frac{\overline{v^*x}}{v^*x}$.
Применяя отражение на первом шаге алгоритма всё работает, но как его задействовать на следующих шагах, мне не совсем понятно. Измененный код представлен ниже:

Используется синтаксис Matlab M

function [ Q, R ] = houseTriangular( A )


[r, ~] = size(A);
Q = eye(r);
k = 1;
for i = 1 : r - 1
    
    
    alpha = norm(A(:,i));%Sum(A(:,k), k);
    e = zeros(1, r);
    e(i) = 1;
    u1 = A(:,k) - (alpha * e');
    v1 = u1 / norm(u1);
    H = eye(r) - ( 1 + conj(v1' * A(:,k))/(v1' * A(:,k))) * (v1 * v1');
    %%%%%%%%%%%%%%%%%%%%%%%
    
%     v = zeros(1,r);
%     v(k) = A(k,k) + sign(A(k,k)) * Sum(A(:,k), k);
%     for j = k + 1 : r
%         v(j) = A(j,k);
%     end
%     H = eye(r) - 2 * ( (v' * v) / (v * v'));
    
    Q = Q * H;
    A = H * A;
    k = k + 1;
    
end
R = A;
end

Подскажите, как же мне получить корректную матрицу Хаусхолдера на последующих шагах после первого?
Пример:
Исходная матрица

A =

 -10.0000 + 4.0000i   4.0000 - 5.0000i
   2.0000 - 1.0000i   8.0000 + 6.0000i
   6.0000 + 0.0000i   0.0000 + 7.0000i

После первого шага:

A =

  12.5300 - 0.0000i  -3.9904 + 7.6616i
   0.0000 - 0.0000i   8.3973 + 4.5919i
   0.0000 + 0.0000i   2.6433 + 4.0973i

После второго шага (ошибка):

A =

   5.2367 - 6.2171i  -0.0000 - 0.0000i
  -3.6274 - 6.9646i  13.7840 - 0.0000i
  -5.0866 + 1.8406i   0.0000 + 0.0000i

[WinterLord1]

Исправил функцию, теперь больше похоже на правду, НО возникают проблемы со знаками. Вот исправленная функция:

Используется синтаксис Matlab M

function [ Q, R ] = houseTriangular( A )
%   QR разложение матрицы A для поля вещественных
%   и комплексных чисел 
%   Входные параметры:  A - симметричная матрица
%   Выходные параметры: Q - унитарная матрица
%                       R - верхняя треугольная матрица

[r, ~] = size(A); 
Q = eye(r);
k = 0;

% для комплексной матрицы
if ~isreal(A)
    for i = 1 : r - 1
        alpha = norm(A(i:r,i));
        e = zeros(1, r);
        e(i) = 1;
        u = A(:,i) - (alpha * e');

        if k ~= 0 % для строк, начиная со второй
            u(1:k) = 0;
        end

        v = u / norm(u);
        H = eye(r) - ( 1 + conj(v' * A(:,i))/(v' * A(:,i))) * (v * v');
        A = H * A;
        Q = Q * H;
        k = k + 1;
    end
R = A;
end

Результат функции:

Цитата:

q =

-0.7981 - 0.3192i -0.0582 - 0.0176i 0.4551 - 0.2244i
0.1773 - 0.0200i 0.7877 - 0.4275i 0.2097 - 0.3478i
0.4496 + 0.1648i -0.4391 - 0.0173i 0.6770 - 0.3453i


r =

12.5300 - 0.0000i -3.9904 + 7.6616i
0.0000 - 0.0000i 10.7413 + 0.0000i
0.0000 - 0.0000i 0.0000 + 0.0000i

Встроенная функция матлаба:

Цитата:

q =

-0.7981 + 0.3192i -0.3036 - 0.2224i 0.3220 - 0.1260i
0.1596 - 0.0798i -0.7472 - 0.4151i -0.4190 + 0.2490i
0.4789 - 0.0000i -0.1779 - 0.3101i 0.8017 - 0.0085i


r =

12.5300 + 0.0000i -3.9904 + 7.6616i
0.0000 + 0.0000i -10.7413 + 0.0000i
0.0000 + 0.0000i 0.0000 + 0.0000i

Как видно, значения разнятся, хоть матрица Q и остается унитарной. Если умножить Q на R:

Цитата:

ans =

-10.0000 - 4.0000i 5.0059 - 5.0293i
2.2220 - 0.2504i 7.9066 - 3.1531i
5.6333 + 2.0654i -7.7732 + 2.6015i

Т.е. $A \ne QR$. Разложение не работает. Проблема осталась, есть какие идеи по этому поводу?

[WinterLord1]

Проблема решена. Собственно, что для комплексных, что для действительных, алгоритм работает одинаково. Моя ошибка была в начальном преобразовании столбца матрицы A в стоку v. Затем, при транспонировании, умножался комплексно-сопряженный столбец на стоку. А должно было быть наоборот.

Используется синтаксис Matlab M

function [ Q, R ] = houseTriangular( A )

[r, ~] = size(A);
Q = eye(r);
k = 1;
for i = 1 : r - 1
    v = zeros(1,r);
    v(k) = A(k,k) + sign(A(k,k)) * Sum(A(:,k), k);
    for j = k + 1 : r
        v(j) = A(j,k);
    end
    H = eye(r) - 2 * ( conj(v' * v) / (v * v')); % дополнительно выполняем комплексное сопряжение 
    Q = Q * H;
    A = H * A;
    k = k + 1;
end
R = A;
end