5.9 QR分解–Gram-Schmidt 分解
最小二乘法需要解方程
A
T
A
x
=
A
T
b
A^TA\mathbf{x} = A^T\mathbf{b}
A T A x = A T b ,需要计算矩阵乘法
A
T
A
A^TA
A T A ,然后再高斯消元法解普通方程,缺点有两个,一是效率比较低,因为需要计算矩阵乘法
A
T
A
A^TA
A T A ,二是不稳定,当矩阵
A
A
A 列向量接近线性相关时,矩阵
A
T
A
A^TA
A T A 列向量会更接近线性相关,更病态,导致高斯消元法极不稳定,为此稳健最小二乘法和正则化方法都是克服这点的。还有一种方法,利用QR分解可以同时克服上面两个缺点,即可以高效率的稳定求出最小二乘解。
第一章指出,线性空间存在标准正交基,任意子空间都存在标准正交基的子集,该子集张成子空间。令标准正交基的子集为矩阵
Q
m
n
,
m
>
n
Q_{mn},m>n
Q m n , m > n ,注意不是方阵,要与正交矩阵
Q
m
m
Q_{mm}
Q m m 进行区分,正交矩阵是方阵。由于矩阵
Q
m
n
Q_{mn}
Q m n 每个列向量都是单位向量且两两正交,所以
Q
m
n
T
Q
m
n
=
E
n
n
Q^T_{mn}Q_{mn} = E_{nn}
Q m n T Q m n = E n n 是单位阵。注意
Q
m
n
Q
m
n
T
≠
E
m
m
Q_{mn}Q^T_{mn} \ne E_{mm}
Q m n Q m n T = E m m ,这也是和正交矩阵的区别,正交矩阵
Q
m
m
Q_{mm}
Q m m 满足:
Q
m
m
T
Q
m
m
=
E
,
Q
m
m
Q
m
m
T
=
E
Q^T_{mm}Q_{mm}=E,Q_{mm}Q^T_{mm}=E
Q m m T Q m m = E , Q m m Q m m T = E 。为了简化符号,本节
Q
Q
Q 不是正交矩阵,而是矩阵
Q
m
n
,
m
>
n
Q_{mn},m>n
Q m n , m > n 。
因为
Q
T
Q
=
E
Q^TQ = E
Q T Q = E ,所以矩阵
Q
Q
Q 的左逆为
Q
T
Q^T
Q T ,投影矩阵为
P
=
Q
Q
T
P=QQ^T
P = Q Q T ,都不涉及逆矩阵,只有矩阵转置,十分简单。矛盾方程
Q
x
=
b
Q\mathbf{x} = \mathbf{b}
Q x = b 的解可以直接得到
x
^
=
Q
T
b
\mathbf{\hat{x}} = Q^T\mathbf{b}
x ^ = Q T b ,这就是最小二乘解。残差向量为
b
−
Q
Q
T
b
=
(
E
−
Q
Q
T
)
b
\mathbf{b}-QQ^T\mathbf{b}=(E-QQ^T)\mathbf{b}
b − Q Q T b = ( E − Q Q T ) b ,矩阵
E
−
Q
Q
T
E-QQ^T
E − Q Q T 称为残差矩阵。
对于列满秩矩阵
A
A
A ,取张成子空间中标准正交向量,组成矩阵
Q
Q
Q ,注意有无穷多种取法。例如三维空间中任意二维子空间即平面,取平面内任意两个垂直的单位向量组成的矩阵即是
Q
Q
Q 。矩阵
A
A
A 任一列向量
a
i
\mathbf{a}_i
a i ,方程
Q
x
=
a
i
Q \mathbf{x}= \mathbf{a}_i
Q x = a i 都存在解
x
=
Q
T
a
i
\mathbf{x} = Q^T\mathbf{a}_i
x = Q T a i ,因为向量
a
i
\mathbf{a}_i
a i 位于子空间内,所以解是精确解,是严格相等的。则矩阵方程
Q
R
=
A
QR=A
Q R = A ,矩阵
R
R
R 每个列向量为对应方程的精确解,所以矩阵
R
R
R 唯一,注意其尺寸为
n
×
n
n \times n
n × n 是方阵。根据矩阵乘积的秩小于每个矩阵的秩,得
n
=
r
a
n
k
A
≤
r
a
n
k
R
n = rank A \le rank R
n = r a n k A ≤ r a n k R ,又因为
r
a
n
k
R
≤
n
rank R \le n
r a n k R ≤ n ,所以
r
a
n
k
R
=
n
rank R = n
r a n k R = n ,即方阵
R
R
R 满秩,是可逆矩阵。方程
A
x
=
b
A\mathbf{x} = \mathbf{b}
A x = b 带入
Q
R
=
A
QR=A
Q R = A 得
Q
R
x
=
b
QR\mathbf{x} = \mathbf{b}
Q R x = b ,所以
R
x
=
Q
T
b
R\mathbf{x} = Q^T\mathbf{b}
R x = Q T b 得
x
^
=
R
−
1
Q
T
b
\mathbf{\hat{x}} = R^{-1}Q^T\mathbf{b}
x ^ = R − 1 Q T b 为最小二乘解。由于需要求逆矩阵
R
−
1
R^{-1}
R − 1 ,效率比较低。由于矩阵
Q
Q
Q 有无穷多种取法,那有没有一种取法,使矩阵
R
R
R 为上三角阵,则逆矩阵
R
−
1
R^{-1}
R − 1 很容易求得?答案是肯定的。
定义 矩阵QR分解 列满秩矩阵
A
A
A 分解为
A
=
Q
R
A=QR
A = Q R ,其中矩阵
Q
Q
Q 每个列向量都是单位向量且两两正交,即满足
Q
T
Q
=
E
Q^TQ=E
Q T Q = E ,矩阵
R
R
R 为上三角阵,对角元素均为正数,可逆。此时左逆为
R
−
1
Q
T
R^{-1}Q^T
R − 1 Q T ,投影矩阵为
Q
Q
T
QQ^T
Q Q T 。
Gram-Schmidt 分解
下面求解矩阵QR分解的矩阵
Q
Q
Q 和矩阵
R
R
R 。因为
A
=
Q
R
A=QR
A = Q R ,则
[
a
1
,
⋯
,
a
n
]
=
[
q
1
,
⋯
,
q
n
]
[
r
1
,
⋯
,
r
n
]
[ \mathbf{a}_1,\cdots,\mathbf{a}_n] = [ \mathbf{q}_1,\cdots,\mathbf{q}_n] [ \mathbf{r}_1,\cdots,\mathbf{r}_n]
[ a 1 , ⋯ , a n ] = [ q 1 , ⋯ , q n ] [ r 1 , ⋯ , r n ] ,所以
a
1
=
[
q
1
,
⋯
,
q
n
]
r
1
\mathbf{a}_1 = [ \mathbf{q}_1,\cdots,\mathbf{q}_n] \mathbf{r}_1
a 1 = [ q 1 , ⋯ , q n ] r 1 ,因为
R
R
R 是上三角阵,所以
r
i
1
=
0
,
i
>
1
r_{i1} = 0, i>1
r i 1 = 0 , i > 1 ,则
a
1
=
q
1
r
11
\mathbf{a}_1 = \mathbf{q}_1 r_{11}
a 1 = q 1 r 1 1 ,因为
q
1
\mathbf{q}_1
q 1 是单位向量,所以对向量
a
1
\mathbf{a}_1
a 1 进行单位化可得
r
11
=
∥
a
1
∥
r_{11} = \|\mathbf{a}_1\|
r 1 1 = ∥ a 1 ∥ ,向量
q
1
=
a
1
/
r
11
\mathbf{q}_1 = \mathbf{a}_1/r_{11}
q 1 = a 1 / r 1 1 。
同理,
a
2
=
[
q
1
,
⋯
,
q
n
]
r
2
\mathbf{a}_2 = [ \mathbf{q}_1,\cdots,\mathbf{q}_n] \mathbf{r}_2
a 2 = [ q 1 , ⋯ , q n ] r 2 ,因为
R
R
R 是上三角阵,所以
r
i
2
=
0
,
i
>
2
r_{i2} = 0, i>2
r i 2 = 0 , i > 2 ,则
a
2
=
q
1
r
12
+
q
2
r
22
\mathbf{a}_2 = \mathbf{q}_1 r_{12} + \mathbf{q}_2 r_{22}
a 2 = q 1 r 1 2 + q 2 r 2 2 ,因为
q
i
\mathbf{q}_i
q i 是单位向量且两两正交,采用向量
q
1
\mathbf{q}_1
q 1 取内积得
q
1
T
a
2
=
q
1
T
q
1
r
12
+
q
1
T
q
2
r
22
\mathbf{q}^T_1\mathbf{a}_2 = \mathbf{q}^T_1\mathbf{q}_1 r_{12} +\mathbf{q}^T_1 \mathbf{q}_2 r_{22}
q 1 T a 2 = q 1 T q 1 r 1 2 + q 1 T q 2 r 2 2 ,因为
q
1
T
q
2
=
0
,
q
1
T
q
1
=
1
\mathbf{q}^T_1 \mathbf{q}_2 = 0, \mathbf{q}^T_1\mathbf{q}_1=1
q 1 T q 2 = 0 , q 1 T q 1 = 1 ,所以
r
12
=
q
1
T
a
2
r_{12} = \mathbf{q}^T_1\mathbf{a}_2
r 1 2 = q 1 T a 2 。所以
a
2
−
q
1
r
12
=
q
2
r
22
\mathbf{a}_2 - \mathbf{q}_1 r_{12} = \mathbf{q}_2 r_{22}
a 2 − q 1 r 1 2 = q 2 r 2 2 得
r
22
=
∥
a
2
−
q
1
r
12
∥
r_{22} = \|\mathbf{a}_2 - \mathbf{q}_1 r_{12} \|
r 2 2 = ∥ a 2 − q 1 r 1 2 ∥ ,
q
2
=
(
a
2
−
q
1
r
12
)
/
r
22
\mathbf{q}_2 = (\mathbf{a}_2 - \mathbf{q}_1 r_{12})/r_{22}
q 2 = ( a 2 − q 1 r 1 2 ) / r 2 2 。
同理对任意列向量
a
i
=
[
q
1
,
⋯
,
q
n
]
r
i
\mathbf{a}_i = [ \mathbf{q}_1,\cdots,\mathbf{q}_n] \mathbf{r}_i
a i = [ q 1 , ⋯ , q n ] r i ,因为
R
R
R 是上三角阵,所以
r
j
i
=
0
,
j
>
i
r_{ji} = 0, j>i
r j i = 0 , j > i ,则
a
i
=
q
1
r
1
i
+
q
2
r
2
i
+
⋯
+
q
i
r
i
i
\mathbf{a}_i = \mathbf{q}_1 r_{1i} + \mathbf{q}_2 r_{2i} + \cdots + \mathbf{q}_i r_{ii}
a i = q 1 r 1 i + q 2 r 2 i + ⋯ + q i r i i ,因为
q
i
\mathbf{q}_i
q i 是单位向量且两两正交,采用向量
q
j
,
j
<
i
\mathbf{q}_j,j < i
q j , j < i 取内积得
q
j
T
a
i
=
q
j
T
q
1
r
1
i
+
q
j
T
q
2
r
2
i
+
⋯
+
q
j
T
q
i
r
i
i
\mathbf{q}^T_j\mathbf{a}_i = \mathbf{q}^T_j\mathbf{q}_1 r_{1i} + \mathbf{q}^T_j\mathbf{q}_2 r_{2i} + \cdots + \mathbf{q}^T_j\mathbf{q}_i r_{ii}
q j T a i = q j T q 1 r 1 i + q j T q 2 r 2 i + ⋯ + q j T q i r i i ,因为
q
j
T
q
k
=
0
(
j
≠
k
)
,
q
j
T
q
j
=
1
\mathbf{q}^T_j \mathbf{q}_k = 0 (j \ne k), \mathbf{q}^T_j\mathbf{q}_j=1
q j T q k = 0 ( j = k ) , q j T q j = 1 ,所以
r
j
i
=
q
j
T
a
i
,
j
<
i
r_{ji} = \mathbf{q}^T_j\mathbf{a}_i, j < i
r j i = q j T a i , j < i 。所以
a
i
−
(
q
1
r
1
i
+
q
2
r
2
i
+
⋯
)
=
q
i
r
i
i
\mathbf{a}_i - (\mathbf{q}_1 r_{1i} + \mathbf{q}_2 r_{2i} + \cdots) = \mathbf{q}_i r_{ii}
a i − ( q 1 r 1 i + q 2 r 2 i + ⋯ ) = q i r i i 得
r
i
i
=
∥
a
i
−
(
q
1
r
1
i
+
q
2
r
2
i
+
⋯
)
∥
r_{ii} = \|\mathbf{a}_i - (\mathbf{q}_1 r_{1i} + \mathbf{q}_2 r_{2i} + \cdots) \|
r i i = ∥ a i − ( q 1 r 1 i + q 2 r 2 i + ⋯ ) ∥ ,
q
i
=
(
a
i
−
(
q
1
r
1
i
+
q
2
r
2
i
+
⋯
)
)
/
r
i
i
\mathbf{q}_i = (\mathbf{a}_i - (\mathbf{q}_1 r_{1i} + \mathbf{q}_2 r_{2i} + \cdots))/r_{ii}
q i = ( a i − ( q 1 r 1 i + q 2 r 2 i + ⋯ ) ) / r i i 。采用向量
q
i
\mathbf{q}_i
q i 取内积得
q
i
T
a
i
=
q
i
T
q
1
r
1
i
+
q
i
T
q
2
r
2
i
+
⋯
+
q
i
T
q
i
r
i
i
\mathbf{q}^T_i\mathbf{a}_i = \mathbf{q}^T_i\mathbf{q}_1 r_{1i} + \mathbf{q}^T_i\mathbf{q}_2 r_{2i} + \cdots + \mathbf{q}^T_i\mathbf{q}_i r_{ii}
q i T a i = q i T q 1 r 1 i + q i T q 2 r 2 i + ⋯ + q i T q i r i i ,可得
r
i
i
=
q
i
T
a
i
r_{ii} = \mathbf{q}^T_i\mathbf{a}_i
r i i = q i T a i ,这个计算方式形式上和
r
j
i
r_{ji}
r j i 统一。
一直计算,直到
i
=
n
i=n
i = n 结束。
该正交化方法就是Gram-Schmidt方法。该方法具有明显的几何图像,向量组
(
q
1
,
⋯
,
q
n
)
(\mathbf{q}_1,\cdots,\mathbf{q}_n)
( q 1 , ⋯ , q n ) 可看作空间的坐标轴,则
r
j
i
=
q
j
T
a
i
,
j
<
i
r_{ji} = \mathbf{q}^T_j\mathbf{a}_i, j < i
r j i = q j T a i , j < i 是列向量
a
i
\mathbf{a}_i
a i 在坐标轴
q
j
T
\mathbf{q}^T_j
q j T 的投影值,即坐标值。
q
i
r
i
i
=
a
i
−
(
q
1
r
1
i
+
q
2
r
2
i
+
⋯
)
\mathbf{q}_i r_{ii} = \mathbf{a}_i - (\mathbf{q}_1 r_{1i} + \mathbf{q}_2 r_{2i} + \cdots)
q i r i i = a i − ( q 1 r 1 i + q 2 r 2 i + ⋯ ) 是列向量
a
i
\mathbf{a}_i
a i 减去各个坐标轴投影分量,或者说
(
q
1
r
1
i
+
q
2
r
2
i
+
⋯
)
(\mathbf{q}_1 r_{1i} + \mathbf{q}_2 r_{2i} + \cdots)
( q 1 r 1 i + q 2 r 2 i + ⋯ ) 是列向量
a
i
\mathbf{a}_i
a i 向坐标轴
(
q
1
,
⋯
,
q
i
−
1
)
(\mathbf{q}_1,\cdots,\mathbf{q}_{i-1})
( q 1 , ⋯ , q i − 1 ) 张成子空间的投影分量,所以
q
i
r
i
i
\mathbf{q}_i r_{ii}
q i r i i 是列向量
a
i
\mathbf{a}_i
a i 垂直于由坐标轴
(
q
1
,
⋯
,
q
i
−
1
)
(\mathbf{q}_1,\cdots,\mathbf{q}_{i-1})
( q 1 , ⋯ , q i − 1 ) 张成子空间的分量,
r
i
i
r_{ii}
r i i 是该分量的长度,
q
i
\mathbf{q}_i
q i 是该分量的单位化向量。如果想象有个超长方体正好包围矩阵
A
A
A 的列向量,则
r
i
i
r_{ii}
r i i 就是长方体的各边长度,
q
i
\mathbf{q}_i
q i 是长方体的各边向量的方向。而且如果长方体有个边长很短,趋近
0
0
0 ,则矩阵
R
R
R 是病态的,从而矩阵
A
A
A 是病态的,方程
A
x
=
b
A\mathbf{x} = \mathbf{b}
A x = b 解不稳定。因为
x
^
=
R
−
1
Q
T
b
\mathbf{\hat{x}} = R^{-1}Q^T\mathbf{b}
x ^ = R − 1 Q T b ,假设
r
i
i
→
0
r_{ii} \to 0
r i i → 0 ,则
x
^
i
=
(
(
Q
T
b
)
i
−
∑
j
=
n
i
+
1
(
(
Q
T
b
)
j
x
^
j
)
)
/
r
i
i
\hat{x}_i = ((Q^T\mathbf{b})_i - \sum^{i+1}_{j=n} ((Q^T\mathbf{b})_j\hat{x}_j))/r_{ii}
x ^ i = ( ( Q T b ) i − ∑ j = n i + 1 ( ( Q T b ) j x ^ j ) ) / r i i ,除以趋于
0
0
0 的数,
x
^
i
\hat{x}_i
x ^ i 会很不稳定。
通过上面计算过程可以发现,矩阵的
Q
R
QR
Q R 分解是唯一的,因为计算每个
r
j
i
,
q
i
r_{ji},\mathbf{q}_i
r j i , q i 都是唯一的。也可以不通过计算过程证明唯一性。假设不唯一,则存在两个分解即
A
=
Q
R
=
Q
′
R
′
A=QR=Q'R'
A = Q R = Q ′ R ′ ,则
B
=
Q
′
T
Q
=
R
′
R
−
1
B=Q'^TQ=R'R^{-1}
B = Q ′ T Q = R ′ R − 1 ,
B
=
Q
′
T
Q
B=Q'^TQ
B = Q ′ T Q 是正交阵,则
B
T
=
B
−
1
B^T=B^{-1}
B T = B − 1 ,因为上三角阵乘积是上三角阵,所以
B
=
R
′
R
−
1
B=R'R^{-1}
B = R ′ R − 1 是上三角阵,上三角阵的逆矩阵也是上三角阵,所以
B
−
1
B^{-1}
B − 1 是上三角阵,则
B
T
B^T
B T 是上三角阵,所以
B
B
B 是下三角阵。这说明
B
B
B 即是上三角阵,又是下三角阵,则只能是对角阵,对角元素等于
r
i
i
/
r
i
i
′
r_{ii}/r'_{ii}
r i i / r i i ′ ,因为
r
i
i
,
r
i
i
′
r_{ii},r'_{ii}
r i i , r i i ′ 都是正数,所以对角元素都是正数。又
B
B
B 是正交矩阵,则只能是单位阵,故
E
=
Q
′
T
Q
=
R
′
R
−
1
E=Q'^TQ=R'R^{-1}
E = Q ′ T Q = R ′ R − 1 ,则
Q
′
=
Q
,
R
′
=
R
Q'=Q,R'=R
Q ′ = Q , R ′ = R ,故分解唯一。
Gram-Schmidt方法的优点是十分直观,几何图像很清晰,缺点是矩阵
Q
Q
Q 不是特别正交,因为数值舍入误差(计算机存储实数时采用有限位数,所以有舍入误差),会导致
q
i
\mathbf{q}_i
q i 偏离理论值,特别的
q
i
=
(
a
i
−
(
q
1
r
1
i
+
q
2
r
2
i
+
⋯
)
)
/
r
i
i
\mathbf{q}_i = (\mathbf{a}_i - (\mathbf{q}_1 r_{1i} + \mathbf{q}_2 r_{2i} + \cdots))/r_{ii}
q i = ( a i − ( q 1 r 1 i + q 2 r 2 i + ⋯ ) ) / r i i ,当
r
i
i
→
0
r_{ii} \to 0
r i i → 0 ,舍入误差导致
q
i
\mathbf{q}_i
q i 出现较大误差,导致
q
j
,
j
>
i
\mathbf{q}_j,j>i
q j , j > i 与
q
i
\mathbf{q}_i
q i 不正交,使矩阵
Q
T
Q
−
E
Q^TQ-E
Q T Q − E 偏离零矩阵较大。Gram-Schmidt分解由于舍入误差导致不稳定,从而导致最优解
x
^
=
R
−
1
Q
T
b
\mathbf{\hat{x}} = R^{-1}Q^T\mathbf{b}
x ^ = R − 1 Q T b 不稳定,特别当矩阵
A
A
A 接近病态时更严重。
当矩阵
A
A
A 是可逆矩阵时,可逆矩阵显然是列满秩矩阵,故存在QR分解,此时
Q
Q
Q 是正交矩阵,投影矩阵为
P
=
Q
Q
T
=
E
P=QQT=E
P = Q Q T = E 是单位矩阵!方程解为
Q
R
x
=
b
,
x
=
R
−
1
Q
T
b
QR\mathbf{x}=\mathbf{b},\mathbf{x}=R^{-1}Q^T\mathbf{b}
Q R x = b , x = R − 1 Q T b ,矩阵
A
A
A 的逆矩阵为
A
−
1
=
R
−
1
Q
T
A^{-1}=R^{-1}Q^T
A − 1 = R − 1 Q T 。
改进 Gram-Schmidt 分解
上面的Gram-Schmidt分解称为经典Gram-Schmidt分解,矩阵
R
R
R 是按列计算的,其实仔细观察计算公式
r
j
i
=
q
j
T
a
i
,
j
<
i
r_{ji} = \mathbf{q}^T_j\mathbf{a}_i, j < i
r j i = q j T a i , j < i ,因为
a
i
\mathbf{a}_i
a i 已知,只要知道了
q
j
T
\mathbf{q}^T_j
q j T ,是可以得到整行
r
j
i
,
i
=
1
,
⋯
,
n
r_{ji},i=1,\cdots,n
r j i , i = 1 , ⋯ , n ,可以按行计算。
改进Gram-Schmidt分解就是矩阵
R
R
R 是按行计算,即先计算
r
11
=
∥
a
1
∥
r_{11} = \|\mathbf{a}_1\|
r 1 1 = ∥ a 1 ∥ ,向量
q
1
=
a
1
/
r
11
\mathbf{q}_1 = \mathbf{a}_1/r_{11}
q 1 = a 1 / r 1 1 ,计算
r
1
i
=
q
1
T
a
i
,
i
=
2
,
⋯
,
n
r_{1i} = \mathbf{q}^T_1\mathbf{a}_i,i=2,\cdots,n
r 1 i = q 1 T a i , i = 2 , ⋯ , n 。
接着计算
r
22
=
∥
a
2
−
q
1
r
12
∥
r_{22} = \|\mathbf{a}_2 - \mathbf{q}_1 r_{12} \|
r 2 2 = ∥ a 2 − q 1 r 1 2 ∥ ,
q
2
=
(
a
2
−
q
1
r
12
)
/
r
22
\mathbf{q}_2 = (\mathbf{a}_2 - \mathbf{q}_1 r_{12})/r_{22}
q 2 = ( a 2 − q 1 r 1 2 ) / r 2 2 ,计算
r
2
i
=
q
2
T
a
i
,
i
=
3
,
⋯
,
n
r_{2i} = \mathbf{q}^T_2\mathbf{a}_i,i=3,\cdots,n
r 2 i = q 2 T a i , i = 3 , ⋯ , n 。
接着计算
r
i
i
=
∥
a
i
−
(
q
1
r
1
i
+
q
2
r
2
i
+
⋯
)
∥
r_{ii} = \|\mathbf{a}_i - (\mathbf{q}_1 r_{1i} + \mathbf{q}_2 r_{2i} + \cdots) \|
r i i = ∥ a i − ( q 1 r 1 i + q 2 r 2 i + ⋯ ) ∥ ,
q
i
=
(
a
i
−
(
q
1
r
1
i
+
q
2
r
2
i
+
⋯
)
)
/
r
i
i
\mathbf{q}_i = (\mathbf{a}_i - (\mathbf{q}_1 r_{1i} + \mathbf{q}_2 r_{2i} + \cdots))/r_{ii}
q i = ( a i − ( q 1 r 1 i + q 2 r 2 i + ⋯ ) ) / r i i ,计算
r
i
j
=
q
i
T
a
j
,
j
=
i
+
1
,
⋯
,
n
r_{ij} = \mathbf{q}^T_i\mathbf{a}_j,j=i+1,\cdots,n
r i j = q i T a j , j = i + 1 , ⋯ , n 。
一直计算,直到
i
=
n
i=n
i = n 结束。
如果这样计算,改进Gram-Schmidt分解和经典Gram-Schmidt分解结果一致,只是计算顺序不同而以,效果不会有任何改进。观察
q
i
=
(
a
i
−
(
q
1
r
1
i
+
q
2
r
2
i
+
⋯
)
)
/
r
i
i
\mathbf{q}_i = (\mathbf{a}_i - (\mathbf{q}_1 r_{1i} + \mathbf{q}_2 r_{2i} + \cdots))/r_{ii}
q i = ( a i − ( q 1 r 1 i + q 2 r 2 i + ⋯ ) ) / r i i ,需要减去
q
j
r
j
i
,
j
<
i
\mathbf{q}_j r_{ji},j < i
q j r j i , j < i ,这些量都是提前计算好的,当它们计算好后,立即减去这些量,就是改进Gram-Schmidt分解。具体为:
即先计算
r
11
=
∥
a
1
∥
r_{11} = \|\mathbf{a}_1\|
r 1 1 = ∥ a 1 ∥ ,向量
q
1
=
a
1
/
r
11
\mathbf{q}_1 = \mathbf{a}_1/r_{11}
q 1 = a 1 / r 1 1 ,计算
r
1
i
=
q
1
T
a
i
,
i
=
2
,
⋯
,
n
r_{1i} = \mathbf{q}^T_1\mathbf{a}_i,i=2,\cdots,n
r 1 i = q 1 T a i , i = 2 , ⋯ , n ,注意此时必须计算
a
i
=
a
i
−
q
1
r
1
i
,
i
=
2
,
⋯
,
n
\mathbf{a}_i = \mathbf{a}_i - \mathbf{q}_1 r_{1i}, i=2,\cdots,n
a i = a i − q 1 r 1 i , i = 2 , ⋯ , n 。令
b
0
=
b
\mathbf{b}_0=\mathbf{b}
b 0 = b ,注意对向量
b
\mathbf{b}
b 也要同样处理,即令
δ
1
=
q
1
T
b
\delta_1 = \mathbf{q}^T_1\mathbf{b}
δ 1 = q 1 T b ,
b
=
b
−
q
1
δ
1
\mathbf{b} = \mathbf{b} - \mathbf{q}_1\delta_1
b = b − q 1 δ 1 。
接着计算
r
22
=
∥
a
2
∥
r_{22} = \|\mathbf{a}_2\|
r 2 2 = ∥ a 2 ∥ ,
q
2
=
a
2
/
r
22
\mathbf{q}_2 = \mathbf{a}_2/r_{22}
q 2 = a 2 / r 2 2 ,计算
r
2
i
=
q
2
T
a
i
,
i
=
3
,
⋯
,
n
r_{2i} = \mathbf{q}^T_2\mathbf{a}_i,i=3,\cdots,n
r 2 i = q 2 T a i , i = 3 , ⋯ , n ,注意此时必须计算
a
i
=
a
i
−
q
2
r
2
i
,
i
=
3
,
⋯
,
n
\mathbf{a}_i = \mathbf{a}_i - \mathbf{q}_2 r_{2i}, i=3,\cdots,n
a i = a i − q 2 r 2 i , i = 3 , ⋯ , n 。注意对向量
b
\mathbf{b}
b 也要同样处理,即令
δ
2
=
q
2
T
b
\delta_2 = \mathbf{q}^T_2\mathbf{b}
δ 2 = q 2 T b ,
b
=
b
−
q
2
δ
2
\mathbf{b} = \mathbf{b} - \mathbf{q}_2 \delta_2
b = b − q 2 δ 2 。
接着计算
r
i
i
=
∥
a
i
∥
r_{ii} = \|\mathbf{a}_i\|
r i i = ∥ a i ∥ ,
q
i
=
a
i
/
r
i
i
\mathbf{q}_i = \mathbf{a}_i/r_{ii}
q i = a i / r i i ,计算
r
i
j
=
q
i
T
a
j
,
j
=
i
+
1
,
⋯
,
n
r_{ij} = \mathbf{q}^T_i\mathbf{a}_j,j=i+1,\cdots,n
r i j = q i T a j , j = i + 1 , ⋯ , n ,注意此时必须计算
a
j
=
a
j
−
q
i
r
i
j
,
j
=
i
+
1
,
⋯
,
n
\mathbf{a}_j = \mathbf{a}_j - \mathbf{q}_i r_{ij}, j=i+1,\cdots,n
a j = a j − q i r i j , j = i + 1 , ⋯ , n 。注意对向量
b
\mathbf{b}
b 也要同样处理,即令
δ
i
=
q
i
T
b
\delta_i = \mathbf{q}^T_i\mathbf{b}
δ i = q i T b ,
b
=
b
−
q
i
δ
i
\mathbf{b} = \mathbf{b} - \mathbf{q}_i\delta_i
b = b − q i δ i 。
一直计算,直到
i
=
n
i=n
i = n 结束。
令向量
d
=
(
δ
1
,
⋯
,
δ
m
)
\mathbf{d} = (\delta_1,\cdots,\delta_m)
d = ( δ 1 , ⋯ , δ m ) ,因为
d
=
Q
T
b
0
\mathbf{d} = Q^T\mathbf{b}_0
d = Q T b 0 ,则最优近似解为
x
^
=
R
−
1
d
\mathbf{\hat{x}} = R^{-1}\mathbf{d}
x ^ = R − 1 d ,即
x
^
i
=
(
δ
i
−
∑
j
=
i
+
1
n
(
δ
j
x
^
j
)
)
/
r
i
i
\hat{x}_i = (\delta_i - \sum^{n}_{j=i+1} (\delta_j\hat{x}_j))/r_{ii}
x ^ i = ( δ i − ∑ j = i + 1 n ( δ j x ^ j ) ) / r i i 。此时残差平方和为
∥
b
∥
2
\|\mathbf{b}\|^2
∥ b ∥ 2 ,因为
b
\mathbf{b}
b 为
b
0
\mathbf{b}_0
b 0 减去
A
A
A 子空间的投影。
改进Gram-Schmidt分解理论上和经典Gram-Schmidt分解结果一致,数学上是一样的。但是改进Gram-Schmidt分解不易受舍入误差影响,获得的正交矩阵
Q
Q
Q 更正交,使矩阵
Q
T
Q
−
E
Q^TQ-E
Q T Q − E 偏离零矩阵较小,最优解
x
^
=
R
−
1
Q
T
b
\mathbf{\hat{x}} = R^{-1}Q^T\mathbf{b}
x ^ = R − 1 Q T b 更稳定。为什么呢?注意到
a
j
=
a
j
−
q
i
r
i
j
,
j
=
i
+
1
,
⋯
,
n
\mathbf{a}_j = \mathbf{a}_j - \mathbf{q}_i r_{ij}, j=i+1,\cdots,n
a j = a j − q i r i j , j = i + 1 , ⋯ , n ,
q
i
r
i
j
\mathbf{q}_i r_{ij}
q i r i j 是向量
a
j
\mathbf{a}_j
a j 位于 $\mathbf{q}_i $ 投影,
a
j
\mathbf{a}_j
a j 减去投影分量后垂直于向量
q
i
\mathbf{q}_i
q i ,所以立即减去投影分量后,新的
a
j
\mathbf{a}_j
a j 都垂直于
q
i
\mathbf{q}_i
q i ,这样它们就是正交的,后面计算的
q
j
,
j
>
i
\mathbf{q}_j,j>i
q j , j > i 都与
q
i
\mathbf{q}_i
q i 正交,
q
i
\mathbf{q}_i
q i 的误差不会传导到后面。
Gram-Schmidt QR分解与高斯消元法对比
当矩阵
A
A
A 为方阵时,QR分解与高斯消元法都能解方程,它们的差别为:首先计算量上,高斯消元法更少,QR分解更多。其次,在数值稳定上,如果矩阵
A
A
A 接近病态时,QR分解得到的最优解要比高斯消元法的好很多。不严谨可以简单比较,因为根据QR分解
q
i
=
a
i
/
r
i
i
\mathbf{q}_i = \mathbf{a}_i/r_{ii}
q i = a i / r i i ,除数是向量
a
i
\mathbf{a}_i
a i 的长度
r
i
i
r_{ii}
r i i ,而高斯消元法的除数是主元,只是向量
a
i
\mathbf{a}_i
a i 的第
i
i
i 个分量
a
i
i
a_{ii}
a i i ,比
r
i
i
r_{ii}
r i i 更小,除数越小,计算越不稳定。最后,高斯消元法当
a
i
i
a_{ii}
a i i 为
0
0
0 时,需要进行行对调操作即选主元,即选择
a
j
i
,
j
=
i
,
⋯
,
m
a_{ji},j=i,\cdots,m
a j i , j = i , ⋯ , m 最大值为主元,然后对调行,可以改善解质量。QR分解也可以选择
∥
a
j
∥
,
j
>
i
\|\mathbf{a}_j\|,j>i
∥ a j ∥ , j > i 最大值为主元,然后对调列,但计算模拟表明解质量的改善很小,或者没有。
总之,如果矩阵
A
A
A 接近病态时,采用改进Gram-Schmidt分解法,否则可以采用高斯消元法(加上必要的选主元)。