宝塔服务器面板，一键全能部署及管理，送你10850元礼包，点我领取

0 解与算符

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

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

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

1 分辨率和协方差

1.1 数据分辨率矩阵

假设：已知一个求解反问题$\mathbf{G}\mathbf{m}=\mathbf{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带入$\mathbf{G}\mathbf{m}=\mathbf{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为数据分辨率矩阵，描述了预测值与观测数据匹配程度有多好
讨论：
$\mathbf{N}=\mathbf{I}$,此时dpre=dobsd^{pre}=d^{obs}dpre=dobs,那么预测误差为0；
$\mathbf{N}\ne \mathbf{I}$, 那么预测误差不为0，此时矩阵结构如下图：

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

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

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

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
式中，$\mathbf{R}$为分辨率矩阵。
如果$\mathbf{R}=\mathbf{I}$，模型参数都是唯一确定的。否则，模型参数的估计实际是真实模型参数的加权平均（如下图）。

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

例子:拉普拉斯变换（Laplace transform）的分辨率：
dc)=∫0∞exp⁡−cz)mz)dz→di=∑j=1Mexp⁡−cizj)mjdc)=\int_{0}^{\infty} \exp -c z) mz) d z \rightarrow d_{i}=\sum_{j=1}^{M} \exp \left-c_{i} z_{j}\right) m_{j} dc)=∫0∞​exp−cz)mz)dz→di​=j=1∑M​exp−ci​zj​)mj​
上式中：
did_idi​, 模型参数mjm_jmj​的加权平均，权重随深度$z$衰减。指数的衰减速率由常数cic_ici​控制，以致较小的$i$对应于深度范围的平均，而较大的$i$对应于较浅深度范围的平均。其浅部模型参数被较好地求解了。如下图，大值（红色）仅出现在主对角线（虚线）的接近顶部（较小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}} [covu​m]=σ−2G−g[covd]G−gT=G−gG−gT
将数据协方差矩阵进行归一化，形成一个单位数据协方差[cov⁡ud]\left[\operatorname{cov}_{u} \mathbf{d}\right][covu​d]，则上式可写为：
[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}} [covu​m]=G−g[covu​d]G−gT

例子：考虑使用一条直线拟合$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] [covu​m]=N∑zi2​−∑zi​)21​[∑zi2​−∑zi​​−∑zi​N​]
下图为使用最小二乘法对不相关的数据（黑圈）拟合为一条直线（红色），此数据具有均匀方差（垂直棒，$1\sigma$置信界限），图中A，B的方程相同，但A图数据没有很好的分开，而B图数据较为分散，因此预测数据的方差（蓝色曲线，$1\sigma$界限）图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} spreadN)=∥N−I∥22​=∑i=1N​∑j=1N​[Nij​−δij​]2spreadR)=∥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 [covu​m])=∥∥∥​[varu​m]1/2∥∥∥​22​=i=1∑M​[covu​m]ii​
式中，平方根是一个元素一个元素施加的。注意，协方差大小的度量没有考虑单位协方差矩阵中偏离主对角线元素的大小。

2.2 展布函数的应用举例

例子，考虑用展布函数去求解超定反问题$\mathbf{G}\mathbf{m}=\mathbf{d}$。
分析$\mathbf{N}$的第$k$行的展布Jk\mathbf{J_k}Jk​:
Jk=∑i=1NNki−δki)2=∑i=1NNki2−2∑i=1NNkiδki+∑i=1Nδki2J_{k}=\sum_{i=1}^{N}\leftN_{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∑N​Nki2​−2i=1∑N​Nki​δ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∑M​Gkj​Gji−g​][p=1∑M​Gkp​Gpi−g​]]​=∂Gqr−g​∂​[i=1∑N​j=1∑M​p=1∑M​Gji−g​Gpi−g​Gkj​Gkp​]=2i=1∑N​j=1∑M​p=1∑M​δjq​δir​Gpi−g​Gkj​Gkp​=2p=1∑M​Gpr−g​Gkq​Gkp​​
第二项：
−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∑N​j=1∑M​Gkj​Gji−g​δki​=−2i=1∑N​j=1∑M​Gkj​δ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∑M​Gkq​Gkp​Gpr−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)α1​spreadN)+α2​spreadR)+α3​size[covu​m])
式中$\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​[covu​d]]=[α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}[covu​d]=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}[covu​d]=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=1Mwi,j)[Rij−δij]2=∑i=1M∑j=1Mwi,j)Rij2\operatorname{spread}\mathbf{R})=\sum_{i=1}^{M} \sum_{j=1}^{M} wi, j)\left[R_{i j}-\delta_{i j}\right]^{2}=\sum_{i=1}^{M} \sum_{j=1}^{M} wi, j) R_{i j}^{2} spreadR)=i=1∑M​j=1∑M​wi,j)[Rij​−δij​]2=i=1∑M​j=1∑M​wi,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∑M​Rij​=[1]i​
这个约束保证了分辨率矩阵对角元素的大小是有限的，并且它的行是施加于真实模型参数的单位平均函数。

将分辨率矩阵的某个行的展布写为JkJ_kJk​，并代入分辨率矩阵的表达式，我们有：
Jk=∑l=1Mwl,k)RklRkl=∑l=1Mwl,k)[∑i=1NGki−gGil][∑j=1NGkj−gGjl]=∑i=1N∑j=1NGki−gGkj−g∑l=1Mwl,k)GilGjl=∑i=1N∑j=1NGki−gGkj−gSijk)\begin{aligned} J_{k} &=\sum_{l=1}^{M} wl, k) R_{k l} R_{k l} \\ &=\sum_{l=1}^{M} wl, 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} wl, 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∑M​wl,k)Rkl​Rkl​=l=1∑M​wl,k)[i=1∑N​Gki−g​Gil​][j=1∑N​Gkj−g​Gjl​]=i=1∑N​j=1∑N​Gki−g​Gkj−g​l=1∑M​wl,k)Gil​Gjl​=i=1∑N​j=1∑N​Gki−g​Gkj−g​Sijk)​​
其中，Sijk)S_{ij}^{k)}Sijk)​定义为：
Sijk)=∑l=1Mwl,k)GilGjlS_{i j}^{k)}=\sum_{l=1}^{M} wl, k) G_{i l} G_{j l} Sijk)​=l=1∑M​wl,k)Gil​Gjl​

约束方程∑jRij=[1]i\sum_j{R_{ij}}=[1]_i∑j​Rij​=[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∑M​Rik​=k=1∑M​[j=1∑N​Gij−g​Gjk​]=j=1∑N​Gij−g​k=1∑M​Gjk​=j=1∑N​Gij−g​uj​
其中，uju_juj​定义为
uj=∑k=1MGjku_{j}=\sum_{k=1}^{M} G_{j k} uj​=k=1∑M​Gjk​

关于广义逆元素最小化JkJ_kJk​，（在给定约束条件下）可以使用拉格朗日乘数来求解。首先定义拉格朗日函数$\mathbf{\Phi }$，表达为
Φ=∑i=1N∑j=1MGki−gGkj−gSijk)+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∑N​j=1∑M​Gki−g​Gkj−g​Sijk)​+2λj=1∑N​Gkj−g​uj​
式中，$2\lambda$为拉格朗日乘数。
对关于广义逆的元素$\Phi$求导，并让结果等于0，得到：
∂Φ∂Gkp−g=2∑i=1NSpik)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∑N​Spik)​Gki−g​+2λup​=0
该方程必须与原来的约束方程同时求解，取出G−g\mathbf{G^{-g}}G−g的第$k$行作为列向量gk)\mathbf{g^{k)}}gk)，并且将Sijk)\mathbf{S_{ij}^{k)}}Sijk)​，我们可以将这些方程写为：
[Sk)uuT0][gk)λ]=[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] [Sk)uT​u0​][gk)λ​]=[O1​]
将设上式中对称矩阵的逆矩阵是存在的，在此将其分割为一个$N\times N$的对称方阵$\mathbf{A}$、向量$\mathbf{b}$和标量$c$。使用上述假设，将对称矩阵左乘其逆矩阵产生单位矩阵：
[AbbTc][Sk)uuT0]=[IO01]=[ASk)+buTAubTSk)+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] [AbT​bc​][Sk)uT​u0​]=[I0​O1​]=[ASk)+buTbTSk)+cuT​AubTu​]
对于未知子矩阵$\mathbf{A}$、向量$\mathbf{b}$和标量$c$可以让子矩阵相等来确定：
ASk)+buT=I\mathbf{A S}^{\mathbf{k})}+\mathbf{b u}^{\mathrm{T}}=\mathbf{I}ASk)+buT=I ==> A=[Sk)]−1[I−buT]\mathbf{A}=\left[\mathbf{S}^{\mathbf{k})}\right]^{-1}\left[\mathbf{I}-\mathbf{b u}^{\mathbf{T}}\right]A=[Sk)]−1[I−buT]
$\mathbf{A}\mathbf{u}=\mathbf{O}$ ==> [Sk)]−1u=buTSk)u\left[\mathbf{S}^{\mathrm{k})}\right]^{-1} \mathbf{u}=\mathbf{b u}^{\mathrm{T}} \mathbf{S}^{\mathrm{k})} \mathbf{u}[Sk)]−1u=buTSk)u 和 b=[Sk)]−1uuT[Sk)]−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[Sk)]−1u[Sk)]−1u​
bTSk)+cuT=0\mathbf{b}^{\mathrm{T}} \mathbf{S}^{\mathrm{k})}+c \mathbf{u}^{\mathrm{T}}=0bTSk)+cuT=0 ==> c=−1uT[Sk)]−1uc=\frac{-1}{\mathbf{u}^{\mathrm{T}}\left[\mathbf{S}^{\mathrm{k})}\right]^{-1} \mathbf{u}}c=uT[Sk)]−1u−1​