矩阵的QR分解以及在最小二乘法中的应用

一、最小二乘法

最小二乘法是一种数学优化方法，通过最小化误差的平方和来拟合数据点。
以线性回归模型为例，如果我们用最小二乘法来求解线性回归的系数，可得：
err(yi−y^)=1n∑i=1n(yi−y^)2=1n∑i=1n(yi−wTxi)2=1n(y−wX)T(y−wX)=1n(yTy−2wXTy+wTXTXw)\begin{aligned} err(y_i-\hat y) &= \frac{1}{n}\sum_{i=1}^n (y_i-\hat y)^2 = \frac{1}{n}\sum_{i=1}^n (y_i-w^Tx_i)^2\\ &= \frac{1}{n}(y-wX)^T(y-wX) \\ &= \frac{1}{n}(y^Ty-2wX^Ty+w^TX^TXw) \end{aligned} err(yi−y^)=n1i=1∑n(yi−y^)2=n1i=1∑n(yi−wTxi)2=n1(y−wX)T(y−wX)=n1(yTy−2wXTy+wTXTXw)
我们要求上式的最小值，要对其求导，然后寻找极小值点。
∂∂werr=1n(2wXTX−2XTy)=0wXTX=XTyw=(XTX)−1XTy\begin{aligned} \frac{\partial}{\partial w}err &= \frac{1}{n}(2wX^TX-2X^Ty) = 0\\ &wX^TX = X^Ty\\ &w=(X^TX)^{-1}X^Ty \end{aligned} ∂w∂err=n1(2wXTX−2XTy)=0wXTX=XTyw=(XTX)−1XTy
由此我们便可以推导出参数的表达式。

二、QR分解

QR分解是把一个矩阵分解为一个正交矩阵和一个上三角矩阵的积。即有实数矩阵A，有A=Q×RA=Q\times RA=Q×R，其中Q为正交矩阵(QT⋅Q=IQ^T\cdot Q=IQT⋅Q=I)，R为上三角矩阵。QR分解常见的算法有Gram–Schmid正交化、Household变换，以及Givens变换。

2.1 Gran-Schmid正交化

设矩阵A=(a1⃗,a2⃗,...,an⃗)A=(\vec{a_1},\vec{a_2},...,\vec{a_n})A=(a1,a2,...,an)，对矩阵A进行Gran-Schmid正交化过程。其中pi⃗\vec{p_i}pi为正交向量，qi⃗\vec{q_i}qi为归一化后的标准正交向量，i=1,2,...,ni=1,2,...,ni=1,2,...,n 。
p1⃗=a1⃗=∥p1⃗∥q1⃗=r11q1⃗p2⃗=a2⃗−a2⃗⋅p1⃗∥p1⃗∥2⋅p1⃗=∥p2⃗∥q2⃗a2⃗=∥p2⃗∥q2⃗+a2⃗⋅p1⃗∥p1⃗∥2⋅p1⃗=∥p2⃗∥q2⃗+a2⃗⋅p1⃗∥p1⃗∥2∥p1⃗∥q1⃗=r21q1⃗+r22q2⃗p3⃗=a3⃗−a3⃗⋅p1⃗∥p1⃗∥2⋅p1⃗−a3⃗⋅p2⃗∥p2⃗∥2⋅p2⃗=∥p3⃗∥q3⃗a3⃗=∥p3⃗∥q3⃗+a3⃗⋅p1⃗∥p1⃗∥2⋅p1⃗+a3⃗⋅p2⃗∥p2⃗∥2⋅p2⃗=r31q1⃗+r32q2⃗+r33q3⃗\begin{aligned} \vec{p_1} &= \vec{a_1} = \lVert \vec{p_1}\rVert\vec{q_1}=r_{11}\vec{q_1}\\ \vec{p_2} &= \vec{a_2} - \frac{\vec{a_2}\cdot \vec{p_1}}{\lVert \vec{p_1}\rVert^2}\cdot \vec{p_1} = \lVert \vec{p_2}\rVert\vec{q_2}\\ \vec{a_2} &= \lVert \vec{p_2}\rVert\vec{q_2} + \frac{\vec{a_2}\cdot \vec{p_1}}{\lVert \vec{p_1}\rVert^2}\cdot \vec{p_1}\\ &= \lVert \vec{p_2}\rVert\vec{q_2} + \frac{\vec{a_2}\cdot \vec{p_1}}{\lVert \vec{p_1}\rVert^2}\lVert \vec{p_1}\rVert \vec{q_1}\\ &= r_{21}\vec{q_1} + r_{22}\vec{q_2}\\ \vec{p_3} &= \vec{a_3} - \frac{\vec{a_3}\cdot \vec{p_1}}{\lVert \vec{p_1}\rVert^2}\cdot \vec{p_1} - \frac{\vec{a_3}\cdot \vec{p_2}}{\lVert \vec{p_2}\rVert^2}\cdot \vec{p_2}\\ &= \lVert \vec{p_3}\rVert\vec{q_3}\\ \vec{a_3} &= \lVert \vec{p_3}\rVert\vec{q_3} + \frac{\vec{a_3}\cdot \vec{p_1}}{\lVert \vec{p_1}\rVert^2}\cdot \vec{p_1} + \frac{\vec{a_3}\cdot \vec{p_2}}{\lVert \vec{p_2}\rVert^2}\cdot \vec{p_2}\\ &= r_{31}\vec{q_1} + r_{32}\vec{q_2} + r_{33}\vec{q_3} \end{aligned} p1p2a2p3a3=a1=∥p1∥q1=r11q1=a2−∥p1∥2a2⋅p1⋅p1=∥p2∥q2=∥p2∥q2+∥p1∥2a2⋅p1⋅p1=∥p2∥q2+∥p1∥2a2⋅p1∥p1∥q1=r21q1+r22q2=a3−∥p1∥2a3⋅p1⋅p1−∥p2∥2a3⋅p2⋅p2=∥p3∥q3=∥p3∥q3+∥p1∥2a3⋅p1⋅p1+∥p2∥2a3⋅p2⋅p2=r31q1+r32q2+r33q3
之后通过分解矩阵AAA:
A=(a1⃗,a2⃗,...,an⃗)=(r11q1⃗,r21q1⃗+r22q2⃗,...,∑i=1nrniqi⃗)=(r11q1⃗,r21q1⃗,r31q1⃗,...,rn1q1⃗)+(0,r22q2⃗,r32q2⃗,...,rn2q2⃗)+......+(0,0,0,...,rnnqn⃗)=(q1⃗,q2⃗,...,qn⃗)⋅[r11r21⋯rn10r22⋯rn2⋮⋮⋱⋮00⋯rnn]=Q⋅R\begin{aligned} A &= (\vec{a_1},\vec{a_2},... ,\vec{a_n})\\ &= (r_{11}\vec{q_1},r_{21}\vec{q_1}+r_{22}\vec{q_2},... ,\sum_{i=1}^n r_{ni}\vec{q_i})\\ &= (r_{11}\vec{q_1},r_{21}\vec{q_1},r_{31}\vec{q_1},...,r_{n1}\vec{q_1})\\ &+ (0,r_{22}\vec{q_2},r_{32}\vec{q_2},...,r_{n2}\vec{q_2})\\ &+...... \\ &+ (0,0,0,...,r_{nn}\vec{q_n})\\ &= (\vec{q_1},\vec{q_2},...,\vec{q_n})\cdot \left[ \begin{matrix} r_{11} & r_{21} & \cdots & r_{n1} \\ 0 & r_{22} & \cdots & r_{n2} \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & r_{nn} \\ \end{matrix} \right]\\ &=Q\cdot R \end{aligned} A=(a1,a2,...,an)=(r11q1,r21q1+r22q2,...,i=1∑nrniqi)=(r11q1,r21q1,r31q1,...,rn1q1)+(0,r22q2,r32q2,...,rn2q2)+......+(0,0,0,...,rnnqn)=(q1,q2,...,qn)⋅⎣⎢⎢⎢⎡r110⋮0r21r22⋮0⋯⋯⋱⋯rn1rn2⋮rnn⎦⎥⎥⎥⎤=Q⋅R

2.2 Householder矩阵与Householder变换

在平面直角坐标系中，将向量a⃗=(c,d)\vec a=(c,d)a=(c,d)作关于x轴的交换，可得到：
b=[c−d]=[100−1][cd]=(I−2[01][01]T)x=Hxb = \left[ \begin{matrix} c\\-d \end{matrix} \right]=\left[ \begin{matrix} 1 & 0\\ 0 & -1 \end{matrix}\right] \left[ \begin{matrix} c\\ d \end{matrix} \right]=(I-2 \left[ \begin{matrix} 0\\ 1 \end{matrix} \right] \left[ \begin{matrix} 0\\ 1 \end{matrix} \right]^T)x=Hx b=[c−d]=[100−1][cd]=(I−2[01][01]T)x=Hx
可将其推广至：

2.2.1 定义

设单位列向量u∈Rnu \in R^nu∈Rn，称H=I−2uuTH=I-2uu^TH=I−2uuT为Householder矩阵(初等反射矩阵)，由Householder矩阵所确定的线性变换(y=Hxy=Hxy=Hx)称为Householder变换。

2.2.2 性质

(1) HT=HH^T=HHT=H（实对称），H−1=HTH^{-1}=H^TH−1=HT（正交）,H2=IH^2=IH2=I（对合），H−1=HH^{-1}=HH−1=H（自逆），det(H)=−1det(H)=-1det(H)=−1
(2) 对于任何非零列向量x∈Rnx\in R^nx∈Rn及任何单位列向量z∈Rnz\in R^nz∈Rn，存在Householder矩阵HHH，使得Hx=∣x∣zHx=|x|zHx=∣x∣z
(3) 初等旋转矩阵(Givens矩阵)是两个初等反射矩阵HHH的乘积
(定理证明请参照参考文献[2])

2.2.3 采用Householder变换的QR分解

A=[b(1)∗]A=\left[\begin{matrix} b^{(1)}& * \end{matrix}\right]A=[b(1)∗]，存在H1H_1H1，使得
H1b(1)=∣b(1)∣e1→H1A=[∣b(1)∣e1nA(1)]H_1b^{(1)} =\left|b^{(1)}\right|e_1\to H_1A=\left[\begin{matrix}\left|b^{(1)}\right|e_1^nA^{(1)}\end{matrix}\right] H1b(1)=∣∣∣b(1)∣∣∣e1→H1A=[∣∣b(1)∣∣e1nA(1)]
A(1)=[b(2)∗]A^{(1)}=\left[\begin{matrix} b^{(2)}& * \end{matrix}\right]A(1)=[b(2)∗]，存在H2H_2H2，使得H2A(1)=[∣b(2)∣e1n−1A(2)]LH_2A^{(1)} =\left[\begin{matrix}\left|b^{(2)}\right|e_1^{n-1}A^{(2)}\end{matrix}\right]LH2A(1)=[∣∣b(2)∣∣e1n−1A(2)]L

A(n−2)=[b(n−1)b(n)]A^{(n-2)}=\left[\begin{matrix} b^{(n-1)}& b^{(n)} \end{matrix}\right]A(n−2)=[b(n−1)b(n)]，存在Hn−1H_{n-1}Hn−1，使得
Hn−1A(n−2)=[an−1,n−1(n−1)an−1,n(n−1)0ann(n−1)]H_{n-1}A^{(n-2)} = \left[ \begin{matrix} a^{(n-1)}_{n-1,n-1} & a^{(n-1)}_{n-1,n}\\ \\ 0 & a_{nn}^{(n-1)} \end{matrix} \right] Hn−1A(n−2)=⎣⎢⎡an−1,n−1(n−1)0an−1,n(n−1)ann(n−1)⎦⎥⎤
令S=[In−200Hn−1][In−300Hn−2]L[I200H3][100H2]H1S=\left[ \begin{matrix} I_{n-2} & 0\\ \\ 0 & H_{n-1} \end{matrix} \right] \left[ \begin{matrix} I_{n-3} & 0\\ \\ 0 & H_{n-2} \end{matrix} \right]L \left[ \begin{matrix} I_{2} & 0\\ \\ 0 & H_{3} \end{matrix} \right] \left[ \begin{matrix} 1 & 0\\ \\ 0 & H_2 \end{matrix} \right] H_1S=⎣⎡In−200Hn−1⎦⎤⎣⎡In−300Hn−2⎦⎤L⎣⎡I200H3⎦⎤⎣⎡100H2⎦⎤H1

Hl+1=In−l−2uuT(u∈Rn−1,uTu=1)H_{l+1}=I_{n-l}-2uu^T(u\in R^{n-1},u^Tu=1)Hl+1=In−l−2uuT(u∈Rn−1,uTu=1)

则

[Il00Hl+1]=[Il00In−l]−2[000uuT]=In−2[0u][0TuT]=In−2vvT\left[ \begin{matrix} I_l & 0\\ \\ 0 & H_{l+1} \end{matrix} \right]= \left[ \begin{matrix} I_l & 0\\ \\ 0 & I_{n-l} \end{matrix} \right]-2 \left[ \begin{matrix} 0 & 0\\ \\ 0 & uu^T \end{matrix} \right]= I_n-2\left[ \begin{matrix} 0\\ \\ u \end{matrix} \right] \left[ \begin{matrix} 0^T & u^T\\ \end{matrix} \right]=I_n-2vv^T⎣⎡Il00Hl+1⎦⎤=⎣⎡Il00In−l⎦⎤−2⎣⎡000uuT⎦⎤=In−2⎣⎡0u⎦⎤[0TuT]=In−2vvT

vTv=[0TuT][0u]=uTu=1v^Tv= \left[ \begin{matrix} 0^T & u^T\\ \end{matrix} \right] \left[ \begin{matrix} 0\\ u \end{matrix} \right]=u^Tu=1vTv=[0TuT][0u]=uTu=1

SA=[a11(1)a12(1)La1n(1)0a22(2)La2n(2)00OM000annn−1]=R，S−1=QSA=\left[ \begin{matrix} a_{11}^{(1)} & a_{12}^{(1)} & L & a_{1n}^{(1)}\\ 0 & a_{22}^{(2)} & L & a_{2n}^{(2)}\\ 0 & 0 & O & M\\ 0 & 0 & 0 & a_{nn}^{n-1} \end{matrix} \right]=R，S^{-1}=QSA=⎣⎢⎢⎡a11(1)000a12(1)a22(2)00LLO0a1n(1)a2n(2)Mannn−1⎦⎥⎥⎤=R，S−1=Q

QQQ为正交矩阵

2.3 Givens矩阵与Givens变换

2.3.1 定义

设实数ccc与实数sss满足c2+s2=1c^2+s^2=1c2+s2=1，称

为Givens矩阵，记作Tij=Tij(c,s)T_{ij}=T_{ij}(c,s)Tij=Tij(c,s)，由Givens矩阵确定的线性变换成为Givens变换(初等旋转变换)。
说明：
(1).实数c2+s2=1c^2+s^2=1c2+s2=1，故存在θ\thetaθ，使c=cosθ，s=sinθc=cos\theta，s=sin\thetac=cosθ，s=sinθ
(2).y=Tijxy=T_{ij}xy=Tijx中TijT_{ij}Tij确定了将向量变成yyy的一种变换，正是Givens变换。二阶情况下，y=[cosθsinθ−sinθcosθ]xy=\left[\begin{matrix}cos\theta&sin\theta\\ -sin\theta&cos\theta\end{matrix}\right]xy=[cosθ−sinθsinθcosθ]x确定的正是平面直角坐标系中绕原点的一个旋转变换(旋转θ\thetaθ度)。

2.3.2 性质

(1) [Tij(c,s)]−1=[Tij(c,s)]T=Tij(c,−s)[T_{ij}(c,s)]^{-1}=[T_{ij}(c,s)]^{T}=T_{ij}(c,-s)[Tij(c,s)]−1=[Tij(c,s)]T=Tij(c,−s)，TijT_{ij}Tij为正交矩阵。det[Tij(c,s)]=1det[T_{ij}(c,s)]=1det[Tij(c,s)]=1
(2) 设x=[a1,a2,⋯,an]T,y=Tijx=[b1,b2,⋯,bn]x=[a_1,a_2,\cdots,a_n]^T,y=T_{ij}x=[b_1,b_2,\cdots,b_n]x=[a1,a2,⋯,an]T,y=Tijx=[b1,b2,⋯,bn]，则有,
{bi=c×ai+s×ajbj=−s×ai+c×ajbk=ak,k≠i,j\begin{cases} b_i = c\times a_i+s\times a_j \\ b_j = -s\times a_i + c \times a_j\\ b_k=a_k,k\neq i,j \end{cases} ⎩⎪⎨⎪⎧bi=c×ai+s×ajbj=−s×ai+c×ajbk=ak,k=i,j
(3) 设x=[a1,a2,⋯,an]T≠0x=[a_1,a_2,\cdots,a_n]^T\neq 0x=[a1,a2,⋯,an]T=0，则存在有限个Givens矩阵的乘积T，使得Tx=∣x∣e1Tx=| x |e_1Tx=∣x∣e1，
说明：(1).∣x∣=∥x∥22=xTx| x |=\sqrt{\lVert x \rVert_2^2}=\sqrt{x^Tx}∣x∣=∥x∥22=xTx
(2).e1=[1,0,0,⋯,0]Te_1=[1,0,0,\cdots,0]^Te1=[1,0,0,⋯,0]T
推论：对于任何非零列向量x∈Rnx \in R^nx∈Rn及任何单位列向量z(∣z∣=1)z(|z|=1)z(∣z∣=1)，均存在着有限个Givens矩阵的乘积T，使Tx=∣x∣zTx=| x |zTx=∣x∣z
(定理证明请参照参考文献[2])

2.3.3 应用Givens旋转的QR分解

可以用下面的这幅图来理解：x是没有改变的元素，m是改变了的元素。每一次箭头都进行了一次Givens旋转。

详细推导参见矩阵QR分解 Givens变换 Household变换

2.4 总结

无论是Gram–Schmid正交化、Household变换，还是Givens变换，对矩阵进行QR分解的思路都是对矩阵进行多次线性变换，直至分解成一个正定矩阵与一个上三角矩阵的乘积。不同的方法有不同的优点和缺点。

三、最小二乘法中的QR分解

已知线性回归的系数向量w=(XTX)−1XTyw=(X^TX)^{-1}X^Tyw=(XTX)−1XTy，我们为什么不直接求解呢？

3.1 条件数

计算机在进行运算的时候，有的时候会因为矩阵自身的特点而产生较大的误差。一个矩阵的条件数是它在计算机计算中的容易程度，条件数较大的时候矩阵计算就比较容易产生误差，此时该矩阵被称为病态矩阵。
条件数计算的公式为：condi(X)=∥X−1∥∥X∥condi(X)=\lVert X^{-1}\rVert \lVert X\rVertcondi(X)=∥X−1∥∥X∥
若∥⋅∥\lVert\cdot\rVert∥⋅∥为2范数，则condi(X)=σmax(X)σmin(X)condi(X)=\frac{\sigma_{max}(X)}{\sigma_{min}(X)}condi(X)=σmin(X)σmax(X)，其中σmax(X)\sigma_{max}(X)σmax(X)与σmin(X)\sigma_{min}(X)σmin(X)分别是矩阵XXX的极大奇异值和极小奇异值(奇异值参考)。
若我们使用最小二乘法推导出的公式计算系数向量，我们可以看到XTXX^TXXTX的条件数：
condi(XTX)=condi(VΣTUTUΣVH)=condi(VΣ2VH)=condi(X)2condi(X^TX)=condi(V\Sigma^TU^TU\Sigma V^H)=condi(V\Sigma ^2V^H)=condi(X)^2 condi(XTX)=condi(VΣTUTUΣVH)=condi(VΣ2VH)=condi(X)2 然而通过使用QR分解，我们可以将条件数降到尽量低。

3.2 最小二乘法与QR分解

我们首先将训练集XXX分解为X=QRX=QRX=QR，再带入上面计算系数矩阵的公式中。
w^∗=(XTX)−1XTyXTXw^∗=XTyRTQTQRw^∗=RTQTyRTRw^∗=RTQTy（QTQ=I）Rw^∗=QTyw^∗=R−1QTy\begin{aligned} \hat w^* &= (X^TX)^{-1}X^Ty \\ X^TX\hat w^* &= X^Ty \\ R^TQ^TQR\hat w^*&= R^TQ^Ty \\ R^TR\hat w^*&= R^TQ^Ty（Q^TQ=I） \\ R\hat w^*&= Q^Ty \\ \hat w^*&= R^{-1}Q^Ty \end{aligned} w^∗XTXw^∗RTQTQRw^∗RTRw^∗Rw^∗w^∗=(XTX)−1XTy=XTy=RTQTy=RTQTy（QTQ=I）=QTy=R−1QTy 我们仅需直到X的QR分解就可以直接计算出www。

四、参考文献

[1]. 用QR分解求最小二乘法的最优闭式解
[2]. 矩阵QR分解 Givens变换 Household变换
[3]. 矩阵的QR分解
[4]. QR decomposition and Givens Rotation
[4]. QR分解与最小二乘
[5]. QR分解之Household变换
[6]. 矩阵论基础知识