Macでcupyをビルドする

えー。また競プロとは全然関係のない記事です。同じことをやっていそうな人が少なくてエラーが出てこなかったのでまとめます。

手順1: CUDA対応MacBook Proを購入する

えー。MacBook Mid 2012-2014 15 inchしか対応していないと思います。しかも15 inchじゃないとだめです。頑張りましょう

手順2: 古いCommandLineToolsを入れる

えー。最近(Xcode 8以降)のだと、ld: library not found for -lstdc++.hみたいなエラーを吐かれるはずなので、2016年ぐらいのCommandLineToolsを入れます。https://developer.apple.com/downloadからダウンロードできます。

手順3: CUDA Toolkit, cudnnenv, nccl, cutensorを入れてビルド

ここはやるだけ。cutensorは即時ダウンロードできないので、リクエストが承認されるまで数日待ちそう。

ビルドはcloneしてきたcupyのディレクトリでpip install -e .をする。エラーメッセージに従ってcudnnとかのパスをちょっとずつ調整していけばよさそう。

自分用Makefileのメモ書き

基本構文

<target>: <files>
    <commands>
<target>:
    <commands>

実行

makeとやると、Makefileに定義された最初のターゲットが実行される。他のターゲットを実行するにはmake clean, make installなど臨機応変にターゲット名を入力。

特別なターゲット名

.PHONY

.PHONY test

とかよくみる。コマンド名とファイル名が重複するとき通常ファイルが優先されて実行されないが.PHONYを書いておくと実行ルーチンが作られるっぽい。

.PRECIOUS

このターゲットは処理が中断してもファイル等が削除されない。

ここが詳しい。 www.ecoop.net

特殊変数

詳しくはhttps://www.gnu.org/software/make/manual/make.html#Automatic-Variables

$@

targetを返す

$<

filesの最初の要素を返す

$^

filesのすべての要素を返す

まだまだあるので適宜上のを見ること

変数・関数

こんなん

cpp_files = $(shell $(FIND) src/ -name "*.cpp" -printf "%P\n")
cxx_obj_files = $(subst .cpp,.o,$(cpp_files))
obj_build_root = $(build_root)/objs
objs = $(addprefix $(obj_build_root)/cxx/,$(cxx_obj_files))

みたいにするだけ。substとかaddprefixとかは後述する特殊関数

特殊関数

多すぎ。ここが天才 qiita.com

自分で使うならともかく読む分には比較的読める

ワイルドカードとしての%

suffixの置換

SRC = main.cpp hello.cpp
OBJ = $(SRC:%.cpp=%.o)

targetとかにも使えるっぽい

$(obj_build_root)/cxx/%.o: src/%.cpp
    $(CXX) $(CXXFLAGS) -MMD -c -o $@ $(filter %.cpp, $^)

など。

リンク

この記事よりもうちょい詳しめ www.jsk.t.u-tokyo.ac.jp

まとめ

Makefile系イキリオタク滅びろ

AtCoderで黄色になりました

伝統に従って記事を書きます。

やったこと

AtCoder, CodeForces, Topcoderのコンテストに(撤退することもあるけどとりあえず)出続ける

短期間で圧倒的精進する人にはやっぱり抜かれてしまうけどこれでレートが落ちることは(赤とかにならない限り)ないかなと思います。

青序盤との差分

数え上げの苦手意識が多少軽減した

いやまだ苦手なんですが。でも、数え上げはわりと1ステップ1ステップごとに考察を詰めやすくて構築みたいに一発天才ゲーではなく素直に考察すれば解けることが多いことに気づきました。橙目指すなら数え上げかなりうまくなる必要がありそうだと思う。

典型DPに慣れた

Earsのような状態変化を記録するタイプのDPは青帯の人にとってはわりとつまるところだと思う。これに慣れたのはでかい。挿入DP, 区間DP,  O(3^n)のDP,  O(N^2)の木DP, 約数の包除原理など。

知識

これは列挙しようと思ったんですがさすがに多すぎて思い出しきれないのでやめました。比較的最近勉強したHeavy-Light Decompositionはあまりに便利なので青の人でも一度勉強してみることをおすすめします。

ゴリ押し用のデータ構造(Implicit Treap: 要素の削除や挿入ができる遅延セグ木のようなもの, や要素の削除や挿入ができる上位k個要素のsumを返すデータ構造など)をいくつか常備していますが、2, 3ヶ月に一度ぐらい役に立つ印象です。

ConvexHullTrickとかFFT, NTTなども使えるようになりました。

あと問題を解くと、こういうふうに帰着しちゃうとヤバいみたいなのがいくつか経験からわかってきて余計な考察を省けたりします。

当面の目標

コンテスト中に1時間ぐらい同じ問題を考察してやっと解けた、みたいな経験が最近少なくて(ICPCとかは実装担当なのでほぼ考察をしていないし)、長時間系のコンテストが苦手になっている気がします。黄色になっちゃうとレート上げチャンスがAGCや企業コンになりますが、前者は長時間系コンテスト(Aはすぐ解けちゃうし、たいていCやBに各1時間以上費やすことが多そう)の代表なので考えることに慣れなきゃいけないなと思います。これなにをすればいいんでしょうかね、やはり精進しかなさそう

ICPC国内予選2019に出場しました

4完24位でした。 f:id:xuzijian629:20190712220639p:plain

感想

この一戦だけに絞って感想を言うなら、「Eをもう少し落ち着いて確実に合わせに行きたかった」となるのですが、さすがに今回のICPCは自分にとって最後のICPCであり、チームGirigiriとはこの一戦のみならず人としての付き合いが長いのでそれほどシンプルにまとめることができません。

やっぱりチームメンバーに一番伝えたい気持ちは感謝の気持ちです。team GIrigiriが競プロモードになったのは去年のICPC国内予選以来ですが、わりとそこから1年間、競プロ熱を一度も切らすことなく切磋琢磨してやってこれたのはチーム独特の雰囲気があったからだと思います。予選に通ることが人生の目標なら別のチームメンバーもありえたかもしれませんが、自分が最大限に成長できたのはやっぱりこのチームだったんじゃないかと思っています。

そんな大好きなチームメートと出た最後のICPCなので(海外Regionalは一応考えてはいますが、まあ実質最後なので)残念だと言って終わらせたくないです。正直レート的に見るとまあ妥当なラインだとは思いますし、ABCDノーペナなのは普通にすごいと思います。Dでぼくが初め考えた嘘解法の穴をいち早く指摘してくれた荻野やDのアルゴリズムを理解して代わりに実装してくれたけんしんはめちゃめちゃ仕事をしたと思うし、中の中の上ぐらいの出来なんじゃないでしょうか。去年は50位だったのでそれにくらべれば大進歩です。

ICPCは終わりましたがGirigiriは研究なりその他の競プロなりでこれからも長い付き合いになりそうなので、2人ともよろしくお願いします。

最後になりましたが、通過した東大チーム・他大チームの方々はアジア地区予選頑張ってください!!!!!

お疲れ様でした!!!

ライブラリ

team Girigiriのライブラリを公開しておきます。ご自由にお使いください。幾何がちょっと弱いので補強したほうがいいです drive.google.com

風景

f:id:xuzijian629:20190712223624j:plain

こんな感じの長机でやっていました。

機械学習におけるカーネルを一言で

2つのデータがどれだけ似ているかを返す関数です。

非負実数を返し、似ていれば値は大きく、似ていないと小さいです。

以上!

カーネルトリックについて

たいてい回帰問題の文脈で現れると思いますが、元のデータ点を直線(もしくは超平面)で回帰することはそう一般にうまくいきません。 そのため、ある関数 \phiによって元のデータを高次元空間にマップしてから、超平面で回帰することを考えます。このときに適切な \phiを知ることは至難の業ですが、実は \phiを直接知らなくとも、さまざまな計算をうまく変形して行うと、 \phi(x), \phi(x')内積 \phi(x)^T \phi(x')の形とは限らず、より一般に \phi(x)^T A\phi(x')の形で現れます)さえわかれば多くの場合で事足りることがわかります。そして、 \phiを定める代わりにこの内積関数をカーネルとして予め定めてあげることで、 \phiを陽に扱わずこの問題が解けることになります。これがカーネルトリックです。

最後に

なんで本の説明はどれを見てもあんなに長々しいんや

機械学習のためのGaussian Processを簡潔に

ガウス過程 (Gaussian Process; GP)を簡潔に説明します。日本語の類書や類ウェブページ?はたくさんありますが以下の差別化を図りました

  • できるだけ数式による説明を心がけます
  • 定義をスムーズにするために関数空間の観点からGPによる回帰を説明し、実は線形モデルの非線形空間への拡張でも同じ式が得られることを示します(多くの説明だと逆?)

GPは○○の性質をもつ!みたいなのをひたすら書いている本を読んだのですが、だいぶ読みにくかったので簡潔さを心がけました。

すべての説明は

http://www.gaussianprocess.org/gpml/chapters/RW.pdf

の第2章に基づいていますが、導入の順序が異なります。

定義

関数 f\colon \mathcal{X} = \mathbb{R}^d \rightarrow \mathbb{R}がGaussian Process  \mathrm{GP}(\mu, k) ( \mu\colon \mathcal{X} \rightarrow \mathbb{R},  k\colon \mathcal{X} \times \mathcal{X} \rightarrow \mathbb{R}_+)に従うとは、任意の有限個の点集合 x_1, \ldots, x_n \in \mathcal{X}について、  f(x_1), \ldots, f(x_n)が多次元正規分布 \mathcal{N}(\mu(\mathbb{x}), k(\mathbb{x}))に従うことをいう。ここで、\mu(\mathbb{x})は平均を表す n次元ベクトルで i番目の要素が \mu(x_i)になってるものであり、 k(\mathbb{x})は共分散を表す n次正方行列で (i, j)要素が k(x_i, x_j)になっているものである。

普通、共分散はデータが与えられて、そこから平均を求めて計算されるものだと思いますが、今の文脈では、共分散関数(カーネルともいう)が予め定められています。つまり、2つの xが与えられるとそれらの値 f(x)は、 f(x)の共分散が共分散関数に従うように分布する、と考えます。

イメージと補足

定義は実はこれだけです。 f(x_1), \ldots, f(x_n)が多次元正規分布に従うというのがあまりにGPの本質なのですが、この多次元正規分布としてよくあるような二次元正規分布のような平面図(もしくは立体図)をイメージしてしまうと回帰にどう役立つんだ???という気持ちになってしまうと思います。いま、観測値がわかっているデータ点が n個あるとすると、この多次元正規分布 n次元になっています。 n次元空間をイメージする代わりに、この n点の近くには存在しやすく、点から遠い場所では存在しにくいような分布を \mathcal{X} = \mathbb{R}^dについて想像してみましょう。簡単な場合として、 \mathcal{X} = \mathbb{R}としますと、以下のような図を想像できるかなと思います(回帰でよく見そうな形をしていますね!)

f:id:xuzijian629:20190621021035p:plain (図は上記pdfから)

ノイズ付きのデータからの推論

 x_1, \ldots, x_n \in \mathcal{X}の各点について、 y_i = f(x_i) + \varepsilonに従うノイズが加わったデータが観測されており、未知のデータ点 x_1^*, \ldots, x_m^*についての予測値 f(x_1^*), \ldots, f(x_m^*)を得たいとします。

データ点 \mathbb{y} = (y_1, \ldots, y_n)^T,  \mathbb{f}^* = (f(x_1^*), \ldots, f(x_m^*))^Tが、GPに従う関数 fによって生成されているとします。いま、 fに関する事前知識がないため、平均が \mathbb{0}の関数を考えます。このとき、 \mathbb{y}, \mathbb{f}^*は次の事前分布に従います。

f:id:xuzijian629:20190621022842p:plain

 \sigma_nはノイズ \varepsilonの分散で、ノイズはガウシアンだとします。ガウス分布の共分散行列は n + m次正方行列で、 (i, j)要素が k(\cdot, \cdot)なのですが、注意すべきは左上の n \times n要素部分で、この部分の対角成分には、ノイズによる分散が入ります。なぜなら、異なるデータ間におけるノイズによる共分散は0ですが、自分自身との間ではそうはいかないからです。右上、左下のブロック行列はそれぞれ、サイズ n \times m,  m \times nの行列です。

さて、これを鬼のように変形すると、 \mathbb{f}^*の分布が次のように求まります(だいぶややこしいですが省略します)。

f:id:xuzijian629:20190622045425p:plain

これで、 \mathbb{f}^*を推定することができました。この推定は、最尤推定と異なり、最頻値(この場合平均に一致)だけではなく、分散も同時に推定可能です。

線形モデルからの拡張に一致すること

GPの説明にはいつもこの話がつきまとうかなと思います。上述した話だけで本来完結するのですが、せっかくなのでこちらも説明しておきます。

線形モデルとは、入力 x \in \mathbb{R}^dに対し、出力 y \in \mathbb{R} f(x) = x^T wによって予測するモデルで、パラメータ wを学習します。

実際のデータはこれに分散 \sigma^2ガウス分布に従うノイズが加わって生成されている(つまり y = f(x) + \varepsilon)とすると、データ \mathbb{y} X = (x_1, \ldots, x_n), wから得られる条件付き確率は

f:id:xuzijian629:20190622050625p:plain

となります。最尤推定だと、これを最大化するようなパラメータ wを求めて終了ですが、いま、Bayesianの考え方をします。すなわち、もともと、パラメータが w \sim \mathcal{N}(0, \Sigma_p)という事前分布に従っていた場合、 \mathbb{y}が得られたあとの wの分布の変化はBayesの法則を使って鬼のように計算をすると

f:id:xuzijian629:20190622051009p:plain

となります(再び省略)。線形モデルで推論するときも、最尤推定とは違い、あらゆる wの可能性について、その重みをとって出力を推定します。このとき出力は

f:id:xuzijian629:20190622051929p:plain

となります。

さて、線形モデルだけではやはり性能に限界があります。というのも、データ点が超平面によって二分されるような場合しかうまく動かないです。そこで、各データ点をある関数 \phi \colon \mathbb{R}^d \rightarrow \mathbb{R}^{d'}によって高次元空間上にマップしてから、その高次元空間上で線形モデルを動かす、という方法が考えられます。すなわち f(x) = \phi(x)^T wによって出力を予想します。ここでの重み wは線形モデルでは d次元でしたが、今は高次元空間の次元 d'となっています。線形モデルのときと全く同様の議論をすることで、推論の出力は

f:id:xuzijian629:20190622051734p:plain

となります。これをさらに、変形すると(難しい)

f:id:xuzijian629:20190622052620p:plain

となりますが、これは、前のセクションの最後で見た、関数空間による定義から導かれた推論の結果に一致しています。関数空間の定義で指定されるべきカーネル関数 k\colon \mathcal{X} \times \mathcal{X} \rightarrow \mathbb{R}_+と、このセクションで指定されるべきマップ関数 \phi \colon \mathbb{R}^d \rightarrow \mathbb{R}^{d'}&  \Sigma_pの間には深いつながりがあります。具体的にはpdfの12ページ。kernel / kernel trickを御覧ください。

まとめ

正直このページでは簡潔に説明しすぎていて(式変形とかも飛ばしているし)、このページだけを読んで理解することは難しいと思います。上述したpdfのchapter 2と並列して読むことをすすめます。1, 2時間程度で理解できるのではないでしょうか。

おもしろい数え上げの問題

知り合いからおもしろい問題を聞いたので解いてみました。おもしろいかは人によると思いますが、数え上げ方がおもしろかったのでまとめます。

タイトルを解法のキーワードにしようと思ったんですが非常にもったいないネタバレをしてしまいそうなので先に問題を紹介します。自分で考えてみたい人は一段落してからまた見に来てください。

問題

黒い石がN個、白い石がM個あり、これをいくつかの(nonemptyな)山に分けます。2つの分け方が異なるとは山のマルチセットが異なることをいいます(山の中での黒い石と白い石の順序は考えます)。何通りの分け方がありますか。ただし答えは非常に大きくなることがあるので1000003で割ったあまりを求めてください。

制約

  • 0 <= N, M <= 50
  • N + M >= 1

もともと多項式時間で数え上げられるかな〜と思ってたところを多項式時間でできたのでオーダー正直なんでもいいやという感じです。ぼくの解法だとNとMが同じ程度として、 O(N^5 \log N)ですが、本気を出したら縮むかもしれません。最大ケースで0.9秒でした。

入力

N M

サンプル

Input

3 0

Output

3

黒い石を0, 白い石を1で表すことにすると、山の分ける方法は[0,0,0], [0,00], [000]の3通りです。[0,00]と[00,0]はマルチセットとして等しくなることに注意してください。

Input

3 3

Output

131

131通りあります。これらには[001,110]や[010,110]が含まれますが、上述したように[001,110]と[110,001]は区別されないことに注意してください。

解法

== ネタバレ回避用スペース ==

== ネタバレ回避用スペース ==

== ネタバレ回避用スペース ==

== ネタバレ回避用スペース ==

== ネタバレ回避用スペース ==

== ネタバレ回避用スペース ==

== ネタバレ回避用スペース ==

== ネタバレ回避用スペース ==

== ネタバレ回避用スペース ==

== ネタバレ回避用スペース ==

== ネタバレ回避用スペース ==

== ネタバレ回避用スペース ==

== ネタバレ回避用スペース ==

== ネタバレ回避用スペース ==

== ネタバレ回避用スペース ==

== ネタバレ回避用スペース ==

== ネタバレ回避用スペース ==

== ネタバレ回避用スペース ==

解法

まず重複なく数え上げることを目指します。こういう数え上げは、(挿入DPみたいに)大きい順に数えるというのが典型な気がするのでそれに従います。すなわち、高さhの山をhが大きい順に確定していきます。これから扱う山の高さをh-1以下に制限することで、高さhの山の集合を独立に考えることができます。

数え上げる関数をdfs(n, m, h, k) := n個の黒い石とm個の白い石があり、高さhの山をk個つくり、残りの山は高さh-1以下にする場合の数、と定義すると、とりあえず答えはdfs(n, m, n + m, 1) + dfs(n, m, n + m, 0)になります。

関数dfsの擬似コードを以下に示します(再帰の終了条件はちょっと複雑なのであえて省いています。詳細は最後のコードを参照)

dfs(n, m, h, k)
  ret = 0
  N = h * k    // 注目する石の数
  for (i = 0; i <= N; i++)
    if (i <= n && N - i <= m)    // 黒い石をi個使用する
      w = N個から黒い石をi個選んで高さhの山をつくる場合の数
      for (d = 0; d <= (n + m - N) / (h - 1); d++)    // 高さh-1の山をdを作る。dは最大で(n+m-N)/(h-1)個
        ret += w * dfs(n - i, m - (N - i), h - 1, d)
  return ret

さて、あとは途中に書いた「N個から黒い石をi個選んで高さhの山をつくる場合の数」を計算できればおしまいです。

再び、これを重複なく数えることを考えます。このままだとちょっと考えにくいので次のように言い換えてみましょう。

問題

n, m, kはこの問題で新しく定義した変数とします。 n (=h)個の石からなる山がm (=N/h)個ある。 このうちk個の石が黒いです。山をマルチセットとして見たとき、異なるマルチセットは何通りあるか。

これを重複なく数える際にも、先程と同じ考え方ができます。すなわち、ある量の大きい順に確定させていくという考え方です。今の場合、山の中の黒い石の数を大きい順に確定させていきながら考えます。くどいようですが、黒い石がj個含まれている山がs個あるとし、残りの山はすべて黒い石がj-1個以下だとすると、このs個の山は独立に考えることができます。挿入DPですね。

計算の流れは先程と同様なので詳細は省きますが、確定させた部分の計算で次のような問題を考えることになります。

問題

n個の石からなる山がs個あり、それぞれの山にはちょうどj個黒い石がある。山をマルチセットとみたときに何通りのマルチセットがあるか

n個のうちj個が黒い石であるような山はnCj通りあります。これらnCj通りから重複をゆるしてs個選ぶ選び方なのでこの場合の数は(nCj)Hs通りになります。 これらは事前計算することでO(1)で求まるので結局併せて以下のようなプログラムで解くことができます。

コード

constexpr int MODULO = 1000003;
constexpr int MAX = 50;

mint fact[MODULO];
mint factinv[MODULO];
mint dp[MAX + 1][MAX + 1][2 * MAX + 1][MAX + 1], dp2[MAX * MAX + 1][MAX + 1][MAX + 1][MAX + 1];
bool visited[MAX + 1][MAX + 1][2 * MAX + 1][MAX + 1], visited2[MAX * MAX + 1][MAX + 1][MAX + 1][MAX + 1];

void init() {
    fact[0] = 1;
    for (int i = 1; i < MODULO; i++) {
        fact[i] = fact[i - 1] * i;
    }
    factinv[MODULO - 1] = fact[MODULO - 1].inv();
    for (int i = MODULO - 2; i >= 0; i--) {
        factinv[i] = factinv[i + 1] * (i + 1);
    }
    assert(factinv[0] == 1);
}

mint ncr(int n, int r) {
    if (n < r) return 0;
    return fact[n] * factinv[r] * factinv[n - r];
}

mint dfs2(int n, int m, int k, int j) {
    if (j == 0) return k == 0;
    if (j * m < k) return 0;
    assert(n <= MAX * MAX && m <= MAX && k <= MAX && j <= MAX);
    if (visited2[n][m][k][j]) return dp2[n][m][k][j];
    mint ret = 0;
    if (k == 0) return 1;
    for (int s = k / j; s >= 0; s--) {
        mint q = ncr(n, j);
        ret += ncr(q.x + s - 1, s) * dfs2(n, m - s, k - j * s, j - 1);
    }
    visited2[n][m][k][j] = true;
    return dp2[n][m][k][j] = ret;
}

mint cnt(int n, int m, int k) {
    assert(k <= n * m);
    return dfs2(n, m, k, min(n, k));
}

mint dfs(int n, int m, int h, int k) {
    if (n + m == 0) return 1;
    if (h == 1) return k == n + m;
    assert(n <= MAX && m <= MAX && h <= 2 * MAX && k <= MAX);
    if (visited[n][m][h][k]) return dp[n][m][h][k];
    int N = h * k;
    mint ret = 0;
    for (int i = 0; i <= N; i++) {
        if (i <= n && N - i <= m) {
            mint w = cnt(h, N / h, i);
            for (int d = 0; d <= (n + m - N) / (h - 1); d++) {
                ret += w * dfs(n - i, m - (N - i), h - 1, d);
            }
        }
    }
    visited[n][m][h][k] = true;
    return dp[n][m][h][k] = ret;
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(0);
    init();

    int n, m;
    cin >> n >> m;
    assert(n <= MAX && m <= MAX);
    cout << dfs(n, m, n + m, 1) + dfs(n, m, n + m, 0) << endl;
}

ダブル挿入DP

挿入DPの中で挿入DPを使う問題を初めて見た。うむ、やっぱあまりに汎用性が高いので僕が知らなかっただけで常識なのではないかと思い始めた。

プロたちによる意見

  • 山を黒と白の個数で適当に順序付ける

dp[a][b][c][d] := 黒a白bの山まで考えたときに合計で黒c白dにする通り数

とすると素直に Ο(N^5)になりそう

(by @noshi91)