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