矢量拟合（2） - Vector Fitting算法原理

在Sanathanan–Koerner算法中：
$$\tilde{H}(s)=\frac{n(s)}{d(s)}=\frac{{\sum }_{n=0}^{\bar{n}}{a}_n{s}^n}{{\sum }_{n=0}^{\bar{n}}{b}_n{s}^n}$$
求解数值逼近时的误差项：
$${\left(e_{SK}^{(i)}\right)}^2=\frac{1}{\bar{k}}\sum _{k=1}^{\bar{k}}{\left|\frac{H_k{\sum }_{n=0}^{\bar{n}}{b}_n^{(i)}(j{\omega }_k{)}^n-{\sum }_{n=0}^{\bar{n}}{a}_n^{(i)}(j{\omega }_k{)}^n}{{\sum }_{n=0}^{\bar{n}}{b}_n^{(i-1)}(j{\omega }_k{)}^n}\right|}^2$$
会包含频率的多次项，求解的最小二乘问题的矩阵将包含 Vandermonde 块，导致恶化的条件，会增加数值的不稳定性，当条件数变大时，计算结果对输入数据的微小变化会变得非常敏感，可能导致数值解不准确甚至不收敛。
为了解决上述的问题，VF算法用部分分式替代之前的传递函数的表达式
$$\tilde{H}(s)=\frac{n(s)}{d(s)}=\frac{{\sum }_{n=0}^{\bar{n}}{a}_n{s}^n}{{\sum }_{n=0}^{\bar{n}}{b}_n{s}^n}$$
将$\tilde{H}(s)$的分子和分母项分别写作：

$$n^{(i)}=c_0^{(i)}+\sum _{n=1}^n\frac{c_n^{(i)}}{s-p_n^{(0)}}\text{\hspace{1em}(1.10)}$$

$$d^{(i)}=1+\sum _{n=1}^n\frac{d_n^{(i)}}{s-p_n^{(0)}}\text{\hspace{1em}(1.11)}$$

其中：

• $n^{(i)}$和$d^{(i)}$ 分别表示分子和分母多项式。
• $c_0^{(i)}$和$c_n^{(i)}$是分子多项式的系数。
• $d_n^{(i)}$是分母多项式的系数。
• $p_n^{(0)}$表示初始极点，与变量$s$一起出现在分式中。

这里提到的初始极点 $p_n^{(0)}\in \mathbb{C}$是一组复数，将在后续讨论它们的选取方式。可以看到，在公式 (1.11) 中，常数项已标准化为 1，这是不失一般性的处理。与 Sanathanan-Koerner 迭代中使用的幂函数基$s^n$相比，分式$\frac{1}{s-p_n^{(0)}}$的变化在 $s$增加时更加受控。如果选择合适的极点$p_n^{(0)}$，则分式形式在频率上的变化更加平缓。此特性可改善条件数，特别是在极点$p_n^{(0)}$彼此不同且分布较开时，效果更显著。
将（1.10）和（1.11）代入之前的误差表达式（1.9）：

$$(e_{SK}^{(i)}{)}^2=\frac{1}{k}\sum _{k=1}^k{\left|\frac{H_k{\sum }_{n=0}^n{b}_n^{(i)}(\text{j}{\omega }_k{)}^n-{\sum }_{n=0}^n{a}_n^{(i)}(\text{j}{\omega }_k{)}^n}{{\left({\sum }_{n=0}^n{b}_n^{(i-1)}(\text{j}{\omega }_k{)}^n\right)}^2}\right|}^2$$
可以得到改写以后的误差表达式为：

$$e_{SK}^{(i)}=\frac{1}{k}\sum _{k=1}^k{\left|\frac{H_k{d}^{(i)}(j{\omega }_k)}{d^{(i-1)}(j{\omega }_k)}-\frac{n^{(i)}(j{\omega }_k)}{d^{(i-1)}(j{\omega }_k)}\right|}^2\text{\hspace{1em}(1.12)}$$

对该公式的分解说明如下：

• $e_{SK}^{(i)}$ 表示第 $i$ 次迭代的误差项。
• $H_k$ 和 $n^{(i)}(j{\omega }_k)$ 可能表示在离散频率 $j{\omega }_k$ 下的频率相关项或传递函数。
• $d^{(i)}(j{\omega }_k)$ 和 $d^{(i-1)}(j{\omega }_k)$ 表示分母中的项。

引入了两个新量，令：

1. $w^{(i)}(s)$：
   $$w^{(i)}(s)=\frac{d^{(i)}(s)}{d^{(i-1)}(s)}\text{\hspace{1em}(1.13)}$$

2. ${\tilde{H}}^{(i)}(s)$：
   $${\tilde{H}}^{(i)}(s)=\frac{n^{(i)}(s)}{d^{(i-1)}(s)}\text{\hspace{1em}(1.14)}$$
   将其代入1.12得到：
   $$(e_{SK}^{(i)}{)}^2=\frac{1}{k}\sum _{k=1}^k{\left|H_k{w}^{(i)}(j{\omega }_k)-{\tilde{H}}^{(i)}(j{\omega }_k)\right|}^2.$$

函数 ${\tilde{H}}^{(i)}(s)$ 可以解释为通过最小化公式 (1.12) 在第 $i$ 次迭代时估计的模型传递函数。具体来说，这个传递函数由当前迭代的分子 $n^{(i)}(s)$（待确定）和前一迭代的分母 $d^{(i-1)}(s)$（已知）构成。由于误差函数的线性化，此近似方法有效地避免了分母中未知量的存在。

函数 $w^{(i)}(s)$ 可以理解为一个频率相关的权重，用来乘以给定的样本 $H_k$。

加权函数$w^{(i)}(s)$有两个主要目的：

• 提供分母 $d^{(i)}(s)$的新估计，
• 补偿通过将${\tilde{H}}^{(i)}(s)$的分母固定为前一次迭代值而引入的近似误差。实际上，权重$w^{(i)}(s)$依赖于新分母估计$d^{(i)}(s)$与前一次估计$d^{(i-1)}(s)$的比值。

将$d^{(i)}$和$d^{(i-1)}$的表达式代入权重$w^{(i)}(s)$的表达式得到：
$$w^{(i)}(s)=\frac{1+{\sum }_{n=1}^N\frac{d_n^{(i)}}{s-p_n^{(0)}}}{1+{\sum }_{n=1}^N\frac{d_n^{(i-1)}}{s-p_n^{(0)}}}\text{\hspace{1em}(1.15a)}$$

1. 分子部分：
   $$1+\sum _{n=1}^N\frac{d_n^{(i)}}{s-p_n^{(0)}}$$
   可以看作一个分式展开的形式，对应于某个分母为 ${\prod }_{n=1}^N(s-p_n^{(0)})$ 的多项式。

2. 分母部分：
   $$1+\sum _{n=1}^N\frac{d_n^{(i-1)}}{s-p_n^{(0)}}$$
   同样可以看作另一个分式展开形式，对应于某个分母为 ${\prod }_{n=1}^N(s-p_n^{(0)})$ 的多项式。

我们需要将分式的加法和乘积展开为分母和分子的一般形式，从而得到目标式的等价关系。

3. 分子部分的通分：
   $$1+\sum _{n=1}^N\frac{d_n^{(i)}}{s-p_n^{(0)}}$$
   通分后分母为 ${\prod }_{n=1}^N(s-p_n^{(0)})$，因此分子为：
   $$\prod _{n=1}^N(s-p_n^{(0)})+\sum _{n=1}^N{d}_n^{(i)}\prod _{m\ne n}(s-p_m^{(0)})$$

4. 分母部分的通分：
   $$1+\sum _{n=1}^N\frac{d_n^{(i-1)}}{s-p_n^{(0)}}$$
   通分后分母同样为 ${\prod }_{n=1}^N(s-p_n^{(0)})$，因此分子为：
   $$\prod _{n=1}^N(s-p_n^{(0)})+\sum _{n=1}^N{d}_n^{(i-1)}\prod _{m\ne n}(s-p_m^{(0)})$$

5. 将上面两式带入原始表达式：
   $$w^{(i)}(s)=\frac{{\prod }_{n=1}^N(s-p_n^{(0)})+{\sum }_{n=1}^N{d}_n^{(i)}{\prod }_{m\ne n}(s-p_m^{(0)})}{{\prod }_{n=1}^N(s-p_n^{(0)})+{\sum }_{n=1}^N{d}_n^{(i-1)}{\prod }_{m\ne n}(s-p_m^{(0)})}$$

6. 分子的多项式形式可以表示为新的形式：
   $$\prod _{n=1}^N(s-p_n^{(i)})$$
   其中 $p_n^{(i)}$ 表示多项式的根，由迭代过程中 $d_n^{(i)}$ 的值确定。

7. 分母的多项式形式可以表示为：
   $$\prod _{n=1}^N(s-p_n^{(i-1)})$$
   其中 $p_n^{(i-1)}$ 表示多项式的根，由迭代过程中 $d_n^{(i-1)}$ 的值确定。

因此，可以直接将分子和分母重新组合为目标形式：
$$w^{(i)}(s)=\frac{\frac{{\prod }_{n=1}^N(s-p_n^{(i)})}{{\prod }_{n=1}^N(s-p_n^{(0)})}}{\frac{{\prod }_{n=1}^N(s-p_n^{(i-1)})}{{\prod }_{n=1}^N(s-p_n^{(0)})}}$$

分子分母的公共项 ${\prod }_{n=1}^N(s-p_n^{(0)})$ 相互抵消，最终得到：
$$w^{(i)}(s)=\frac{{\prod }_{n=1}^N(s-p_n^{(i)})}{{\prod }_{n=1}^N(s-p_n^{(i-1)})}\text{\hspace{1em}(1.15b)}$$

这个表达式也可以写作极点和残差的形式：
$$w^{(i)}(s)=1+\sum _{n=1}^N{w}_n^{(i)}\frac{1}{s-p_n^{(i-1)}}\text{\hspace{1em}(1.15c)}$$

在上述表达式中，$w^{(i)}(s)$ 和 $w_n^{(i)}$ 表示不同的量，具体区别如下：

1. $w^{(i)}(s)$：
   $w^{(i)}(s)$ 是一个关于变量 $s$ 的函数。它表示在第 $i$ 次迭代时的整个权重函数，用于在模型拟合中调整传递函数的误差。$w^{(i)}(s)$ 是一个复变量的函数，其形式为：

   $$w^{(i)}(s)=1+\sum _{n=1}^N{w}_n^{(i)}\frac{1}{s-p_n^{(i-1)}}$$

   这个表达式表明 $w^{(i)}(s)$ 是通过每个 $w_n^{(i)}$ 和分量 $\frac{1}{s-p_n^{(i-1)}}$ 的线性组合构成的。

2. $w_n^{(i)}$：
   $w_n^{(i)}$ 是 $w^{(i)}(s)$ 中的单个权重系数，不依赖于 $s$。它表示每个分量 $\frac{1}{s-p_n^{(i-1)}}$ 的权重。换句话说，$w_n^{(i)}$ 控制了第 $i$ 次迭代中第 $n$ 个极点 $p_n^{(i-1)}$ 对权重函数 $w^{(i)}(s)$ 的影响大小。

总结来说，$w^{(i)}(s)$ 是完整的权重函数，而 $w_n^{(i)}$ 是该权重函数中各个项的权重系数。

在公式 (1.15b) 中，$p_n^{(i)}$是 $d^{(i)}(s)$的零点，因此也是${\tilde{H}}^{(i+1)}(s)$的极点。我们将 $w^{(i)}(s)$ 表达为一组新的极点 $p_n^{(i-1)}$ 的形式，这些极点会在每次迭代中发生变化，如 (1.15c) 的表达式所示。

对 ${\tilde{H}}^{(i)}(s)$ 也可以进行相同的操作，最终得到：
$${\tilde{H}}^{(i)}(s)=\frac{c_0^{(i)}+{\sum }_{n=1}^{\bar{n}}\frac{c_n^{(i)}}{s-p_n^{(0)}}}{1+{\sum }_{n=1}^{\bar{n}}\frac{d_n^{(i-1)}}{s-p_n^{(0)}}}=r_0^{(i)}+\sum _{n=1}^{\bar{n}}\frac{r_n^{(i)}}{s-p_n^{(i-1)}}.\text{\hspace{1em}(1.16)}$$

将新的权重函数的表达式和传递函数的表达式代入，可以得到优化后的误差表达式误差函数 $(e_{SK}^{(i)}{)}^2$ 表示为：

$$(e_{SK}^{(i)}{)}^2=\frac{1}{k}\sum _{k=1}^K{\left|H_k\left(1+\sum _{n=1}^n\frac{w_n^{(i)}}{\text{j}{\omega }_k-p_n^{(i-1)}}\right)-\left(r_0^{(i)}+\sum _{n=1}^n\frac{r_n^{(i)}}{\text{j}{\omega }_k-p_n^{(i-1)}}\right)\right|}^2(1.17)$$

这是矢量拟合（VF）中实际使用的误差函数，用于将模型拟合到给定的样本数据上。公式 (1.17) 和 (1.9) 的主要区别在于线性化误差如何通过迭代加权逐步收敛到公式 (1.6)。

在公式 (1.12) 中，权重 $\frac{1}{d^{(i-1)}(j{\omega }_k)}$ 被显式地应用，但这种方法会降低数值条件的稳定性。

而在公式 (1.17) 中，相同的权重通过在每次迭代中重新定位极点 $p_n^{(i-1)}$ 的方式被隐式地应用。

一旦公式 (1.17) 被最小化后，可以通过$d^{(i)}(s)$ 的零点来找到下一次迭代的更新极点 $p_n^{(i)}$，正如从公式 (1.15c) 表达式中可以看出的那样。可以证明，这些零点可以通过一个矩阵的特征值来计算：

{ p n ( i ) } = eig ( A ( i − 1 ) − b w ( c w ( i ) ) T ) , (1.18) \{p_n^{(i)}\} = \text{eig}(A^{(i-1)} - b_w (c_w^{(i)})^T), \tag{1.18} {pn(i)​}=eig(A(i−1)−bw​(cw(i)​)T),(1.18)
其中：

• A ( i − 1 ) = diag { p 1 ( i − 1 ) , … , p n ( i − 1 ) } A^{(i-1)} = \text{diag}\{p_1^{(i-1)}, \dots, p_n^{(i-1)}\} A(i−1)=diag{p1(i−1)​,…,pn(i−1)​}是一个对角矩阵，由前一次迭代的极点 $p_n^{(i-1)}$ 构成；
• $b_w$ 是一个$n\times 1$ 的全 1 向量；
• $(c_w^{(i)}{)}^T=[w_1^{(i)},\dots ,w_n^{(i)}]$ 是权重系数的转置向量。

在收敛时，$p_n^{(i-1)}\to p_n^{(i)}$，它们成为所得到模型$H^{(i)}(s)$的极点。当这种情况发生时，根据公式 (1.15c) 表达式，可以看到 $w^{(i)}\to 1$。同时，线性化误差 (1.12) 逐渐收敛到 (1.6)，达到了目标。