Search on the blog

2016年7月17日日曜日

最急降下法とニュートン法の比較

autogradを使って、最急降下法(steepest descent method)とニュートン法(Newton's method)の比較をしてみた。

最急降下法のおさらい
  1. 適当な初期点xを選ぶ
  2. xでの勾配ベクトルを求めて勾配ベクトル方向(最も急な方向)に点を移動させる
  3. 収束するまで2.を繰り返す
ニュートン法のおさらい
関数fの点xnのまわりでのテイラー展開は以下のようになる。


ただし、3次以上の項は無視している。

最適化問題で求めたいのはfの停留点である。
ということで、上のテイラー展開を使って微分が零になる点を近似的に計算してみる。


なので、fをxで微分したものと、fをΔxで微分したものは同じ。
よって、

となる。この値が0になるようなΔxを求めると、

である。

つまりこの方向に進めば、fの微分が0になりそうな場所に向かって、xを逐次的に更新することができる。

ニュートン法のアルゴリズムを以下に示す。
  1. 適当な初期点xを選ぶ
  2. 上記の手順でx = x + Δxと更新する
  3. 収束するまで2.を繰り返す
実験
autogradを使って、最急降下法とニュートン法を実装する。
必要なライブラリのimport。
import autograd
import matplotlib.pyplot as plt
import numpy as np

最急降下法の実装。
グラフに描画したいので、xの最適値だけじゃなくて、xの更新履歴も返すようにしている。
def steepest_descent(f, x0, alpha=0.01, eps=1e-9, maxitr=500):
  x = x0
  history = [np.array(x)]
  for _ in range(maxitr):
    g = autograd.grad(f)
    if (np.linalg.norm(g(x)) < eps):
      break
    x -= alpha * g(x)
    history.append(np.array(x))
  return x, history

ニュートン法の実装。
def newton(f, x0, alpha=0.01, eps=1e-9, maxitr=500):
  x = x0
  history = [np.array(x)]
  for _ in range(maxitr):
    g = autograd.grad(f)
    H = autograd.hessian(f)
    if (np.linalg.norm(g(x)) < eps):
      break
    x += alpha * np.linalg.solve(H(x), -g(x))
    history.append(np.array(x))
  return x, history

軸がずれた楕円形の関数を使って、両手法を比較してみる。

def f(x, y):
  s = 2*x + y
  t =  -x + y
  return 4 * s**2 + 3 * t**2

# draw contour of f(x, y)
delta = 0.01
x = np.arange(-2.2, 0.2, delta)
y = np.arange(-2.2, 0.2, delta)
X, Y = np.meshgrid(x, y)
Z = f(X, Y)
plt.contour(X, Y, Z, 80, colors='black')

# execute steepest descent
x0 = np.array([-1.5, -2.0])
wrapper = lambda xs: f(*xs)
optimal, history = steepest_descent(wrapper, x0, alpha=0.001)
for x in history:
  plt.plot(x[0], x[1], 'ro', markeredgewidth=0.0)

# execute newton method
x0 = np.array([-1.5, -2.0])
wrapper = lambda xs: f(*xs)
optimal, history = newton(wrapper, x0)
for x in history:
  plt.plot(x[0], x[1], 'bo', markeredgewidth=0.0)
  
plt.show()

結果は以下のとおり。
最急降下法(赤点)の場合は、等高線と直交する方向にxが更新されているのが分かる。
それに対して、ニュートン法(青点)の場合は、近似的に計算した関数の停留点に向かって真直ぐxが更新されているのが分かる。

Codeforces Round #121 Div1 C. Fools and Roads

問題
n個の町が道路で結ばれている。グラフ(町, 道路)は木構造になっている。
以下の情報がk個与えられる。

u v

これは、ノードuとノードv間の単純路を人が移動したことを意味する。
各道路について、その道路を通った人の数を求めよ。

解法
u vという情報に対して、以下の処理を実行すればよい。
  • uからrootまでのすべての枝に1を足す
  • vからrootまでのすべての枝に1を足す
  • wからrootまでのすべての枝から2を引く
ただし、wはLCA(u, v)とする。
毎回計算していたら遅いので、情報をすべて集めたあとに、葉から根方向にDPしながら和を計算すればよい。

実装
LCAのライブラリを作っておけば、張ってちょっと足して終わり。
#include <algorithm>
#include <bitset>
#include <cassert>
#include <climits>
#include <cmath>
#include <cstdio>
#include <cstring>
#include <iomanip>
#include <iostream>
#include <list>
#include <map>
#include <numeric>
#include <queue>
#include <set>
#include <sstream>
#include <stack>
#include <string>
#include <vector>
#include <unordered_set>
#include <unordered_map>

using namespace std;

#define REP(i,n) for(int i=0; i<(int)(n); i++)
#define FOR(i,b,e) for (int i=(int)(b); i<(int)(e); i++)
#define ALL(x) (x).begin(), (x).end()

int n;
vector<int> edges[100000];
int par[100000];
int depth[100000];
long long cost[100000];
long long r[100000];
map<pair<int, int>, int> etoi;

class LCA {
  int V, logV;
  vector<int> depth;
  vector<vector<int> > parent;
  
  void build() {
    for (int k = 0; k + 1 < logV; k++) {
      for (int v = 0; v < V; v++) {
        if (parent[k][v] < 0) parent[k+1][v] = -1;
        else parent[k+1][v] = parent[k][parent[k][v]];
      }
    }
  }
public:
  // V: maximum number of nodes
  LCA(int V) {
    this->V = V;
    logV = 0;
    while (V > (1LL<<logV)) logV++;
    this->depth = vector<int>(V);
    this->parent = vector<vector<int> >(logV, vector<int>(V));
  }
  
  // N: number of nodes
  // p: parent of nodes (p[root] = -1)
  // d: depth of nodes (d[root = 0])
  void init(int N, int p[], int d[]) {
    for (int i = 0; i < N; i++) {
      parent[0][i] = p[i];
      depth[i] = d[i];
    }
    this->build();
  }
  
  // p: parent of nodes (p[root] = -1)
  // d: depth of nodes (d[root = 0])
  void init(vector<int> p, vector<int> d) {
    int N = d.size();
    for (int i = 0; i < N; i++) {
      parent[0][i] = p[i];
      depth[i] = d[i];
    }
    this->build();
  }
  
  int query(int u, int v) {
    if (depth[u] > depth[v]) swap(u, v);
    for (int k = 0; k < logV; k++) {
      if ((depth[v] - depth[u]) >> k & 1)
        v = parent[k][v];
    }
    if (u == v) return u;
    
    for (int k = logV-1; k >= 0; k--) {
      if (parent[k][u] != parent[k][v]) {
        u = parent[k][u];
        v = parent[k][v];
      }
    }
    return parent[0][u];
  }
};

void rec(int v, int p = -1) {
  for (auto w: edges[v]) {
    if (w != p) {
      rec(w, v);
      r[etoi[make_pair(v, w)]] += cost[w];
      cost[v] += cost[w];
    }
  }
}

void dfs(int v, int p = -1, int d = 0) {
  par[v] = p;
  depth[v] = d;
  for (auto &w: edges[v]) {
    if (w != p)
      dfs(w, v, d+1);
  }
}

int main() {
  scanf(" %d", &n);
  REP (i, n-1) {
    int x, y;
    scanf(" %d %d", &x, &y);
    --x, --y;
    etoi[make_pair(x, y)] = i;
    etoi[make_pair(y, x)] = i;
    edges[x].push_back(y);
    edges[y].push_back(x);
  }

  dfs(0);

  LCA lca(100000);
  lca.init(n, par, depth);

  int k;
  scanf(" %d", &k);
  REP (i, k) {
    int x, y;
    scanf(" %d %d", &x, &y);
    --x, --y;
    ++cost[x];
    ++cost[y];
    int z = lca.query(x, y);
    cost[z] -= 2;
  }

  rec(0);
  REP (i, n-1)
    cout << r[i] << " ";
  cout << endl;

  return 0;
}

便利なPythonライブラリ(8)autograd

その名の通り自動微分してくれるライブラリです。

インストール
pipでインストール出来ます。

サンプル(1)
f(x) = 3x2という関数の一次微分および二次微分を計算してみます。

from autograd import grad
import matplotlib.pyplot as plt
import numpy as np

def f(x):
  return 3 * x**2

df = grad(f)
d2f = grad(df)

x = np.linspace(-3, 3, 100)
plt.plot(x, f(x), label='f')
plt.plot(x, np.vectorize(df)(x), label='df')
plt.plot(x, np.vectorize(d2f)(x), label='d2f')
plt.legend()
plt.show()

以下のようなグラフが表示されます。
サンプル(2)
もちろん、多変数関数の微分も可能です。
ナブラとヘッセ行列を求めてみます。
簡単のため、二次形式の関数をサンプルとして使います。

from autograd import grad as nabla
from autograd.convenience_wrappers import hessian
import numpy as np

def f(x):
    Q = [[3, 2],
        [2, 7]]
    return 1./2 * np.dot(x, np.dot(Q, x))

x = np.array([1.0, 2.0])
print (nabla(f)(x))
print (hessian(f)(x))

出力結果を見ると、正しく微分できていることが確認できます。
$ python sample_autograd.py
[  7.  16.]
[[ 3.  2.]
 [ 2.  7.]]

2016年7月16日土曜日

SRM 657 Div2 1000 PolynomialRemainder

問題
a, b, cが与えられる。
a x^2 + b x + c = 0 (mod 1,000,000,000)
となるようなxを求めよ。

解法
10^9で割って0になるので、
2^9で割って0、5^9で割って0となるようなxを求めて、中国剰余定理をすればよい。

実装
using namespace std;

#define ALL(x) (x).begin(), (x).end()
#define EACH(itr,c) for(__typeof((c).begin()) itr=(c).begin(); itr!=(c).end(); itr++)  
#define FOR(i,b,e) for (int i=(int)(b); i<(int)(e); i++)
#define MP(x,y) make_pair(x,y)
#define REP(i,n) for(int i=0; i<(int)(n); i++)

long long extgcd(long long a, long long b, long long &x, long long &y) {
  long long d = a;
  if (b != 0) {
    d = extgcd(b, a%b, y, x);
    y -= a/b * x;
  } else {
    x = 1;
    y = 0;
  }
  return d;
}

pair<long long, long long> crm(long long a, long long n, long long b, long long m) {
  a %= n;
  b %= m;
  long long p, q;
  extgcd(n, m, p, q);
  long long mod = n * m;
  long long x = (a * q  * m) % mod + (b * p * n) % mod;
  x = (x % mod + mod) % mod;
  return make_pair(x, mod);
}

long long calc(long long a, long long b, long long c, long long x, long long mod) {
  long long ret = ((a % mod * x) % mod) * x % mod;
  ret = (ret + b * x) % mod;
  ret = (ret + c) % mod;
  return ret;
}

class PolynomialRemainder {
public:
  int findRoot(int a, int b, int c) {
    long long s, t;

    int n = 1;
    REP (i, 9) n *= 2;
    for (s = 0; s < n; s++) {
      if (calc(a,b,c,s,n) == 0)
        break;
    }
    if (s == n)
      return -1;

    int m = 1;
    REP (i, 9) m *= 5;
    for (t = 0; t < m; t++) {
      if (calc(a,b,c,t,m) == 0)
        break;
    }
    if (t == m)
      return -1;

    cout << s << " " << t << endl;
    auto ret = crm(s, n, t, m);
    return ret.first;
  }
};

2016年7月3日日曜日

SRM 653 Div2 250 CountryGroup

問題概要
椅子にn人の人が座っている。同じ国から来た人は必ず隣に座っているらしい。
左から順に、あなたの国からは何人の人がここにいますか?と聞いていく。
これに対する回答がn個与えられるので、矛盾しないかどうかと矛盾しない場合は何カ国の人がここにいるかを求めよ。

解法
stackを使うと綺麗に書ける。
using namespace std;

#define ALL(x) (x).begin(), (x).end()
#define EACH(itr,c) for(__typeof((c).begin()) itr=(c).begin(); itr!=(c).end(); itr++)  
#define FOR(i,b,e) for (int i=(int)(b); i<(int)(e); i++)
#define MP(x,y) make_pair(x,y)
#define REP(i,n) for(int i=0; i<(int)(n); i++)

class CountryGroup {
  public:
  int solve(vector<int> a) {
    int r = 0;
    while (a.size()) {
      int x = a.back();
      REP (i, x) {
        if (a.empty() || a.back() != x)
          return -1;
        a.pop_back();
      }
      ++r;
    }
    return r;
  }
};

2016年7月2日土曜日

WSGI(2)

前回の続き。
uwsgiを起動するときのオプションは、設定ファイルに保存しておくことができる。

kenjih$ cat uwsgi.ini 
[uwsgi]
base = /Users/kenjih/tmp/python/wsgi
app = myflaskapp
module = %(app)
pythonpath = %(base)
callable = app
master = true
processes = 4
threads = 2
socket = 127.0.0.1:3031
daemonize=/var/log/uwsgi/%(app).log
touch-logreopen = %(base)/.logreopen_trigger
touch-reload = %(base)/.reload_trigger
harakiri = 60

例えば上のような設定ファイルを作っておくと、
$ uwsgi --ini uwsgi.ini
のようにuwsgiを起動できる。

daemonizeを指定することで、uwsgiをデーモンとして起動できる。引数には標準出力/標準エラー出力の出力先ファイルを指定できる。

touch-logreopenには、ログファイルのロテートのtriggerとなるファイル名を指定できる。これは試してみたがうまく動かなかった。これはlogrotateなどでロテートしたときに、ファイルがロテートされたのでファイルをreopenしてくれとuWSGI側に合図を送るためのファイル。指定したファイルをtouchすればログがロテートされるというものではないので注意。logrotateでロテートしたい場合は、postrotateスクリプトで指定したファイルをtouchするというような使い方をする。

touch-reloadには、リソースのデプロイを行うためのtriggerファイルを指定できる。指定したファイルをtouchすれば、リソースをホットデプロイしてくれる。

harakiriは指定した秒数以上かかるタスクを強制終了するという設定。

WSGI

WSGIとは?
PythonのアプリケーションサーバとWebサーバを接続するための標準インターフェース定義。ウィズギーと読む。

昔々、PythonでWebアプリケーションをつくる際、
使うフレームワークによって使えるWebサーバが決まっていたり、
逆に使うWebサーバによって使えるフレームワークが決まっていたりした。

この問題を解決するために、WSGIという標準インターフェースが定義された。

uWSGIとは?
WSGIで定義されたインターフェースを実装したもの。

インストール
uwsgiをインストール。pipで入る。
pip install uwsgi

Webサーバとアプリサーバの連携
Flaskで作ったアプリをnginxで動かしてみる。
公式ページにサンプルがあったのでそれを使ってみる。

以下をmyflaskapp.pyというファイル名で保存しておく。
from flask import Flask

app = Flask(__name__)

@app.route('/')
def index():
    return "I am app 1"

以下でuwsgiを起動。
uwsgi --socket 127.0.0.1:3031 --wsgi-file myflaskapp.py --callable app --processes 4 --threads 2

nginxの設定ファイルを以下のように変更する。
$ diff -u nginx.conf.default nginx.conf
--- nginx.conf.default 2016-07-02 17:45:31.000000000 +0900
+++ nginx.conf 2016-07-02 18:02:36.000000000 +0900
@@ -33,7 +33,7 @@
     #gzip  on;
 
     server {
-        listen       8080;
+        listen       80;
         server_name  localhost;
 
         #charset koi8-r;
@@ -41,8 +41,8 @@
         #access_log  logs/host.access.log  main;
 
         location / {
-            root   html;
-            index  index.html index.htm;
+            include uwsgi_params; 
+            uwsgi_pass 127.0.0.1:3031;
         }

nginx起動。
$ sudo nginx

80番ポートにアクセスすると、flaskのアプリがみれた。