Householder 变换
反射矩阵
对任意模为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
简单证明:
- H T = ( I − 2 u u T ) T = I − 2 ( u T ) T u T = I − 2 u u T = H H^T=(I-2uu^T)^T=I-2(u^T)^Tu^T=I-2uu^T=H HT=(I−2uuT)T=I−2(uT)TuT=I−2uuT=H
- H ∗ H = ( I − 2 u u T ) ∗ ( I − 2 u u T ) = I − 4 u u T + 4 u ( u T u ) u T = I H*H=(I-2uu^T)*(I-2uu^T)=I-4uu^T+4u(u^Tu)u^T=I H∗H=(I−2uuT)∗(I−2uuT)=I−4uuT+4u(uTu)uT=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 , y ∈ R x ,y\in R x,y∈R,则存在反射变幻 H H H,使 H x = y Hx=y Hx=y
对
∀
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开始:
计算 ρ = − s g n ( a 11 ) ∑ i = 1 m a i 1 \rho=-sgn(a_{11})\overset{m}{\underset{i=1}{\sum}}a_{i1} ρ=−sgn(a11)i=1∑mai1,其中 s g n sgn sgn为取符号操作, − s g n ( a 11 ) -sgn(a_{11}) −sgn(a11)是为了保证 ∣ u 1 ∣ = ∣ a 11 − ρ ∣ |u_1|=|a_{11}-\rho| ∣u1∣=∣a11−ρ∣尽可能地大,避免数值精度损失,于是变换后的 A ^ 1 = [ ρ , 0 , … , 0 ] T \hat A_1=[\rho,0,\dots,0]^T A^1=[ρ,0,…,0]T。
计算 u 1 u_1 u1( x − y x-y x−y),由式(5)可知,这里 x x x指的是变换前的 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, y y y指的是变换后的 A ^ 1 = [ ρ , 0 , … , 0 ] T \hat A_1=[\rho,0,\dots,0]^T A^1=[ρ,0,…,0]T。显然, u 1 ( 1 ) = a 11 − ρ , u ( n ) = A 1 ( n ) , n > 1 u_1(1)=a_{11}-\rho,u(n)=A_1(n),n>1 u1(1)=a11−ρ,u(n)=A1(n),n>1
计算缩放因子 β 1 \beta_1 β1, β 1 = ∣ ∣ x − y ∣ ∣ = ∣ ∣ u ∣ ∣ = − 1 ρ u 1 ( 1 ) \beta_1=||x-y||=||u||=\frac{-1}{\rho u_1(1)} β1=∣∣x−y∣∣=∣∣u∣∣=ρu1(1)−1
于是:
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
在实际工程计算中,后续继续增加待定参数时,无需重复计算 u u u和 β \beta β,提高计算效率。
具体代码实现如下:
#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,β,第二种情况则为存储,根据具体需求进行选择。
注意:
- s的计算策略不同(
s1=sqrt(s1)
或者s1 = -sqrt(s1) * sign(1, src[n][n])
),最后展示结果略有差异,后者稳定性更佳,通过矩阵的基本原理可知,这几种结果均是等价的- 实际工程计算中,s为0时将会影响后续对矩阵进行求逆,需要谨慎