終点枕崎

にゃーんo(^・x・^)o

バッチ正規化層でのgradの逆伝播

「バッチ正規化 逆伝播」で検索すると、計算グラフを用いた説明ばかり出てきます。数式での証明も欲しいものです。そこで、この記事は、数式を使って導出することを目標とします。一応Chainerの実装 batch_normalization.py で確認しましたが、間違いや誤植などあれば教えていただけると幸いです。

バッチ正規化は、レイヤー間を流れるデータの分布が、適度な広がりと偏りを持つように正規化することを目的とします。具体的には、各データ要素に対し、ミニバッチ全体の平均と分散を求め、データの平均が0、分散が1になるように調整します。

データ要素を表す添字は省略
\(N\) : バッチサイズ
\(x_i\) : 入力。iはミニバッチのi番目のデータであることを表す添字。
\(y_i\) : 出力。iはミニバッチのi番目のデータであることを表す添字。
\( \epsilon \) : 定数。小さい。

\begin{eqnarray}
\overline{x} & = & \frac{1}{N} \sum_{k = 1}^{N} x_k \\
s^2 & = & \frac{1}{N} \sum_{k = 1}^{N} (x_i - \overline{x})^2 \\
\hat{x_i} & = & \frac{x_i - \overline{x}}{ \sqrt{s^2 + \epsilon} } \\
y_i & = & \gamma \hat{x_i} + \beta
\end{eqnarray}

逆伝播では、\( \partial E / \partial {\bf y} \) を受け取って、最終的に\( \partial E / \partial {\bf x} \) を計算します。ただし、誤差関数を \(E\) としました。

バッチ正規化が学習するパラメータは、\( \gamma \) と \( \beta \) です。なので、まずは \( \gamma \) と \( \beta \) についての偏微分を求めます。\( \gamma \) と \( \beta \) はどちらも、変化すると任意の \( y_i \) が変化するため、連鎖律から

\begin{eqnarray}
\frac{\partial E}{\partial \gamma} & = & \sum_{k} \frac{\partial E}{\partial y_k} \frac{\partial y_k}{\partial \gamma} = \sum_{k} \frac{\partial E}{\partial y_k} \hat{x_k} \\
\frac{\partial E}{\partial \beta} & = & \sum_{k} \frac{\partial E}{\partial y_k} \frac{\partial y_k}{\partial \beta} = \sum_{k} \frac{\partial E}{\partial y_k}
\end{eqnarray}

と計算できます。
\( y_i \) と \( \hat{x_i} \) に関する偏微分の関係は、(4)から次のようになります。

\begin{equation}
\frac{\partial E}{\partial \hat{x_i}} = \gamma \frac{\partial E}{\partial y_i}
\end{equation}

あとは、 \( \hat{x_i} \) と \( x_i \) に関する偏微分の関係を求めればOKです。

\begin{equation}
\frac{\partial E}{\partial x_i} = \sum_{k} \frac{\partial E}{\partial \hat{x_k} } \frac{\partial \hat{x_k} }{\partial x_i} = \frac{\partial E}{\partial \hat{x_i} } \frac{\partial \hat{x_i} }{\partial x_i} + \sum_{k \neq i} \frac{\partial E}{\partial \hat{x_k} } \frac{\partial \hat{x_k} }{\partial x_i}
\end{equation}

ここで、\( \frac{\partial}{\partial x_i} (\sqrt{s^2 + \epsilon}) \) をあらかじめ計算しておきます。

\begin{equation}
\frac{\partial }{\partial x_i} (\sqrt{s^2 + \epsilon}) = \frac{s}{\sqrt{s^2 + \epsilon}} \frac{\partial s}{\partial x_i} \\
\end{equation}
(2)を微分して
\begin{eqnarray}
2s \frac{\partial s}{\partial x_i} & = & \frac{1}{N} \biggl[ 2(x_i - \overline{x}) \biggl(1 - \frac{1}{N} \biggr) + \sum_{k \neq i} 2(x_k - \overline{x}) \biggl(-\frac{1}{N} \biggr) \biggr] \nonumber \\
& = & \frac{1}{N} \bigg[ 2(x_i - \overline{x}) - 2\sum_k \frac{x_k - \overline{x}}{N} \bigg] \nonumber \\
& = & \frac{2}{N} (x_i - \overline{x}) \nonumber
\end{eqnarray}
より
$$
\frac{\partial s}{\partial x_i} = \frac{1}{N} \frac{x_i - \overline{x}}{s}
$$
(3)と(9)を使って
\begin{equation}
\frac{\partial}{\partial x_i} (\sqrt{s^2 + \epsilon}) = \frac{\hat{x_i}}{N}
\end{equation}
次に、(3)を微分し、(8)の右辺第1項の \( \partial \hat{x_i} / \partial x_i \)を求めます。途中で(10)を用います。
\begin{eqnarray}
\frac{\partial \hat{x_i}}{\partial x_i} & = & \frac{1}{s^2 + \epsilon} \bigg[ \bigg( 1 - \frac{1}{N} \bigg) \sqrt{s^2 + \epsilon} - (x_i - \overline{x}) \frac{\hat{x_i}}{N} \bigg] \nonumber \\
& = & \frac{1}{\sqrt{s^2 + \epsilon}} \bigg[ \bigg( 1 - \frac{1}{N} \bigg) - \frac{\hat{x_i}^2}{N} \bigg]
\end{eqnarray}
次に、\( k \neq i \) のもと、(3)を微分し、(8)の右辺第2項の \( \partial \hat{x_k} / \partial x_i \)を求めます。
\begin{eqnarray}
\frac{\partial \hat{x_k}}{\partial x_i} & = & \frac{1}{s^2 + \epsilon} \bigg[ - \frac{1}{N} \sqrt{s^2 + \epsilon} - (x_k - \overline{x}) \frac{\hat{x_i}}{N} \bigg] \nonumber \\
& = & \frac{1}{\sqrt{s^2 + \epsilon}} \bigg[ -\frac{1}{N} - \frac{\hat{x_i} \hat{x_k}}{N} \bigg]
\end{eqnarray}

(11)と(12)を(8)に代入し整理すると、
\begin{eqnarray}
\frac{\partial E}{\partial x_i} & = & \frac{1}{\sqrt{s^2 + \epsilon}} \Bigg[ \frac{\partial E}{\partial \hat{x_i}} \bigg[ \bigg( 1 - \frac{1}{N} \bigg) - \frac{\hat{x_i}^2}{N} \bigg]
+ \sum_{k \neq i} \frac{\partial E}{\partial \hat{x_k} } \bigg[ -\frac{1}{N} - \frac{\hat{x_i} \hat{x_k}}{N} \bigg] \Bigg] \nonumber \\
& = & \frac{1}{\sqrt{s^2 + \epsilon}} \Bigg[ \frac{\partial E}{\partial \hat{x_i}} - \frac{1}{N} \bigg[ \sum_k \frac{\partial E}{\partial \hat{x_k}} + \hat{x_i} \sum_k \frac{\partial E}{\partial \hat{x_k}} \hat{x_k} \bigg] \Bigg] \nonumber \\
& = & \frac{\gamma}{\sqrt{s^2 + \epsilon}} \Bigg[ \frac{\partial E}{\partial y_i} - \frac{1}{N} \bigg[ \sum_k \frac{\partial E}{\partial y_k} + \hat{x_i} \sum_k \frac{\partial E}{\partial y_k} \hat{x_k} \bigg] \Bigg] \nonumber
\end{eqnarray}
ただし途中で(7)を使いました。最後に、(5)と(6)を代入すると、
\begin{equation}
\frac{\partial E}{\partial x_i} = \frac{\gamma}{\sqrt{s^2 + \epsilon}} \Bigg[ \frac{\partial E}{\partial y_i} - \frac{1}{N} \bigg[ \frac{\partial E}{\partial \beta} + \hat{x_i} \frac{\partial E}{\partial \gamma} \bigg] \Bigg]
\end{equation}

と比較的簡単な式で表せました。実際、Chainerもこの式で導出しています。まとめると、バッチ正規化レイヤーでのbackwardの演算は、
1. (5), (6)式で \( \partial E / \partial \gamma \) と \( \partial E / \partial \beta \) を求める
2. (13)式で \( \partial E / \partial x_i \) を求める
という流れになります。実装は省略します。

畳み込みニューラルネットで手書き文字認識をやってみる(MNIST)

ここ最近やったことを書き留めておきます。

前回からニューラルネットの実装にいくつか昨日の追加や変更を行ないました。

追加点

  • レイヤー
    • 畳み込みレイヤー
    • プーリングレイヤー
      • max pooling layer
      • average pooling layer
    • 多変数関数用活性化レイヤー(中間層にsoftmaxなどが利用できるようにするため)
    • ドロップアウトレイヤー
    • 誤差関数レイヤー
  • Optimizer
    • 確率的勾配法
    • モメンタム付き確率的勾配法
    • AdaGrad
    • RMSprop
    • AdaDelta
    • Adam
  • 重み減衰

RMSpropは、有名な手法ですが論文化されておらず、大学の講義で発表したのみだそうです。面白いですね。

変更点

前回の設計が適当すぎたので、設計を大幅に見直しました。最終的には、ゼロから作るDeepLearningの実装に似たものになりました。実はゼロ作の実装はChainerの実装とそっくりで、今回はChainerも参考に設計しています。

  • 線形変換層と活性化層の分離
  • 各層をクラス化
  • レイヤー間を流れるデータの次元の順番を変更

実装で突っかかったところ

畳み込み層はim2colで実装しました。im2colは、出力層の各ピクセルについて、フィルターの適用範囲のピクセルを重複ありで一覧のように並べ換える変換です。

f:id:uchifuji:20180219015316p:plain

その逆伝播の実装で、im2colの逆のようなことをするcol2imを使うのですが、ここで重要なのは、実はcol2imはim2colの逆写像ではないということです。col2imは偏微分を伝播するのだから当たり前ですけどね。はじめこれに気づかずに実装しましたが、im2colは1変数対1変数の変換ではないので、偏微分の連鎖律を適用しなければいけません。(実際は、間違った方法でもなぜか上手く動きました。不思議です。)

f:id:uchifuji:20180219021020p:plain
正しくはこう。

まだやってないこと

BatchNormalizationの逆伝播の計算が途中なので、計算が完了し次第実装するつもりです.

MNISTに再チャレンジ

CNNが構築できるようになったので、再度MNIST(手書き文字認識)に挑戦してみます。
ネットワークの構成はTensorFlowのMNISTチュートリアルとほぼ同じです。OptimizerはAdamを選択、学習率は0.001、{\beta_1}{\beta_2}{\epsilon}は推奨値ををそのまま設定しました。バッチサイズを50とし、5epoch分学習させます。

f:id:uchifuji:20180219050120p:plain

学習結果は以下になります。

  • 正解率(全データ)
テストデータ 訓練データ
99.14% 99.55%
  • 正解率の推移(先頭1000データ)

f:id:uchifuji:20180219050402p:plain

さすがAdamといった感じで、学習速度も正確性も良好でした。ちなみにMNIST世界一位は99.79%で、DropConnectというドロップアウトを重みに適応した手法を用いてるそうです。

最後に、MNISTの学習に使用したソースコードを載せておきます。作ったライブラリの実装はgithubに上げておいたのでそちらをご覧ください。

github.com


参考文献

  • 深層学習 - 岡谷貴之
  • ゼロから作るDeepLearning - 斎藤 康毅

深層学習

お久しぶりです。

最近、岡谷貴之先生の「深層学習」を読みました。楽しいです。

深層学習 (機械学習プロフェッショナルシリーズ)

深層学習 (機械学習プロフェッショナルシリーズ)

大学1年次でやった微積線形代数がバンバンたくさん出てきて、やっと今年やった数学のありがたみが少しわかった気がします。

理論を学んだら今度は実装をしようと「ゼロから作るDeepLearning」を少し読んだわけですが、どうしても横ベクトルが受け付けられず、結局自前で実装することに決めました。

今回は定番の手書き文字認識をやります。なのでソースコード自体はスッキリしていると思います。今度は畳み込み層にも対応させたいな…

以下詳細について書きますが、初学者ゆえ、間違いやお見苦しい点が多々あるかもしれません。もしお気づきになられたら指摘していただけると幸いです。

学習の大まかな流れとしては(「深層学習」 p50 - p52 参照)、まずは順伝播。各層の総入力Uと出力Zはresultとして記憶しておいて、次の誤差逆伝播法で誤差関数の重さについての勾配を求めるのに使います。

デルタは逐次的に求め、重さに関しての偏微分を全て求めて、勾配法で重みを学習していきます。

出力層のデルタは、出力層の活性化関数と誤差関数の種類によって計算方法が違うので、そこだけ個別に対応します。

手書き文字認識では、MNISTというデータセットが著名らしいので、今回はそれを使用。28 * 28 のモノクロ画像と数字 0〜9 を対応させたデータセットです。

今回、学習に用いたニューラルネットは3層からなり、以下のような構成です。

入力 第1層 第2層 第3層
ユニット数 28 * 28 200 40 10
活性化関数 - ランプ関数 ランプ関数 ソフトマックス関数

分類問題なので、誤差関数は交差エントロピーを使用します。学習係数は学習するにつれ線形に減少させます。

ソースコード

  • test.py
import numpy as np
import neitwork as nw
import copy
import matplotlib.pyplot as plt
import pickle

def read_image(path):
	with open(path, mode = "rb") as f:
		src = f.read()
		num = int.from_bytes(src[4 : 8], "big")
		img_size = int.from_bytes(src[8 : 12], "big") * int.from_bytes(src[12 : 16], "big")
		img = [ x / 255.0 for x in src[16 : 16 + num * img_size]]
		return np.array(img).reshape(-1, img_size).T

def read_label(path):
	with open(path, mode = "rb") as f:
		src = f.read()
		num = int.from_bytes(src[4 : 8], "big")
		L = np.zeros((10, num))
		for i in range(num):
			L[src[8 + i]][i] = 1.0
		return L

def accuracy(N, X, T):
	Y = nw.forward_all(N, X)
	Y = np.argmax(Y, axis = 0)
	T = np.argmax(T, axis = 0)
	return np.sum(Y == T) / X.shape[1]

print("loading...", flush = True)
# 各データセットのバイナリファイルのパス
Xtr = read_image("train-images-idx3-ubyte")
Ttr = read_label("train-labels-idx1-ubyte")
Xts = read_image("t10k-images-idx3-ubyte")
Tts = read_label("t10k-labels-idx1-ubyte")
print("completed.", flush = True)

iter_num = 100
call_learn_num = 50
batch_size = 100
tr_size = Xtr.shape[1]
init_learning_rate = 1.0
goal_learning_rate = 1e-2
learning_rate = lambda j : init_learning_rate - (j / call_learn_num) * (init_learning_rate - goal_learning_rate)

N = nw.make_network(unit_num = [28 * 28, 200, 40, 10], func = [nw.ramp, nw.ramp, nw.soft_max])

train_accuracy_list = []
test_accuracy_list = []

#100イテレーションごとに計測
for j in range(call_learn_num):
	N = nw.learn(N, Xtr, Ttr, nw.cross_entropy_error, lambda Y, T : Y - T, iter_num, batch_size, learning_rate(j))
	train_accuracy = accuracy(N, Xtr, Ttr)
	test_accuracy = accuracy(N, Xts, Tts)
	train_accuracy_list.append(train_accuracy)
	test_accuracy_list.append(test_accuracy)
	print(train_accuracy, test_accuracy, flush = True)

epoch_num = iter_num * batch_size * call_learn_num / tr_size
x_axis = np.linspace(epoch_num / call_learn_num, epoch_num, call_learn_num)
plt.ylim([0, 1])
plt.plot(x_axis, np.array(train_accuracy_list), label = "train accuracy")
plt.plot(x_axis, np.array(test_accuracy_list), label = "test accuracy")
plt.legend()
plt.xlabel("epochs")
plt.ylabel("accuracy")
plt.savefig("1.png")

with open("save.pkl", mode = "wb") as f:
	pickle.dump(N, f)
  • neitwork.py
import numpy as np
import copy

class layer:
	def __init__(self, shape):
		self.W = np.empty(shape)
		self.b = np.empty(shape[0])
		self.f = ident

class forward_result:
	def __init__(self, shape):
		self.U = np.empty(shape)
		self.Z = np.empty(shape)

def make_network(unit_num, func, init_std = 0.01):
	ret = []
	layer_num = len(unit_num) - 1
	for i in range(layer_num):
		l = layer((unit_num[i + 1], unit_num[i]))
		l.W = np.random.standard_normal(l.W.shape) * init_std
		l.b = np.zeros(l.b.shape)
		l.f = func[i]
		ret.append(l)

	return ret

def forward_layer(l, X):
	ret = forward_result((l.W.shape[0], X.shape[1]))
	ret.U = np.dot(l.W, X) + np.dot(np.c_[l.b], np.ones(X.shape[1]).reshape(1, X.shape[1]))
	ret.Z = l.f(ret.U)

	return ret

def forward_all(network, X):
	for l in network:
		X = forward_layer(l, X).Z

	return X

def forward_all_with_result(network, X):
	ret = []

	I = forward_result(X.shape)
	I.Z = X
	ret.append(I)

	for l in network:
		ret.append(forward_layer(l, ret[-1].Z))

	return ret

def grad_bp(network, result, delta):
	layer_num = len(network)
	N = result[0].U.shape[1]

	ret = []

	Z = result[-2].Z
	dL = layer(network[-1].W.shape)
	dL.W = np.dot(delta, Z.T) / N
	dL.b = np.dot(delta, np.c_[np.ones(N)]).reshape(delta.shape[0]) / N
	ret.append(dL)

	for i in reversed(range(layer_num - 1)): #今 0-origin で i 層目のデルタを計算する
		W = network[i + 1].W # i + 1 層目のW
		U = result[i + 1].U # i 層目のU
		df = derivative_dict[network[i].f] # i 層目の活性化関数の導関数
		delta = df(U) * np.dot(W.T, delta)
		Z = result[i].Z # i - 1 層目のZ
		dL = layer(network[i].W.shape)
		dL.W = np.dot(delta, Z.T) / N
		dL.b = np.dot(delta, np.c_[np.ones(N)]).reshape(delta.shape[0]) / N
		ret.insert(0, dL)

	return ret

def learn(network, X, T, loss_function, last_delta, iter_num, batch_size, learning_rate):
	N = X.shape[1] # データ数
	layer_num = len(network)

	network = copy.deepcopy(network)

	for i in range(iter_num):
		batch_mask = np.random.choice(N, batch_size)
		x_batch = X.T[batch_mask].T
		t_batch = T.T[batch_mask].T

		result = forward_all_with_result(network, x_batch)
		y_batch = result[-1].Z
		grad = grad_bp(network, result, last_delta(y_batch, t_batch))

		for l in range(layer_num):
			network[l].W -= learning_rate * grad[l].W
			network[l].b -= learning_rate * grad[l].b

	return network

def logistic_sigmoid(x):
	return 1 / (1 + np.exp(-x))
def logistic_sigmoid_d(x):
	return 1 / (np.exp(0.5 * x) + np.exp(-0.5 * x)) ** 2

def tanh(x):
	return np.tanh(x)
def tanh_d(x):
	return 1 / np.cosh(x) ** 2

def ramp(x):
	return np.maximum(0, x)
def ramp_d(x):
	return 1 * (x > 0)

def ident(x):
	return x
def ident_d(x):
	return np.ones_like(x)

def soft_max(a):
	c = np.max(a, axis = 0, keepdims = True)
	exp_a = np.exp(a - c)
	sum_exp_a = np.sum(exp_a, axis = 0, keepdims = True)
	return exp_a / sum_exp_a

def mean_squared_error(y, t):
	return 0.5 * np.sum((y - t) ** 2, axis = 0, keepdims = True)

def cross_entropy_error(y, t):
	delta = 1e-7
	return -np.sum(t * np.log(y + delta), axis = 0, keepdims = True)

derivative_dict = { logistic_sigmoid : logistic_sigmoid_d, tanh : tanh_d, ramp : ramp_d, ident : ident_d }

活性化関数とその1階導関数の対応づけをどのようにするか悩んだのですが、もっといい方法ありそうです。

以下学習結果。
f:id:uchifuji:20180212165521p:plain
テストデータでの正解率は 97.84% 。若干過適合気味なのかな?これ以上時間をかけてもテストデータに対しての正確性は上がらないみたいです。畳み込みニューラルネットなどを使うと99%の大台に乗せることも可能らしいです。ロマンがありますね。

さて、次はいつになるかわかりませんが、畳み込みを実装して画像認識にも挑戦できたらいいなーって思っています。

終わり

参考文献

  • 「深層学習」岡谷 貴之
  • 「ゼロから作るDeepLearning」斎藤 康毅

ABCに初参加した - 僕と競プロ

11/4の AtCoder Beginner Contest (ABC)に参加してきました。ABCはAtCoderが主催するプログラミングコンテストで、100分で4問、主に競プロ初心者を対象としています。

A問題

マス目の数が 2 x 3 で固定だったのでナイーブに実装しました。

#include <cstdio>
#include <cstdlib>

#define rep(i, a, b) for(int i = (a); i < (b); i++)
#define repr(i, a, b) for(int i = (a); i >= (b); i--)
#define range(x) begin(x), end(x)

using namespace std;

typedef long long ll;
typedef unsigned long long ull;
typedef unsigned int uint;

const int iinf = 1 << 29;
const ll linf = 1l << 61;

int main(int argc, char* argv[])
{
	char a[2][4];
	rep(i, 0, 2)
	{
		scanf("%s", a[i]);
	}
	if(a[0][0] == a[1][2] && a[1][1] == a[0][1] && a[0][2] == a[1][0])
	{
		printf("YES\n");
	}
	else
	{
		printf("NO\n");
	}
	return 0;
}

B問題

pow 使ってときました。模範解答は線形探索してたので、この解き方はちょっとずるかったかも。usingしてるのに、std::powとか頭悪いことしてます。

#include <cstdio>
#include <cmath>

using namespace std;

typedef long long ll;
typedef unsigned long long ull;
typedef unsigned int uint;

int main(int argc, char* argv[])
{
	ll N;
	scanf("%lld", &N);
	ll s = std::pow(N, 0.5);
	printf("%lld", s * s);
	return 0;
}

C問題

サイズの小さい順にソートして、累積和 + 動的計画法 でときました。ソートを除けばO(n)。
解説は中段を走査 + 上段・下段を2分探索で解いてました。よく考えたら祭壇の段数は3段で固定なので、別に難しいこと考える必要なかったんですね...。頭が硬いです。

#define rep(i, a, b) for(int i = (a); i < (b); i++)
#define repr(i, a, b) for(int i = (a); i >= (b); i--)
#define range(x) begin(x), end(x)

using namespace std;

typedef long long ll;
typedef unsigned long long ull;
typedef unsigned int uint;

const int iinf = 1 << 29;
const ll linf = 1l << 61;

int N;
int size[3][100800];
// i + 1 段目 j 番目のパーツを一番上にして、i + 1 段目、...、3段目のピースを使ってできる祭壇の種類の数を今 comb[i][j] とする。
// ただしcomb[i][0] = 0 (i = 0, 1, 2), comb[2][j] = 1 (j = 1, 2, ..., N)
// sum[i][j] = comb[i][0] + comb[i][1] + ... + comb[i][j] と定義する。
// size[i][j] < size[i + 1][k] を満たす最小のkをk(i, j)とすると、(満たすkがない時k(i, j) = N + 1)
// comb[i][j] = comb[i + 1][k(i, j)] + comb[i + 1][k(i, j) + 1] + ... + comb[i + 1][N] = sum[i + 1][N] - sum[i + 1][k(i, j) - 1]
// sum[i][j] = sum[i][j - 1] + comb[i][j]
// が成り立ち、結果
// sum[i][j] = sum[i][j - 1] + (sum[i + 1][N] - sum[i + 1][k(i, j) - 1])
// この漸化式を解いていけば良い(sum[0][N]が答えになる)
ll sum[3][100800];

int main(int argc, char* argv[])
{
	scanf("%d", &N);
	rep(j, 0, 3)
	{
		size[j][0] = 0;
		for(int i = 1; i <= N; i++)
		{
			scanf("%d", &size[j][i]);
		}
		sort(size[j] + 1, size[j] + N + 1);
	}

	for(int i = 2; i >= 0; i--)
	{
		int k = 1;
		sum[i][0] = 0;
		for(int j = 1; j <= N; j++)
		{
			if(i == 2)
			{
				sum[i][j] = j;
			}
			else
			{
				while(k <= N && size[i][j] >= size[i + 1][k]) k++;
				sum[i][j] = sum[i][j - 1] + (sum[i + 1][N] - sum[i + 1][k - 1]);
			}
		}
	}
	printf("%lld\n", sum[0][N]);

	return 0;
}

sumとsizeの添え字を1つずらしたせいで色々バグらせましたがなんとかAC。ただ、この解法なら任意の段数でも効率よく解けるはず。

D問題

意味不明でした。解説を見ても意味不明。グラフの問題に帰結させてる?01BFSってなんですか。

感想

楽しかったですが、C問題で、中段を中心に探索するのを思いつけなかったのが最大の心残りです。まだまだ頭が硬いですね...。

僕と競プロ

美談でもなく、競プロとの出会い、JOIが僕の人生に大きな影響を与えました。いずれ話したいなと思っていたことでしたが、今回この場を借りて語ろうと思います。



話は2年前のちょうどこの頃、高校2年の秋まで遡ります。地方の自称進学校でゴミみたいな成績に甘んじていた頃のことです。

中学生の時趣味でCを書いていたので、その惰性で高校でも細々とプログラミングをしていました。そのことを知った情報科目の教師が、JOI予選に出ることを勧めてきました。面白そうだったので参加することにしましたが、競技プログラミングというという単語を初めて聞いたレベルだったので全く期待はしていませんでした。

しかしなぜだか予選を突破してしまい、本選に出ることになりました。当時とても喜びましたが、その後すぐに、本選には進学校のパソコン部で鍛え上げられた強者がわらわら居るということを知って、戦意喪失。その代わり、プロと一度だけでも話してみたい、そう思うようになりました。

そして本選当日。案の定1完で敗退。夜の自己紹介で、約90人の前で「蟻本昨日買いました(本当)」という盛大なギャグをかましたところ「ひっでえなあww」という正直なコメントをもらったのをよく覚えています。そして、その懇親会で出会った情報のプロたちに圧倒されました。自分の遥か上には、頭が良くてかつ魅力的な人がたくさん居るということを、この時初めて知ったわけです。



しかしここで一つ、心の変化がありました。

「勉強は面白い。この人たちと一緒に勉強したい。いい大学に行きたい。」そう思うようになりました。それが、僕が受験勉強を続ける大きなモチベーションになったということは言うまでもありません。裏を返せば、もしもこの人たちに出会わなければ、僕は今の大学に進学していなかっただろうし、大学で勉強などほとんどしなかったことでしょう。



そしてもう一つ、考え方の変化がありました。

それまで僕は、プログラミングができることをあまり人に知られたくありませんでした。そう思っていた理由は今考えるとよくわかりませんが、多分変な奴だと思われたくなかったんだと思います。なので、学校でJOIについて表彰されることになった時、卑屈にも当時の僕はそれを拒否しました。当時の僕にとっては、表彰は全校生徒の前で腹踊りをすることに等しかったのです。しかしながら、先生に「いいから、いいから」と押し切られ、表彰されることになりました。

当日、なぜか自分はトリ。すでに帰りたさ度MAX。そのうち、とうとう僕の番が回ってきました。穴があったら入りたい。

校長「右は何たらかんたら〜〜〜代読。おめでとう。」
校長「みなさん何の事だか分からないと思うので、ちょっと情報オリンピックについて簡単に説明してください。」 え、説明しろなんて聞いてないんだが
僕 「えー...、与えられた問題をプログラムを組んで解く能力を競う、えっと...まあニッチな大会です。」 ちょっとウケた。よかった。
校長「えーみなさんよく分からないと思いますが、彼は灘や開成の生徒たちに混じって戦ってきました。これは非常に名誉なことです。どうぞ拍手を送ってください。」

友人曰く、僕の表彰が一番拍手が大きかったそうです。思ってもみませんでした。みんな「?」みたいな顔をしながらまばらな拍手が起こるのが関の山だと思っていました。

その後、友達やクラスメイト、先生から賞賛を受けて、「プログラミングは胸を張って言える特技の一つなんだ」と気づかされました。自分の好きな事、得意なことを知ってもらうのに、後ろめたさを感じる必要はないのです。それ以来、プログラミングに限らず、自分の特技や趣味を知ってもらうことに抵抗感を持たなくなりました。この考え方の変化は、僕にとって大きなプラスとなりました。それまでと比べて、胸のつかえが取れた気がしました。

本選を通して、僕が競プロ初心者以上の何かになれたかと言えばなれなかったのですが、その代わりに大きな刺激を受けることができました。本選に出場できたことを本当に幸福に思います。



現在、僕は大学で競技プログラミングの自主ゼミに参加し、競プロを一から学んでいます。他にも数学の自主ゼミにも楽しく参加させてもらっています。一緒に勉強できる仲間がいるということは素晴らしいことです。しかし、2年前もし考え方を変えていなかったら、一緒に学ぶ仲間を作ろうともせず、一人細々と趣味に興じていたことでしょう。自分よりも遥か上の存在を知って、刺激を受ける。これが何事においても向上の近道となるわけですね。

区間スケジューリング問題

蟻本p43。社畜の問題です。

問題

{N} 個の仕事 {w_i\ (i=1, 2, \cdots, N)} がある。各仕事 {w_i} の開始時刻は {s_i} 、終了時刻は {t_i} である。あなたは各仕事について参加するか参加しないかを選ばなければならない。ある仕事に参加するならば、その仕事に始めから終わりまで参加しなくてはならない。すなわち、他に参加する仕事と時間か重なってはならない。開始・終了の瞬間だけが重なるのも許されない。

参加できる仕事の数の最大値を求めよ。

制約

  • {1 \leq N \leq 10^5 }
  • {1 \leq s_i < t_i \leq 10^9 }

例えば次のような仕事の集合が入力として与えられた時、
f:id:uchifuji:20171024025451p:plain
最適なスケジュール(のうちの一つ)は赤の仕事。つまり解は4。
f:id:uchifuji:20171024025445p:plain

頭の悪い人

f:id:uchifuji:20171024030424j:plain
計算量のオーダーが { O(2^n) }なのでだめ。

頭のいい人

f:id:uchifuji:20171024030920j:plain
終了時刻の早い順に貪欲に仕事を選んでいくけば最適解が求まります。計算量は平均で {O(n \log n)}。というかソートアルゴリズムに依ります。

開始時刻の早い順に貪欲に選ぶのは不正解です。
次のような仕事の場合、間違った解を導きだします。
f:id:uchifuji:20171024031738p:plain

そうパッと思いつく解法ではないですね。蟻本に証明の概略が載っていましたがいまいちピンとこないので自分で考えます。

証明

アルゴリズム {A} : 選べる仕事(今まで選んだ仕事と重なっていない仕事)の中で終了時刻が最も早いものを選ぶ。これを繰り返す。

与えられた {N} 個の仕事の集合 {W = \{ w_1, w_2, \cdots, w_N \} } に対してアルゴリズム {A} が導きだしたスケジュールに含まれる仕事の個数を {M} 、スケジュールに含まれる仕事の集合を {W' = \{ w_{k_1}, w_{k_2}, \cdots, w_{k_M} \} } とする。ただし { t_{k_1} < t_{k_2} < \cdots < t_{k_M} }{ a = \min\{ s_1, s_2, \cdots, s_N \}, b = \max\{ t_1, t_2, \cdots, t_N \} } とする。

ある仕事 {w_i \in W} があって、適当な仕事 {w_{k_j} \in W' }{W'} から選べば { t_{k_j - 1} < s_i < t_i < t_{k_j} } となる...(1) 、又は { a \leq s_i
 < t_i < t_{k_1} \lor t_{k_M} < s_i < t_i \leq b } ...(2) となると仮定する。

(1)の場合
アルゴリズム {A}{ w_{k_{j-1}} } が選ばれた時点で、開始時刻が { t_{k_{j-1}} } よりも遅い仕事を選ぶことができることは保証される。よって {w_i} は選ぶことができる仕事である。{ t_i < t_{k_j} } であるから、この時アルゴリズム {A} は仕事 {w_{k_j}} ではなく仕事 {w_{i}} を選んでいるはずである。

(2)の場合
(1)と同様の理由で、{w_i}アルゴリズム {A} によって選ばれているはずだが実際には選ばれていない。

いずれにせよ矛盾が発生するから、仮定が間違っていることがわかる。

これを言い換えると、任意の仕事 {w_i \in W} に対し、{ s_i \leq t_{k_j} \leq t_i } となるような仕事 {w_{k_j} \in W'} が存在する。

任意の合法な(時間のかぶる仕事が存在しない)仕事スケジュールを考える。その仕事スケジュールに含まれる仕事の個数を {M''} とし、このスケジュールに含まれる仕事の集合を { W'' = \{ w_{l_1}, w_{l_2}, \cdots, w_{l_{M''}} \} } とする。{ i = 1, 2, \cdots, M'' } に対し { s_{l_i} \leq t_{k_j} \leq t_{l_i} } を満たす {j} のうちの一つを {u_i} とおく。({1 \leq u_i \leq M}) このスケジュールは合法であるから、{u_1, u_2, \cdots, u_{M''}} は互いに異なる。鳩の巣原理から、{M'' \leq M} である。

より、任意の合法なスケジュールの仕事の個数は {M} を超えない。よってアルゴリズム {A} は正しい。

簡単なことを言っているのに

改めて読むと物凄く分かりにくい気がします。説明が下手なのか、はたまたアルゴリズムを数学の言葉で説明すること自体が難しいのか。図を使って説明した方がよかったと今更ながら思います。

急性腸炎にかかった

急に寒くなったと思えば真夏の天気になったりと、体に優しくないこの頃。風邪が流行っていますね。僕の周りでも数人がダウンしている模様です。

そんな中、空気を読まずに腸炎にかかりました。俗に胃腸風邪とか、ノロウイルスとか言われるやつです。胃腸炎という病気は胃腸に炎症が起こる病気全般を指すっていう中々にガバガバ定義で、原因は大きく分けて2つ、細菌が原因のものとウイルスが原因のものがあるそうです。細菌が引き起こす胃腸炎で、食事が原因になるものは食中毒って言われます。O-157とか有名。原因となるウイルスもいくつか種類があって、ノロウイルスアデノウイルス、他にもたくさんあります。

ウイルスごとの特効薬とかないんで、どのウイルスであっても治療法は変わらないらしいです。なのでどのウイルスが原因か特定する必要はないと。詳しい検査もされませんでした。整腸剤を飲みながら自然治癒を待つしかないらしいです。

同じ病気にかかった人の参考になればということで、病状の経過を記したいと思います。

10/6(金)

サークルに行きました。いつもよりも疲労感が溜まって「ちょっといつもと違うな」と感じました。

10/7(土)

朝起きたら鈍い腹痛が。眠気と倦怠感のせいで、なかなかベッドから起きられませんでした。この時点ではただの腹痛だろうと思っていたのでバイトに行きましたが、いざ働き出すと腹痛と倦怠感で死にそうになりました。4時間働いたのち、上司に言って早めに上がらせてもらいましたが、家に帰るのも精一杯な状態。

ここでカレーを食べるという悪手を打ちました。あとで知りましたが、カレーは胃腸に負担がかかるので胃腸炎の時は食べてはいけないらしいです。

家について秒で眠りに落ちて、3時間後に起床。熱を測ると39.7度。これはまずいことになったと思い、この時間でもやってる病院を探して診察を受けました。
しかし症状を話しても「よくわかんない」と言われる始末。解熱剤と整腸剤を処方されました。

10/8(日)

解熱剤が効いたのか38.5度まで熱が下がりました。何食べても下痢で出て辛いので、食事の量を絞り始めました。夜中に腹痛で叩き起こされた上に吐き気まで出てきました。

10/9(月)

日中は平熱まで下がりましたが下痢が続きます。もう出すものもないはずなのに腹痛と水っぽい下痢だけは続くんですよね。夕方から微熱。

10/10(火)

症状変わらず。夜に吐き気が出ました。絶食を始めました。

10/11(水)

近くの総合病院の外来に行って、血液検査とエコー検査、尿検査を受けました。血液検査の結果は、CRP定量が3.5(基準0.3)と高値。これは炎症反応の度合いを示す度数で、この数値が高いと体のどこかで炎症が起きていることがわかります。「小腸のリンパ節が腫れています。ウイルス性の腸炎でしょう。他の臓器に異常はありません。」との診断でした。土曜日に行った病院でもらった整腸剤が効いてる気がしないので、違う種類の整腸剤を処方してもらいました。

夜、とうとう絶食に耐えられなくなり天一のチャーハンを食べたら案の定死にました。

10/12(木)

ようやく熱が下がりました。下痢も一番ひどかった時期に比べれば少なくなりましたね。

10/13(金)

やっと下痢が止まりました。大学にも少し行きました。全快です。

本来2、3日で治る病気らしいんですけど、多分バイトとカレーのせいで引きずった気がします。バイトやめたいです。

std::arrayの初期化

次のコードをclangで-Wallでコンパイルすると警告が出ます。gccでは出ないです。

clangのバージョン : Apple LLVM version 9.0.0 (clang-900.0.37)

[test1.cpp]

#include <array>
#include <iostream>

int main(int argc, char* argv[])
{
	std::array<int, 4> ar = { 0, 1, 2, 3 };

	for(int i : ar)
	{
		std::cout << i << std::endl;
	}

	return 0;
}

[コンパイル]


clang++ test1.cpp -std=c++11 -Wall -o test1.out

test1.cpp:6:28: warning: suggest braces around initialization of subobject [-Wmissing-braces]
std::array ar = { 0, 1, 2, 3 };
^~~~~~~~~~
{ }
1 warning generated.

中括弧の数が足りないと言っています。次のコードでは警告は出ません。

#include <array>
#include <iostream>

int main(int argc, char* argv[])
{
	std::array<int, 4> ar = { { 0, 1, 2, 3 } };

	for(int i : ar)
	{
		std::cout << i << std::endl;
	}

	return 0;
}

std::arrayの実装を読むと、かなり大雑把に言えば次のようになっています。

template<typename _Tp, std::size_t _Nm>
struct array
{
    _Tp _M_elems[Nm];
    //以下メンバ関数
    //......
};

std::arrayには、publicなメンバ変数である_Tp型配列以外のメンバ変数がありません。さらにaggregateの条件を満たしているため、通常の配列と同様の初期化の仕方が可能になります。ただし、std::arrayは配列そのものではなくて配列を唯一メンバ変数に持つ構造体なので、括弧が二重に必要です。

つまり、一重括弧の方が文法的にグレーなのであって、二重括弧の方が本来正しい記述の仕方というわけです。しかしstd::arrayは_M_elemsよりも前にメンバ変数を持たないため、一重括弧でも意図した通りに動くというわけです。

結論としては、二重括弧で書いた方がいいと思います。

余談ですが、std::arrayでも生配列と同じように0初期化が可能です。

std::array<int, 4> ar = {};
std::array<int, 4> ar = { 0 };

clangでは、前者だとワーニングは出ませんが後者だと出ます。

二次元配列でも同じように初期化が可能です。{}内の要素が足りない場合、
・数値なら0で初期化
・ポインタならnullptrで初期化
・aggregateならこのルールを再帰的に適用
というルールで初期化を行うからです。次のコードも全て0で初期化できますね。

#include <array>
#include <iostream>

struct vec
{
	int x;
	int y;
};

int main(int argc, char* argv[])
{
	std::array<std::array<vec, 4>, 4> ar = {};

	for(const auto& line : ar)
	{
		for(const auto& cell : line)
		{
			std::cout << "(" << cell.x << ", " << cell.y << "), ";
		}
		std::cout << std::endl;
	}

	return 0;
}

[実行結果]
(0, 0), (0, 0), (0, 0), (0, 0),
(0, 0), (0, 0), (0, 0), (0, 0),
(0, 0), (0, 0), (0, 0), (0, 0),
(0, 0), (0, 0), (0, 0), (0, 0),