{(x_i, f(x_i)} が等間隔に並ぶような点列 {x_i} を求める

divfunc.py

import numpy as np

from crfmnes import CRFMNES

np.set_printoptions(formatter={'float_kind': lambda f: f"{f:.2f}"})

def divide(f, x0, xn, n):
    '''
    divide: p_i = (x_i, f(x_i)) (i = 0, 1, ..., n) について、
       x_0 = x0, x_n = xn, (x_{i-1}, f_{x_{i-1}}) = (x_i, f_{x_i}) 
       となる分割を求める。
    [x0, xn]: 区間
    n: 区間の区分数
    '''

    xs = (xn - x0) * np.arange(0, n+1) / n + x0
    xs = xs[1:-1]

    def sq_distance(x0, x1):
        return np.square(x1 - x0) + np.square(f(x1) - f(x0))

    def dists(xs):
        xs = xs.reshape(-1, 1)
        d = np.empty((n,1))
        d[0]    = sq_distance(xs[0, 0], x0)
        d[1:-1] = sq_distance(xs[1:],   xs[:-1])
        d[-1]   = sq_distance(xn,       xs[-1, 0])
        return np.sqrt(d)

    def objective(xs):
        return np.var(dists(xs))

    print('xs:', xs)
    print('distances:', dists(xs).reshape(-1,))
    print()

    optimizer = CRFMNES(np.size(xs), objective, xs.reshape(-1, 1), 1, 16)
    xs, f_best = optimizer.optimize(1000)
    xs = xs.reshape(-1,)

    print('xs:', xs)
    print('distances:', dists(xs).reshape(-1,))

print('f(x) = x*x')
divide(lambda x: x * x, -3, 3, 10)

print('f(x) = x*x*x')
divide(lambda x: x * x * x, -3, 3, 10)

考察

• 目的関数を np.var(dists(xs)) とするより np.sum(dists(xs)) とする方が素直だと思ったのだけど関数の性質によってはおかしなことに

  • 1点がx0となり、残りがすべてxnになる
    になるのが残念。区間 R について max f(R) > |R| なときにこの問題が発生する。

• もともとは dists はユークリッド距離ではなく、その二乗を返していた。それでも大丈夫かと思っていたのだけど、結果は滅茶苦茶だった。