機械学習（4）確率的勾配降下法

前回とほとんど同じなのですが、

$\vec{w}^{(t+1)} = \vec{w}^{(t)} - \eta^{(t)}\nabla\ell(\vec{x}^{(t)}, y^{(t)}, \vec{w}^{(t)})$

ηという学習率を掛けるところが違います。これは時間の関数になっていて、t = 1, 2, …に対し、恐らく単調減少になっているほうが収束が速くなると思われます。

η= 1

次元	1000サンプル	10000サンプル	100000サンプル
2	0.0433	0.0147	0.0064
10	0.1718	0.0800	0.0369
100	0.4307	0.1888	0.0865

$\eta=\frac{1}{\sqrt{t}}$

次元	1000サンプル	10000サンプル	100000サンプル
2	0.0149	0.0043	0.0013
10	0.0535	0.0148	0.0046
100	0.3538	0.0488	0.0134

$\eta=t^{-0.6}$

次元	1000サンプル	10000サンプル	100000サンプル
2	0.0093	0.0020	0.0005
10	0.0441	0.0078	0.0018
100	0.4771	0.0527	0.0062

$\eta=t^{-0.7}$

次元	1000サンプル	10000サンプル	100000サンプル
2	0.0095	0.0011	0.0002
10	0.0701	0.0054	0.0007
100	0.6732	0.1601	0.0106

このあたりに最適な関数があるようです。
コードは並列化しました。でも思ったほど速くならない。

from itertools import combinations, islice, izip
from math import sqrt, log, cos, acos, pi
import random
import sys
import time
import multiprocessing as mp

def average(v):
    n, s = reduce(lambda (x, y), e: (x + 1, y + e), v, (0, 0))
    return s / n

# inner product
def ip(xs, ys):
    return sum(x * y for x, y in izip(xs, ys))

def Box_Muller():
    r1 = random.random()
    r2 = random.random()
    return sqrt(-2.0 * log(r1)) * cos(2.0 * pi * r2)

def random_hyper_sphere(D):
    xs = [ Box_Muller() for _ in range(D) ]
    r = sqrt(sum(x * x for x in xs))
    return [ x / r for x in xs ]

def make_hyper_plane(D):
    return random_hyper_sphere(D) + [ 0.0 ]

def gen_points(D):
    while True:
        yield [ random.uniform(-1.0, 1.0) for _ in range(D) ]

def perceptron(D, N):
    W = make_hyper_plane(D)
    
    w = [ 0.0 ] * (D + 1)
    for t, x in izip(xrange(1, N + 1), gen_points(D)):
        y = 1 if ip(W, x) > 0 else -1
        y1 = 1 if ip(w[:-1], x) + w[-1] > 0 else -1
        if y * y1 <= 0:
            eta = 1.0 / sqrt(t)
            w = ([ v + eta * y * x1 for v, x1 in izip(w, x) ] +
                                                [ w[-1] + eta * y ])
    
    mag = sqrt(sum(v * v for v in w[:-1]))
    w1 = [ v / mag for v in w ]
    return acos(sum(v * V for v, V in izip(w1, W)))

def mean_diff(D, N, M):
    pool = mp.Pool(3)
    diffs = pool.map(lambda k: perceptron(D, N), xrange(M))
    return average(diffs)

if __name__ == '__main__':
    t0 = time.clock()
    D = int(sys.argv[1])    # dimension
    N = int(sys.argv[2])    # number of samples
    M = int(sys.argv[3])    # number of trial
    print "%.4f" % mean_diff(D, N, M)
    print time.clock() - t0