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');