Tensorflowで数千パラメータの最小二乗法fittingに挑む その2

fitting ステージ2

今回もまだ、Tensorflow fitting のための初期値推算としての fitting を行う。

前回得られた結果 R5 を RDD からオンメモリの numpy array に変換する。 R5 の構造は

((47, 61), 41800.269085989807), ((47, 62), 41800.351512553934), ((47, 63), 41806.114347614603)]

のようなものだった。

これを

>>> R6=np.array(list(map(lambda l: (l[0][0], l[0][1], l[1]), R5.collect())))
>>> print(R6)
>>> [[  0.00000000e+00   0.00000000e+00   4.40547959e+04]
 [  0.00000000e+00   1.00000000e+00   4.40470300e+04]
 [  0.00000000e+00   2.00000000e+00   4.40379737e+04]
 ..., 
 [  4.70000000e+01   6.10000000e+01   4.18002691e+04]
 [  4.70000000e+01   6.20000000e+01   4.18003515e+04]
 [  4.70000000e+01   6.30000000e+01   4.18061143e+04]]

と、mapで並べ直したものを list 化し（python3の都合でiterableで返って来るため） np.array 化する。R6を転置すれば　px, py, fitting で得られた Lorentz関数 peak 中心になる。

[ここで後述するように、どうしても fitting が合わないので調べてみるとこのモデルでは合わせられないデータを使っていたことが判明してデータを作り直している。前処理の方法は変わらず、データの形状も変わっていない。
また、scipy でも検証計算をやるために以下の DeltaY の演算も R6 の作成の中に入れた。]

>>> R6=np.array(list(map(lambda l: (l[0][0], l[0][1]*dy, l[1]), R5)))
>>> print(R6)
>>> [[  0.00000000e+00   1.87500000e-02   3.33514943e+04]
 [  0.00000000e+00   2.81250000e-02   3.33532640e+04]
 [  0.00000000e+00   3.75000000e-02   3.33575154e+04]
 ..., 
 [  7.10000000e+01   5.53125000e-01   3.43357848e+04]
 [  7.10000000e+01   5.62500000e-01   3.43314231e+04]
 [  7.10000000e+01   5.71875000e-01   3.43252143e+04]]

>>> print(len(R6))
>>> 4253

ここから面倒が続く。元データの px, py は整数、すなわちインデックス値だった。これを実座標に変換する。ところが、 Y については固定値 DeltaY をかけた値で固定されるが X に関しては fitting パラメータを含む関数で変換してやる必要がある。以下にその変換関数の一つを載せる。 gap が fitting パラメータで、 fitting の際に毎回変換する必要があるため X については事前の変換は行わず fitting の内部で変換し、一方 Y については上に示したように事前に変換している。

def posX1(px, deltaX, gap):
    return px*deltaX + (np.trunc(px/8)*gap)

さらに、空間座標 px, py から、ある係数を計算する関数がある。

def transFunc1(Params, data):
    """
    @param Params [X, Y, al, bl, Inc_x, Inc_y, Inc_z]
    @param data [[pixel x][pixel y]]
    @returns (al+sqrt(r^2+bl^2)) * sin(theta)
    """
    lx = data[0]-Params[0]
    ly = data[1]-Params[1]
    r2 = lx*lx + ly*ly

    r0 = np.sqrt(r2+Params[3]*Params[3])
    I0 = np.sqrt(Params[4]*Params[4] + Params[5]*Params[5] + Params[6]*Params[6])

    I1 = lx * Params[4]
    I2 = ly * Params[5]

    I3 = I1 + I2 + (-Params[3])*Params[6]
    sin_theta = np.sqrt((1-I3 / (I0*r0))/2)

    C1 = (Params[2] + r0) * sin_theta

    return C1

変換パラメータが7変数あり、 numpy array で受け取った px, py から、係数ベクトルを計算する関数である。この係数とある数の積が上述の Lorentz 関数のピーク中心と一致するように最小二乗法 fitting を行うのだが、このある数自身も正しい値が不明なので fitting で求めてしまう。変換パラメータは7つのうち5変数が fitting 対象で、2変数は今回固定である。

これらを使って残差関数を定義する。この残差関数は fitting パラメータ、実データ、固定パラメータの3引数を取る。

def estimate(Params, data, P0):
    """
    estimation of foot of perpendicular
    @param Params [X, Y, Inc_x, Inc_y, Inc_z, pvar, gap]
    @param data   [pixel_X, pixel_Y, peak_center]
    @param P0     [al, bl, posX, deltaX]
    """

    C0 = transFunc1((Params[0],Params[1], P0[0], P0[1], Params[2], Params[3], Params[4]), (P0[2](data[0], P0[3], Params[6]), data[1]))

    C1 = (C0*Params[5])-data[2]
    return C1

Paramsが fitting パラメータ、dataが px, py, peak中心のベクトルで、残りがP0に押し込んである。P0の前2つのパラメータは transFunc1 の変換パラメータに挿入される。 P0 の第3パラメータには X の座標変換関数 posX1() を渡し、posX1() の引数の片方は　Params から取り、もう一つは P0 の第4引数から取る。

これを Tensorflow に移植して解かせる。

ステージ2 fitting with TF

x, y, dat は上の R6 から取り出して feed_dict を使って流し込む。

    R6X = R6.T[0]
    R6Y = R6.T[1]
    R6P = R6.T[2]

    p_x=tf.placeholder(tf.float32)
    p_y=tf.placeholder(tf.float32)
    p_dat=tf.placeholder(tf.float32)

関数も移植してゆく。gap は fitting パラメータなので　tf.Variable とし、 dx, dy は python の通常の変数である。 trunc がないので round を使っている。

    gap = tf.Variable(0.0, tf.float32)
    posX = p_x * dx + tf.round(p_x/8-0.5) * gap

truncFunc1() は名前空間なしで全部フラットに展開される。al, blは通常のpython変数でfittingに対しては固定パラメータである。

    p_1 = tf.Variable(0.0, tf.float32)
    p_2 = tf.Variable(0.0, tf.float32)
    p_3 = tf.Variable(0.0, tf.float32)
    p_4 = tf.Variable(0.0, tf.float32)
    p_5 = tf.Variable(1.0, tf.float32)

    lx = posX - p_1
    ly = p_y  - p_2
    r2 = lx*lx + ly*ly

    r0 = tf.sqrt(r2+bl*bl)
    I0 = tf.sqrt(p_3*p_3 + p_4*p_4 + p_5*p_5)

    I1 = lx * p_3
    I2 = ly * p_4

    I3 = I1 + I2 + (-bl)*p_5
    sin_theta = tf.sqrt((1-I3 / (I0*r0))/2)

    C1 = (al + r0) * sin_theta

残差二乗和まで求める。残差二乗和 RSS(Residual Sum of Squares) は Est と p_dat の差から計算する。p_dat は fitting する相手となる実測値 R6.T の第３項で Est が推算値である。

    p_6 = tf.Variable(505.56, tf.float32)
    Est = C1 * p_6
    rss = tf.reduce_mean(tf.square(Est - p_dat))

    optimizer = tf.train.GradientDescentOptimizer(0.0000001)
    lsq = optimizer.minimize(rss) 
    sess.run(tf.initialize_all_variables())

学習係数が極端に小さな値になっているが、このオーダーでなければ即RSSがNANとなって収束させられなかったためである。

    print(sess.run(Est, feed_dict={p_x: R6X, p_y: R6Y, p_dat: R6P}))
    print(sess.run(Est-p_dat, feed_dict={p_x: R6X, p_y: R6Y, p_dat: R6P}))
    print(sess.run(rss, feed_dict={p_x: R6X, p_y: R6Y, p_dat: R6P}))

    for step in range(1601):
        sess.run(lsq, feed_dict={p_x: R6X, p_y: R6Y, p_dat: R6P})
        if step % 20 == 0:
            print(step+1, sess.run(rss, feed_dict={p_x: R6X, p_y: R6Y, p_dat: R6P}),
                  sess.run(p_1), sess.run(p_2), sess.run(p_3), sess.run(p_4), sess.run(p_5),
                  sess.run(p_6), sess.run(gap))

    print(sess.run(sin_theta, feed_dict={p_x: R6X, p_y: R6Y, p_dat: R6P}))
    print(sess.run(Est, feed_dict={p_x: R6X, p_y: R6Y, p_dat: R6P}))
    print(sess.run(Est-p_dat, feed_dict={p_x: R6X, p_y: R6Y, p_dat: R6P}))

これで計算してみると、最後の方の出力は以下のようになった。

1501 8.00984e+06 0.521582 0.295551 0.0 0.0 1.0 578.134 0.01
1521 7.82816e+06 0.521913 0.295554 0.0 0.0 1.0 578.737 0.01
1541 7.65061e+06 0.522246 0.295557 0.0 0.0 1.0 579.333 0.01
1561 7.47715e+06 0.522584 0.29556 0.0 0.0 1.0 579.922 0.01
1581 7.30765e+06 0.522924 0.295563 0.0 0.0 1.0 580.504 0.01
1601 7.14205e+06 0.523269 0.295566 0.0 0.0 1.0 581.08 0.01
[ 0.98966622  0.98980838  0.98994589 ...,  0.99174905  0.99160904
  0.9914642 ]
[ 31103.36914062  31107.1328125   31110.77734375 ...,  31158.5703125
  31154.85546875  31151.015625  ]
[-2248.12695312 -2246.1328125  -2246.73828125 ..., -3177.21484375
 -3176.56640625 -3174.19921875]

これでも gap および p_3, p_4, p_5 は tf.constant() として固定したもので、これらの値はほぼこの値になるはずにも関わらず、そのほかの変数については全く正しく収束していない。
また、1600stepも進めている割には残差二乗和は巨大な数のままであり、しかも此の期に及んで未だ収束しきっている様子でもない。

あまりに酷いので検証用に別に scipy で以下のように fitting を行ってみると

    r=leastsq(estimate, (0.0, 0.0,  0.0, 0.0, 1.0, 505.56, 0.0),
              #factor=0.1,
              full_output=True,
              args=(R6.T, (al, bl, posX1, dx)))
    print(r[0])
    print(r[3])
    print("num of calls", r[2]["nfev"])

    I0=np.array((r[0][2:5]))
    I0 /= np.linalg.norm(I0)
    print("estimated center:", r[0][0:2])
    print("estimated vector:", I0.tolist())
    print("d=", r[0][5]/505.56)
    print("gap=", r[0][6])

以下のような結果となって

[  1.17016513e+00   3.14490697e-01   1.30752216e+00   1.11290020e+00
   1.70130624e+07   6.36972788e+02   1.02428484e-02]
Both actual and predicted relative reductions in the sum of squares
  are at most 0.000000
num of calls 236
estimated center: [1.1701651346484467, 0.31449069722719936]
estimated vector: [7.685401522791991e-08, 6.541445481803861e-08, 0.999999999999995]
d= 1.25993509795
gap= 0.0102428484062

妥当な値が得られることが確認できた。

とりあえずここまでのまとめ

Tensorflow の optimizer に使われている最急降下法は、勾配が浅いところでは　1step に変化させるパラメータの刻みが小さくなる。 scipy (Levenerg-Marquardt 法)で 236 step で収束しているところ、 tensorflow では 1600 step でもまだ収束する見込みが見られない。
また、最急降下法では local minimum にトラップされる恐れが大きいが、今回の結果でも scipy の収束値とは異なる場所で解の探索を行っている様子が確認できた。
tensorflow を最小二乗法 solver として使うには恐ろしく高いハードルがある。

収束はしないのかというと、step数を増やしてやれば

4801 238972.0 0.644014 0.29685 0.0 0.0 1.0 622.827 0.01
9601 23594.8 0.902129 0.301772 0.0 0.0 1.0 631.677 0.01
14401 13046.4 1.01098 0.305508 0.0 0.0 1.0 633.85 0.01

と scipy の値に近づきつつはある。ただし 4変数は固定である。全部有効にすると以下のようになって

1601 15889.9 8.16842 -12.9467 11.9864 -15.8953 12.4997 507.482 -0.500613
4801 15870.9 8.16576 -12.9602 12.0444 -15.8461 12.5089 507.482 -0.498695
9601 15843.3 8.16513 -12.9785 12.1298 -15.7728 12.5226 507.482 -0.495877
14401 15816.6 8.16513 -12.9969 12.2122 -15.6996 12.5357 507.482 -0.493181

local minimum にトラップされていることがわかる。