Household 变换原理与代码设计

Household 变换

文章目录

Household 变换

反射矩阵

对任意模为1的向量 u ∈ R u \in R u∈R,有矩阵 H = I − 2 u u T H=I-2uu^T H=I−2uuT,并满足：

H = H T H=H^T H=HT
H ∗ H = I H*H=I H∗H=I

记超平面 S S S（过原点以 u u u为法向）：
S = { x ∣ u T x = 0 , ∀ x ∈ R n } S=\{x|u^Tx=0,\forall x \in R^n\} S={x∣uTx=0,∀x∈Rn}
任一向量 z ∈ R n z \in R^n z∈Rn可以在子空间 S S S和 s p a n { u } span\{u\} span{u}做正交分解：
z = x + y , x ∈ S , y ∈ α ⋅ u , α ∈ R z=x+y,x \in S, y \in \alpha \cdot u,\alpha \in R z=x+y,x∈S,y∈α⋅u,α∈R
将矩阵 H H H作用到 z z z上，于是：
H z = H x + H y = ( I − 2 u u T ) x + ( I − 2 u u T ) y = x + y − 2 u u T y = x + y − 2 α u ( u T u ) = x − y \begin{aligned} Hz&=Hx+Hy \\ &=(I-2uu^T)x+(I-2uu^T)y \\ &=x+y-2uu^Ty \\&=x+y-2\alpha u(u^T u) \\&=x - y \end{aligned} Hz=Hx+Hy=(I−2uuT)x+(I−2uuT)y=x+y−2uuTy=x+y−2αu(uTu)=x−y
通过矩阵 H H H的变换，将 x + y x+y x+y以 S S S为镜面反射到 x − y x-y x−y ，因此矩阵 H H H叫做反射矩阵，而这一变换过程叫做Householder变换。

上三角化处理

对 ∀ x ∈ R n , ∣ ∣ x ∣ ∣ = ρ \forall x \in R^n, ||x||=\rho ∀x∈Rn,∣∣x∣∣=ρ，可以通过Householder变幻H得到：
H x = [ ± ρ , 0 , . . . , 0 ] T Hx=[\pm \rho,0,...,0]^T Hx=[±ρ,0,...,0]T
构建反射矩阵H时，需满足 u u u的模为1，因此取：
u = x − y ∣ ∣ x − y ∣ ∣ u=\frac{x-y}{||x-y||} u=∣∣x−y∣∣x−y
对于给定方阵，可通过反射变幻逐列进行上三角化操作：
[ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ] ⟶ H 1 = I − 2 u 1 u 1 T [ ∗ ∗ ∗ ∗ 0 ∗ ∗ ∗ 0 ∗ ∗ ∗ 0 ∗ ∗ ∗ ] ⟶ H 2 = I − 2 u 2 u 2 T [ ∗ ∗ ∗ ∗ 0 ∗ ∗ ∗ 0 0 ∗ ∗ 0 0 ∗ ∗ ] ⟶ H 3 = I − 2 u 3 u 3 T [ ∗ ∗ ∗ ∗ 0 ∗ ∗ ∗ 0 0 ∗ ∗ 0 0 0 ∗ ] \begin{bmatrix} *&*&*&*\\ *&*&*&*\\ *&*&*&*\\ *&*&*&*\\ \end{bmatrix} \overset{H_1=I-2u_1u_1^T}{\longrightarrow} \begin{bmatrix} *&*&*&*\\ 0&*&*&*\\ 0&*&*&*\\ 0&*&*&*\\ \end{bmatrix} \overset{H_2=I-2u_2u_2^T}{\longrightarrow} \begin{bmatrix} *&*&*&*\\ 0&*&*&*\\ 0&0&*&*\\ 0&0&*&*\\ \end{bmatrix} \overset{H_3=I-2u_3u_3^T}{\longrightarrow} \begin{bmatrix} *&*&*&*\\ 0&*&*&*\\ 0&0&*&*\\ 0&0&0&*\\ \end{bmatrix} ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗ ⟶H1=I−2u1u1T ∗000∗∗∗∗∗∗∗∗∗∗∗∗ ⟶H2=I−2u2u2T ∗000∗∗00∗∗∗∗∗∗∗∗ ⟶H3=I−2u3u3T ∗000∗∗00∗∗∗0∗∗∗∗
通过一系列 H H H作用到 A A A，可将其逐渐化为上三角矩阵，于是上述过程可描述为：
H n − 1 H 2 H 1 ⋅ A = R , P A = R , A = P T R = Q R H_{n-1}H_2H_1\cdot A=R, PA=R,A=P^TR=QR Hn−1H2H1⋅A=R,PA=R,A=PTR=QR
注：这里的 H n H_n Hn实际上为 d i a g ( I , H n ) diag(I,H_n) diag(I,Hn),仅为表达公式的合理性。在实际代码实现过程中，我们并不会将其补全再计算。

逐列上三角化操作流程

例如，对于 m × n m\times n m×n维矩阵 A m × n ( m ≥ n ) A_{m\times n}(m\ge n) Am×n(m≥n)进行上三角化操作，从第一列 A 1 = [ a 11 a 21 ⋯ a m 1 ] T A_1=[a_{11}\space a_{21}\cdots \space a_{m1} ]^T A1=[a11 a21⋯ am1]T开始：

于是：
T u 1 A 1 = [ ρ , 0 , … , 0 ] T T_{u_1}A_1=[\rho,0,\dots,0]^T Tu1A1=[ρ,0,…,0]T
上述公式即完成了对第一列元素的更新，主对角线元素以下皆为0。对于第 2 , 3 , … , n 2,3,\dots,n 2,3,…,n列元素的更新，可按照下式进行：
T u j A j = A j − β ( A j T u ) u T_{u_j}A_j=A_j-\beta(A^T_ju)u TujAj=Aj−β(AjTu)u
当n列元素均更新完毕后，至此， H 1 ⋅ A = R 1 H_1\cdot A=R_1 H1⋅A=R1已计算完毕，已完成等式（6）中第一步至第二步的转化。仿照上述步骤，逐列对矩阵 A m × n A_{m\times n} Am×n进行上三角化操作。当矩阵加入新的参数时，为了避免重复计算，需要记录 u 1 , u 2 , … , u n u_1,u2,\dots,u_n u1,u2,…,un和 β 1 , β 2 , … , β n \beta_1,\beta_2,\dots,\beta_n β1,β2,…,βn。

C语言代码实现

仿照上述公式实现下C语言代码，在工程应用中，为了节省空间，可以将需要记录的 u 1 , u 2 , … , u n u_1,u2,\dots,u_n u1,u2,…,un和 β 1 , β 2 , … , β n \beta_1,\beta_2,\dots,\beta_n β1,β2,…,βn记录在上三角化之后的矩阵中。若是如此，则需要确保矩阵空间开辟时，行数要增加2。

例如：对矩阵 A 4 × 4 A_{4\times 4} A4×4进行上三角化操作（实际开辟维度为 6 × 4 6\times4 6×4），当第一列完成上三角化后：
[ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ] ⟶ H 1 = I − 2 u 1 u 1 T [ ∗ ∗ ∗ ∗ 0 ∗ ∗ ∗ 0 ∗ ∗ ∗ 0 ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ] ⟶ A i 1 = u i , ( i > 1 ) [ ∗ ∗ ∗ ∗ u 1 ∗ ∗ ∗ u 2 ∗ ∗ ∗ u 3 ∗ ∗ ∗ u 4 ∗ ∗ ∗ β 1 ∗ ∗ ∗ ] \begin{bmatrix} *&*&*&*\\ *&*&*&*\\ *&*&*&*\\ *&*&*&*\\ *&*&*&*\\ *&*&*&*\\ \end{bmatrix} \overset{H_1=I-2u_1u_1^T}{\longrightarrow} \begin{bmatrix} *&*&*&*\\ 0&*&*&*\\ 0&*&*&*\\ 0&*&*&*\\ *&*&*&*\\ *&*&*&*\\ \end{bmatrix} \overset{A_{i1}=u_i,(i>1)}{\longrightarrow} \begin{bmatrix} *&*&*&*\\ u_{1}&*&*&*\\ u_{2}&*&*&*\\ u_{3}&*&*&*\\ u_{4}&*&*&*\\ \beta_1&*&*&*\\ \end{bmatrix} ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗ ⟶H1=I−2u1u1T ∗000∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗ ⟶Ai1=ui,(i>1) ∗u1u2u3u4β1∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
当所有列均上三角完成后，矩阵A为：
[ ∗ ∗ ∗ ∗ u 11 ∗ ∗ ∗ u 12 u 21 ∗ ∗ u 13 u 22 u 31 ∗ u 14 u 23 u 32 u 41 β 1 β 2 β 3 β 4 ] \begin{bmatrix} *&*&*&*\\ u_{11}&*&*&*\\ u_{12}&u_{21}&*&*\\ u_{13}&u_{22}&u_{31}&*\\ u_{14}&u_{23}&u_{32}&u_{41}\\ \beta_1&\beta_2&\beta_3&\beta_4\\ \end{bmatrix} ∗u11u12u13u14β1∗∗u21u22u23β2∗∗∗u31u32β3∗∗∗∗u41β4

具体代码实现如下：

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include <stdbool.h>

#define MAX 1024
#define ROW 4
#define COL 4

// 返回x的绝对值，y的符号
double sign(double x, double y)
{
    return fabs(x) * (y >= 0 ? 1 : -1);
}
// 获取每列元素，u,s,beta
void get_ele(double **src, int n, int row, double u[], double *s, double *beta)
{
    int i, j;
    double s1;
    *beta = *s = s1 = 0.0;
    for (i = n; i < row; i++)
    {
        s1 += src[i][n] * src[i][n];
    }
    s1 = -sqrt(s1) * sign(1, src[n][n]);
    if (s1 == 0.0)
    {
        printf("*****WARRING(get_ele):s equal 0.0, column: %d\n", n + 1);
        return;
    }
    u[0] = src[n][n] - s1;
    if (u[0] == 0.0)
    {
        printf("*****WARRING(get_ele):u[0] equal 0.0, column: %d\n", n + 1);
        return;
    }
    for (i = 1; i < row - n; i++)
    {
        u[i] = src[n + i][n];
    }
    *beta = 1.0 / s1 / u[0];
    *s = s1;
}
// householder变换
void householder(double **src, int row, int col, bool save)
{
    int i, j, k, flag = 1, tmp;
    double s, beta, gama;
    double u[MAX] = {0};

    if (src == NULL)
    {
        printf("*****ERROR(householder):src Matrix is NULL, exit!\n");
        exit(-1);
    }
    for (i = 0; i < col; i++)
    {
        // check elemental first
        flag = 1;
        tmp = 0;
        for (j = i + 1; j < row; j++)
        {
            if (src[j][i] == 0)
            {
                tmp++;
            }
        }
        if (tmp == row)
        {
            flag = 0;
            if (save)
            {
                src[j + 1][i] = 0;
                src[j = 2][i] = 0;
            }
        }
        if (flag)
        {
            get_ele(src, i, col, u, &s, &beta);
            if (s == 0.0)
                continue;
            src[i][i] = s;
            for (j = i + 1; j < col; j++)
            {
                gama = 0.0;
                for (k = i; k < row; k++)
                {
                    gama += u[k - i] * src[k][j];
                }
                if (gama == 0.0)
                    continue;
                gama *= beta;
                for (k = i; k < row; k++)
                {
                    src[k][j] += gama * u[k - i];
                }
            }
        }
        if (save)
        {
            for (j = i + 1; j < row + 1; j++)
            {
                src[j][i] = u[j - i - 1];
            }
            src[j][i] = beta;
        }
        else
        {
            for (j = i + 1; j < row; j++)
            {
                src[j][i] = 0.0;
            }
        }
    }
}

测试用例：

void test1(bool save)
{
    int i, j;
    double **arr = (double **)malloc(sizeof(double *) * (ROW + 2));
    for (i = 0; i < ROW + 2; i++)
    {
        arr[i] = (double *)malloc(sizeof(double) * COL);
        memset(arr[i], 0, sizeof(double) * COL);
    }
    arr[0][0] = 1;
    arr[1][0] = 1;
    arr[2][0] = 1;
    arr[3][0] = 1;

    arr[0][1] = 2;
    arr[1][1] = 0;
    arr[2][1] = 0;
    arr[3][1] = 2;

    arr[0][2] = 0;
    arr[1][2] = 3;
    arr[2][2] = 3;
    arr[3][2] = 0;

    arr[0][3] = 1;
    arr[1][3] = 1;
    arr[2][3] = 2;
    arr[3][3] = 2;

    householder(arr, ROW, COL, save);
	
    
    for (i = 0; i < ROW + 2; i++)
    {
        for (j = 0; j < COL; j++)
        {
            printf("%5.2lf  ", arr[i][j]);
        }
        putchar('\n');
    }
}

int main(int argc, char *argv[])
{
    test1(false);
    test1(true);
    return 0;
}

结果：

第一种情况为不存储 u , β u,\beta u,β，第二种情况则为存储，根据具体需求进行选择。

Household 变换原理与代码设计-LMLPHP

是王久久阿