8.2 雷米兹算法

雷米兹介绍

  切比雪夫只是提出理论，但是没给出计算方法。乌克兰数学家雷米兹Evgeny Yakovlevich Remez提出了计算切比雪夫近似值的算法。

  雷米兹对切比雪夫近似值提供一个操作简单的算法，这是个迭代过程。
  N次迭代，每次两步。对于这种迭代过程，是没法预测迭代多少步，类似的迭代过程还有线性代数中的Gram-Schmidt分解。但是他们有个共同的特点，就是数学家们证明了迭代过程的收敛性，假如不收敛，迭代再多次也是无益的，就比如Hessenberg矩阵的QR分解，迭代无数次还是Hessenberg矩阵，因为不收敛。雷米兹算法的收敛性，我这里不提供证明，因为太难了，哈哈，我也看不懂。

算法第一步——解方程

  假设是求sin(x)在(-π/16, π/16)的九阶切比雪夫近似值。这个手算是不太可能的，会比较累。九阶多项式，有10个我们要求的未知数，其形式如下：
$$a{x}^9+b{x}^8+c{x}^7+d{x}^6+e{x}^5+f{x}^4+g{x}^3+h{x}^2+ix+j$$
  我们把未知数有a到j十个字母来表示。N次迭代，按下面步骤来：
  在区间内，找任意11(根据切比雪夫的N+2定义，$N=9,9+2=11$)个点。
  第一步 求出系数，这个系数要求在11个点（N+2）上，切换正负号之后，错误幅值相等。
  将这11个点的值代入公式：
$$\begin{gathered} a{x}_0^9+b{x}_0^8+c{x}_0^7+d{x}_0^6+e{x}_0^5+f{x}_0^4+g{x}_0^3+h{x}_0^2+i{x}_0+j=sin(x_0)-e \\ a{x}_1^9+b{x}_1^8+c{x}_1^7+d{x}_1^6+e{x}_1^5+f{x}_1^4+g{x}_1^3+h{x}_1^2+i{x}_1+j=sin(x_1)+e \\ ⋮ \\ a{x}_{10}^9+b{x}_{10}^8+c{x}_{10}^7+d{x}_{10}^6+e{x}_{10}^5+f{x}_{10}^4+g{x}_{10}^3+h{x}_{10}^2+i{x}_{10}+j=sin(x_{1}0)-e \end{gathered}$$
  e是误差。所以误差e也是个未知数。这11个等式，计算出来的a~j和误差e肯定不是最终结果。这个计算过程，扩展到任意阶是这样的公式：
$$\begin{gathered} c_0+c_1{h}_i+\cdot \cdot \cdot +c_n{h}_i^n+(-1{)}^{i}e=f(x_i) \\ i=0,1,2,...,n+1 \end{gathered}$$
  公式中$h_i^n$是点的n次方,c是要求的系数。上面的例子，这一步我们得到了一个11元一次方程组，这个时候用高斯消元、LU分解或LUP分解可以解出10个系数和一个误差e。
  但是第一步的11个点，并不恰好是最大误差的11个点。他们的误差值e，也不是最小误差值。

算法第二步——重新选点

  把第一步得到的方程。利用二分法，把第一步的11个点，然后加上左边界和右边界重新划分为12个区域。
  在这12个区域，找最大误差或最小误差，就是说找误差绝对值最大的点。误差值按以下公式计算：
$$e(x)=sin(x)-p(x)$$
  在误差绝对值最大的点中去掉一个最小值，产生一个新的集合。
  一直循环迭代，直到第一步计算出的误差值e与第二步的误差最大值相等。
  比如三阶切比雪夫多项式的第二步就是下图这个样子:

  图中的函数是误差函数，图上*号是第一步的五个点，可以看出一正一负的，绝对值相等。这些区间里找到了最大值的点$x_0^{\ast },x_1^{\ast },\cdots ,x_{n+1}^{\ast }$，也就是图上x坐标的圆圈。但是这样产生的点要比需要的点多一个的，所以必须去掉一个最大误差点。图中有一个区域，是误差绝对值最小的区域，我勾出来了，是所有区域里最大误差绝对值最小的区域，它的最大误差点也就是就是被剔除的点：

Java实现

雷米兹的核心代码如下：

/**
 * @param n 几次多项式
 */
public static BigDecimal[] remez(final int n,
                         final BigDecimal min,
                         final BigDecimal max,
                         final Function<BigDecimal, BigDecimal> target) {
    Random random = new Random();
    BigDecimal step1Error = null;
    BigDecimal step2Error = null;
    final BigDecimal[] horners = new BigDecimal[n + 1];
    while ((step1Error == null && step2Error == null) || step1Error.subtract(step2Error).abs().compareTo(MAX_ERROR) > 0) {
        // step1
        // 区间内区N个任意值
        // 0 和 1 之间 产生 3和6之间 6-3 =3.0 min + 3.0 * r
        BigDecimal[] randomValues = initRandomValue(n, min, max, random);
        // 产生一个n+2阶一次方程组
        final EquationGroup equationGroup = getEquationGroup(n, target, randomValues);
        final BigDecimal[] resolve = equationGroup.resolve();
        // 第一步最大误差
        step1Error = resolve[resolve.length - 1];
        // 将第一步得到11个值，利用二分法，形成10个分界点，这10个分界点与左右边界重新形成11个区域
        Region[] regions = getRegion(randomValues, min, max);
        // 在每个区域里找最大误差值或最小误差值
        // 第一个区间是 sin(x) - p(x) 正数，也就是最大值，第二个区间是最小值
        System.arraycopy(resolve, 0, horners,0, n + 1);
        final BigDecimal[] errors = new BigDecimal[n + 2];
        for (int i = 0; i < regions.length; i++) {
            Region region = regions[i];
            // 找出这些点
            final Error maxErrorPoint = getMaxErrorPoint(region, target, horners);
            randomValues[i] = maxErrorPoint.getPoint();
            // 然后求出最大误差
            errors[i] = maxErrorPoint.getError();
        }
        step2Error = Arrays.stream(errors).max(BigDecimal::compareTo).orElse(BigDecimal.ONE);

    }

    return horners;
}

Git地址

https://e.coding.net/buildt/learn/remez.git