0 解与算符

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

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

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

1 分辨率和协方差

1.1 数据分辨率矩阵

假设：已知一个求解反问题Gm=d\mathbf{Gm=d}Gm=d的广义逆，并求解出模型参数估计mest=G−gd\mathbf{m^{est}=G^{-g}d}mest=G−gd
那么，这个模型参数的估计值mest\mathbf{m^{est}}mest对实测数据dobs\mathbf{d^{obs}}dobs的拟合有多好呢？（即预测数据有多接近观测数据？）
通过将mest\mathbf{m^{est}}mest带入Gm=d\mathbf{Gm=d}Gm=d中，可得到：
dpre =Gmest =G[G−gdobs ]=[GG−g]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 }} dpre =Gmest =G[G−gdobs ]=[GG−g]dobs =Ndobs
上式中：
上表pre,obs分别表示预测和观测；
N=GG−g\mathbf{N=GG^{-g}}N=GG−g为数据分辨率矩阵，描述了预测值与观测数据匹配程度有多好
讨论：
N=I\mathbf{N=I}N=I,此时dpre=dobsd^{pre}=d^{obs}dpre=dobs,那么预测误差为0；
N≠I\mathbf{N \ne I}N=I, 那么预测误差不为0，此时矩阵结构如下图：

数据分辨率矩阵N的行，描述了相邻实测数据可以多大程度上独立地预测或求解。如在第iii行，这一行的元素的值为：
[…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.[…0000.10.80.1000….]
那么，第iii个数据：
dipre=∑j=1NNijdjobs=0.1di−1obs+0.8diobs+0.1di+1obsd_{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}} dipre=j=1∑NNijdjobs=0.1di−1obs+0.8diobs+0.1di+1obs
可以看出，dipred_{i}^{\mathrm{pre}}dipre是三个相邻观测数据的加权平均。上图中，如果行元素有单个尖锐的极大值，且它以主对角线为中心，那么数据可以被很好的恢复。若，图中的峰值非常宽，那么数据恢复较差。

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

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

1.2 模型分辨率矩阵

数据分辨率描述了数据是否可以被独立地预测或拟合。对于模型参数要探索的问题是，模型参数的一个特定估计mestm^{est}mest有多接近真实的模型参数？
假设：存在一个真实的模型参数mtrue\mathbf{m^{true}}mtrue,则Gmtrue=dobs\mathbf{Gm^{true}=d^{obs}}Gmtrue=dobs,
那么，模型参数的一个特定估计mestm^{est}mest有多接近真实的模型参数？
将Gmtrue=dobs\mathbf{Gm^{true}=d^{obs}}Gmtrue=dobs，代入mest=G−gdobs\mathbf{m^{est}}=\mathbf{G^{-g}d^{obs}}mest=G−gdobs，得到
mest=G−gdobs=G−g[Gmtrue]=[G−gG]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}} mest=G−gdobs=G−g[Gmtrue]=[G−gG]mtrue=Rmtrue
式中，R\mathbf{R}R为分辨率矩阵。
如果R=I\mathbf{R=I}R=I，模型参数都是唯一确定的。否则，模型参数的估计实际是真实模型参数的加权平均（如下图）。

同数据分辨率矩阵一样，它是数据核和先验信息的函数。可作为实验设计的另一个重要工具。

例子:拉普拉斯变换（Laplace transform）的分辨率：
d(c)=∫0∞exp⁡(−cz)m(z)dz→di=∑j=1Mexp⁡(−cizj)mjd(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} d(c)=∫0∞exp(−cz)m(z)dz→di=j=1∑Mexp(−cizj)mj
上式中：
did_idi, 模型参数mjm_jmj的加权平均，权重随深度zzz衰减。指数的衰减速率由常数cic_ici控制，以致较小的iii对应于深度范围的平均，而较大的iii对应于较浅深度范围的平均。其浅部模型参数被较好地求解了。如下图，大值（红色）仅出现在主对角线（虚线）的接近顶部（较小z值），表明在较大z值处分辨率较差。

1.3 单位协方差矩阵

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

例子：考虑使用一条直线拟合(z,d)(z,d)(z,d)数据的问题。截距m1m_1m1和斜率m2m_2m2的单位协方差矩阵为：
[cov⁡um]=1N∑zi2−(∑zi)2[∑zi2−∑zi−∑ziN]\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] [covum]=N∑zi2−(∑zi)21[∑zi2−∑zi−∑ziN]
下图为使用最小二乘法对不相关的数据（黑圈）拟合为一条直线（红色），此数据具有均匀方差（垂直棒，1σ1\sigma1σ置信界限），图中A，B的方程相同，但A图数据没有很好的分开，而B图数据较为分散，因此预测数据的方差（蓝色曲线，1σ1\sigma1σ界限）图A较大，而图B较小。

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

2.1 狄利克雷展布函数

对于分辨率和协方差矩阵，当它们为单位矩阵时为最佳，那么可以通过对分辨率矩阵中的非对角元素的大小和展布情况来作为分辨率的度量：
spread⁡(N)=∥N−I∥22=∑i=1N∑j=1N[Nij−δij]2spread⁡(R)=∥R−I∥22=∑i=1M∑j=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} spread(N)=∥N−I∥22=∑i=1N∑j=1N[Nij−δij]2spread(R)=∥R−I∥22=∑i=1M∑j=1M[Rij−δij]2
上述这种基于分辨率矩阵与单位矩阵之间差值的L2L_2L2范数定义的度量，可称之为狄利克雷展布函数（Dirichlet spread functions）。

单位协方差的度量为：
size ([cov⁡um])=∥[var⁡um]1/2∥22=∑i=1M[cov⁡um]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} size ([covum])=∥∥∥[varum]1/2∥∥∥22=i=1∑M[covum]ii
式中，平方根是一个元素一个元素施加的。注意，协方差大小的度量没有考虑单位协方差矩阵中偏离主对角线元素的大小。

2.2 展布函数的应用举例

例子，考虑用展布函数去求解超定反问题Gm=d\mathbf{Gm=d}Gm=d。
分析N\mathbf{N}N的第kkk行的展布Jk\mathbf{J_k}Jk:
Jk=∑i=1N(Nki−δki)2=∑i=1NNki2−2∑i=1NNkiδki+∑i=1Nδki2J_{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=i=1∑N(Nki−δki)2=i=1∑NNki2−2i=1∑NNkiδki+i=1∑Nδki2
对Jk\mathbf{J_k}Jk最小化，∂Jk∂Gqr−g=0\frac{\partial J_{k}}{\partial G_{q r}^{-g}}=0∂Gqr−g∂Jk=0
对JkJ_kJk的三项分别求导，
第一项为：
∂∂Gqr−g[∑i=1N[∑j=1MGkjGji−g][∑p=1MGkpGpi−g]]=∂∂Gqr−g[∑i=1N∑j=1M∑p=1MGji−gGpi−gGkjGkp]=2∑i=1N∑j=1M∑p=1MδjqδirGpi−gGkjGkp=2∑p=1MGpr−gGkqGkp\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} ∂Gqr−g∂[i=1∑N[j=1∑MGkjGji−g][p=1∑MGkpGpi−g]]=∂Gqr−g∂[i=1∑Nj=1∑Mp=1∑MGji−gGpi−gGkjGkp]=2i=1∑Nj=1∑Mp=1∑MδjqδirGpi−gGkjGkp=2p=1∑MGpr−gGkqGkp
第二项：
−2∂∂Gqr−g∑i=1N∑j=1MGkjGji−gδki=−2∑i=1N∑j=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} −2∂Gqr−g∂i=1∑Nj=1∑MGkjGji−gδki=−2i=1∑Nj=1∑MGkjδjqδirδki=−2Gkqδkr
第三项是0，整理可得：
∑p=1MGkqGkpGpr−g=Gkqδkr\sum_{p=1}^{M} G_{k q} G_{k p} G_{p r}^{-g}=G_{k q} \delta_{k r} p=1∑MGkqGkpGpr−g=Gkqδkr
对所有的k求和，并转化为矩阵形式，即得到：
GTGG−g=GT\mathbf{G}^{\mathrm{T}} \mathbf{G} \mathbf{G}^{-\mathrm{g}}=\mathbf{G}^{\mathrm{T}} GTGG−g=GT
由于GTG\mathbf{G^TG}GTG是方阵，两边同时乘以它的逆，即可得到广义逆
G−g=[GTG]−1GT\mathbf{G}^{-\mathrm{g}}=\left[\mathbf{G}^{\mathrm{T}} \mathbf{G}\right]^{-1} \mathbf{G}^{\mathrm{T}} G−g=[GTG]−1GT
最小二乘广义逆可以被看作是某种逆，它既可以被认为是最小化预测误差，也可以被认为是最小化数据分辨率矩阵的狄利克雷展布。
同理，对于纯欠定情况，最小长度广义逆为：
G−g=GT[GGT]−1\mathbf{G}^{-\mathbf{g}}=\mathbf{G}^{\mathrm{T}}\left[\mathbf{G} \mathbf{G}^{\mathrm{T}}\right]^{-1} G−g=GT[GGT]−1
既可以看作是最小化解的长度，也可以看作是最小化模型分辨率矩阵展布。

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

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

3 旁瓣和巴克斯-吉尔伯特展布函数

用狄利克雷展布函数来计算广义逆，通常会产生旁瓣（sidelobes），即在分辨率矩阵中远离主对角线处存在大幅值区域。 为此希望找到一种不产生旁瓣的广义逆，即使代价是加宽主对角线附近非零元素带的宽度，因为具有这种分辨率矩阵的解是可解释的，并被解释为物理上相邻的模型参数的局部平均。

巴克斯-吉尔伯特展布函数（Backus-Gilbert spreading function; Backus and Gilbert,1967,1968）,由下式给出：
spread⁡(R)=∑i=1M∑j=1Mw(i,j)[Rij−δij]2=∑i=1M∑j=1Mw(i,j)Rij2\operatorname{spread}(\mathbf{R})=\sum_{i=1}^{M} \sum_{j=1}^{M} w(i, j)\left[R_{i j}-\delta_{i j}\right]^{2}=\sum_{i=1}^{M} \sum_{j=1}^{M} w(i, j) R_{i j}^{2} spread(R)=i=1∑Mj=1∑Mw(i,j)[Rij−δij]2=i=1∑Mj=1∑Mw(i,j)Rij2
对于数据分辨率展布，相似表达式同样成立。
有了新的展布函数，那么就可以推导新的广义逆。新的广义逆的旁瓣将比基于狄利克雷展布函数的广义逆更小，同时，它们在以用其他标准判定时有时会较差。

3.1 欠定问题的巴克斯-吉尔伯特广义逆

我们谋求特定的广义逆G−g\mathbf{G^{-g}}G−g，它最小化了模型分辨率的巴克斯-吉尔伯特展布。由于没有给予模型分辨率的对角元素任何权重，我们需要得到的模型分辨率矩阵满足如下方程：
∑j=1MRij=[1]i\sum_{j=1}^{M} R_{i j}=[1]_{i} j=1∑MRij=[1]i
这个约束保证了分辨率矩阵对角元素的大小是有限的，并且它的行是施加于真实模型参数的单位平均函数。

将分辨率矩阵的某个行的展布写为JkJ_kJk，并代入分辨率矩阵的表达式，我们有：
Jk=∑l=1Mw(l,k)RklRkl=∑l=1Mw(l,k)[∑i=1NGki−gGil][∑j=1NGkj−gGjl]=∑i=1N∑j=1NGki−gGkj−g∑l=1Mw(l,k)GilGjl=∑i=1N∑j=1NGki−gGkj−gSij(k)\begin{aligned} J_{k} &=\sum_{l=1}^{M} w(l, k) R_{k l} R_{k l} \\ &=\sum_{l=1}^{M} w(l, k)\left[\sum_{i=1}^{N} G_{k i}^{-g} G_{i l}\right]\left[\sum_{j=1}^{N} G_{k j}^{-g} G_{j l}\right] \\ &=\sum_{i=1}^{N} \sum_{j=1}^{N} G_{k i}^{-g} G_{k j}^{-g} \sum_{l=1}^{M} w(l, k) G_{i l} G_{j l} \\ &=\sum_{i=1}^{N} \sum_{j=1}^{N} G_{k i}^{-g} G_{k j}^{-g} S_{i j}^{(k)} \end{aligned} Jk=l=1∑Mw(l,k)RklRkl=l=1∑Mw(l,k)[i=1∑NGki−gGil][j=1∑NGkj−gGjl]=i=1∑Nj=1∑NGki−gGkj−gl=1∑Mw(l,k)GilGjl=i=1∑Nj=1∑NGki−gGkj−gSij(k)
其中，Sij(k)S_{ij}^{(k)}Sij(k)定义为：
Sij(k)=∑l=1Mw(l,k)GilGjlS_{i j}^{(k)}=\sum_{l=1}^{M} w(l, k) G_{i l} G_{j l} Sij(k)=l=1∑Mw(l,k)GilGjl

约束方程∑jRij=[1]i\sum_j{R_{ij}}=[1]_i∑jRij=[1]i的左边也可以写为广义逆的表达式，即
∑k=1MRik=∑k=1M[∑j=1NGij−gGjk]=∑j=1NGij−g∑k=1MGjk=∑j=1NGij−guj\sum_{k=1}^{M} R_{i k}=\sum_{k=1}^{M}\left[\sum_{j=1}^{N} G_{i j}^{-g} G_{j k}\right]=\sum_{j=1}^{N} G_{i j}^{-g} \sum_{k=1}^{M} G_{j k}=\sum_{j=1}^{N} G_{i j}^{-g} u_{j} k=1∑MRik=k=1∑M[j=1∑NGij−gGjk]=j=1∑NGij−gk=1∑MGjk=j=1∑NGij−guj
其中，uju_juj定义为
uj=∑k=1MGjku_{j}=\sum_{k=1}^{M} G_{j k} uj=k=1∑MGjk

关于广义逆元素最小化JkJ_kJk，（在给定约束条件下）可以使用拉格朗日乘数来求解。首先定义拉格朗日函数Φ\mathbf{\Phi}Φ，表达为
Φ=∑i=1N∑j=1MGki−gGkj−gSij(k)+2λ∑j=1NGkj−guj\Phi=\sum_{i=1}^{N} \sum_{j=1}^{M} G_{k i}^{-g} G_{k j}^{-g} S_{i j}^{(k)}+2 \lambda \sum_{j=1}^{N} G_{k j}^{-g} u_{j} Φ=i=1∑Nj=1∑MGki−gGkj−gSij(k)+2λj=1∑NGkj−guj
式中，2λ2\lambda2λ为拉格朗日乘数。
对关于广义逆的元素Φ\PhiΦ求导，并让结果等于0，得到：
∂Φ∂Gkp−g=2∑i=1NSpi(k)Gki−g+2λup=0\frac{\partial \Phi}{\partial G_{k p}^{-g}}=2 \sum_{i=1}^{N} S_{p i}^{(k)} G_{k i}^{-g}+2 \lambda u_{p}=0 ∂Gkp−g∂Φ=2i=1∑NSpi(k)Gki−g+2λup=0
该方程必须与原来的约束方程同时求解，取出G−g\mathbf{G^{-g}}G−g的第kkk行作为列向量g(k)\mathbf{g^{(k)}}g(k)，并且将Sij(k)\mathbf{S_{ij}^{(k)}}Sij(k)，我们可以将这些方程写为：
[S(k)uuT0][g(k)λ]=[O1]\left[\begin{array}{ll} \mathbf{S}^{(k)} & \mathbf{u} \\ \mathbf{u}^{\mathrm{T}} & 0 \end{array}\right]\left[\begin{array}{l} \mathbf{g}^{(\mathrm{k})} \\ \lambda \end{array}\right]=\left[\begin{array}{l} \mathbf{O} \\ 1 \end{array}\right] [S(k)uTu0][g(k)λ]=[O1]
将设上式中对称矩阵的逆矩阵是存在的，在此将其分割为一个N×NN \times NN×N的对称方阵A\mathbf{A}A、向量b\mathbf{b}b和标量ccc。使用上述假设，将对称矩阵左乘其逆矩阵产生单位矩阵：
[AbbTc][S(k)uuT0]=[IO01]=[AS(k)+buTAubTS(k)+cuTbTu]\left[\begin{array}{ll} \mathbf{A} & \mathbf{b} \\ \mathbf{b}^{\mathrm{T}} & c \end{array}\right]\left[\begin{array}{ll} \mathbf{S}^{(\mathrm{k})} & \mathbf{u} \\ \mathbf{u}^{\mathrm{T}} & 0 \end{array}\right]=\left[\begin{array}{ll} \mathbf{I} & \mathbf{O} \\ 0 & 1 \end{array}\right]=\left[\begin{array}{ll} \mathbf{A} \mathbf{S}^{(\mathrm{k})}+\mathbf{b} \mathbf{u}^{\mathrm{T}} & \mathbf{A} \mathbf{u} \\ \mathbf{b}^{\mathrm{T}} \mathbf{S}^{(\mathrm{k})}+\mathbf{c u}^{\mathrm{T}} & \mathbf{b}^{\mathrm{T}} \mathbf{u} \end{array}\right] [AbTbc][S(k)uTu0]=[I0O1]=[AS(k)+buTbTS(k)+cuTAubTu]
对于未知子矩阵A\mathbf{A}A、向量b\mathbf{b}b和标量ccc可以让子矩阵相等来确定：
AS(k)+buT=I\mathbf{A S}^{(\mathbf{k})}+\mathbf{b u}^{\mathrm{T}}=\mathbf{I}AS(k)+buT=I ==> A=[S(k)]−1[I−buT]\mathbf{A}=\left[\mathbf{S}^{(\mathbf{k})}\right]^{-1}\left[\mathbf{I}-\mathbf{b u}^{\mathbf{T}}\right]A=[S(k)]−1[I−buT]
Au=O\mathbf{A} \mathbf{u}=\mathbf{O}Au=O ==> [S(k)]−1u=buTS(k)u\left[\mathbf{S}^{(\mathrm{k})}\right]^{-1} \mathbf{u}=\mathbf{b u}^{\mathrm{T}} \mathbf{S}^{(\mathrm{k})} \mathbf{u}[S(k)]−1u=buTS(k)u 和 b=[S(k)]−1uuT[S(k)]−1u\mathbf{b}=\frac{\left[\mathbf{S}^{(\mathbf{k})}\right]^{-1} \mathbf{u}}{\mathbf{u}^{\mathrm{T}}\left[\mathbf{S}^{(\mathbf{k})}\right]^{-1} \mathbf{u}}b=uT[S(k)]−1u[S(k)]−1u
bTS(k)+cuT=0\mathbf{b}^{\mathrm{T}} \mathbf{S}^{(\mathrm{k})}+c \mathbf{u}^{\mathrm{T}}=0bTS(k)+cuT=0 ==> c=−1uT[S(k)]−1uc=\frac{-1}{\mathbf{u}^{\mathrm{T}}\left[\mathbf{S}^{(\mathrm{k})}\right]^{-1} \mathbf{u}}c=uT[S(k)]−1u−1

将方程左乘对称逆矩阵得到g(k)=b\mathbf{g^{(k)}=b}g(k)=b和λ=c\mathbf{\lambda}=cλ=c，由此，广义逆写为求和形式后为
Gkl−g=∑i=1Nui[(S(k))−1]il∑i=1N∑j=1Nui[(S(k))−1]ijujG_{k l}^{-g}=\frac{\sum_{i=1}^{N} u_{i}\left[\left(\mathbf{S}^{(k)}\right)^{-1}\right]_{i l}}{\sum_{i=1}^{N} \sum_{j=1}^{N} u_{i}\left[\left(\mathbf{S}^{(k)}\right)^{-1}\right]_{i j} u_{j}} Gkl−g=∑i=1N∑j=1Nui[(S(k))−1]ijuj∑i=1Nui[(S(k))−1]il
这个广义逆的巴克斯-吉尔伯特公式可类比于最小长度解。

例子，下图为拉普拉斯变换问题的狄利克雷解和巴克斯-吉尔伯特解，图a为真实模型（红色）包含了一系列的尖刺。使用巴克斯-吉尔伯特展布函数估计的模型（蓝色）非常光滑，并且光滑的宽度随zzz的增加而增大。图b为对应的模型分辨率矩阵R\mathbf{R}R。图c和d与上述相同，但对应狄利克雷展布函数。注意，巴克斯-吉尔伯特分辨率矩阵具有较低强度的旁瓣，却有较宽的中心条带。

巴克斯-吉尔伯特解使两个解中较光滑的那个，且所对应的模型分辨率矩阵由一个沿主对角线的条带构成。狄利克雷解拥有更多的细节，但也拥有更多的假象（如在z≈3z \approx 3z≈3处的负值）。它们与对应的模型分辨率矩阵中的大振幅旁瓣有关联。

3.2 引入协方差大小

αspread⁡(R)+(1−α)size⁡([cov⁡um])=α∑i=1M∑j=1Mw(i,j)Rij2+(1−α)∑i=1M[cov⁡um]ii\alpha \operatorname{spread}(\mathbf{R})+(1-\alpha) \operatorname{size}\left(\left[\operatorname{cov}_{\mathbf{u}} \mathbf{m}\right]\right)=\alpha \sum_{i=1}^{M} \sum_{j=1}^{M} w(i, j) R_{i j}^{2}+(1-\alpha) \sum_{i=1}^{M}\left[\operatorname{cov}_{u} \mathbf{m}\right]_{i i} αspread(R)+(1−α)size([covum])=αi=1∑Mj=1∑Mw(i,j)Rij2+(1−α)i=1∑M[covum]ii
那么，第kkk行的度量Jk′J_{k}^{\prime}Jk′是
Jk′=α∑l=1Mw(k,l)Rkl2+(1−α)[covum]kk=α∑i=1N∑j=1NGki−gGkj−g[Sij]k+(1−α)∑i=1N∑j=1NGki−gGkj−g[covud]ij=∑i=1N∑j=1NGki−gGkj−gSij′(k)\begin{aligned} J_{k}^{\prime} &=\alpha \sum_{l=1}^{M} w(k, l) R_{k l}^{2}+(1-\alpha)\left[\mathrm{cov}_{\mathrm{u}} \mathbf{m}\right]_{k k} \\ &=\alpha \sum_{i=1}^{N} \sum_{j=1}^{N} G_{k i}^{-g} G_{k j}^{-g}\left[S_{i j}\right]_{k}+(1-\alpha) \sum_{i=1}^{N} \sum_{j=1}^{N} G_{k i}^{-g} G_{k j}^{-g}\left[\mathrm{cov}_{\mathrm{u}} \mathbf{d}\right]_{i j} \\ &=\sum_{i=1}^{N} \sum_{j=1}^{N} G_{k i}^{-g} G_{k j}^{-g} S_{i j}^{\prime(k)} \end{aligned} Jk′=αl=1∑Mw(k,l)Rkl2+(1−α)[covum]kk=αi=1∑Nj=1∑NGki−gGkj−g[Sij]k+(1−α)i=1∑Nj=1∑NGki−gGkj−g[covud]ij=i=1∑Nj=1∑NGki−gGkj−gSij′(k)
其中，量Sij′(k)S_{ij}^{\prime(k)}Sij′(k)由如下方程定义：
Sij′(k)=αSij(k)+(1−α)[cov⁡ud]ijS_{i j}^{\prime(k)}=\alpha S_{i j}^{(k)}+(1-\alpha)\left[\operatorname{cov}_{\mathrm{u}} \mathbf{d}\right]_{i j} Sij′(k)=αSij(k)+(1−α)[covud]ij
最后，广义逆仅需要将前面结果中的Sij(k)S_{ij}^{(k)}Sij(k)替换Sij′(k)S_{ij}^{\prime(k)}Sij′(k):
Gkl−g=∑i=1Nui[(S′(k))−1]il∑i=1N∑j=1Nui[(S′(k))−1]ijujG_{k l}^{-g}=\frac{\sum_{i=1}^{N} u_{i}\left[\left(\mathbf{S}^{\prime(k)}\right)^{-1}\right]_{i l}}{\sum_{i=1}^{N} \sum_{j=1}^{N} u_{i}\left[\left(\mathbf{S}^{\prime(k)}\right)^{-1}\right]_{i j} u_{j}} Gkl−g=∑i=1N∑j=1Nui[(S′(k))−1]ijuj∑i=1Nui[(S′(k))−1]il
这个广义逆是类比于阻尼最小长度解的巴克斯-吉尔伯特广义逆。

4 分辨率和协方差的折中

如医学层析成像问题中的X射线不透明度（如下图），试图确定一系列的模型参数。

如果离散化模型的网络非常小（图B），那么X射线将难以把每个方格都采样到，这个问题是欠定的。如果我们要确定每一个方格的不透明度，因为很少有方格被多条X射线穿过，以致于几乎不会发生误差被平均掉的情况，那么不透明度的估计将具有相当大的方差。但由于格子较小，那么可以检测到很小的特诊。若方格较大（图A），大部分方格可以被多条X射线穿过，可以平均掉误差（噪声），大方差可以降低。但由于格子较大，小的特征不能被检测到，因此X射线不透明度的分辨率变得较差。

上如的例子，描述了模型分辨率与方差大小之间的一种重要折中。一个量的降低只能以另一个量的增加为代价。通过折中思想，最小化分辨率展布和方差大小的加权求和：
αspread⁡(R)+(1−α)size⁡([cov⁡um])\alpha \operatorname{spread}(\mathbf{R})+(1-\alpha) \operatorname{size}\left(\left[\operatorname{cov}_{\mathbf{u}} \mathbf{m}\right]\right) αspread(R)+(1−α)size([covum])
若α\alphaα设定在1附近，那么广义逆的模型分辨率矩阵将拥有小的展布，但模型参数将拥有大方差。若α\alphaα设定为接近0，那么模型参数将具有较小的方差，但模型分辨率矩阵将拥有较大的展布。

线性高斯反问题--广义逆相关推荐

线性高斯反问题的解--最大似然法
0 一组测量的均值假设,在一个实验中对一模型参数m1m_1m1在有噪声的情况下测量了N次,每次都采集到一个独立数据did_idi.那么,如果随机变量是高斯型的,那么它们的联合概率密度函数能够以方 ...
求解线性最小二乘问题的奇异值分解及广义逆法的C++实现
求解线性最小二乘问题的广义逆法的C++实现 1,功能 2,方法说明 3,函数语句与形参说明第一步,求对系数矩阵进行奇异值分解(muav函数) #include "stdlib.h" ...
求非线性方程组的最小二乘解的广义逆法C实现
求非线性方程组的最小二乘解的广义逆法 #include "math.h"#include "stdlib.h"#include "6gmiv.c&qu ...
矩阵分析与应用（一）——矩阵基础知识、广义逆
文章目录前言部分符号约定关于矩阵理论的碎碎念一些基础知识与本门课知识串讲矩阵奇异与线性无关向量空间.内积.范数行列式.特征值.迹逆.广义逆矩阵方程.向量化.Kronecker积向量 ...
机器人学中的状态估计学习笔记（二）第三章线性高斯系统的状态估计
机器人学中的状态估计学习笔记(二)第三章线性高斯系统的状态估计 3.1 离散时间的批量估计问题 3.1.1 问题定义 3.1.2 最大后验估计 3.1.3 贝叶斯推断 3.1.4 存在性.唯一性与能观 ...
广义逆高斯分布及推广(To be continued)
1.资料广义逆高斯分布(Generalized Inverse Gaussian Distribution)及修正贝塞尔函数广义逆高斯分布及其特列
Moore-Penrose广义逆：可解决MATLAB报错“矩阵接近奇异值，或者缩放错误。结果可能不准确”
上一篇博文讲到:<方程AX=b的解的讨论(特解.通解.零空间向量等概念)及其MATLAB实现>,程序中用到的是mldivide或者A\b的方法(二者相同)来解方程. 但实际上运行过程中我们 ...
MATLAB中：左右除法、逆inv、广义逆pinv的区别
以下是通过实验得出的一些结论: 左除行相等,右除列相等.只要满足此条件便可运算,且左.右除意义不相同! 逆inv()仅针对"非奇异方阵|A|≠0",使得A^(-1) A=A ...
最优化理论与方法(part2)--矩阵求逆和广义逆
学习笔记,仅供参考,有错必纠文章目录最优化理论与方法矩阵求逆和广义逆定理 1.2.3 (von Neumann引理) 定理 1.2.3' 子空间最优化理论与方法矩阵求逆和广义逆定理 1. ...

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