线性高斯反问题--广义逆

0 解与算符

对于线性反问题Gm=d\mathbf{Gm=d}的方法,长度方法基于分析解的两个属性:预测误差和解的简单程度。而广义逆的方法,则希望通过研究算符矩阵来获得更多关于反问题的属性

假设求解反问题Gm=d\mathbf{Gm = d},存在一个广义逆可求解模型估计,mest=Ggd\mathbf{m^{est}=G^{-g}d},这里的GgG^{-g}为广义逆。

通过广义逆,可以求取:
1)数据分辨率矩阵
2)模型分辨率矩阵
3)单位协方差矩阵

1 分辨率和协方差

1.1 数据分辨率矩阵

假设:已知一个求解反问题Gm=d\mathbf{Gm=d}的广义逆,并求解出模型参数估计mest=Ggd\mathbf{m^{est}=G^{-g}d}
那么,这个模型参数的估计值mest\mathbf{m^{est}}对实测数据dobs\mathbf{d^{obs}}的拟合有多好呢?(即预测数据有多接近观测数据?)
通过将mest\mathbf{m^{est}}带入Gm=d\mathbf{Gm=d}中,可得到:
dpre =Gmest =G[Ggdobs ]=[GGg]dobs =Ndobs  \mathbf{d}^{\text {pre }}=\mathbf{G m}^{\text {est }}=\mathbf{G}\left[\mathbf{G}^{-\mathbf{g}} \mathbf{d}^{\text {obs }}\right]=\left[\mathbf{G} \mathbf{G}^{-\mathbf{g}}\right] \mathbf{d}^{\text {obs }}=\mathbf{N} \mathbf{d}^{\text {obs }}
上式中:
上表pre,obs分别表示预测和观测;
N=GGg\mathbf{N=GG^{-g}}为数据分辨率矩阵,描述了预测值与观测数据匹配程度有多好
讨论:
N=I\mathbf{N=I},此时dpre=dobsd^{pre}=d^{obs},那么预测误差为0;
NI\mathbf{N \ne I}, 那么预测误差不为0,此时矩阵结构如下图:
线性高斯反问题--广义逆
数据分辨率矩阵N的行,描述了相邻实测数据可以多大程度上独立地预测或求解。如在第ii行,这一行的元素的值为:
[0000.10.80.1000.]\left[\begin{array}{lllllllll} \ldots & 0 & 0 & 0 & 0.1 & 0.8 & 0.1 & 0 & 0 & 0 \ldots . \end{array}]\right.
那么,第ii个数据:
dipre=j=1NNijdjobs=0.1di1obs+0.8diobs+0.1di+1obs d_{i}^{\mathrm{pre}}=\sum_{j=1}^{N} N_{i j} d_{j}^{\mathrm{obs}}=0.1 d_{i-1}^{\mathrm{obs}}+0.8 d_{i}^{\mathrm{obs}}+0.1 d_{i+1}^{\mathrm{obs}}
可以看出,dipred_{i}^{\mathrm{pre}}是三个相邻观测数据的加权平均。上图中,如果行元素有单个尖锐的极大值,且它以主对角线为中心,那么数据可以被很好的恢复。若,图中的峰值非常宽,那么数据恢复较差。

需要指出的是,数据分辨率矩阵不是数据的函数,它仅是数据核G\mathbf{G}和施加于反问题的先验信息的函数。。因此,可以在没有实际进行实验的情况下计算和研究数据分辨率矩阵,是实验设计过程中的有用工具。

例子:考虑对数据点(z,d)(z,d)拟合一条直线的问题。如下图,用一条直线拟合100个沿zz轴等间距分布数据的情况,大值(红色)仅出现在主对角线(虚线)的两头,表明在zz取中间值时分辨率较差

线性高斯反问题--广义逆

1.2 模型分辨率矩阵

数据分辨率描述了数据是否可以被独立地预测或拟合。对于模型参数要探索的问题是,模型参数的一个特定估计mestm^{est}有多接近真实的模型参数?
假设:存在一个真实的模型参数mtrue\mathbf{m^{true}},则Gmtrue=dobs\mathbf{Gm^{true}=d^{obs}},
那么,模型参数的一个特定估计mestm^{est}有多接近真实的模型参数?
Gmtrue=dobs\mathbf{Gm^{true}=d^{obs}},代入mest=Ggdobs\mathbf{m^{est}}=\mathbf{G^{-g}d^{obs}},得到
mest=Ggdobs=Gg[Gmtrue]=[GgG]mtrue=Rmtrue \mathbf{m}^{\mathrm{est}}=\mathbf{G}^{-\mathrm{g}} \mathbf{d}^{\mathrm{obs}}=\mathbf{G}^{-\mathrm{g}}\left[\mathbf{G} \mathbf{m}^{\mathrm{true}}\right]=\left[\mathbf{G}^{-\mathrm{g}} \mathbf{G}\right] \mathbf{m}^{\mathrm{true}}=\mathbf{R} \mathbf{m}^{\mathrm{true}}
式中,R\mathbf{R}为分辨率矩阵。
如果R=I\mathbf{R=I},模型参数都是唯一确定的。否则,模型参数的估计实际是真实模型参数的加权平均(如下图)。
线性高斯反问题--广义逆
同数据分辨率矩阵一样,它是数据核和先验信息的函数。可作为实验设计的另一个重要工具。

例子:拉普拉斯变换(Laplace transform)的分辨率:
d(c)=0exp(cz)m(z)dzdi=j=1Mexp(cizj)mj d(c)=\int_{0}^{\infty} \exp (-c z) m(z) d z \rightarrow d_{i}=\sum_{j=1}^{M} \exp \left(-c_{i} z_{j}\right) m_{j}
上式中:
did_i, 模型参数mjm_j的加权平均,权重随深度zz衰减。指数的衰减速率由常数cic_i控制,以致较小的ii对应于深度范围的平均,而较大的ii对应于较浅深度范围的平均。其浅部模型参数被较好地求解了。如下图,大值(红色)仅出现在主对角线(虚线)的接近顶部(较小z值),表明在较大z值处分辨率较差。
线性高斯反问题--广义逆

1.3 单位协方差矩阵

模型参数的协方差依赖于:数据的协方差和 误差从数据映射至模型参数的方式。用单位协方差矩阵,可以描述发生在映射过程中的误差幅度
假设,数据是不相关的,且具有均匀的方差σ2\sigma^2,那么单位协方差为:
[covum]=σ2Gg[covd]GgT=GgGgT \left[\operatorname{cov}_{u} \mathbf{m}\right]=\sigma^{-2} \mathbf{G}^{-\mathbf{g}}[\operatorname{cov} \mathbf{d}] \mathbf{G}^{-\mathbf{g} \mathbf{T}}=\mathbf{G}^{-\mathbf{g}} \mathbf{G}^{-\mathbf{g} \mathbf{T}}
将数据协方差矩阵进行归一化,形成一个单位数据协方差[covud]\left[\operatorname{cov}_{u} \mathbf{d}\right],则上式可写为:
[covum]=Gg[covud]GgT \left[\mathrm{cov}_{\mathrm{u}} \mathbf{m}\right]=\mathbf{G}^{-\mathrm{g}}\left[\operatorname{cov}_{\mathrm{u}} \mathbf{d}\right] \mathbf{G}^{-\mathrm{g} \mathrm{T}}

例子:考虑使用一条直线拟合(z,d)(z,d)数据的问题。截距m1m_1和斜率m2m_2的单位协方差矩阵为:
[covum]=1Nzi2(zi)2[zi2ziziN] \left[\operatorname{cov}_{u} \mathbf{m}\right]=\frac{1}{N \sum z_{i}^{2}-\left(\sum z_{i}\right)^{2}}\left[\begin{array}{cc} \sum z_{i}^{2} & -\sum z_{i} \\ -\sum z_{i} & N \end{array}\right]
下图为使用最小二乘法对不相关的数据(黑圈)拟合为一条直线(红色),此数据具有均匀方差(垂直棒,1σ1\sigma置信界限),图中A,B的方程相同,但A图数据没有很好的分开,而B图数据较为分散,因此预测数据的方差(蓝色曲线,1σ1\sigma界限)图A较大,而图B较小。
线性高斯反问题--广义逆

2 反映分辨率和协方差好坏的度量

2.1 狄利克雷展布函数

对于分辨率和协方差矩阵,当它们为单位矩阵时为最佳那么可以通过对分辨率矩阵中的非对角元素的大小和展布情况来作为分辨率的度量
spread(N)=NI22=i=1Nj=1N[Nijδij]2spread(R)=RI22=i=1Mj=1M[Rijδij]2 \begin{array}{l} \operatorname{spread}(\mathbf{N})=\|\mathbf{N}-\mathbf{I}\|_{2}^{2}=\sum_{i=1}^{N} \sum_{j=1}^{N}\left[N_{i j}-\delta_{i j}\right]^{2} \\ \operatorname{spread}(\mathbf{R})=\|\mathbf{R}-\mathbf{I}\|_{2}^{2}=\sum_{i=1}^{M} \sum_{j=1}^{M}\left[R_{i j}-\delta_{i j}\right]^{2} \end{array}
上述这种基于分辨率矩阵与单位矩阵之间差值的L2L_2范数定义的度量,可称之为狄利克雷展布函数(Dirichlet spread functions)。

单位协方差的度量为:
 size ([covum])=[varum]1/222=i=1M[covum]ii \text { size }\left(\left[\operatorname{cov}_{u} \mathbf{m}\right]\right)=\left\|\left[\operatorname{var}_{u} \mathbf{m}\right]^{1 / 2}\right\|_{2}^{2}=\sum_{i=1}^{M}\left[\operatorname{cov}_{u} \mathbf{m}\right]_{i i}
式中,平方根是一个元素一个元素施加的。注意,协方差大小的度量没有考虑单位协方差矩阵中偏离主对角线元素的大小。

2.2 展布函数的应用举例

例子,考虑用展布函数去求解超定反问题Gm=d\mathbf{Gm=d}
分析N\mathbf{N}的第kk行的展布Jk\mathbf{J_k}:
Jk=i=1N(Nkiδki)2=i=1NNki22i=1NNkiδki+i=1Nδki2 J_{k}=\sum_{i=1}^{N}\left(N_{k i}-\delta_{k i}\right)^{2}=\sum_{i=1}^{N} N_{k i}^{2}-2 \sum_{i=1}^{N} N_{k i} \delta_{k i}+\sum_{i=1}^{N} \delta_{k i}^{2}
Jk\mathbf{J_k}最小化,JkGqrg=0\frac{\partial J_{k}}{\partial G_{q r}^{-g}}=0
JkJ_k的三项分别求导,
第一项为:
Gqrg[i=1N[j=1MGkjGjig][p=1MGkpGpig]]=Gqrg[i=1Nj=1Mp=1MGjigGpigGkjGkp]=2i=1Nj=1Mp=1MδjqδirGpigGkjGkp=2p=1MGprgGkqGkp \begin{aligned} \frac{\partial}{\partial G_{q r}^{-g}}\left[\sum_{i=1}^{N}\left[\sum_{j=1}^{M} G_{k j} G_{j i}^{-g}\right]\left[\sum_{p=1}^{M} G_{k p} G_{p i}^{-g}\right]\right] &=\frac{\partial}{\partial G_{q r}^{-g}}\left[\sum_{i=1}^{N} \sum_{j=1}^{M} \sum_{p=1}^{M} G_{j i}^{-g} G_{p i}^{-g} G_{k j} G_{k p}\right] \\ &=2 \sum_{i=1}^{N} \sum_{j=1}^{M} \sum_{p=1}^{M} \delta_{j q} \delta_{i r} G_{p i}^{-g} G_{k j} G_{k p} \\ &=2 \sum_{p=1}^{M} G_{p r}^{-g} G_{k q} G_{k p} \end{aligned}
第二项:
2Gqrgi=1Nj=1MGkjGjigδki=2i=1Nj=1MGkjδjqδirδki=2Gkqδkr -2 \frac{\partial}{\partial G_{q r}^{-g}} \sum_{i=1}^{N} \sum_{j=1}^{M} G_{k j} G_{j i}^{-g} \delta_{k i}=-2 \sum_{i=1}^{N} \sum_{j=1}^{M} G_{k j} \delta_{j q} \delta_{i r} \delta_{k i}=-2 G_{k q} \delta_{k r}
第三项是0,整理可得:
p=1MGkqGkpGprg=Gkqδkr \sum_{p=1}^{M} G_{k q} G_{k p} G_{p r}^{-g}=G_{k q} \delta_{k r}
对所有的k求和,并转化为矩阵形式,即得到:
GTGGg=GT \mathbf{G}^{\mathrm{T}} \mathbf{G} \mathbf{G}^{-\mathrm{g}}=\mathbf{G}^{\mathrm{T}}
由于GTG\mathbf{G^TG}是方阵,两边同时乘以它的逆,即可得到广义逆
Gg=[GTG]1GT \mathbf{G}^{-\mathrm{g}}=\left[\mathbf{G}^{\mathrm{T}} \mathbf{G}\right]^{-1} \mathbf{G}^{\mathrm{T}}
最小二乘广义逆可以被看作是某种逆,它既可以被认为是最小化预测误差,也可以被认为是最小化数据分辨率矩阵的狄利克雷展布
同理,对于纯欠定情况,最小长度广义逆为:
Gg=GT[GGT]1 \mathbf{G}^{-\mathbf{g}}=\mathbf{G}^{\mathrm{T}}\left[\mathbf{G} \mathbf{G}^{\mathrm{T}}\right]^{-1}
既可以看作是最小化 解的长度,也可以看作是最小化模型分辨率矩阵展布。

2.3 一般情况的狄利克雷展布展布函数

通过折中思想(正则化),可以通过最小化狄利克雷展布函数和协方差大小的加权求和来定义一个有意义的解:
最小化:α1spread(N)+α2spread(R)+α3size([covum])\alpha_{1} \operatorname{spread}(\mathbf{N})+\alpha_{2} \operatorname{spread}(\mathbf{R})+\alpha_{3} \operatorname{size}\left(\left[\operatorname{cov}_{\mathbf{u}} \mathbf{m}\right]\right)
式中α\alpha为任意的权重因子。
通过对上式求解(略去求解过程),所得结果为一个广义逆方程:
α1[GTG]Gg+Gg[α2[GGT]+α3[covud]]=[α1+α2]GT \alpha_{1}\left[\mathbf{G}^{\mathrm{T}} \mathbf{G}\right] \mathbf{G}^{-\mathrm{g}}+\mathbf{G}^{-\mathrm{g}}\left[\alpha_{2}\left[\mathbf{G} \mathbf{G}^{\mathrm{T}}\right]+\alpha_{3}\left[\operatorname{cov}_{\mathrm{u}} \mathbf{d}\right]\right]=\left[\alpha_{1}+\alpha_{2}\right] \mathbf{G}^{\mathrm{T}}
这种形式的方程被称为西尔维斯特方程(Sylvester equation)。当指定权重因子的大小,可以求出该方程的显式解。
α1=1,α2=α3=0\alpha_1=1,\alpha_2=\alpha_3=0时,对应最小二乘解
α1=0,α2=1,α3=0\alpha_1=0,\alpha_2=1,\alpha_3=0时,对应最小长度解
α1=1,α2=0,α3=ε2\alpha_1=1,\alpha_2=0,\alpha_3=\varepsilon^2,且[covud]=I[cov_u\mathbf{d}]=\mathbf{I},那么广义逆为:
Gg=[GTG+ε2I]1GT \mathbf{G}^{-g}=\left[\mathbf{G}^{\mathrm{T}} \mathbf{G}+\varepsilon^{2} \mathbf{I}\right]^{-1} \mathbf{G}^{\mathrm{T}}
此为阻尼最小二乘的广义逆
α1=0,α2=1,α3=ε2\alpha_1=0,\alpha_2=1,\alpha_3=\varepsilon^2,且[covud]=I[cov_u\mathbf{d}]=\mathbf{I}时:
Gg=GT[GGT+ε2I]1 \mathbf{G}^{-g}=\mathbf{G}^{\mathrm{T}}\left[\mathbf{G} \mathbf{G}^{\mathrm{T}}+\varepsilon^{2} \mathbf{I}\right]^{-1}
为阻尼最小长度的广义逆