范德蒙德和Teoplitz方程组的解法

简单介绍

工程中的很多实际问题的处理，比如说图像处理的某些情况，最后往往归结为比较容易处理的Vandermonde方程组和Teoplitz方程组的求解问题，因此对于这类特殊问题的快速求解方法就显得尤为重要。

所谓Vandermonde方程组，指的是如下方程组：
Vz=bVz=bVz=b
其中V是Vandermonde矩阵， V=V(x0,…,xn)=(11⋯1x0x1⋯xn⋮⋮⋮x0nx1n⋯xnn)V = V(x_0,\dots,x_n) = \left( {\begin{array}{ccc} 1&1& \cdots &1\\ {{x_0}}&{{x_1}}& \cdots &{{x_n}}\\ \vdots & \vdots &{}& \vdots \\ {x_0^n}&{x_1^n}& \cdots &{x_n^n} \end{array}} \right)V=V(x0,…,xn)=⎝⎜⎜⎜⎛1x0⋮x0n1x1⋮x1n⋯⋯⋯1xn⋮xnn⎠⎟⎟⎟⎞ 而Toeplitz方程组指的是： Tnx=bT_nx=bTnx=b

其中TnT_nTn是Toeplitz矩阵， Tn=(1γ1⋯γn−1γ11⋯⋮⋮⋮⋱γ1γn−1⋯γ11)T_n = \left( {\begin{array}{ccc} 1&{{\gamma _1}}& \cdots &{{\gamma _{n - 1}}}\\ {{\gamma _1}}&1& \cdots & \vdots \\ \vdots & \vdots & \ddots &{{\gamma _1}}\\ {{\gamma _{n - 1}}}& \cdots &{{\gamma _1}}&1 \end{array}} \right)Tn=⎝⎜⎜⎜⎜⎛1γ1⋮γn−1γ11⋮⋯⋯⋯⋱γ1γn−1⋮γ11⎠⎟⎟⎟⎟⎞

一般矩阵的求解有各种分解，那么对于这类特殊的方程组，有么有什么特殊的解法呢？下面做一个简单的介绍。

基本思想

Vandermonde方程组

假设我们已经得到了范德蒙德逆转置的的一个分解V−T=ULV^{-T} = ULV−T=UL，那么V−1=LTUTV^-1 = L^TU^TV−1=LTUT，可以得到方程组的解为z=V−1b=LTUTbz=V^-1b=L^TU^Tbz=V−1b=LTUTb。
所以我们主要把力量放在的到范德蒙德逆转置的一个分解。

首先考虑这样的一个方程组求解:
VTv=yV^Tv = yVTv=y 这里v=(a0,a1,…an)Tv = (a_0,a_1,\dots a_n)^Tv=(a0,a1,…an)T.

注意这里的VTV^TVT是带转置的，并不是原来的Vandermonde方程组的求解问题了。因为范德蒙德矩阵的特殊性，不难想到，这其实是以vvv为稀系数的多项式在(xi,yi)(x_i,y_i)(xi,yi)点的插值问题，如果用cic_ici表示牛顿插值表达式的系数，其实就是kkk阶差商，我们可以通过差商之间的递推关系，找到cic_ici
和yiy_iyi的一个关系，进而，通过将牛顿插值多项式以某种巧妙的方式展开，找到cic_ici和插值多项式系数aia_iai之间的关系，换言之，就是找到了vvv和yyy之间的一个关系。这个过程比较繁琐，但是细细品味，也是很有意思的。这里考虑到期末考试将至，就不再细述了。给出一个结果，同过一步步推导，最后可以得到vvv和yyy之间的形如下面这样的一个关系：

v=ULyv = ULyv=ULy

结合前式，可以得到V−Ty=ULyV^-Ty =ULyV−Ty=ULy，有yyy的任意性，我们就得到了V−TV^-TV−T的一个形如V−T=ULV^-T=ULV−T=UL的分解，这就得到了我们想要的。

Teoplitz方程组

在说Teoplitz方程组前，可以先提一下Yule-Walker方程组，Yule-Walker方程组指的是Teoplitz方程组b=−(γ1,…,γn−1,γn)Tb = -(\gamma_1,\dots,\gamma_{n-1},\gamma_n)^Tb=−(γ1,…,γn−1,γn)T（γn\gamma_nγn自由）的特殊情况。简单地说Yule-Walker方程组的求解，是通过矩阵分块的方法，找到yky_kyk和yk+1y_{k+1}yk+1之间的关系，进而可以从一阶解y1y_1y1出发，递推得到方程组的解，这一个过程所需的运算量为3n2/23n^2/23n2/2。类似于Yule-Walker方程组的求解过程，我们也可以递推地得到一般Teoplitz方程组的解，由于时间关系，就不再细述。

程序和结果

下面直接给出程序和简单的算例和结果，我这里用的环境是VS2013：

Vandermonde方程组

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define N 6int main(int argc, char** argv) {double b[N] = { 1, 2, 3, 4, 5, 6 };double x[N] = { 5, 6, 7, 8, 9, 10 };int i, k, k0;int n = N - 1;printf("要求解的方程组为：\n x = \n");for (i = 0; i < N; i++){printf("%f ", x[i]);}printf("\n b = \n");for (i = 0; i < N; i++){printf("%f\n ", b[i]);}for (k = 0; k <= n - 1; k++){for (i = n; i > k; i--){b[i] = b[i] - x[k] * b[i - 1];}}for (k0 = k; k >= 0; k--){for (i = k + 1; i <= n; i++){b[i] = b[i] / (x[i] - x[i - k - 1]);}for (i = k; i <= n - 1; i++){b[i] = b[i] - b[i + 1];}}printf("\nthe solution is :\n");for (i = 0; i <= n; i++){printf("%f \n", b[i]);}getchar();return 0;
}

Toeplitz方程组

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define N 3int main(int argc, char** argv){double gamma[N - 1] = { 0.3, 0.6 };double b[N] = { 5, 6, 7 };int n = N;int i, k;double s = 0;double zeta[N], xi[N], delta, alpha, mu, v[N], sigma[N];zeta[0] = -gamma[0];xi[0] = b[0];delta = 1.0;alpha = -gamma[0];k = 1;printf("gamma （from 1 to %d） is\n", N - 1);for (i = 0; i <= N - 2; i++){printf("%f ", gamma[i]);}printf("\n\nb  is \n");for (i = 0; i <= N - 1; i++){printf("%f ", b[i]);}while (1)//while(k<n-1){delta = (1.0 - alpha*alpha)*delta;for (i = 1; i <= k; i++) {s = s + gamma[i - 1] * xi[k - i];}mu = (b[k] - s) / delta;s = 0.0;for (i = 1; i <= k; i++) {v[i - 1] = xi[i - 1] + mu*zeta[k - i];}for (i = 1; i <= k; i++) {xi[i - 1] = v[i - 1];}xi[k] = mu;if (k >= n - 1) break;for (i = 1; i <= k; i++){s = s + zeta[k - i] * gamma[i - 1];}alpha = -(gamma[k] + s) / delta;s = 0.0;for (i = 1; i <= k; i++) {sigma[i - 1] = zeta[i - 1] + alpha*zeta[k - i];}for (i = 1; i <= k; i++) {zeta[i - 1] = sigma[i - 1];}zeta[k] = alpha;k = k + 1;}printf("\n\nsolution:\n");for (i = 0; i <= n - 1; i++) {printf("%lf\n", xi[i]);}getchar();return 0;
}