矩阵分析与应用-4.7-QR分解及其应用-Section2

前言

本文学习过程来源是《矩阵分析与应用-张贤达》一书. 可以通过 z-lib 下载.

一、采用 $\mathrm{G}\mathrm{i}\mathrm{v}\mathrm{e}\mathrm{n}\mathrm{s}$ 旋转的 $\mathrm{Q}\mathrm{R}$ 分解

$\mathrm{G}\mathrm{i}\mathrm{v}\mathrm{e}\mathrm{n}\mathrm{s}$ 旋转也可以用来计算 $\mathrm{Q}\mathrm{R}$ 分解. 这里以 $4\times 3$ 矩阵为例, 说明 $\mathrm{G}\mathrm{i}\mathrm{v}\mathrm{e}\mathrm{n}\mathrm{s}\mathrm{Q}\mathrm{R}$ 分解的思想.

$$\begin{array}{rl}& \left[\begin{array}{ccc}\times & \times & \times \\ \times & \times & \times \\ \otimes & \times & \times \\ \otimes & \times & \times \end{array}\right]\stackrel{G(3,4)}{⟶}\left[\begin{array}{ccc}\times & \times & \times \\ \otimes & \times & \times \\ \otimes & \times & \times \\ 0& \times & \times \end{array}\right]\stackrel{G(2,3)}{⟶}\left[\begin{array}{ccc}\otimes & \times & \times \\ \otimes & \times & \times \\ 0& \times & \times \\ 0& \times & \times \end{array}\right]\stackrel{G(1,2)}{⟶}\left[\begin{array}{ccc}\times & \times & \times \\ 0& \times & \times \\ 0& \otimes & \times \\ 0& \otimes & \times \end{array}\right]\stackrel{G(3,4)}{⟶}\\ & \left[\begin{array}{ccc}\times & \times & \times \\ 0& \otimes & \times \\ 0& \otimes & \times \\ 0& 0& \times \end{array}\right]\stackrel{G(2,3)}{⟶}\left[\begin{array}{ccc}\times & \times & \times \\ 0& \times & \times \\ 0& 0& \otimes \\ 0& 0& \otimes \end{array}\right]\stackrel{G(3,4)}{⟶}\left[\begin{array}{ccc}\times & \times & \times \\ 0& \times & \times \\ 0& 0& \times \\ 0& 0& 0\end{array}\right]\end{array}$$

其中 $\otimes$ 代表用 $\mathrm{G}\mathrm{i}\mathrm{v}\mathrm{e}\mathrm{n}\mathrm{s}$ 旋转进行变换的元素. 变换过程就是乘以箭头上用 $G(i,j)$ 表示的 $\mathrm{G}\mathrm{i}\mathrm{v}\mathrm{e}\mathrm{n}\mathrm{s}$ 矩阵.

从上述说明中易得出结论: 如果令 $G_j$ 代表约化过程中的第 $j$ 次 $\mathrm{G}\mathrm{i}\mathrm{v}\mathrm{e}\mathrm{n}\mathrm{s}$ 旋转, 则 $Q^{\mathrm{T}}A=R$ 是上三角矩阵, 其中, $Q=G_t{G}_{t-1}\cdots G_1$, 而 $t$ 是总的旋转次数.

归根到底还是为了解方程, 不论是有解还是最小二乘法, $\mathrm{Q}\mathrm{R}$ 分解都是一个不错的选择.

二、基于 $\mathrm{Q}\mathrm{R}$ 分解的参数估计问题

系统辨识问题的提法是: 已知系统输入 $x(k)$ 和输出观测值 $y(k)$, 其中, $k=1,2,\cdots ,n$ 估计系统参数向量 $\theta$. 在时变系统的辨识中, 则要求在已估计 $n$ 时刻的系统参数向量 ${\theta }_n$ 的情况下, 使用增加的 $x(n+1),y(n+1)$ 值, 通过简单的运算, 递推出 $n+1$ 时刻的系统参数向量 ${\theta }_{n+1}$. $n$ 时刻的系统辨识问题可以化为最小二乘问题.

看起来有点像预测方面的问题.

$$\underset{{\theta }_n}{\min }\Vert A_n{\theta }_n-y_n{\Vert }_2^2\text{\hspace{1em}(1)}$$

求解, 并且其解由 “法方程”

$$A_n^{\mathrm{T}}{A}_n{\theta }_n=A_n^{\mathrm{T}}{y}_n或者R_{xx}{\theta }_n=r_n\text{\hspace{1em}(2)}$$

确定. 式中, $R_{xx}=A_n^{\mathrm{T}}{A}_n$ 代表系统输入 $x(k)$ 的协方差矩阵, $r_n=A_n^{\mathrm{T}}{y}_n$.

之间求解式 (2) 的方法叫做协方差方法.

引理 1: 若 $A_n=Q_n\left[\begin{array}{c}{R}_n\\ O\end{array}\right],Q_n^{\mathrm{T}}{y}_n=\left[\begin{array}{c}{\bar{y}}_n\\ {\tilde{y}}_n\end{array}\right]$, 其中, $Q_n$ 是正交矩阵, $R_n$ 是上三角矩阵. 故有

$$\begin{array}{rlr}{\theta }_{n+1}& =\underset{\theta }{\mathrm{argmin}}\Vert A_{n+1}\theta -y_{n+1}{\Vert }_2^2& =\underset{\theta }{\mathrm{argmin}}\Vert \left[\begin{array}{c}\lambda R_n\\ x_{n+1}^{\mathrm{T}}\end{array}\right]\theta -\left[\begin{array}{c}\lambda {\bar{y}}_n\\ y(n+1)\end{array}\right]{\Vert }_2^2\end{array}\text{\hspace{1em}(3)}$$

算法 1: 系统参数的自适应估计算法

Step 1 : 对矩阵 $\bar{R}=\left[\begin{array}{c}\lambda R_n\\ x_{n+1}^{\mathrm{T}}\end{array}\right]$ 进行 $\mathrm{Q}\mathrm{R}$ 分解, 得到

$$Q_{n+1}^{\mathrm{T}}\bar{R}=Q_{n+1}^{\mathrm{T}}\left[\begin{array}{c}\lambda R_n\\ x_{n+1}^{\mathrm{T}}\end{array}\right]=\left[\begin{array}{c}{R}_{n+1}\\ O\end{array}\right]\text{\hspace{1em}(4)}$$

式子中, $Q_{n+1}$ 是 $(n+1)\times (n+1)$ 正交矩阵, $R_{n+1}$ 为 $(p+1)\times (p+1)$ 上三角矩阵, 且 $O$ 是 $(n-p)\times (p+1)$ 零矩阵.

Step 2 : 进行分块运算

$$Q_{n+1}^{\mathrm{T}}{y}_{n+1}=\left[\begin{array}{c}{\bar{y}}_{n+1}\\ {\tilde{y}}_{n+1}\end{array}\right]$$

其中, ${\bar{y}}_{n+1}$ 为 $(p+1)\times 1$ 向量, ${\tilde{y}}_{n+1}$ 为 $(n-p)\times 1$ 向量

Step 3 : 求解三角矩阵方程 $R_{n+1}{\theta }_{n+1}={\bar{y}}_{n+1}$ 得到 ${\theta }_{n+1}$

三、基于 $\mathrm{H}\mathrm{o}\mathrm{u}\mathrm{s}\mathrm{e}\mathrm{h}\mathrm{o}\mathrm{l}\mathrm{d}\mathrm{e}\mathrm{r}$ 变换的快速时变参数估计

考查 $n\times (p+1)$ 矩阵

$$A_n=\left[\begin{array}{cccc}{a}_{11}& a_{12}& \cdots & a_{1,p+1}\\ a_{21}& a_{22}& \cdots & a_{2,p+1}\\ ⋮& ⋮& & ⋮\\ a_{n1}& a_{n2}& \cdots & a_{n,p+1}\end{array}\right]$$

的 $\mathrm{H}\mathrm{o}\mathrm{u}\mathrm{s}\mathrm{e}\mathrm{h}\mathrm{o}\mathrm{l}\mathrm{d}\mathrm{e}\mathrm{r}\mathrm{Q}\mathrm{R}$ 分解, 即

$$H_n{A}_n=\left[\begin{array}{cccc}{a}_{11}^{\ast }& a_{12}^{\ast }& \cdots & a_{1,p+1}^{\ast }\\ 0& a_{22}^{\ast }& \cdots & a_{2,p+1}^{\ast }\\ ⋮& ⋮& & ⋮\\ 0& 0& \cdots & a_{p+1,p+1}^{\ast }\\ 0& 0& \cdots & 0\\ ⋮& ⋮& & ⋮\\ 0& 0& \cdots & 0\end{array}\right]\text{\hspace{1em}(5)}$$

为了得到上述 $\mathrm{Q}\mathrm{R}$ 分解, 应该选择 $H_n$ 为 $p$ 个 $\mathrm{H}\mathrm{o}\mathrm{u}\mathrm{s}\mathrm{e}\mathrm{h}\mathrm{o}\mathrm{l}\mathrm{d}\mathrm{e}\mathrm{r}$ 变换矩阵之积, 即

$$H_n=H_n(p)H_n(p-1)\cdots H_n(1)\text{\hspace{1em}(6)}$$

式中

$$H_n(j)=I-u_j{u}_j^{\mathrm{T}}/{\sigma }_j,j=1,2,\cdots ,p\text{\hspace{1em}(7)}$$

是对矩阵 $A_n^{(j)}=H_n(j-1)H_n(2)H_n(1)A_n$ 第 $j$ 列向量 $[a_{1j}^{(j)},a_{2j}^{(j)},\cdots ,a_{nj}^{(j)}{]}^{\mathrm{T}}$ 进行的 $\mathrm{H}\mathrm{o}\mathrm{u}\mathrm{s}\mathrm{e}\mathrm{h}\mathrm{o}\mathrm{l}\mathrm{d}\mathrm{e}\mathrm{r}$ 变换矩阵, 其参数选择方法为

$$\begin{array}{rlr}{\alpha }_j& =\sqrt{\sum _{i=j}^n[a_{ij}^{(j)}{]}^2}& \\ {\sigma }_j& ={\alpha }_j({\alpha }_j+|a_{jj}^{(j)}|)& \\ u_j(i)& =\{\begin{array}{cc}0& i<j\\ a_{jj}^{(j)}+\mathrm{s}\mathrm{g}\mathrm{n}(a_{jj}^{(j)}){\alpha }_j& j=i\\ a_{ij}^{(j)}& i>j\end{array}& \end{array}\},j=1,2,\cdots ,p\text{\hspace{1em}(8)}$$

其中

$$A_n^{(j+1)}=A_n^{(j)}-u_j{q}_j^{\mathrm{T}}\text{\hspace{1em}(9)}$$

并且

$$q_j^{\mathrm{T}}=u_j^{\mathrm{T}}{A}_n^{(j)}/{\sigma }_j\text{\hspace{1em}(10)}$$

算法 2: 基于 $\mathrm{H}\mathrm{o}\mathrm{u}\mathrm{s}\mathrm{e}\mathrm{h}\mathrm{o}\mathrm{l}\mathrm{d}\mathrm{e}\mathrm{r}QR$ 分解的快速自适应参数估计算法

四、基于 $\mathrm{G}\mathrm{i}\mathrm{v}\mathrm{e}\mathrm{n}\mathrm{s}$ 旋转的时变参数估计

递推求解 ${\sigma }_n$ 的变换量 ${\delta }_n$, 而不是之间递推求 ${\sigma }_{n+1}$ 本身. 其式子应该为

$${\sigma }_{n+1}={\sigma }_n+{\delta }_n\text{\hspace{1em}(11)}$$

问题的关键就在于更新 ${\delta }_n$

假定正交矩阵 $\tilde{Q}$ 为已知, 满足

$$\tilde{Q}\left[\begin{array}{c}\lambda R_n\\ x_{n+1}^{\mathrm{T}}\end{array}\right]=\left[\begin{array}{c}{R}_{n+1}\\ O\end{array}\right]\text{\hspace{1em}(12)}$$

化简得到

$${\delta }_n=\underset{{\delta }_n}{\mathrm{argmin}}\Vert \left[\begin{array}{c}{R}_{n+1}\\ O\end{array}\right]{\delta }_n-\tilde{Q}\left[\begin{array}{c}0\\ u(n+1)\end{array}\right]\Vert \text{\hspace{1em}(13)}$$

式中, $u(n+1)=y(n+1)-x_{n+1}^{\mathrm{T}}{\theta }_n$. 因此, ${\delta }_n$ 可以从三角矩阵方程

$$R_{n+1}{\delta }_n={\bar{y}}_{n+1}\text{\hspace{1em}(14)}$$

解出, 其中, ${\bar{y}}_{k+1}$ 满足

$$\left[\begin{array}{c}{\bar{y}}_{n+1}\\ r(n+1)\end{array}\right]=\tilde{Q}\left[\begin{array}{c}0\\ u(n+1)\end{array}\right]\text{\hspace{1em}(15)}$$

为求出 $\tilde{Q}$, 需要对增广矩阵

$$\left[\begin{array}{cc}\lambda R_n& 0\\ x_{n+1}^{\mathrm{T}}& u(n+1)\end{array}\right]\text{\hspace{1em}(16)}$$

( ! ! ! 在这个地方存疑, 不能很好的理解这个增广矩阵的由来 )

执行所需要的清零. 综上所述, 每一步递推更新需要的步骤如下

(1) 计算预测误差 $y_{k+1}-{\varphi }_{k+1}^{\mathrm{T}}{\theta }_k$

(2) 形成式子 (16) 中的 $(n+1)\times (n+1)$ 矩阵

(3) 利用一系列 $\mathrm{G}\mathrm{i}\mathrm{v}\mathrm{e}\mathrm{n}\mathrm{s}$ 旋转将上述矩阵最底一行的左边 $n$ 个元素扫除为零

(4) 解上三角矩阵方程得到 ${\delta }_k$