最小二成问题的QR分解法

最小二成问题的QR分解法

一、QR分解法

$$\begin{gathered} \Vert Ax-b{\Vert }_2=\underset{z\in C^n}{\min }\Vert Az-b{\Vert }_2 \\ \Vert Q^T(Ax-b){\Vert }_2^2=\Vert Rx-f_1{\Vert }_2^2+\Vert f_2{\Vert }_2^2 \\ 其中f=Q^{T}b=\left[\begin{array}{c}{Q}_1^T\\ Q_2^T\end{array}\right]b=\left[\begin{array}{c}{f}_1\\ f_2\end{array}\right] \\ Rx=f_1 \end{gathered}$$
对$\hat{A}=[A,b]$进行QR分解
$$\begin{gathered} H_n\dots H_1[A,b]=[R,\hat{b}] \\ \Rightarrow Q^T[A,b]=[R,\hat{b}] \\ A=QR,\hat{b}=Q^{T}b,Q=H_1\dots H_n \end{gathered}$$
所以对于列满秩矩阵$A\in R^{m\ast n}$，LS问题QR分解方法将问题转化为求解
$$Rx=f_1$$
可先进行QR分解，得到此方程，再用回代法求解$x_{LS}$

二、求解秩亏最小二乘问题的数值解法

考虑 $\min \Vert Ax-b{\Vert }_2$，其中$A\in R_r^{m\ast n}(m\ge n)$,$rank(A)=r<n$。此时最小二乘问题有无穷个解。

1. 列主元QR分解法

A是秩亏损时，QR分解不一定能给出列空间R(A)的一组标准正交基，但是我们可以通过计算列置换之后的矩阵$\hat{A}=AP$的QR分解来求解来求解，即AP=QR，其中P是置换矩阵
设$b=b_1+b_2,b_1\in R(A),b_2\in R^{\perp }(A)$
可以证明：$Ax=b_1$
现在假定$R(A)=R(Q_1)$其中$Q_1\in R^{m\ast r}$且$Q^{T}Q=I_r$，即$Q_1$构成R(A)的一组标准正交基，则存在矩阵$S\in R^{r\ast n}$和向量$h\in R^r$，使得
$$\begin{gathered} A=Q_{1}S,b_1=Q_{1}h \\ 又\because Ax=b_1 \\ \therefore Sx=h \end{gathered}$$
显然方程组总是有解的，进一步可知
$$\begin{gathered} S=Q_1^{T}A, \\ h=Q_1^T{b}_1=Q_1^T(b-b_2)=Q_1^{T}b \\ {b}_2\in R^{\perp }(Q_1) \end{gathered}$$
由于rank(A)=r < n,一般不能直接产生R(A)的一组标准正交基，但是如果对A进行适当的排列使其前r列线性无关，然后再进行QR分解，则仍然可产生R(A)的一组标准正交基
设有分解式
$$AP=Q\left[\begin{array}{cc}{R}_{11}& R_{12}\\ 0& 0\end{array}\right]$$
其中P为置换矩阵，Q为正交矩阵，$R_{11}$为非奇异上三角阵，则Q的前r列就是R(A)的一组标准正交基
一旦求出分解式，就有
$$\begin{gathered} (Q^{T}AP)(P^{T}x)=Q^T{b}_1 \\ \left[\begin{array}{cc}{R}_{11}& R_{12}\\ 0& 0\end{array}\right]\left[\begin{array}{c}w\\ z\end{array}\right]=Q^T{b}_1=\left[\begin{array}{c}{Q}_1^T\\ Q_2^T\end{array}\right]b_1=\left[\begin{array}{c}h\\ g\end{array}\right] \\ \Rightarrow \left[\begin{array}{c}{R}_{11}w+R_{12}z\\ 0\end{array}\right]=\left[\begin{array}{c}h\\ g\end{array}\right] \\ \therefore R_{11}w+R_{12}z=h,g=0 \\ w=R_{11}^{-1}(h-R_{12}z),g=0 \end{gathered}$$

$$x=P\left[\begin{array}{c}w\\ z\end{array}\right]=P\left[\begin{array}{c}{R}_{11}^{-1}(h-R_{12}z)\\ z\end{array}\right]=P\left[\begin{array}{c}{R}_{11}^{-1}h\\ 0\end{array}\right]+\left[\begin{array}{c}-R_{11}^{-1}{R}_{12}\\ I_{n-r}\end{array}\right]z$$

$$令x_b=P\left[\begin{array}{c}{R}_{11}^{-1}h\\ 0\end{array}\right]$$
$x$就是最小二乘基本解
$$A{x}_b=\left[\begin{array}{c}{b}_1\\ 0\end{array}\right]$$

求解P

考虑分解式(8)的具体实施过程。假设对对某一正整数k，已求k-1个Householder变换$H_1,H_2,\cdots ,H_{k-1}$和k-1个初等变换，$P_1,P_2,\cdots ,P_{k-1}$使得

$R_{k-1}=H_{k-1}\cdots H_2{H}_{1}A{P}_1{P}_2\cdots P_{k-1}=\left[\begin{array}{cc}{R}_{11}^{(k-1)}& R_{12}^{(k-1)}\\ 0& R_{22}^{(k-1)}\end{array}\right]$，其中$R_{11}^{(k-1)} $为非奇异上三角矩阵，记$R_{22}^{(k-1)}=[V_k{(k-1)},V_{k+1}^{(k-1)}，\cdots ,V_n^{(k-1)}] $，其中$V_i^{{(k-1)}$表示$R_{22}}{(k-1)}$的第（i-k+1)列，下一步是确定指标P($k\le P \le n$)满足

∥ V p ( k − 1 ) ∥ 2 = max ⁡ { ∥ V k ( k − 1 ) ∥ 2 , ∥ V k + 1 ( k − 1 ) ∥ 2 , ⋯   , ∥ V n ( k − 1 ) ∥ 2 } \lVert V_p^{(k-1)}\rVert_2=\max \{\lVert V_k^{(k-1)}\rVert_2,\lVert V_{k+1}^{(k-1)}\rVert_2,\cdots,\lVert V_ n^{(k-1)}\rVert_2\} ∥Vp(k−1)​∥2​=max{∥Vk(k−1)​∥2​,∥Vk+1(k−1)​∥2​,⋯,∥Vn(k−1)​∥2​}

若$\Vert V_p^{(k-1)}{\Vert }_2=0$，停止计算。否则取$P_k$为第k列与第p列交换的初等列变换矩阵，并确定一个Householder变换，${\hat{H}}_k\in R^{(m-k+1)\times (m-k+1)}$使得${\hat{H}}_k{V}_p^{(k-1)}={\gamma }_{kk}{e}_1$

令$H_k=diag(I_{k-1},{\hat{H}}_k)$，则有$R_k=\left[\begin{array}{cc}{R}_{11}^{(k)}& R_{12}^{(k)}\\ 0& R_{22}^{(k)}\end{array}\right]$，其中$R_{11}^{(k-1)} $为非奇异上三角矩阵$R_k=H_k\begin{bmatrix}R_{11}^{(k-1)} & R_{12}^{(k-1)} \0 & R_{22}^{(k-1)}\end{bmatrix}P_k$

$R_{12}$约化为0矩阵，$z_1,z_2,\cdots ,z_r$和置换，$\hat{P}$，使得$[R_{11},\cdots ,R_{1r}]z_1\cdots z_r\hat{P}=[T,0]$，T为上三角矩阵，令

$Z=P{z}_1\cdots z_r$，则有$Q^{T}AZ=\left[\begin{array}{cc}T& 0\\ 0& 0\end{array}\right]$，由此可得$x_{LS}=z\left[\begin{array}{c}{T}^{-1}C\\ 0\end{array}\right]$，C为由$Q^{T}b$前r个分量构成的r维向量

算法

秩亏最小二乘问题的列主元QR分解法

• Step1：输入矩阵A，计算$P(j):=j，{\sigma }_j=A(:j{)}^{T}A(:j)$，置k=1
• Step2： 确定r使满足${\sigma }_r:{\max }_{k\le j\le n}{\sigma }_j$，若${\sigma }_r=0$停算
• Step3：交换$P(r)$与$P(k)$，${\sigma }_r$与${\sigma }_k$以及A(:,r)与A(:,k)
• Step4：确定一个Householder矩阵${\hat{H}}_k$使得${\hat{H}}_{r}A(k:m,k)={\gamma }_{kk}{e}_1$并计算$[A,b]:=diag(I_{k-1},{\hat{H}}_{k-1})[A,b],{\sigma }_j:={\sigma }_j-a_{kj}^2$
• Step5：置k:=k+1，转步骤2

变换技巧
$$\begin{gathered} \Vert Q^{T}x{\Vert }_2^2=\Vert [\alpha ,z{]}^T{\Vert }_2^2={\alpha }^2+\Vert z{\Vert }_2^2 \\ \Vert z{\Vert }_2^2=\Vert x{\Vert }_2^2-{\alpha }^2 \\ \therefore \Vert x^{(j)}{\Vert }_2^2=\Vert x^{(j-1)}{\Vert }_2^2-{\alpha }_{kj}^2 \end{gathered}$$

二、奇异值分解法

设$A\in R^{m\times n}(m>n)$的奇异值分解为$A=U\Sigma V^T$,$\sigma =\left[\begin{array}{cc}{\Sigma }_r& 0\\ 0& 0\end{array}\right]$，其中$U=[u_1,\cdots ,u_m],V=[v_1,\cdots ,v_n],{\Sigma }_r=diag({\sigma }_1,\cdots ,{\sigma }_r)$，则有
$$\begin{gathered} x_{LS}=A^{+}b=V\left[\begin{array}{cc}{\Sigma }_r^{-1}& 0\\ 0& 0\end{array}\right]U^T \\ =\sum _{i=1}^r\frac{u_i^{T}b}{{\sigma }_i}{v}_i \end{gathered}$$