clear; clc;
% 1. 输入待分解的矩阵A
A = [1 2 3 4;
4 5 6 5;
7 8 9 6;
1 3 5 7];
fprintf('待进行QR分解的原始矩阵A:\n');
disp(A); % 显示原始矩阵,便于后续结果对比
[m,n] = size(A); % 获取矩阵A的行列数
Q = eye(n); % 初始化正交矩阵Q为n阶单位矩阵
% 循环构造Householder矩阵,逐步将A化为上三角矩阵
for j = 1 : n - 1
sigma = sign(A(j,j)) * norm(A(j:n,j), 2);
u = A(j:n,j) + sigma * eye(n - j + 1, 1);
beta = sigma * (sigma + A(j,j));
I = eye(n - j + 1);
Hj = I - (1 / beta) * (u * u');
Hk = blkdiag(eye(j - 1), Hj);
Q = Q * Hk;
A = Hk * A;
end
R = A; % 循环结束后,A已变为上三角矩阵R
fprintf('\n==================== QR分解结果 ====================\n');
fprintf('1. 正交矩阵Q:\n');
disp(Q); % 显示正交矩阵Q
fprintf('\n2. 上三角矩阵R:\n');
disp(R); % 显示上三角矩阵R
% 验证分解正确性(Q*R应接近原始矩阵,误差为浮点计算误差)
fprintf('\n==================== 结果验证 ====================\n');
fprintf('Q*R的计算结果(应与原始矩阵A近似):\n');
disp(Q * R);
fprintf('\n说明:Q*R与原始矩阵的微小差异为浮点计算误差,可忽略。\n');