householder方法化为Hessenberg矩阵(附Matlab、C代码)

1. Householder 变换¹

Householder 变换也称为初等反射变换，下面先定义 Householder 矩阵。用 Householder 矩阵左乘一个向量（或矩阵），即实现 Householder 变换。设向量 $\mathbf{w}\in {\mathbf{R}}^n$ 并且 ${\mathbf{w}}^{\mathbf{T}}\mathbf{w}=1$, 称矩阵 $\mathbf{H}(\mathbf{w})=\mathbf{I}-2\mathbf{w}{\mathbf{w}}^{\mathrm{T}}$​ 为 Householder 矩阵
$$\mathbf{H}(\mathbf{w})=\left[\begin{array}{cccc}1-2w_1^2& -2w_1{w}_2& \cdots & -2w_1{w}_n\\ -2w_2{w}_1& 1-2w_2^2& \cdots & -2w_2{w}_n\\ ⋮& ⋮& \ddots & ⋮\\ -2w_n{w}_1& -2w_n{w}_2& \cdots & 1-2w_n^2\end{array}\right],w=(w_1,w_2,\cdot \cdot \cdot ,w_n{)}^{\mathrm{T}}$$
定理：设 $\mathbf{x}={\left[x_1,x_2,\cdots ,x_n\right]}^{\mathrm{T}}\ne 0$ , 则存在 Householder 矩阵 $\mathbf{H}$ , 使 $\mathbf{Hx}=-\sigma {\mathbf{e}}_1$ , 其中,

$$\sigma =\mathrm{sign}\left(x_1\right)\parallel \mathbf{x}{\parallel }_2,{\mathbf{e}}_1={[1,0,\cdots ,0]}^{\mathrm{T}},\mathrm{sign}\left(x_1\right)=\{\begin{array}{ll}1,& x_1\ge 0\\ -1,& x_1<0\end{array}$$

构造满足定理要求的 Householder 矩阵时，可取向量
$$w=(x+\sigma e_1)/\parallel x+\sigma e_1{\parallel }_2$$

$$Hx=\left(\mathbf{I}-2\frac{\mathbf{v}{\mathbf{v}}^{\mathrm{T}}}{{\mathbf{v}}^{\mathrm{T}}\mathbf{v}}\right)\mathbf{x}=\mathbf{x}-2\frac{{\mathbf{v}}^{\mathrm{T}}\mathbf{x}}{{\mathbf{v}}^{\mathrm{T}}\mathbf{v}}\mathbf{v}$$

其中 $\mathbf{v}=(x+\sigma e_1)$​，只需计算向量 v 与 x 的内积，而不需要计算矩阵与向量的乘法

2. 将矩阵化简为上 Hessenberg 矩阵 ²

基本思路是: 先通过 Householder 变换将一般的矩阵 $\mathbf{A}$ 正交相似变换为上 Hessenberg 矩阵, 然后用 QR 算法求上 Hessenberg 矩阵的特征值, 便得到原矩阵的特征值。上 Hessenberg 矩阵形如
$$\left[\begin{array}{cccc}\times & \times & \times & \times \\ \times & \times & \times & \times \\ & \times & \times & \times \\ & & \times & \times \end{array}\right]$$
下面考虑用初等反射阵来正交相似约化一般矩阵和对称矩阵. 设
$$\mathbf{A}=\left(\begin{array}{cccc}{a}_{11}& a_{12}& \cdots & a_{1n}\\ a_{21}& a_{22}& \cdots & a_{2n}\\ ⋮& ⋮& & ⋮\\ a_{n1}& a_{n2}& \cdots & a_{nn}\end{array}\right)\equiv \left(\begin{array}{cc}{a}_{11}& {\mathbf{A}}_{12}^{(1)}\\ {\mathbf{a}}_{21}^{(1)}& {\mathbf{A}}_{22}^{(1)}\end{array}\right),$$
步 1 不妨设 ${\mathbf{a}}_{21}^{(1)}\ne 0$ , 否则这一步不需约化, 选择初等反射阵 ${\mathbf{R}}_1$ 使 ${\mathbf{R}}_1{\mathbf{a}}_{21}^{(1)}=$ $-{\sigma }_1{e}_1$ , 其中
$$\{\begin{array}{l}{\sigma }_1=\mathrm{sgn}\left(a_{21}\right){\left(\sum _{i=2}^n{a}_{i1}^2\right)}^{\frac{1}{2}},\\ {\mathbf{u}}_1={\mathbf{a}}_{21}^{(1)}+{\sigma }_1{\mathbf{e}}_1,\\ {\rho }_1=\frac{1}{2}{\parallel {\mathbf{u}}_1\parallel }_2^2={\sigma }_1\left({\sigma }_1+a_{21}\right),\\ {\mathbf{R}}_1=\mathbf{I}-{\rho }_1^{-1}{\mathbf{u}}_1{\mathbf{u}}_1^{\mathrm{T}}.\end{array}$$
令 ${\mathbf{U}}_1=\left(\begin{array}{cc}\mathbf{I}& 0\\ 0& R_1\end{array}\right)$ , 则
$${\mathbf{A}}_2={\mathbf{U}}_1{\mathbf{A}}_1{\mathbf{U}}_1=\left(\begin{array}{cc}{a}_{11}& {\mathbf{A}}_{21}^{(1)}{\mathbf{R}}_1\\ {\mathbf{R}}_1{\mathbf{a}}_{21}^{(1)}& {\mathbf{R}}_1{\mathbf{A}}_{22}^{(1)}{\mathbf{R}}_1\end{array}\right)\equiv \left(\begin{array}{ccc}{\mathbf{A}}_{11}^{(2)}& {\mathbf{a}}_{12}^{(2)}& {\mathbf{A}}_{13}^{(2)}\\ \mathbf{O}& {\mathbf{a}}_{22}^{(2)}& {\mathbf{A}}_{23}^{(2)}\end{array}\right),$$
其中
$${\mathbf{A}}_{11}^{(2)}\in {\mathbf{R}}^{2\times 1},{\mathbf{a}}_{22}^{(2)}\in {\mathbf{R}}^{n-2},{\mathbf{A}}_{23}^{(2)}\in {\mathbf{R}}^{(n-2)\times (n-2)}.$$
步 $\mathbf{k}$ 设对 $\mathbf{A}$ 已进行了第 $k-1$ 步正交相似约化, 即 ${\mathbf{A}}_k$ 有形式

$${\mathbf{A}}_k={\mathbf{U}}_{k-1}{\mathbf{A}}_{k-1}{\mathbf{U}}_{k-1}=\left[\begin{array}{ccccccccc}{a}_{11}& a_{12}^{(2)}& \cdots & a_{1k}^{(k)}& a_{1,k+1}^{(k)}& \cdots & a_{1n}^{(k)}& & \\ -{\sigma }_1& a_{22}^{(2)}& \cdots & a_{2k}^{(k)}& a_{2,k+1}^{(k)}& \cdots & a_{2n}^{(k)}& & \\ & \ddots & \ddots & \ddots & ⋮& ⋮& ⋮& & \\ & & \ddots & -{\sigma }_{k-1}& a_{kk}^{(k)}& a_{k,k+1}^{(k)}& \cdots & a_{k,n}^{(k)}& \\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & & \\ & & & & a_{k+1,k}^{(k)}& a_{k+1,k+1}^{(k)}& \cdots & a_{k+1,n}^{(k)}& \\ & & & & ⋮& ⋮& ⋮& & \\ & & & & & a_{nk}^{(k)}& a_{n,k+1}^{(k)}& \cdots & a_{nn}^{(k)}\end{array}\right]$$

$$\equiv \left(\begin{array}{ccc}{\mathbf{A}}_{11}^{(k)}& {\mathbf{a}}_{12}^{(k)}& {\mathbf{A}}_{13}^{(k)}\\ \mathbf{O}& {\mathbf{a}}_{22}^{(k)}& {\mathbf{A}}_{23}^{(k)}\end{array}\right),$$

其中
$${\mathbf{A}}_{11}^{(k)}\in {\mathbf{R}}^{k\times (k-1)},{\mathbf{a}}_{22}^{(k)}\in {\mathbf{R}}^{n-k},{\mathbf{A}}_{23}^{(k)}\in {\mathbf{R}}^{(n-k)\times (n-k)}.$$

设 ${\mathbf{a}}_{22}^{(k)}\ne 0$ , 选择初等反射阵 ${\mathbf{R}}_k$ , 使 ${\mathbf{R}}_k{\mathbf{a}}_{22}^{(k)}=-{\sigma }_k{\mathbf{e}}_1$ , 其中

$$\{\begin{array}{l}{\sigma }_k=\mathrm{sgn}\left(a_{k+1,k}^{(k)}\right){\left(\sum _{i=k+1}^n{a}_{ik}^2\right)}^{\frac{1}{2}},\\ {\mathbf{u}}_k={\mathbf{a}}_{22}^{(k)}+{\sigma }_k{\mathbf{e}}_1,\\ {\rho }_k=\frac{1}{2}{\parallel {\mathbf{u}}_k\parallel }_2^2={\sigma }_k\left({\sigma }_k+a_{k+1,n}^{(k)}\right),\\ {\mathbf{R}}_k=\mathbf{I}-{\rho }_k^{-1}{\mathbf{u}}_k{\mathbf{u}}_k^{\mathrm{T}}.\end{array}$$

设 ${\mathbf{U}}_k=\left(\begin{array}{cc}\mathbf{I}& \mathbf{O}\\ \mathbf{O}& {\mathbf{R}}_k\end{array}\right)$ , 则

$${\mathbf{A}}_{k+1}={\mathbf{U}}_k{\mathbf{A}}_k{\mathbf{U}}_k=\left(\begin{array}{ccc}{\mathbf{A}}_{11}^{(k)}& {\mathbf{a}}_{12}^{(k)}& {\mathbf{A}}_{13}^{(k)}{\mathbf{R}}_k\\ \mathbf{O}& {\mathbf{R}}_k{\mathbf{a}}_{22}^{(k)}& {\mathbf{R}}_k{\mathbf{A}}_{23}^{(k)}{\mathbf{R}}_k\end{array}\right)=\left(\begin{array}{ccc}{\mathbf{A}}_{11}^{(k)}& {\mathbf{a}}_{12}^{(k)}& {\mathbf{A}}_{13}^{(k)}{\mathbf{R}}_k\\ \mathbf{O}& -{\sigma }_k{\mathbf{e}}_1& {\mathbf{R}}_k{\mathbf{A}}_{23}^{(k)}{\mathbf{R}}_k\end{array}\right).$$

由上式知, ${\mathbf{A}}_{k+1}$ 的左上角 $k+1$ 阶子阵为上 Hessenberg 阵, 从而约化又进 了一步, 重复这过程, 直到

$${\mathbf{U}}_{n-2}\cdots {\mathbf{U}}_2{\mathbf{U}}_1\mathbf{A}{\mathbf{U}}_1{\mathbf{U}}_2\cdots {\mathbf{U}}_{n-2}=\left(\begin{array}{cccccc}{a}_{11}& \times & \times & \cdots & \times & \\ -{\sigma }_1& a_{22}^{(2)}& \times & \cdots & \times & \\ & -{\sigma }_2& a_{33}^{(3)}& \ddots & ⋮& \\ & & \ddots & \ddots & \times & \\ & & & \ddots & \times & \\ & & & & -{\sigma }_{n-1}& a_{nn}^{(n-1)}\end{array}\right)={\mathbf{A}}_{n-1}.$$

定理：如果 $A\in {\mathbf{R}}^{n\times n}$ 为对称阵,则存在初等反射阵 ${\mathbf{U}}_1,{\mathbf{U}}_2,\cdots ,{\mathbf{U}}_{n-2}$ ,使

$${\mathbf{U}}_{n-2}\cdots {\mathbf{U}}_2{\mathbf{U}}_1\mathbf{A}{\mathbf{U}}_1{\mathbf{U}}_2\cdots {\mathbf{U}}_{n-2}={\mathbf{A}}_{n-1}=\left(\begin{array}{ccccc}{c}_1& b_1& & & \\ b_1& c_2& b_2& & \\ & \ddots & \ddots & \ddots & \\ & & b_{n-2}& c_{n-1}& b_{n-1}\\ & & & b_{n-1}& c_n\end{array}\right)\equiv \mathbf{C}.$$

由上面讨论可知,在由 ${\mathbf{A}}_k\to {\mathbf{A}}_{k+1}={\mathbf{U}}_k{\mathbf{A}}_k{\mathbf{U}}_k$ 一步计算过程中,只需计算 ${\mathbf{R}}_k$ 和 ${\mathbf{R}}_k{\mathbf{A}}_{22}^{(k)}{\mathbf{R}}_k$ . 由于 $\mathbf{A}$ 的对称性,故只需计算 ${\mathbf{R}}_k{\mathbf{A}}_{23}^{(k)}{\mathbf{R}}_k$ 的对角线下面的元素. 注意到

$${\mathbf{R}}_k{\mathbf{A}}_{23}^{(k)}{\mathbf{R}}_k=\left(\mathbf{I}-{\rho }_k^{-1}{\mathbf{u}}_k{\mathbf{u}}_k^{\mathrm{T}}\right)\left({\mathbf{A}}_{23}^{(k)}-{\rho }_k^{-1}{\mathbf{A}}_{23}^{(k)}{\mathbf{u}}_k{\mathbf{u}}_k^{\mathrm{T}}\right),$$
引进记号
$${\mathbf{r}}_k={\rho }_k^{-1}{\mathbf{A}}_{23}^{(k)}{\mathbf{u}}_k,{\mathbf{t}}_k={\mathbf{r}}_k-\frac{{\rho }_k^{-1}}{2}\left({\mathbf{u}}_k^{\mathrm{T}}{\mathbf{r}}_k\right){\mathbf{u}}_k,$$

则 ${\mathbf{R}}_k{\mathbf{A}}_{23}^{(k)}{\mathbf{R}}_k={\mathbf{A}}_{23}^{(k)}-{\mathbf{u}}_k{\mathbf{t}}_k^{\mathrm{T}}-{\mathbf{t}}_k{\mathbf{u}}_k^{\mathrm{T}}(i=k+1,\cdots ,n;j=k+1,\cdots ,i).$

一般矩阵的变换Matlab代码³如下

function [H,P]=hessenb(A)
%使用Householder变换将矩阵A变为Hessenberg矩阵
%[H,P] H为变换后的Hessenberg矩阵，P为变换矩阵
[n,n]=size(A);
E=eye(n);
P1=E;
for k=1:n-2
    s=-sign(A(k+1,k))*norm(A(k+1:n,k));
    W(1:k)=zeros(1,k);
    W(k+1)=(A(k+1,k)+s);
    W(k+2:n)=A(k+2:n,k)';
    if norm(W)~=0
        W=W/norm(W);
    end
    P=E-2*W'*W;
    A=P*A*P;
    P1=P1*P;
end
H=A;
P=P1;

针对实对称矩阵的Matlab代码如下

function [Q,A]=rsmhessenb(A)
%实对称矩阵化为三对角矩阵
n=length(A(:,1));
Q=eye(n);
for i=2:n
    sigma=sign(A(i,i-1))*sqrt(sum(A(i:n,i-1).^2));
    rou=sigma*(A(i,i-1)+sigma);
    A(i,i-1)=A(i,i-1)+sigma;
    R=eye(n-i+1)-1.0/rou*A(i:n,i-1)*A(i:n,i-1)';
    A(i,i-1)=-sigma;
    A(i+1:n,i-1)=0;
    A(i-1,i:n)=A(i:n,i-1);%A^T=A
    A(i:n,i:n)=R*A(i:n,i:n)*R;
    m=length(R(:,1));
    U=[eye(n-m,n-m),zeros(n-m,m);zeros(m,n-m),R];
    Q=Q*U;
end

将三对角化的数学过程转为C语言代码实现

void rsmhessenb(double A[], int n, double Q[],double u[],double t[]) {
    //输入实对称矩阵A->三对角矩阵，输出变换矩阵Q，主对角线u，次对角线t
    int i, j, k;
    double sigma, rou,d;//记录求和值    
    // 初始化Q为单位矩阵
    for (i = 0; i < n; i++) {
        for (j = 0; j < n; j++) {
            Q[i * n + j] = (i == j) ? 1.0 : 0.0;
        }
    }
    for (i = 1; i < n; i++){
        //初始化向量t、u
        for(j=0;j<n;j++){
            t[j]=0.0;
            u[j]=0.0;
        }
        d=0.0;//记录矩阵内积值，每次循环清零
        // 计算sigma
        sigma = 0.0;
        for (j = i; j < n; j++) {
            sigma += A[j * n + i-1] * A[j * n + i-1];
        }
        sigma =(A[i*n+i-1] >= 0) ? sqrt(sigma): -sqrt(sigma);
        // 计算ρ
        rou = sigma * (A[i * n + i - 1] + sigma);
        //构建向量u
        u[i]=sigma+A[i * n + i - 1];
        for(j=i+1;j<n;j++){
            u[j] = A[j*n+i-1];
        }
        //t=A*u/ρ，作为中间过程的r
        for(j=i;j<n;j++){
            for(k=i;k<n;k++){
               t[j]+=A[j*n+k]*u[k]/rou;
        }
        }
        //d=(u^T)*t/(2ρ)
        for(j=i;j<n;j++){
            d+=u[j]*t[j]/2/rou;
        }
        //计算真实的t
        for(j=i;j<n;j++){
            t[j]=t[j]-u[j]*d;
        }
        //计算R*A*R=A-u*(t^T)-t*(u^T)，更新下方矩阵
        for(j=i;j<n;j++){
            for(k=j;k<n;k++){
                A[j*n+k]=A[j*n+k]-u[j]*t[k]-u[k]*t[j];
                if(k!=j){
                    A[k*n+j]=A[j*n+k];//更新对角元素 
                }
                }
            if(j+1<n){
                A[(j+1)*n+i-1]=0.0;//更新0元素
                A[(i-1)*n+j+1]=0.0;
            }
        }
        A[i*n+i-1]=-sigma;//更新次对角线
        A[(i-1)*n+i]=-sigma;
        sigma=0.0;//清零用于变换矩阵的计算
        for(j=0;j<n;j++){
            sigma+=u[j]*u[j];
        }
        //清零用于存储Q*u
        for(j=0;j<n;j++){
            t[j]=0.0;
        }
        //计算Q*u
        for(j=0;j<n;j++){
            for(k=0;k<n;k++){
                t[j]+=Q[j*n+k]*u[k];
            }
        }
        //计算Q*变换矩阵=Q-2*u*(u^T)/sigma        
        for(j=0;j<n;j++){
            for(k=0;k<n;k++){
                Q[j*n+k]=Q[j*n+k]-t[j]*u[k]*2/sigma;
            }
        }
    }
    //返回对角元素
    for (i = 0; i < n; i++) {
        u[i]=A[i*n+i];
        if((i+1)<n){
            t[i]=A[(i+1)*n+i];
        }
    }
}

3. 结果验证

​ 一般矩阵的Hessenberg化验证

图1. 一般矩阵变换结果

​ 实对称矩阵的hessenberg化验证

图2. 实对称测试矩阵

图3. 实对称矩阵变换结果

​ C语言实对称hessenberg化验证，测试矩阵与Matlab保持一致

图4. C代码运行结果

​ 由上图可知，C代码输出结果与Matlab结果相同。

4. 参考文献

---

1. 喻文健主编. 数值分析与算法 第3版[M]. 北京：清华大学出版社, 2020.03 ↩︎

2. 李庆扬，王能超，易大义. 数值分析 第5版[M]. 武汉：华中科技大学出版社, 2018.04. ↩︎

3. 夏省祥，于正文著. 常用数值算法及其MATLAB实现[M]. 北京：清华大学出版社, 2014.04. ↩︎