Givens矩阵
Givens矩阵有三个参数组成,
i
,
j
,
θ
i,j,\theta
i,j,θ.它长这样:
G
(
i
,
j
,
θ
)
=
(
I
0
0
0
0
0
cos
(
θ
)
0
sin
(
θ
)
0
0
0
I
0
0
0
−
sin
(
θ
)
0
cos
(
θ
)
0
0
0
0
0
I
)
G(i,j,\theta)=\begin{pmatrix} I & 0 & 0 & 0 & 0\\ 0 & \cos(\theta) & 0 & \sin(\theta) & 0 \\ 0 & 0 & I & 0 & 0\\ 0 & -\sin(\theta) & 0 & \cos(\theta) & 0\\ 0 & 0 & 0 & 0 & I \end{pmatrix}\\
G(i,j,θ)=
I00000cos(θ)0−sin(θ)000I000sin(θ)0cos(θ)00000I
如果角度不好计算也可以把
cos
θ
\cos \theta
cosθ和
sin
θ
\sin \theta
sinθ改成
c
,
s
c,s
c,s,于是Given矩阵就由四个参数定义:
G
(
i
,
j
,
c
,
s
)
=
(
I
0
0
0
0
0
c
0
s
0
0
0
I
0
0
0
−
s
0
c
0
0
0
0
0
I
)
G(i,j,c,s)=\begin{pmatrix} I & 0 & 0 & 0 & 0\\ 0 & c & 0 & s & 0 \\ 0 & 0 & I & 0 & 0\\ 0 & -s & 0 & c & 0\\ 0 & 0 & 0 & 0 & I \end{pmatrix}\\
G(i,j,c,s)=
I00000c0−s000I000s0c00000I
整个矩阵除了对角线和
a
i
i
,
a
i
j
,
a
j
i
,
a
j
j
a_{ii},a_{ij},a_{ji},a_{jj}
aii,aij,aji,ajj,其余位置都是0。对角线上除了第
i
i
i行和第
j
j
j行,其余位置都是1.
Givens矩阵和Householder差不多,都可以把向量变成第一个坐标为它的长度,其余坐标变成0.不同的是这样变,只需要一个Householder矩阵,而Givens变换需要多个Givens矩阵。
旋转到标准基
回忆利用Householder变换进行QR分解就是将向量投影到标准基上。Givens旋转也是一样,不过旋转比较麻烦。把任意向量x旋转到标准基
e
1
e_1
e1上,这个旋转变换是
n
n
n个Givens矩阵的乘积。乘积的每一项定义如下:
T
1
i
=
G
(
1
,
i
,
∑
j
=
1
i
−
1
x
j
2
∑
j
=
1
i
x
j
2
,
a
i
∑
j
=
1
i
x
j
2
)
T_{1i}=G(1,i,\frac{\sqrt{\sum_{j=1}^{i-1}x_j^2}}{\sqrt{\sum_{j=1}^{i}x_j^2}},\frac{a_i}{\sqrt{\sum_{j=1}^{i}x_j^2}})
T1i=G(1,i,∑j=1ixj2∑j=1i−1xj2,∑j=1ixj2ai)
把所有的这些Givens矩阵倒序乘起来得到一个矩阵:
T
=
∏
i
=
n
2
T
1
i
T=\prod_{i=n}^2T_{1i}
T=i=n∏2T1i
这个矩阵乘以向量会把向量旋转到
e
1
e_1
e1方向,也就是第一个坐标为向量的长度,其余坐标为0.说得这么多不如举个例子:
x
=
(
3
4
3
4
5
)
x=\begin{pmatrix}3\\ 4\\ 3\ 4\\ 5\\ \end{pmatrix}
x=
343 45
分解步骤:
T
12
=
(
0.6
0.8
0
0
0
−
0.8
0.6
0
0
0
0
0
1
0
0
0
0
0
1
0
0
0
0
0
1
)
T
13
=
(
0.857
0
0.514
0
0
0
1
0
0
0
−
0.514
0
0.857
0
0
0
0
0
1
0
0
0
0
0
1
)
T
14
=
(
0.825
0
0
0.566
0
0
1
0
0
0
0
0
1
0
0
−
0.566
0
0
0.825
0
0
0
0
0
1
)
T
15
=
(
0.816
0
0
0
0.577
0
1
0
0
0
0
0
1
0
0
0
0
0
1
0
−
0.577
0
0
0
0.816
)
T
=
T
15
T
14
T
13
T
12
=
(
0.346
0.462
0.346
0.462
0.577
−
0.8
0.6
0
0
0
−
0.309
−
0.412
0.857
0
0
−
0.291
−
0.388
−
0.291
0.825
0
−
0.245
−
0.327
−
0.245
−
0.327
0.816
)
T
x
=
(
8.66
0
0
0
0
)
T_{12}=\begin{pmatrix}0.6 & 0.8 & 0 & 0 & 0\\ -0.8 & 0.6 & 0 & 0 & 0\\ 0 & 0 & 1 & 0 & 0\\ 0 & 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 0 & 1\\ \end{pmatrix}\\ T_{13}=\begin{pmatrix}0.857 & 0 & 0.514 & 0 & 0\\ 0 & 1 & 0 & 0 & 0\\ -0.514 & 0 & 0.857 & 0 & 0\\ 0 & 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 0 & 1\\ \end{pmatrix}\\ T_{14}=\begin{pmatrix}0.825 & 0 & 0 & 0.566 & 0\\ 0 & 1 & 0 & 0 & 0\\ 0 & 0 & 1 & 0 & 0\\ -0.566 & 0 & 0 & 0.825 & 0\\ 0 & 0 & 0 & 0 & 1\\ \end{pmatrix}\\ T_{15}= \begin{pmatrix}0.816 & 0 & 0 & 0 & 0.577\\ 0 & 1 & 0 & 0 & 0\\ 0 & 0 & 1 & 0 & 0\\ 0 & 0 & 0 & 1 & 0\\ -0.577 & 0 & 0 & 0 & 0.816\\ \end{pmatrix}\\ T=T_{15}T_{14}T_{13}T_{12}=\begin{pmatrix}0.346 & 0.462 & 0.346 & 0.462 & 0.577\\ -0.8 & 0.6 & 0 & 0 & 0\\ -0.309 & -0.412 & 0.857 & 0 & 0\\ -0.291 & -0.388 & -0.291 & 0.825 & 0\\ -0.245 & -0.327 & -0.245 & -0.327 & 0.816\\ \end{pmatrix}\\ Tx=\begin{pmatrix}8.66\\ 0\\ 0\\ 0\\ 0\\ \end{pmatrix}
T12=
0.6−0.80000.80.6000001000001000001
T13=
0.8570−0.51400010000.51400.857000001000001
T14=
0.82500−0.566001000001000.566000.825000001
T15=
0.816000−0.5770100000100000100.5770000.816
T=T15T14T13T12=
0.346−0.8−0.309−0.291−0.2450.4620.6−0.412−0.388−0.3270.34600.857−0.291−0.2450.462000.825−0.3270.5770000.816
Tx=
8.660000
QR分解步骤
思想和Householder变换一样,按列循环,对于第一列找到一个
T
1
T_1
T1把矩阵的第一列变成上三角矩阵的第一列。对于第二列,就用同样的办法,生成Givens矩阵,连乘起来变成变换阵,把第二行以下的数字变成0,但是这个时候Givens矩阵的
c
,
s
c,s
c,s参数的计算方式要变了,变成这个样子:
T
k
i
=
G
(
k
,
i
,
∑
j
=
k
i
−
1
x
j
2
∑
j
=
k
i
x
j
2
,
a
i
∑
j
=
k
i
x
j
2
)
T_{ki}=G(k,i,\frac{\sqrt{\sum_{j=k}^{i-1}x_j^2}}{\sqrt{\sum_{j=k}^{i}x_j^2}},\frac{a_i}{\sqrt{\sum_{j=k}^{i}x_j^2}})
Tki=G(k,i,∑j=kixj2∑j=ki−1xj2,∑j=kixj2ai)
还是以这个向量为例子,通过一系列Givens矩阵可以将其第二行以下变成0:
T
23
=
(
1
0
0
0
0
0
0.8
0.6
0
0
0
−
0.6
0.8
0
0
0
0
0
1
0
0
0
0
0
1
)
T
24
=
(
1
0
0
0
0
0
0.781
0
0.625
0
0
0
1
0
0
0
−
0.625
0
0.781
0
0
0
0
0
1
)
T
25
=
(
1
0
0
0
0
0
0.788
0
0
0.615
0
0
1
0
0
0
0
0
1
0
0
−
0.615
0
0
0.788
)
T
2
=
T
25
T
24
T
23
=
(
1
0
0
0
0
0
0.492
0.369
0.492
0.615
0
−
0.6
0.8
0
0
0
−
0.5
−
0.375
0.781
0
0
−
0.384
−
0.288
−
0.384
0.788
)
T
2
x
=
(
3
8.124
0
0
0
)
T_{23}=\begin{pmatrix}1 & 0 & 0 & 0 & 0\\ 0 & 0.8 & 0.6 & 0 & 0\\ 0 & -0.6 & 0.8 & 0 & 0\\ 0 & 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 0 & 1\\ \end{pmatrix}\\ T_{24}=\begin{pmatrix}1 & 0 & 0 & 0 & 0\\ 0 & 0.781 & 0 & 0.625 & 0\\ 0 & 0 & 1 & 0 & 0\\ 0 & -0.625 & 0 & 0.781 & 0\\ 0 & 0 & 0 & 0 & 1\\ \end{pmatrix}\\ T_{25}= \begin{pmatrix}1 & 0 & 0 & 0 & 0\\ 0 & 0.788 & 0 & 0 & 0.615\\ 0 & 0 & 1 & 0 & 0\\ 0 & 0 & 0 & 1 & 0\\ 0 & -0.615 & 0 & 0 & 0.788\\ \end{pmatrix}\\ T_2=T_{25}T_{24}T_{23}=\begin{pmatrix}1 & 0 & 0 & 0 & 0\\ 0 & 0.492 & 0.369 & 0.492 & 0.615\\ 0 & -0.6 & 0.8 & 0 & 0\\ 0 & -0.5 & -0.375 & 0.781 & 0\\ 0 & -0.384 & -0.288 & -0.384 & 0.788\\ \end{pmatrix}\\ T_2x=\begin{pmatrix}3\\ 8.124\\ 0\\ 0\\ 0\\ \end{pmatrix}
T23=
1000000.8−0.60000.60.8000001000001
T24=
1000000.7810−0.62500010000.62500.781000001
T25=
1000000.78800−0.615001000001000.615000.788
T2=T25T24T23=
1000000.492−0.6−0.5−0.38400.3690.8−0.375−0.28800.49200.781−0.38400.615000.788
T2x=
38.124000
所以完整的QR分解就是这样按列迭代。但是要注意的是这一系列迭代是这样的过程:
T
n
−
1
⋯
T
1
A
=
R
A
=
T
1
−
1
T
2
−
1
⋯
T
n
−
1
−
1
R
=
Q
R
T_{n-1}\cdots T_1A=R\\ A=T_1^{-1}T_2^{-1}\cdots T_{n-1}^{-1}R=QR
Tn−1⋯T1A=RA=T1−1T2−1⋯Tn−1−1R=QR
Givens矩阵是正交阵,他们的乘积也是正交阵,所以有:
T
i
−
1
=
T
i
T
T_i^{-1}=T_i^T
Ti−1=TiT
所以有:
Q
=
T
1
T
T
2
T
⋯
T
n
−
1
T
Q=T_1^TT_2^T\cdots T_{n-1}^T
Q=T1TT2T⋯Tn−1T
Python代码
经过上述讲解,代码写起来就容易了,最终代码如下:
@staticmethod
def givens_matrix(n, i, j, c, s):
array = Matrix.unit_matrix(n)
array[i][i] = c
array[j][j] = c
array[i][j] = -s
array[j][i] = s
return Matrix(array)
@staticmethod
def get_givens_matrix(vector, line, i):
sum_vector = 0
for j in range(line, i):
sum_vector += vector[j] * vector[j]
numerator = math.sqrt(sum_vector)
vector_i_pow = vector[i] * vector[i]
denominator = math.sqrt(sum_vector + vector_i_pow)
alpha = numerator / denominator
beta = math.sqrt(vector_i_pow) / denominator
return Matrix.givens_matrix(len(vector), line, i, alpha, beta)
@staticmethod
def get_givens_rotator(vector, line):
rotator = Matrix(Matrix.unit_matrix(len(vector)))
n = len(vector)
for i in range(n - 1, line, -1):
givens_matrix = Matrix.get_givens_matrix(vector, line, i)
rotator = rotator * givens_matrix
return rotator
def givens_qr(self):
n = len(self.__vectors)
r = self
q = Matrix(Matrix.unit_matrix(n))
for i in range(n-1):
vector = r.__vectors[i]
givens = Matrix.get_givens_rotator(vector, i)
# 正交矩阵的逆矩阵就是它的转置矩阵
q = q * givens.transpose_matrix()
r = givens * r
return q, r