Fisher-Score Algorithm of WTLS Estimation for PEIV Model

PEIV模型^[7]可表示为：

式中，L是n×1观测向量；X是t×1未知参数向量；I_n是n×n单位矩阵；Δ是n×1随机误差向量；h是nt×1已知向量；B是nt×s已知结构矩阵；a是系数矩阵A中不同随机元素所组成的s×1已知向量，s表示系数矩阵A中不同随机元素的个数；a为向量a的真值，e =a－a，A=invvec(h+Ba)，invvec是一个将nt×1向量变换为n×t矩阵的函数^[17]。

其中，σ²为未知的单位权方差。

基于WTLS原理，构造如下目标函数:

式(3)分别对a和X求偏导，并令其为零，可得未知参数向量X和向量a的估值分别为：

其中，$\mathit{\boldsymbol{h = }}\left[\begin{gathered} {\mathit{\boldsymbol{h}}_1} \hfill \\ {\mathit{\boldsymbol{h}}_2} \hfill \\ \vdots \hfill \\ {\mathit{\boldsymbol{h}}_t} \hfill \\ \end{gathered} \right], \mathit{\boldsymbol{B}} = \left[\begin{gathered} {\mathit{\boldsymbol{B}}_1} \hfill \\ {\mathit{\boldsymbol{B}}_2} \hfill \\ \vdots \hfill \\ {\mathit{\boldsymbol{B}}_t} \hfill \\ \end{gathered} \right]$, N_h、N_B、N_hB和N_Bh的第(i, j)个元素分别为：N_h(i, j)=h_i^TQ_Δ^－1h_j, ${\mathit{\boldsymbol{N}}_\mathit{\boldsymbol{B}}}\left( {i, j} \right) = {\hat {\bar a}^{\text{T}}}\mathit{\boldsymbol{B}}_i^{\text{T}}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}^{-1}{\mathit{\boldsymbol{B}}_j}\mathit{\boldsymbol{\hat {\bar a}}}$, ${\mathit{\boldsymbol{N}}_{\mathit{\boldsymbol{hB}}}}\left( {i, j} \right) = \mathit{\boldsymbol{h}}_i^{\text{T}}Q_\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}^{-1}{\mathit{\boldsymbol{B}}_j}\mathit{\boldsymbol{\hat {\bar a}}}$, ${\mathit{\boldsymbol{N}}_{\mathit{\boldsymbol{Bh}}}}\left( {i, j} \right) = {\mathit{\boldsymbol{\hat {\bar a}}}^{\text{T}}}\mathit{\boldsymbol{B}}_i^{\text{T}}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}^{-1}{\mathit{\boldsymbol{h}}_j}$, i, j=1, 2…t；${\mathit{\boldsymbol{S}}_\mathit{\boldsymbol{X}}} = \left( {{{\mathit{\boldsymbol{\hat X}}}^{\text{T}}} \otimes {\mathit{\boldsymbol{I}}_n}} \right)\mathit{\boldsymbol{B = }}\sum\limits_{i = 1}^t {{\mathit{\boldsymbol{B}}_i}{{\mathit{\boldsymbol{\hat X}}}_i}} $, ${\mathit{\boldsymbol{u}}_\mathit{\boldsymbol{h}}} = \left[\begin{gathered} \mathit{\boldsymbol{h}}_1^{\text{T}} \hfill \\ \mathit{\boldsymbol{h}}_2^{\text{T}} \hfill \\ \vdots \hfill \\ \mathit{\boldsymbol{h}}_t^{\text{T}} \hfill \\ \end{gathered} \right]\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}^{ -1}\mathit{\boldsymbol{L}}$，${\mathit{\boldsymbol{u}}_\mathit{\boldsymbol{B}}} = \left[\begin{gathered} {{\mathit{\boldsymbol{\hat {\bar a}}}}^{\text{T}}}\mathit{\boldsymbol{B}}_1^{\text{T}} \hfill \\ {{\mathit{\boldsymbol{\hat {\bar a}}}}^{\text{T}}}\mathit{\boldsymbol{B}}_2^{\text{T}} \hfill \\ \vdots \hfill \\ {{\mathit{\boldsymbol{\hat {\bar a}}}}^{\text{T}}}\mathit{\boldsymbol{B}}_t^{\text{T}} \hfill \\ \end{gathered} \right]\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}^{ -1}\mathit{\boldsymbol{L}}$。

从式(4)和式(5)可以发现，上述公式是交替迭代计算$\mathit{\boldsymbol{\hat X}}$和$\mathit{\boldsymbol{\hat {\bar a}}}$，较易影响计算效率。在统计领域，Fisher-Score方法是一种计算非线性模型参数估计的迭代方法，具有快速收敛性和重要应用价值^[8]。鉴于此，本文充分利用观测值误差和系数矩阵误差的统计性质，提出PEIV模型WTLS估计的Fisher-Score方法。

对式(1)作如下变形：

或

显然上式消掉了待估参数$\mathit{\boldsymbol{\bar a}}$，将其体现在系数矩阵A中，故下面仅讨论估计未知参数X即可。结合观测误差和系数矩阵误差的随机性质，有

利用加权残差平方和最小原理以及式(7)和式(8)，可构造如下非线性目标函数：

这里协方差矩阵$\mathit{\boldsymbol{ \boldsymbol{\tilde \varSigma} }}$可看作总误差的协方差矩阵。将(10)式的两边关于X求偏导，得

其中，X_i为X的第i个分量。

应用矩阵微分等式^[17]：

式中，l表示x的函数，则式(11)中$\frac{{\partial {{\mathit{\boldsymbol{ \boldsymbol{\tilde \varSigma} }}}^{-1}}}}{{\partial {\mathit{\boldsymbol{X}}_j}}}$的具体表达式为：

其中，${\mathit{\boldsymbol{\tilde h}}_j}$为第j个元素为1的单位向量，即${\mathit{\boldsymbol{\tilde h}}_j}$=[0…1…0]^T。注意到式(14)：

令${\mathit{\boldsymbol{S}}_{{B_j}}} = \left( {{\mathit{\boldsymbol{X}}^{\text{T}}} \otimes {I_n}} \right)\mathit{\boldsymbol{B}}{\mathit{\boldsymbol{Q}}_{{e_a}}}{\mathit{\boldsymbol{B}}^{\text{T}}}\left( {{{\mathit{\boldsymbol{\tilde h}}}_j} \otimes {I_n}} \right)$，则式(11)可表示为：

将目标函数Ω对X求二阶偏导数，是在式(15)的基础上再对X求一阶偏导数。为简便起见，先对单个分量求导，设定A_i为A的第i列元素所组成的向量，则有：

对式(16)中的第一项变形，得式(17)。而式(16)中的第二项的具体求解如式(18)所示。

将式(17)、式(18)代入式(16)中，并进一步求其期望，得:

其中，tr[]表示矩阵的迹。令${g_{ij}} = {\text{tr}}\left[{{{\mathit{\boldsymbol{ \boldsymbol{\tilde \varSigma} }}}^{-1}}\left( {{\mathit{\boldsymbol{S}}_{{B_j}}}{{\mathit{\boldsymbol{ \boldsymbol{\tilde \varSigma} }}}^{-1}}\;{\mathit{\boldsymbol{S}}_{{B_i}}} + \mathit{\boldsymbol{S}}_{{B_j}}^{\text{T}}{{\mathit{\boldsymbol{ \boldsymbol{\tilde \varSigma} }}}^{-1}}\;{\mathit{\boldsymbol{S}}_{{B_i}}} + {\mathit{\boldsymbol{S}}_{{B_i}}}{{\mathit{\boldsymbol{ \boldsymbol{\tilde \varSigma} }}}^{ - 1}}\;{\mathit{\boldsymbol{S}}_{{B_j}}} + {\mathit{\boldsymbol{S}}_{{B_i}}}{{\mathit{\boldsymbol{ \boldsymbol{\tilde \varSigma} }}}^{ - 1}}\;\mathit{\boldsymbol{S}}_{{B_j}}^{\text{T}} - \left( {\mathit{\boldsymbol{\tilde h}}_j^{\text{T}} \otimes {\mathit{\boldsymbol{I}}_n}} \right)\mathit{\boldsymbol{B}}{\mathit{\boldsymbol{Q}}_{{e_a}}}{\mathit{\boldsymbol{B}}^{\text{T}}}\left( {{{\mathit{\boldsymbol{\tilde h}}}_i} \otimes {\mathit{\boldsymbol{I}}_n}} \right)} \right)} \right]$，G=[g_ij]，于是目标函数Ω关于X的二阶偏导数的期望为:

按照Fisher-Score算法的原理，可得计算PEIV模型参数WTLS估计的迭代公式为:

而计算PEIV模型WTLS估计的Fisher-Score方法的具体实施步骤设计如下。

1) 给定初始值

2) 对于任意i，计算S^T_Bj和G；

3) 计算未知参数估值

4) 如果$\left\| {{{\mathit{\boldsymbol{\hat X}}}^{\left( {i + 1} \right)}}-{{\mathit{\boldsymbol{\hat X}}}^{\left( i \right)}}} \right\| < \delta $，则迭代结束；否则转到步骤2)。

算法的具体计算流程如图 1所示。

Figure 1.  Flow Chart of Computing the WTLS Estimator for PEIV Model

线性回归模型式为:

式中，(x, y)为点坐标；e_x和e_y为对应的误差；ξ_i(i=1, 2)为待估的拟合参数。采用文献[1]中的实验数据进行计算分析。在上述线性模型中，系数矩阵中含常数列[1…1…1]^T，故含非随机元素。

以Think pad E430(Intel Core i5 CPU, 2.50GHz)个人笔记本电脑和Matlab 7.0软件为实验平台，设定迭代终止条件中的阈值为δ=1×10^－5，运用文献[7]中算法和本文算法可计算得到相应的拟合参数估值，其结果列于表 1。为了便于比较，将文献[1]的计算结果也列于表 1。

从表 1的结果可知，本文算法得到的拟合参数估值与文献[1]的结果更加接近。主要原因在于文献[7]中算法过早收敛，导致不能得到更加真实和可靠的结果。为了消除随机因素的影响，将文献[7]算法和Fisher-Score算法在相同条件下重复计算了200次，然后取平均结果(平均计算时间和平均迭代次数)列于表 2。

表 2的结果表明，本文算法仅6次迭代就满足收敛条件，具备快速收敛性，且低于文献[7]算法的181次。就平均计算时间而言，Fisher-Score算法的计算效率得到显著提升。

残差平方和是评价算法优劣性的重要指标，故这里给出了文献[7]中算法和本文算法的残差平方和随迭代次数的变化趋势，如图 2所示。鉴于文献[7]算法的第一次迭代的残差平方和为LS残差平方和，其值相对较大，故没有给予展示。从图 2中展示结果可知，Fisher-Score方法的初始残差平方和较小，并且在1次迭代后，变化趋于平缓；而文献[7]算法的初始残差平方和较大，且需多次迭代才趋于平缓。这一结果充分体现了本文算法在计算效率方面的优越性。

Figure 2.  Sum of Squared Residuals with Iterations for Reference [7] and Fisher-Score Algorithm in Linear Fitting

考虑两坐标系的n个公共点(X_i, Y_i, Z_i)和(X′_i, Y′_i, Z′_i), i=1, 2…n(n≥3)，其中(X_i, Y_i, Z_i)和(X′_i, Y′_i, Z′_i)分别为同一点在原坐标系和新坐标系的坐标，则3D坐标转换模型^[18]表示为：

为了对一般的WTLS估计方法进行分析，构造如下权阵：

然后模拟服从正态分布N(0, σ²P^－1)的随机误差，将其分别加入到5个点在两套坐标系中的坐标值中，得到新的观测值。显然在上述模型中，含非随机元素0和1以及不同位置含相同的随机元素如X。

该算例的实验平台同§3.1。表 3给出了由文献[7]算法和Fisher-Score算法得到的变换参数估值。表 3的结果说明由两种算法所得到的参数估值近似相等。为了消除随机因素的影响，在相同条件下，在计算机上运用上述算法计算坐标变换估值50次，然后统计其平均计算时间以及平均迭代次数，结果如表 4所示。

从表 4可以看出，文献[7]算法一次运算达到收敛需要迭代375次，而本文算法仅3次，具备快速收敛性。而两种算法50次计算的平均计算时间分别为1.150 66 s和0.151 62 s，即本文算法的计算时间仅为文献[7]算法的13%，计算效率得到了极大提升。文献[7]中算法和Fisher-Score算法的残差平方和随迭代次数的变化趋势如图 3所示，结论与线性拟合算例亦相同。

Figure 3.  Sum of Squared Residuals with Iterations for Reference [7] and Fisher-Score Algorithm in 3D Coordinate Transformation

1) 顾及系数矩阵含非随机元素和不同位置含相同随机元素的结构化特征，PEIV模型较一般的EIV模型，理论上更为严密。

2) 考虑到现有计算PEIV模型WTLS估计的方法需多次迭代，进而影响计算效率。为此，本文结合PEIV模型的特征和随机误差的统计性质，构造了仅基于未知参数的非线性模型，利用加权残差平方和最小原理，推导了相应的计算公式，并设计了WTLS估计的Fisher-Score算法。

3) 以线性拟合和三维坐标转换为例进行了数值实验，从平均迭代次数、平均计算时间和两种算法的残差平方和随迭代次数的变化图可知, Fisher-score算法迭代次数较少，计算效率得到了较大提升。因此，实际应用中应利用迭代次数少的Fisher-Score算法。