矩阵分析与应用-4.7-QR分解及其应用-Section2
前言
本文学习过程来源是《矩阵分析与应用-张贤达》一书. 可以通过 z-lib 下载.
一、采用 \(\mathrm{Givens}\) 旋转的 \(\mathrm{QR}\) 分解
\(\mathrm{Givens}\) 旋转也可以用来计算 \(\mathrm{QR}\) 分解. 这里以 \(4 \times 3\) 矩阵为例, 说明 \(\mathrm{Givens \ QR}\) 分解的思想.
\[ \begin{aligned} & \begin{bmatrix} \times & \times & \times \\ \times & \times & \times \\ \otimes & \times & \times\\ \otimes & \times & \times \end{bmatrix} \overset{G(3,4)}{\longrightarrow} \begin{bmatrix} \times & \times & \times \\ \otimes & \times & \times \\ \otimes & \times & \times\\ 0 & \times & \times \end{bmatrix} \overset{G(2,3)}{\longrightarrow} \begin{bmatrix} \otimes & \times & \times \\ \otimes & \times & \times \\ 0 & \times & \times\\ 0 & \times & \times \end{bmatrix} \overset{G(1,2)}{\longrightarrow} \begin{bmatrix} \times & \times & \times \\ 0 & \times & \times \\ 0 & \otimes & \times\\ 0 & \otimes & \times \end{bmatrix} \overset{G(3,4)}{\longrightarrow} \\ & \begin{bmatrix} \times & \times & \times \\ 0 & \otimes & \times \\ 0 & \otimes & \times\\ 0 & 0 & \times \end{bmatrix} \overset{G(2,3)}{\longrightarrow} \begin{bmatrix} \times & \times & \times \\ 0 & \times & \times \\ 0 & 0 & \otimes\\ 0 & 0 & \otimes \end{bmatrix} \overset{G(3,4)}{\longrightarrow} \begin{bmatrix} \times & \times & \times \\ 0 & \times & \times \\ 0 & 0 & \times\\ 0 & 0 & 0 \end{bmatrix} \end{aligned} \]
其中 \(\otimes\) 代表用 \(\mathrm{Givens}\) 旋转进行变换的元素. 变换过程就是乘以箭头上用 \(G(i,j)\) 表示的 \(\mathrm{Givens}\) 矩阵.
从上述说明中易得出结论: 如果令 \(G_j\) 代表约化过程中的第 \(j\) 次 \(\mathrm{Givens}\) 旋转, 则 \(Q^{\mathrm{T}}A=R\) 是上三角矩阵, 其中, \(Q = G_tG_{t-1} \cdots G_1\), 而 \(t\) 是总的旋转次数.
归根到底还是为了解方程, 不论是有解还是最小二乘法, \(\mathrm{QR}\) 分解都是一个不错的选择.
二、基于 \(\mathrm{QR}\) 分解的参数估计问题
系统辨识问题的提法是: 已知系统输入 \(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\) 时刻的系统辨识问题可以化为最小二乘问题.
看起来有点像预测方面的问题.
\[ \min_{\theta_n} \lVert A_n\theta_n - y_n \rVert^2_2 \tag{1} \]
求解, 并且其解由 "法方程"
\[ A_n^{\mathrm{T}}A_n\theta_n = A_n^{\mathrm{T}}y_n \ 或者 \ R_{xx}\theta_n = r_n \tag{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 \begin{bmatrix} R_n\\ O\end{bmatrix}, Q_n^{\mathrm{T}}y_n= \begin{bmatrix} \bar{y}_n\\ \tilde{y}_n \end{bmatrix}\), 其中, \(Q_n\) 是正交矩阵, \(R_n\) 是上三角矩阵. 故有
\[ \begin{aligned} \theta_{n+1} &= argmin_{\theta} \lVert A_{n+1}\theta - y_{n+1} \rVert^2_2 &= argmin_{\theta} \bigg \lVert \begin{bmatrix} \lambda R_n\\ x_{n+1}^{\mathrm{T}}\end{bmatrix}\theta - \begin{bmatrix} \lambda \bar{y}_n\\ y(n+1)\end{bmatrix} \bigg \rVert^2_2 \end{aligned} \tag{3} \]
算法 1: 系统参数的自适应估计算法
Step 1 : 对矩阵 \(\bar{R} = \begin{bmatrix} \lambda R_n\\ x_{n+1}^{\mathrm{T}}\end{bmatrix}\) 进行 \(\mathrm{QR}\) 分解, 得到
\[ Q_{n+1}^{\mathrm{T}} \bar{R} = Q_{n+1}^{\mathrm{T}}\begin{bmatrix} \lambda R_n\\ x_{n+1}^{\mathrm{T}}\end{bmatrix} = \begin{bmatrix} R_{n+1}\\ O\end{bmatrix} \tag{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} = \begin{bmatrix} \bar{y}_{n+1} \\ \tilde{y}_{n+1}\end{bmatrix} \]
其中, \(\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{Householder}\) 变换的快速时变参数估计
考查 \(n \times (p+1)\) 矩阵
\[ A_n = \begin{bmatrix} a_{11}& a_{12}& \cdots& a_{1,p+1}\\ a_{21}& a_{22}& \cdots& a_{2,p+1}\\ \vdots& \vdots& &\vdots \\ a_{n1}& a_{n2}& \cdots& a_{n,p+1} \end{bmatrix} \]
的 \(\mathrm{Householder \ QR}\) 分解, 即
\[ H_nA_n = \begin{bmatrix} a_{11}^*& a_{12}^*& \cdots& a_{1,p+1}^*\\ 0& a_{22}^*& \cdots& a_{2,p+1}^*\\ \vdots& \vdots& &\vdots \\ 0 & 0 & \cdots & a_{p+1,p+1}^*\\ 0 & 0 & \cdots & 0\\ \vdots& \vdots& &\vdots \\ 0 & 0 & \cdots & 0\\ \end{bmatrix} \tag{5} \]
为了得到上述 \(\mathrm{QR}\) 分解, 应该选择 \(H_n\) 为 \(p\) 个 \(\mathrm{Householder}\) 变换矩阵之积, 即
\[ H_n = H_n(p)H_n(p-1)\cdots H_n(1) \tag{6} \]
式中
\[ H_n(j) = I - u_ju_j^{\mathrm{T}}/\sigma_j, \quad j = 1,2,\cdots,p \tag{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{Householder}\) 变换矩阵, 其参数选择方法为
\[ \left.\begin{aligned} \alpha_j &= \sqrt{\sum_{i=j}^{n}[a_{ij}^{(j)}]^2} & \quad\\ \sigma_j &= \alpha_j(\alpha_j + |a_{jj}^{(j)}|) \\ u_j(i) &= \left\{\begin{matrix} 0& \ i < j \\ a_{jj}^{(j)} + \mathrm{sgn}(a_{jj}^{(j)})\alpha_j& \ j = i \\ a_{ij}^{(j)} & i > j \end{matrix}\right. \end{aligned}\right\}, \qquad j = 1,2,\cdots,p \tag{8} \]
其中
\[ A_n^{(j+1)} = A_n^{(j)} - u_jq_j^{\mathrm{T}} \tag{9} \]
并且
\[ q_j^{\mathrm{T}} = u_j^{\mathrm{T}}A_n^{(j)}/\sigma_j \tag{10} \]
算法 2: 基于 \(\mathrm{Householder} \ QR\) 分解的快速自适应参数估计算法
四、基于 \(\mathrm{Givens}\) 旋转的时变参数估计
递推求解 \(\sigma_n\) 的变换量 \(\delta_n\), 而不是之间递推求 \(\sigma_{n+1}\) 本身. 其式子应该为
\[ \sigma_{n+1} = \sigma_{n} + \delta_n \tag{11} \]
问题的关键就在于更新 \(\delta_n\)
假定正交矩阵 \(\tilde{Q}\) 为已知, 满足
\[ \tilde{Q} \begin{bmatrix} \lambda R_n\\ x_{n+1}^{\mathrm{T}}\end{bmatrix} = \begin{bmatrix} R_{n+1}\\ O \end{bmatrix} \tag{12} \]
化简得到
\[ \delta_n = argmin_{\delta_n} \bigg \lVert \begin{bmatrix} R_{n+1}\\ O \end{bmatrix}\delta_n - \tilde{Q} \begin{bmatrix} 0\\ u(n+1) \end{bmatrix} \bigg \rVert \tag{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} \tag{14} \]
解出, 其中, \(\bar{y}_{k+1}\) 满足
\[ \begin{bmatrix} \bar{y}_{n+1}\\ r(n+1) \end{bmatrix} = \tilde{Q} \begin{bmatrix} 0\\ u(n+1) \end{bmatrix} \tag{15} \]
为求出 \(\tilde{Q}\), 需要对增广矩阵
\[ \begin{bmatrix} \lambda R_n & 0\\ x_{n+1}^{\mathrm{T}} & u(n+1) \end{bmatrix} \tag{16} \]
( ! ! ! 在这个地方存疑, 不能很好的理解这个增广矩阵的由来 )
执行所需要的清零. 综上所述, 每一步递推更新需要的步骤如下
计算预测误差 \(y_{k+1} - \phi_{k+1}^{\mathrm{T}}\theta_k\)
形成式子 (16) 中的 \((n+1) \times (n+1)\) 矩阵
利用一系列 \(\mathrm{Givens}\) 旋转将上述矩阵最底一行的左边 \(n\) 个元素扫除为零
解上三角矩阵方程得到 \(\delta_k\)