% Householder约化脚本（一般实矩阵→上海森伯格矩阵/对称矩阵→对称三对角矩阵）
% 功能：通过Householder变换实现矩阵约化，输出约化后矩阵H、正交矩阵UL和UR（满足UL=UR'，A≈UL*H*UR'）
% 运行方式：直接运行此m脚本，无需命令行输入或调用其他函数
clear; clc;  % 清空工作区变量与命令行窗口，避免历史数据干扰


%% 1. 输入待约化的矩阵A
% 示例1：非对称实矩阵（约化后为上海森伯格矩阵）
A = [1 2 3 4; 
     4 5 6 5; 
     7 8 9 6;
     1 3 5 7];  
% 示例2：对称实矩阵（约化后为对称三对角矩阵，可替换上方A）
% A = [17 3 4; 3 1 12; 4 12 8];

fprintf('待约化的原始矩阵A：\n');
disp(A);  % 显示原始矩阵，便于后续结果对比
[m,n] = size(A);  % 获取矩阵A的行列数


%% 2. Householder约化核心逻辑（整合原H_householder函数核心代码）
UR = eye(n);  % 初始化正交矩阵UR为n阶单位矩阵（累积Householder变换）

% 循环构造Householder矩阵，逐步完成约化
for j = 1 : n - 2 
    % 步骤1：计算当前列的关键参数sigma（保证数值稳定性，基于第j+1行元素符号）
    sigma = sign(A(j+1,j)) * norm(A(j+1:n,j), 2);
    % 步骤2：构造Householder向量u（目标是消去第j列下方非次对角线元素）
    u = A(j+1:n,j) + sigma * eye(n - j, 1);
    % 步骤3：计算归一化系数beta，避免除法运算误差
    beta = sigma * (sigma + A(j+1,j));
    % 步骤4：构造子空间Householder矩阵Rj（作用于第j+1行及以下子空间）
    I = eye(n - j);
    Rj = I - (1 / beta) * (u * u');
    % 步骤5：构造完整Householder矩阵Uk（将Rj嵌入单位矩阵对应位置，不改变前j行/列）
    Uk = blkdiag(eye(j), Rj);
    
    % 步骤6：更新矩阵A和正交矩阵UR
    A = Uk * A * Uk;    % 对A执行正交相似变换，逐步约化
    UR = UR * Uk;       % 累积Householder变换，UR最终为正交矩阵
end 

% 约化结果赋值（与原函数输出一致）
H = A;         % H为约化后的矩阵（上海森伯格/对称三对角矩阵）
UL = UR';      % UL为UR的转置（因UR是正交矩阵，UL=UR逆）


%% 3. 输出约化结果
fprintf('\n==================== Householder约化结果 ====================\n');
% 说明约化后矩阵类型
if isequal(A, A')  % 判断矩阵是否对称（对称则H为对称三对角矩阵）
    fprintf('1. 约化后的对称三对角矩阵H：\n');
else                % 非对称则H为上海森伯格矩阵
    fprintf('1. 约化后的上海森伯格矩阵H：\n');
end
disp(H);

fprintf('\n2. 正交矩阵UL（UR的转置）：\n');
disp(UL);
fprintf('\n3. 正交矩阵UR（累积Householder变换）：\n');
disp(UR);


%% 4. 结果验证（验证UL*H*UR' ≈ 原始矩阵A，确认变换正确性）
fprintf('\n==================== 结果验证 ====================\n');
recovered_A = UL * H * UR';  % 由约化结果反推原始矩阵
fprintf('UL*H*UR'' 的计算结果（应与原始矩阵A近似）：\n');
disp(recovered_A);
% 计算误差（ Frobenius范数，量化差异大小）
error = norm(recovered_A - A, 'fro');
fprintf('\n约化反推误差（Frobenius范数）：%.6e\n', error);
fprintf('说明：误差为浮点计算误差，数值极小时表明约化正确。\n');