技術書典15に申し込みます

技術書典15にサークルで申し込みをします

見出しの通りです。「マジックタンノリターンズ~タンノはもはやユビキタス~」というサークル名で 技術書典15に申し込みます。以下がサークルの公式Twitterアカウントです。
twitter.com

メンバー

ハンドルネーム (Twitterアカウント)

Yoka (Yoka (@Kalmia24879903) / X)

人参食べたい (人参食べたい (@A_kyoutoyenamta) / X)

初学者リョウ (最強つちだズ (@syogakusyaryou) / X)

ゆざきりつ (今回は不参加)

技術書典とは

技術書典とは、自分で書いた技術書を頒布・販売するイベントです。言うなれば、技術書版コミケです。オンラインとオフラインの両方の形式で開催され、オフライン会場であれば今年は2023/11/12 (日) に池袋サンシャインシティ 展示ホールD(文化会館ビル2F)で開催されます。楽しそうですね。オンラインの場合は、2023/11/11 (土) 〜2023/11/25(日)の期間内であればいつでもやっており、技術書典 :技術書のオンラインマーケット開催中から電子書籍の形式でいつでも購入できます。

techbookfest.org

頒布予定の技術書の内容

タンノ電子計算vol.2~手持ちの測定器で電波天文~ (人参食べたい)

こちらは電子工作関連の書籍です。著者本人が書いた説明文は以下

小惑星、地球外生命体や地球類似惑星の探査、あるいは宇宙の成り立ちや構成元素、銀河など宇宙についての研究対象は多岐にわたる。電波天文はこうした研究対象に対する手段として存在する。

本書では簡単な電波天文の概要を述べたあと、歴史から現代の技術までをざっくりさらう。その後、天文学に対する基礎知識として距離の単位や等級、赤方偏移などの説明をする。次に電波天文のキモと言える受信機について、量子数、バンド理論を通して素子の性質から定性的に説明している。その後、距離の計測方法から一般的な観測機の仕組み、フーリエ解析などを説明し、手持ちの道具で太陽電波の観測を試みる。とても真面目な一冊である。

全体的に専門的な知識は要らないように構成されているが、簡単な三角関数についての知識があるとフーリエ解析の部分などは読みやすいだろう。また、専門家にとってはバンド理論の説明などの量子力学的な部分については物足りないかもしれない。

また、最終章は筆者の趣味であるトランシーバ製作のために前回技術書の売上で購入したスペクトラムアナライザを利用している。測定機やアンテナは有り合わせのものであり、そこまで精度良く観測出来ておらず定量的な考察はできていない。GHzレベルの観測ができればもっと精度良く出来うるので機器が揃えば追記するが現状はおまけ程度に考えていただいても構わない。

本書では電子工学専攻のアマチュア無線家がそういった観点から電波天文を語っているので最新技術や歴史については他書にゆずるとする。



初学者がC#マンデルブロ集合を描画するまで (初学者リョウ)

次はCG系の書籍です。技術書典14で頒布予定だったのですが、間に合わなっかったので今回の技術書典で頒布します。
マンデルブロ集合とはフラクタル図形の一種です。フラクタル図形の厳密な定義は難しいですが、ざっくりと説明すると、ある描画ルールを再帰的に適用してできる図形の総称です。植物などを見ていると、まずは枝が2つに分かれて、さらにその2つの枝が2つに分かれて・・・ という再帰的な構造を有していることに気づきます。あれらも一種のフラクタル図形です。マンデルブロ集合の詳しい説明は、本人の書籍に譲りますが(自分もよく分からないので)、とても面白い見た目の図形ですのでネットで画像検索してみてはいかがでしょう。ちなみにYokaのツイッタープロ画フラクタル図形です(隙自語)。

フーリエ変換学習ノート (Yoka)

この本は筆者がフーリエ変換について学んで得た知識をまとめた本です。内容ははっきりとは決まっていませんが、最終的には簡単な画像処理をできる限りライブラリに頼らずに実装することを目標としています。

当選しますように

以上、サークルが頒布予定の書籍の概要です。頒布予定の書籍の一覧から分かる通り、メンバーそれぞれが興味を持ったことを書籍にしているので、分野は多岐にわたります。
メンバーの平均年齢は低く、怖い人もいないので、興味があれば気軽に話しかけてみてください(オフライン会場に当選した場合)。当選しますように・・・

技術書典14に申し込みます

技術書典14にサークルで申し込みをします

見出しの通りです。「マジックタンノリターンズ~タンノはもうメイドじゃない!!~」というサークル名で 技術書典14に申し込みます。以下がサークルの公式Twitterアカウントです(まだ作ったばかり)。
twitter.com

技術書典とは

技術書典とは、自分で書いた技術書を頒布・販売するイベントです。言うなれば、技術書版コミケです。オンラインとオフラインの両方の形式で開催され、オフライン会場であれば今年は2023/05/21 (日) に池袋サンシャインシティ 展示ホールD(文化会館ビル2F)で開催されます。楽しそうですね。オンラインの場合は、2023/05/20 (土) 〜2023/06/04(日)の期間内であればいつでもやっており、技術書典 :技術書のオンラインマーケット開催中から電子書籍の形式でいつでも購入できます。

techbookfest.org

頒布予定の技術書の内容

** モンテカルロ木探索を用いてリバーシAIを作ろう(Yoka)
別記事でも紹介したことがある通り、私はC#モンテカルロ木探索(MCTS)を用いたリバーシAIを開発していました(現在はC++で書き直し中)。その過程で得たMCTSに関する知識のうち、初歩的な部分をまとめた書籍です(予定)。また、ただMCTSについてまとめるだけではなく、実際にPythonを用いて読者にリバーシAIを実装してもらう形式で進めていきます。ゲームAIに関する初歩的な知識がなくとも、基本的な部分(minimax法や\alpha \beta法など)から解説していきますので、ぜひ気軽に手に取ってみてください(技術書典14に当選したらの話ですが)。

2023/04/24追記: Yokaの書籍に関しては内容に大きな変更があります。内容としましては、GUID Partition Tabelに関するものになります。

プログラムからパーティションを弄ろう(Yoka)

この書籍では、GPT(GUID Partition Table)の基本について説明し、実際にC言語を用いて生のディスクデータを書き換えて、パーティション構成を変えてみようという内容です。

50MHz DSBトランシーバー製作記(人参食べたい)

こちらは電子工作関連の書籍です。他メンバーが頒布予定の書籍ですので、詳しいことは説明できませんが、いくつか製作中の写真をもらったので掲載します。

※画像は開発中のものです

初学者がC#マンデルブロ集合を描画するまで(初学者リョウ)

次はCG系の書籍です。マンデルブロ集合とはフラクタル図形の一種です。フラクタル図形の厳密な定義は難しいですが、ざっくりと説明すると、ある描画ルールを再帰的に適用してできる図形の総称です。植物などを見ていると、まずは枝が2つに分かれて、さらにその2つの枝が2つに分かれて・・・ という再帰的な構造を有していることに気づきます。あれらも一種のフラクタル図形です。マンデルブロ集合の詳しい説明は、本人の書籍に譲りますが(自分もよく分からないので)、とても面白い見た目の図形ですのでネットで画像検索してみてはいかがでしょう。ちなみにYokaのツイッタープロ画フラクタル図形です(隙自語)。

困ったらドキュメントを読めの話(ゆざきりつ)

タイトルの通りです。皆さんも困ったときは、適当なサイトのコードをコピペして場当たり的に問題を解決しようとせず、使用言語や使用ライブラリの公式ドキュメントを当たってみましょう。おそらくそんな話です。

当選しますように

以上、サークルが頒布予定の書籍の概要です。サークルとは言っていますが、みんなが自分の興味のあることに取り組む自由なサークルです。怖い人もいないので、(もし出展できたら)気軽に立ち寄ってください。当選しますように・・・

C#で行列と列ベクトルの加算を実装した話

概要

C# + Math.NETでニューラルネットワークを実装しようと思ったけど, NumPyでいうブロードキャストに該当する機能が存在しない

ブロードキャストは行列と列ベクトルを加算するときなどに必要

自分で実装しよう

ブロードキャストとは

NumPyに実装されている便利な機能。線形代数においてMxNの行列AとM次元列ベクトルBについて, A+Bという演算は未定義。でもNumPyではエラーにならない。なぜか?NumPyではブロードキャストによってベクトルBをN個並べてMxNの行列に変換した後に行列Aに加算される*1。だから正しく計算できる。

ブロードキャストの使いどころ

例えばNNでは, WX + Bという計算を頻繁に行う。ここでWはNNの重み, XはNNへの入力をバッチの要素ごとに各列に並べたもの, そしてBはNNのバイアスベクトルを各列に並べたもの。しかし, 実際のところBの各列は同じ列ベクトルであるため, 実装上はわざわざバイアスベクトルを各列に並べるのはメモリの無駄になる。どうせならバイアスベクトルを\boldsymbol{b}として, WX + \boldsymbol{b}と書けた方が楽。このようなときにNumPyのブロードキャストが役に立つ。

Math.NETにブロードキャストが見当たらない

C#用(厳密には.NET用)の数学ライブラリとしてMath.NETがある。Math.NETは内部でIntel MKLやOpenBLASなどの数値計算ライブラリを用いているため, 行列積の計算などが非常に高速。しかし, ざっと調べたところ, Math.NETにはNumPyのブロードキャストに当たる機能が存在しなかった(自分が見逃しているだけかもしれないけど)。

自力で実装

以下, 行列とベクトルの各要素の型はfloat(32bit浮動小数点数)であることを前提に進める。
AddColumnVector(DenseMatrix left, DenseVector right, DenseMatrix sum)というメソッドを実装することにした。このメソッドは行列leftと列ベクトルrightを加算した結果を行列sumに格納する。色々と試行錯誤したのだがそれを書くには余白が足りなすぎるので割愛。結論から示すと以下のような実装になった(例外処理は省略)。

static void AddColumnVector(DenseMatrix left, DenseVector right, DenseMatrix sum)
{
    var leftArray = left.AsColumnMajorArray();
    var rightArray = right.AsArray();
    var sumArray = sum.AsColumnMajorArray();
    var (rowNum, colNum) = (sum.RowCount, sum.ColumnCount);

    unsafe
    {
        const int UNROLL_STEP_NUM = 8;
        fixed (float* leftPtr = leftArray)
        fixed (float* rightPtr = rightArray)
        fixed (float* sumPtr = sumArray)
        {
            for (var j = 0; j < colNum; j++)
            {
                var l = leftPtr + j * rowNum;
                var s = sumPtr + j * rowNum;
                for (var i = 0; i < rowNum; i += UNROLL_STEP_NUM)
                {
                    var leftVec = Avx.LoadVector256(l + i);
                    var rightVec = Avx.LoadVector256(rightPtr + i);
                    var sumVec = Avx.Add(leftVec, rightVec);
                    Avx.Store(s + i, sumVec);
                }
            }

            // 端数処理
            for (var j = 0; j < colNum; j++)
            {
                var l = leftPtr + j * rowNum;
                var s = sumPtr + j * rowNum;
                for (var i = UNROLL_STEP_NUM * (rowNum / UNROLL_STEP_NUM); i < rowNum; i++)
                    s[i] = l[i] + rightPtr[i];
            }
        }
    }
}

C#ではAVX組み込み関数を用いることができるので, それを用いて8要素同時に加算を行っている。列ベクトルの次元(rowNum)を8で割った余りは単純に2重ループを回して計算している(//端数処理 以降)。

実行速度

測定方法

1024x1024の行列に1024次元の列ベクトルを加算する処理を1000回行う。行列とベクトルの各要素は0.0~1.0の一様乱数で初期化する。
また、比較対象としてNumPyでも同様の処理を行う。
ソースコードは"付録"の項を参照。

マシンスペック
CPU Intel(R) Core(TM) i7-10510U
メモリ DDR3 16GB
NumPyの情報

バージョン: 1.23.3
NumPyの設定によって結果が変わると思うので, numpy.__config__.show()の結果を載せておく。

>>> np.__config__.show()
openblas64__info:
library_dirs = ['D:\\a\\numpy\\numpy\\build\\openblas64__info']
libraries = ['openblas64__info']
language = f77
define_macros = [('HAVE_CBLAS', None), ('BLAS_SYMBOL_SUFFIX', '64_'), ('HAVE_BLAS_ILP64', None)]
blas_ilp64_opt_info:
library_dirs = ['D:\\a\\numpy\\numpy\\build\\openblas64__info']
libraries = ['openblas64__info']
language = f77
define_macros = [('HAVE_CBLAS', None), ('BLAS_SYMBOL_SUFFIX', '64_'), ('HAVE_BLAS_ILP64', None)]
openblas64__lapack_info:
library_dirs = ['D:\\a\\numpy\\numpy\\build\\openblas64__lapack_info']
libraries = ['openblas64__lapack_info']
language = f77
lapack_ilp64_opt_info:
library_dirs = ['D:\\a\\numpy\\numpy\\build\\openblas64__lapack_info']
libraries = ['openblas64__lapack_info']
language = f77
define_macros = [('HAVE_CBLAS', None), ('BLAS_SYMBOL_SUFFIX', '64_'), ('HAVE_BLAS_ILP64', None), ('HAVE_LAPACKE', None)]
Supported SIMD extensions in this NumPy install:
baseline = SSE,SSE2,SSE3
found = SSSE3,SSE41,POPCNT,SSE42,AVX,F16C,FMA3,AVX2
not found = AVX512F,AVX512CD,AVX512_SKX,AVX512_CLX,AVX512_CNL,AVX512_ICL

※(2022/09/13 追記) NumPyをAVX2まで対応させたものに変更しました。以下が新たなnumpy.__config__.show()の結果です。

>>> np.__config__.show()
blas_armpl_info:
NOT AVAILABLE
blas_mkl_info:
NOT AVAILABLE
blis_info:
NOT AVAILABLE
openblas_info:
NOT AVAILABLE
accelerate_info:
NOT AVAILABLE
atlas_3_10_blas_threads_info:
NOT AVAILABLE
atlas_3_10_blas_info:
NOT AVAILABLE
atlas_blas_threads_info:
NOT AVAILABLE
atlas_blas_info:
NOT AVAILABLE
blas_info:
NOT AVAILABLE
blas_src_info:
NOT AVAILABLE
blas_opt_info:
NOT AVAILABLE
lapack_armpl_info:
NOT AVAILABLE
lapack_mkl_info:
NOT AVAILABLE
openblas_lapack_info:
NOT AVAILABLE
openblas_clapack_info:
NOT AVAILABLE
flame_info:
NOT AVAILABLE
atlas_3_10_threads_info:
NOT AVAILABLE
atlas_3_10_info:
NOT AVAILABLE
atlas_threads_info:
NOT AVAILABLE
atlas_info:
NOT AVAILABLE
lapack_info:
NOT AVAILABLE
lapack_src_info:
NOT AVAILABLE
lapack_opt_info:
NOT AVAILABLE
numpy_linalg_lapack_lite:
language = c
define_macros = [('HAVE_BLAS_ILP64', None), ('BLAS_SYMBOL_SUFFIX', '64_')]
Supported SIMD extensions in this NumPy install:
baseline = SSE,SSE2,SSE3,SSSE3,SSE41,POPCNT,SSE42,AVX,F16C,FMA3,AVX2
found =
not found = XOP,FMA4,AVX512F,AVX512CD,AVX512_SKX,AVX512_CLX,AVX512_CNL,AVX512_ICL

測定結果

平均実行速度[ms]
マイコー 0.473
NumPy(SSE3まで対応)のadd関数 1.233

※(2022/09/13 追記) この実験で用いたNumPyはSSE3までにしか対応していないため, NumPyにとってかなり不利な比較になっていました。後日, AVXに対応したものを用意して再実験した結果を掲載する予定です。

※(2022/09/13 追記) 新たな実験結果を追加しました。NumPyをAVXに対応させても, 特に大きな速度差はありませんでした。

平均実行速度[ms]
マイコー 0.473
NumPyのadd関数 1.248

NumPyよりも2倍以上高速に動作した。

まとめ

C#を用いて行列と列ベクトルの加算を実装した。実は行列と行ベクトルの加算も実装したが, 内容はほとんど同じなので割愛。自分は高速化に関して素人なので, もっとよい実装方法があったらコメントで指摘してくださるとうれしいです。

付録

C#ベンチマークコード
using System;
using System.Diagnostics;
using System.Runtime.Intrinsics.X86;

using MathNet.Numerics.LinearAlgebra.Single;

class Program
{
    const int DIM = 1024;
    const int SAMPLE_NUM = 1000;

    public static void Main()
    {
        Benchmark(AddColumnVector);
    }

    static void Benchmark(Action<DenseMatrix, DenseVector, DenseMatrix> f)
    {
        var result = new DenseMatrix(DIM, DIM);

        var sw = new Stopwatch();
        for (var i = 0; i < SAMPLE_NUM; i++)
        {
            var left = InitMatrixWithRand(DIM, DIM);
            var right = DenseVector.OfArray(InitMatrixWithRand(DIM, 1).AsColumnMajorArray());
            sw.Start();
            f(left, right, result);
            sw.Stop();
        }
        Console.WriteLine($"Ellapsed = {(float)sw.ElapsedMilliseconds / SAMPLE_NUM}[ms]");
    }

    static DenseMatrix InitMatrixWithRand(int rowNum, int colNum)
    {
        var ret = new DenseMatrix(rowNum, colNum);
        var array = ret.AsColumnMajorArray();
        for (var i = 0; i < array.Length; i++)
            array[i] = Random.Shared.NextSingle();
        return ret;
    }

    static void AddColumnVector(DenseMatrix left, DenseVector right, DenseMatrix sum)
    {
        var leftArray = left.AsColumnMajorArray();
        var rightArray = right.AsArray();
        var sumArray = sum.AsColumnMajorArray();
        var (rowNum, colNum) = (sum.RowCount, sum.ColumnCount);

        unsafe
        {
            const int UNROLL_STEP_NUM = 8;
            fixed (float* leftPtr = leftArray)
            fixed (float* rightPtr = rightArray)
            fixed (float* sumPtr = sumArray)
            {
                for (var j = 0; j < colNum; j++)
                {
                    var l = leftPtr + j * rowNum;
                    var s = sumPtr + j * rowNum;
                    for (var i = 0; i < rowNum; i += UNROLL_STEP_NUM)
                    {
                        var leftVec = Avx.LoadVector256(l + i);
                        var rightVec = Avx.LoadVector256(rightPtr + i);
                        var sumVec = Avx.Add(leftVec, rightVec);
                        Avx.Store(s + i, sumVec);
                    }
                }

                // 端数処理
                for (var j = 0; j < colNum; j++)
                {
                    var l = leftPtr + j * rowNum;
                    var s = sumPtr + j * rowNum;
                    for (var i = UNROLL_STEP_NUM * (rowNum / UNROLL_STEP_NUM); i < rowNum; i++)
                        s[i] = l[i] + rightPtr[i];
                }
            }
        }
    }
}
Pythonベンチマークコード
import time
import numpy as np

DIM = 1024
SAMPLE_NUM = 1000

ellapsed = 0
s = np.zeros((DIM, DIM))
for _ in range(SAMPLE_NUM):
    left = np.random.rand(DIM, DIM)
    right = np.random.rand(DIM, 1)
    start = time.perf_counter()
    np.add(left, right, s)
    end = time.perf_counter()
    ellapsed += end - start
    
print(f"Ellapsed = {(ellapsed / SAMPLE_NUM) * 1000}[ms]")
    

*1:もちろん, わざわざMxN行列用のメモリを確保してベクトルBをそこに書き込んでいるわけではなく, 実装上は普通に行列Aの各列にベクトルBを加算している

自作リバーシAI Kalmiaの仕組み

はじめに

Kalmiaはモンテカルロ木探索(MCTS)ベースのリバーシの思考エンジンです。思考エンジンと呼ぶと馴染みが無いかもしれませんが、誤解を恐れずに端的に述べればいわゆるAIです。KalmiaのMCTSではロールアウトの代わりにEdax*1で用いられているような静的評価関数を用いていることが特徴です。強さはトッププログラムの足元にも及びませんが、おそらく人間のトッププレイヤーやそれ以上の強さであると推測しています。
この記事ではまずKalmiaで用いている技術の基本的な説明から始め、その後にKalmiaの仕組みについて説明します。ですので既にMCTSリバーシの評価関数に関する知識を持っている方は後半からいきなり読んでもらって構いません。 分かりにくい部分や不適切な部分があればコメント欄もしくはTwitter(以下リンク)で指摘をお願いします。

twitter.com

本記事の構成

前半: Kalmiaで用いている技術の基本的な説明
後半: Kalmiaの仕組み

モンテカルロ木探索

モンテカルロ木探索とは"選択"、"シミュレーション"、"バックアップ"、"展開"の4つのフェーズを繰り返しながら探索木を成長させていくアルゴリズムです。元々は囲碁プログラムで使われ始めたアルゴリズムであり、2006年にRémi Coulom氏がCrasy Stoneという囲碁プログラムに採用してブレイクスルーとなりました。Crasy Stoneで採用されたMCTSを簡単に説明すると、乱数による対局シミュレーション(モンテカルロ法)の勝敗で局面を評価する方法と木探索を組み合わせたものです。ただし、MiniMax探索のようにあらゆる局面をしらみつぶしに探索するのではなく、有望そうな手から優先的にシミュレーションを行って木を成長(先を読む)させます。MCTSで肝となるのは"有望そうな手"を選択する方法です。Crasy Stoneでは独自の手法で有望そうな手を選択しています*2。しかし、後に発表された論文“Bandit based Monte-Carlo Planning”*3で紹介された、UCB1方策(後述)によって手を選択する手法が後の囲碁プログラムでは一般的になります。この手法はUCB1方策を木探索に適用することからUCB applied to trees(UCT)と呼ばれます。現在では単にMCTSといった場合、このUCT(もしくはその亜種)のことを指すことが多いです。実際にKalmiaで用いているMCTSもUCTです。UCTについては後に詳しく述べます。

UCB1方策

UCB1(Upper Confidence Bound 1)方策とは、多腕バンディット問題を解くための方策のうちの1つです。多腕バンディット問題とは以下のような問題です。

"当たりの出る確率がそれぞれ異なるスロットマシーンがX台ある。任意のスロットマシーンに合計N回挑戦できるとき、報酬をできる限り最大化するにはどのような方策でスロットマシーンを回せばよいだろうか。ただし、スロットマシーンの当たりが出る確率は公開されておらず、同じマシーンに複数回挑戦しても良いものとする。"
※ここでいう報酬は"当たりが出た回数"という意味で用いています。以降は"当たりが出た回数"と書くと文が長くなるので"報酬"と記述します。

この問題に対する最も簡単な方策は、それぞれのスロットマシーンに一定の回数だけ挑戦して報酬の期待値を推測する方法です。そして残りの回数は報酬の期待値が高そうなスロットマシーンだけを回し続けます。具体例で説明すると、3台のスロットマシーンに計50回挑戦できるとき、それぞれのスロットマシーンを10回ずつ回して当たりが出た回数をカウントし、もっとも当たりが出たスロットマシーンを残りの20回だけ回すという方策です。しかしこの方策では、一定回数は全てのマシーンを回す必要があり非常に非効率です。そこで今まで回したスロットマシーンの結果を活用しつつ、まだ十分に回せていないマシーンも回すという方策が提案されます。それがUCB1方策です。UCB1方策では以下のUCB1式で求めたスコアが最も大きいスロットマシーンを選択して回していきます。

 u_i = \frac{w_i}{n_i} + \sqrt{\frac{2\log{N}}{n_i}}\tag{1}
 N = \sum_{i}^{}{n_i}

u_i: スロットマシーンiのUCBスコア
w_i: スロットマシーンiの累積報酬(今まで当たりが出た回数)
n_i: スロットマシーンiを回した回数
N: 各スロットマシーンを回した回数の総和

この式の第1項は今までの報酬の標本平均となっていることが分かります。 第2項を無視すれば、このスコアは報酬の標本平均が大きければ大きいほどそのスロットマシーンが選ばれやすくなるということを意味しています。今までの報酬の標本平均が高いのであれば、報酬の期待値も高いだろうと推測できるので、これは直感にも合うと思います。では第2項は一体何でしょう。対数や平方根を無視すれば、スロットマシーンiを回した回数の逆数になっています。これの意味することとはつまり、まだ十分に回していないスロットマシーンも選ばれるということです。第1項と第2項をまとめれば、UCB1方策とは「報酬の期待値が高そうなスロットマシーンを選びつつ、まだ十分に回せていないスロットマシーンも回す」という方策です。

UCB1式についてもう少し詳しい説明をすると、UCB1式と報酬の期待値E[X_i]との間には以下の関係が成り立ちます。

P(E[X_i] \leq \frac{w_i}{n_i} + \sqrt{\frac{2\log{N}}{n_i}}) = P(E[X_i] - \frac{w_i}{n_i} \leq \sqrt{\frac{2\log{N}}{n_i}}) \geq 1 - \frac{1}{N^4}\tag{2}

X_i: スロットマシーンiが当たりのときに1、そうではないときに0をとる確率変数

上式からUCB1式の第2項は、報酬の期待値と報酬の標本平均との差についての、信頼度(1-\frac{1}{N^4})以上における信頼区間の上限(Upper Confidence Bound)を表していることが分かります。つまりUCB1方策では、標本平均から推測される報酬の期待値と真の報酬の期待値との乖離が大きいと思われるスロットマシーンも選んでいるのです。また、信頼度にNが含まれていることから分かる通り、全体の試行回数が大きくなるにつれ、信頼度が1に近づくので、UCB1式の第2項の効果は薄れます。
より詳しい説明については、他の方が書いた以下の記事が詳しいです。

blog.monochromegane.com

UCT

UCTとはUCB applied to treesの略称であり、先ほどのUCB1方策を木探索に適用した手法です。MCTSでは有望な手を選択しながら探索木を成長させていくと述べましたが、この有望な手を選択するフェーズでUCB1方策を用います。UCTの具体的な流れを説明します。

選択フェーズ

各子ノードの訪問回数(ノードの選択回数)と対局シミュレーションの勝数、および兄弟ノードの訪問回数の総和からUCB1式に基づいてスコアを算出します。そしてそのスコアが最も高い手を選択します。この操作を末端ノード(葉ノード)に到達するまで続けます。なお、UCTで用いるUCB1式では、以下のように事前に決めた定数を第2項に作用させることが多いです。このようにすることで第2項の強さを変更することができます。c=\sqrt{2}とすれば、式(1)と同じです。

 u(s,a) = \frac{w(s,a)}{n(s,a)} + c\sqrt{\frac{\log{\sum_{b}{n(s,b)}}}{n(s,a)}}\tag{3}

s: 現在の局面(状態)
a: 着手(行動)
u(s,a): 局面sにおける, 着手aのUCBスコア
w(s,a): 局面sで着手aを行って得られた累積報酬(末端ノードで行われた対局シミュレーションの勝敗のこと. 勝利: +1, 敗北: 0, 引き分け: +0.5)
n(s,a): 局面sで着手aを行うノードに訪問した回数
\sum_{b}{n(s,b)}: 兄弟ノードの訪問回数の総和
 c: 第2項の強さを決める定数. 大きくすればするほど、探索が不十分なノードをより優先するようになるが,その代わり木の深さは浅くなる.

なお、KalmiaではAlphaZero*4で用いられている方法で定数部分が探索の進み具合に応じて変化するようにしています。詳しくは後に述べます。

(図1) 選択フェーズでUCBが最大の辺(緑)を選択していく様子
シミュレーションフェーズ

末端ノードに到達したら対局シミュレーションを1回行い、その勝敗を末端ノードに記録します(図2)。
Kalmiaのシミュレーションフェーズでは対局シミュレーションを行っていないのですが、それについては後に述べます。

(図2) シミュレーションフェーズで白の勝ちだった場合
バックアップフェーズ

末端ノードのシミュレーション結果をそのノードに至るまでに経由したノードに伝えます。例えば、末端ノードでのシミュレーション結果が白の勝ちであった場合(図3)、白の着手に関するノードでは累積報酬と訪問回数に1を加算し、黒の着手に関するノードでは訪問回数にのみ加算を行います。白の勝ちは黒にとっての負けだからです。また、末端ノードでのシミュレーション結果が引き分けであった場合は、どちらの着手のノードであるかに関わらず、訪問回数には1を加算、累積報酬には0.5を加算します。当たり前のことですが、相手が引き分けであれば自分も引き分けだからです。

(図3) 対局シミュレーションの結果が白の勝ちであった場合のバックアップ
展開フェーズ

末端ノードのシミュレーション回数があらかじめ決められた閾値に達した場合は、そこからさらにノードを展開します。もちろん、末端ノードの状態が終局であれば展開を行いません。

(図4) N(s2, a21)が展開閾値に達した場合の展開フェーズ

パターンの重み付き線形和を用いた評価関数

ゲームにおいて現在の局面の有利不利を判断することはとても重要なことです。リバーシほどの広い探索空間を持つゲームであれば、初期局面から終局まで完全に読み切ることは不可能だからです。そこで、ある程度の深さで探索を打ち切り、数手先の局面の有利不利を判断する必要があります。その役割を担うのが評価関数です。
リバーシAIの分野では1997年に当時のトッププレイヤーに6戦全勝したLogistello*5というプログラムが存在し、そのプログラムで用いられた石の並びのパターンの重み付き線形和による評価方法が現在の主流です。リバーシプログラムの中ではトップレベルの強さであるEdaxやZebra*6などのプログラムでもLogistelloの評価関数の亜種が用いられています。Logistelloの評価関数ではまず局面をいくつかの要素に分割し、それぞれの要素の石の並びのパターン(特徴)*7に点数を付け、最後にその点数を足し合わせます(式(4))。

v(s) = \boldsymbol{\phi}(s)\cdot\boldsymbol{w} = \sum_{i}^{}{\phi_i(s) w_i}\tag{4}

 s: 局面
 v(s): 局面sの価値(評価値)
 \boldsymbol{\phi}(s): 局面sの特徴ベクトル(各要素は0か1をとる. 特徴iが出現した時はi番目の要素が1,そうでなければ0となる. )
 \boldsymbol{w}: 重みベクトル(各要素はそれぞれの特徴の点数)

より詳しい実装については以下の記事が参考になります。

qiita.com

評価関数の最適化

以下、石の並びのパターンのことを"特徴"、特徴の点数のことを"重み"と呼びます。
リバーシの評価関数では、局面の各特徴に重みをつけて足し合わせると説明しました。では、それぞれの特徴の重みはどのように決めればよいのでしょうか。もしあなたがとても強いリバーシのプレイヤーあったとしても、全ての特徴の重みを手作業で決めることは不可能に近いです。なぜなら、リバーシで出現しうる特徴は20万以上に及ぶからです。さらにリバーシでは、ゲームの進行度によって評価関数を使い分けるので、実際は数百万個にも及ぶ重みを決める必要があります。そこで、それぞれの特徴の重みは、実際に打たれた棋譜から勾配降下法で最適化していきます。勾配降下法では評価関数によって予測されたゲームの結果(石差や勝敗など)と、実際の棋譜における結果との差が小さくなるように点数を調整していきます。ここでは評価関数を用いて最終的な石差を予測したい場合について説明します。まずは以下のような損失関数を考えます。

 E(\boldsymbol{w};\boldsymbol{D})=\frac{1}{2}\sum_{s \in D}^{}{(v^*(s) - v(s))^2}\tag{5}

 \boldsymbol{D}: 棋譜に含まれる局面の集合
 v^*(s): 局面sの最終石差
 v(s): 評価関数によって求めた局面sの最終石差の予測値
 \boldsymbol{w}: 評価関数v(s)の重みベクトル

単純に、実際の最終石差と評価関数による予測値との差についての二乗和を取っているだけです。\frac{1}{2}を掛けている理由は、後に微分する時に定数が消えるからです。
最適な評価関数を得るには、損失関数 E(\boldsymbol{w};\boldsymbol{D})の値が最小となるような \boldsymbol{w}を求めればいいわけですから、まずは \boldsymbol{w}の各要素で損失関数を偏微分します。 \boldsymbol{w}の任意の要素w_jで損失関数 E(\boldsymbol{w};\boldsymbol{D})偏微分すると、以下のようになります。

 \frac{\partial E(\boldsymbol{w};\boldsymbol{D})}{\partial w_j}=\sum_{s \in D}^{}{(v^*(s) - v(s))(-\frac{\partial v(s)}{\partial w_j})}

ここでv(s) = \boldsymbol{\phi}(s)\cdot\boldsymbol{w} = \sum_{i}^{}{\phi_i(s) w_i}より、

\frac{\partial v(s)}{\partial w_j}=\phi_j(s)

\phi_j(s)は特徴jが出現していれば1、そうでなければ0であるから

 \frac{\partial E(\boldsymbol{w};\boldsymbol{D})}{\partial w_j}=
\begin{cases}
  \sum_{s \in D}^{}{(v(s) - v^*(s))} &  (特徴jが局面sに含まれる)\\
  0 &  (それ以外)
  \end{cases}\tag{6}

このように偏微分した結果はとてもシンプルになりました。重みベクトルの各要素で損失関数を偏微分した値は、その点における損失関数の勾配(傾き)を表しています。故に重みベクトルの各要素の値を、勾配を降下する方向に動かすことで、損失関数の値が小さくなることが期待できます。勾配降下法では、重みで偏微分 -> あらかじめ決めた学習率\etaの分だけ重みを変化 -> もう一度重みで偏微分 -> ・・・ という操作を繰り返しながら損失関数の値を小さくしていきます(式(7))

 {w'}_j = w_j - \eta \frac{\partial E(\boldsymbol{w};\boldsymbol{D})}{\partial w_j}\tag{7}

Kalmiaの仕組みへ

以上、Kalmiaで用いている技術について概ね解説しました。ここからは実際のKalmiaの仕組みについて解説していきます。

KalmiaのMCTS

KalmiaのMCTSは選択フェーズにUCBを用いたUCTを用いています。ただし、シミュレーションフェーズと選択フェーズに若干の工夫があります。

シミュレーションフェーズ

通常のMCTSでは、このフェーズで対局シミュレーションを1回だけ行い、その勝敗を記録します。しかし、対局シミュレーションを行う必要があったのは、そもそも囲碁の評価関数を作るのが困難であったことが理由です。それならば、既にLogistelloの評価関数のような計算が高速かつ、それなりに精度の高いものが存在するリバーシでは不要です。そこでKalmiaでは、対局シミュレーションの代わりに、勝率を予測する価値関数(評価関数)v(s)を用いています。この価値関数v(s)は、ある局面sに対して決まった値を出力するので、通常のMCTSのように展開閾値を決める必要はありません。末端ノードの評価が1回終わったらすぐにノードを展開します。すなわち展開閾値は常に1です。

選択フェーズ

選択フェーズではUCBスコアがもっとも大きい子ノードを選択しますが、以下のような例外を設けています。

(1) 勝ちに至る子ノードが存在する場合
UCBスコアに関わらずそのノードを選ぶ。

(2) 負けに至る子ノードが存在する場合
UCBスコアに関わらずそのノードは選ばない。

(3) 引き分けに至る子ノードと負けに至る子ノードのみが存在する場合
引き分けを選ぶ(負けるよりはマシなので)。

また、選択フェーズにおいて、勝ちに至る子ノードを1つでも見つけた場合は、その親ノードも勝ちになるので、親ノードに"勝ち"とマークします。逆に子ノードが全て負けの時は親ノードも負けになるので、親ノードに"負け"とマークします。引き分けと負けの子ノードしか存在しない場合は、親ノードに"引き分け"とマークします。このように、確定したゲームの勝敗を親ノードへ伝播することで、終盤では木の末端でAND/OR木のような挙動をとります。この方法はチェスプログラムのLeela Chess Zero*8MCTSを参考にしました。

KalmiaのUCB式

KalmiaではUCTの説明の際に紹介した式(4)を少し改変しています。その理由は、式(4)をそのまま用いた場合、特定の手を極端に深く読んでしまう現象が発生したからです。具体的には、他の候補手は12手程度先しか読んでいないのにもかかわらず、ある1手は50手ぐらい先を読んでしまうような現象が発生しました。定数cの値を何回か変更してもこの傾向は改善しなかったため、AlphaZeroで用いられていたPUCT式を参考に以下のような改変を加えました。

 u(s,a) = \frac{w(s,a)}{n(s,a)} + c(s)\sqrt{\frac{\log{\sum_{b}{n(s,b)}}}{n(s,a)}}\tag{8}
c(s) = \log{\frac{1 + N(s) + C_{base}}{C_{base}} + C_{init}}
N(s) = \sum_{b}{n(s,b)}

 C_{base}, C_{init}: 事前に決めるハイパーパラメータ

このようにすることで、親ノードの訪問回数によってc(s)が動的に変化します。具体的には、親ノードの訪問回数が少ないうちは価値の高い子ノードが比較的優先され、親ノードの訪問回数が多くなるにつれ、訪問回数が少ない(探索が十分でない)ノードも優先されます。KalmiaではC_{base} = 19652, C_{init} =0.35 と設定していますが、何回か動かして適当に決めた値なので最適な値とは言えません。いつかベイズ最適化などでより良い値にチューニングする予定です。

式(8)は\sum_{b}{n(s,b)} = 0(どの子ノードも未訪問)のときと、n(s,a)=0のときに未定義です。Kalmiaでは、ルート直下の未訪問ノードに関しては勝ちであると仮定し、それ以外の未訪問ノードは親ノードの価値で初期化します。未訪問ノードを親ノードの価値で初期化することで、親ノードよりも価値が高そうな子ノードが見つかれば、しばらくの間はそのノードが選ばれ続けるようになります。

MCTSの並列化

Kalmiaでは、MCTSの "選択"、"シミュレーション"、"バックアップ"、"展開" の一連の操作をCPUの論理コア数の分だけ並列に行います。ただし、単純に並列化してしまうと複数のスレッドが同じノードを選んでしまって効率が悪くなるので、あるスレッドが探索中のノードにはvirtual lossというペナルティをノードの訪問回数に加えます。こうすることで一時的にノードの価値が下がり、他のスレッドからは選ばれづらくなります。Kalmiaではvirtual lossを3に設定しています。

終盤探索

リバーシでは、空きマス数がある程度少なくなれば終局まで読み切ることが可能です。Kalmiaでは、終盤21手は必勝読み(最終的な勝敗を読む)、終盤20手は石差が最大となるような手を読みます*9。終盤探索では特に珍しいことは行っておらず、Alpha-Beta探索で読み切っています。ただし、局面の合流に関しては置換表で管理しています。Move orderingは単純な速さ優先なので非常に重い終盤ソルバーとなっています。終盤探索に関してはつい最近実装した機能なので今後改良していく予定です。

価値関数(評価関数)

KalmiaのMCTSの説明において対局シミュレーションの代わりに価値関数v(s)の値を利用すると述べました。ここではv(s)の構造について説明します。
価値関数v(s)は、局面からEdaxと同様の特徴を抽出し、それらの重み付き線形和で局面の価値を算出します。ただし、Kalmiaは最終石差を予測するEdaxとは異なり勝率を予測するので、重み付き線形和を標準シグモイド関数(式(9))に通して0.0~1.0の値に変換しています。

 f(x) = \frac{1}{1 + e^{-x}}\tag{9}

v(s) = f(\boldsymbol{\phi}(s)\cdot\boldsymbol{w}) = f(\sum_{i}^{}{\phi_i(s) w_i})\tag{10}

 s: 局面
 v(s): 局面sの価値(予測勝率)
 \boldsymbol{\phi}(s): 局面sの特徴ベクトル(各要素は0か1をとる. 特徴iが出現した時はi番目の要素が1,そうでなければ0となる. )
 \boldsymbol{w}: 重みベクトル(各要素はそれぞれの特徴の点数)

1つの価値関数のみでは精度があまり良くないので、棋譜を4手ごとに分割した後、それぞれについて価値関数を最適化します。つまり価値関数は15個存在します。

価値関数の最適化

価値関数の最適化にはリバーシ棋譜用いており、勝ちを1、負けを0、引き分けを0.5として価値関数の予測値と棋譜の勝敗との誤差が小さくなるように最適化します。「評価関数の最適化」の項で紹介した損失関数E(\boldsymbol{w};\boldsymbol{D})では二乗和を用いましたが、Kalmiaの価値関数では勝率を予測するので、損失関数には2値クロスエントロピーを用います。クロスエントロピーとは2つの確率分布 p, qがどれくらい似ているかを表す指標です(式(11))。

E(\boldsymbol{w};\boldsymbol{D}) = -v^{*}(s)log(v(s)) - (1 - v^{*}(s))log(1- v(s)) \tag{11}

 \boldsymbol{D}: 棋譜に含まれる局面の集合
v^{*}: 棋譜の勝敗(勝ち: 1, 負け: 0, 引き分け: 0.5)
v(s): 価値関数による予測勝率
 \boldsymbol{w}: 評価関数v(s)の重みベクトル

この損失関数E(\boldsymbol{w};\boldsymbol{D})の値をできる限り小さくするような \boldsymbol{w}を求めればよいので、まずは \boldsymbol{w}の任意の要素w_jで損失関数 E(\boldsymbol{w};\boldsymbol{D})偏微分します(式(12))。

 \frac{\partial E(\boldsymbol{w};\boldsymbol{D})}{\partial w_j} = \sum_{s \in D}^{}{(- \frac{v^{*}(s)}{v(s)} \cdot \frac{\partial{v(s)}}{\partial{w_j}} + \frac{1 - v^{*}(s)}{1 - v(s)} \cdot \frac{\partial{v(s)}}{\partial{w_j}})} \tag{12}

標準シグモイド関数導関数(式(13))より、\frac{\partial{v(s)}}{\partial{w_j}}は以下の通り(式(14))

f'(x) = (1 - f(x))f(x) \tag{13}

\frac{\partial v(s)}{\partial w_j} = (1- v(s))v(s) \frac{\partial  (\sum_{i}^{}{\phi_i(s) w_i})}{\partial w_j} \tag{14}

よって損失関数 E(\boldsymbol{w};\boldsymbol{D}) w_jにおける勾配は式(14)を式(12)に代入することで得られます(式(15))。

\frac{\partial E(\boldsymbol{w};\boldsymbol{D})}{\partial w_j}  = (v(s) - v^{*}(s))\frac{\partial (\sum_{i}^{}{\phi_i(s) w_i})}{\partial w_j} \tag{15}

ここで\phi_i(s)は特徴iが出現していれば1、そうでなければ0であるから、最終的に式(15)は以下の形になります(式(16))。

 \frac{\partial E(\boldsymbol{w};\boldsymbol{D})}{\partial w_j}=
\begin{cases}
  \sum_{s \in D}^{}{(v(s) - v^*(s))} &  (特徴jが局面sに含まれる)\\
  0 &  (それ以外)
  \end{cases}\tag{16}

結局のところ、「評価関数の最適化」で紹介した式(6)と全く同じ形になります。あとはこの勾配の式を利用して勾配降下法で最適化するだけです。

Kalmiaの欠点

ここではKalmiaの抱えている欠点について紹介します。

石差にこだわらない

Kalmiaでは勝率を予測しながら手を決めるので石差にこだわらない打ち方をします。それ故に、勝勢になるといきなり手が緩みます。具体的には、大差で自分が勝てるような局面でも2石差で勝つような手を選択することがあるのです。逆に劣勢になった際にはとにかく勝つために無理な手を選択し始めます。本当は僅差で負ける局面なのに、大差で負けてしまうこともしばしばあります。しかし、終盤完全読みを実装したお陰でこの傾向は従来よりはマシになりました。

メモリ使用量が多い

MCTSアルゴリズムの特性上、一度訪問したノードは探索中ずっとメモリに残ります。それ故に長時間思考を行うとあっというまにメモリを食い潰します。

読み抜けを起こしやすい

MCTSでは探索が進むにつれ有望でないノードが選ばれづらくなります。これは前向き枝刈りのような効果を見せる反面、これが原因で有利な手を見逃してしまうことがあります。特にMCTSでは、ほとんどの手順では不利な手であるが、ある特定の手順では有利であるような手に弱いです。

Kalmiaの今後の改善点

速さを犠牲に評価関数の精度を上げる

Kalmiaの欠点で述べたように、現状ではメモリをすぐに食い潰してしまうので、敢えて価値関数の精度を上げて重くすることによりメモリの使用率を抑えようと考えています。新たな価値関数としては以下のようなものを考えています。

・末端ノードで浅いAlpha-Beta探索を行う
現状のKalmiaでは末端ノードに到達した際に、価値関数の出力をそのままノードに記録していますが、これを浅いAlpha-Beta探索によって求めた価値で記録するようにします。このようにすることで純粋な価値関数の出力よりも高い精度が期待できます。

非線形モデルなどのより自由度の高いモデルを用いる
現在の価値関数では単純な線形モデルを用いていますが、これをニューラルネットワークのような非線形モデルに変更します。モデルの自由度が上がれば必然的に予測精度も向上するため大きく強さに影響することが期待できます。その反面、予測精度の向上によるメリットが低速化のデメリットよりも下回れば弱体化することになります。また、一般的に自由度の高いモデルは過適合を起こしやすく最適化が難しいです。

まとめと今後の予定

ここまでKalmiaで用いている技術の基本的な説明とKalmiaの仕組みについて説明しました。今後は「Kalmiaの今後の改善点」の項で挙げた改善を施す予定ですが、その前にKalmia用のGUIアプリ、もしくはWebアプリを開発しようかと考えています。というのも、現在のKalmiaはGoGuiというGUIアプリとGTP(Go Text Protocol)というプロトコルで通信することで人間との対局を実現しており、それらのソフトウェアの導入が少々面倒なのです。これでは多くの方に遊んでもらうことができないため、より遊びやすい形で配布したいと思っています。
Kalmiaのソースコードは以下のリポジトリで公開中です(開発言語はC#)。何か改善点やバグなどがあれば遠慮なく指摘してもらって構いません。ライセンスはGPL3です。

github.com

*1:https://github.com/abulmo/edax-reversi

*2:Rémi Coulom (2007). “Efficient Selectivity and Backup Operators in Monte-Carlo Tree Search”. Computers and Games, 5th International Conference, CG 2006, Turin, Italy, May 29–31, 2006.

*3:http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.102.1296

*4:https://www.science.org/doi/suppl/10.1126/science.aar6404/suppl_file/aar6404-silver-sm.pdf

*5:https://skatgame.net/mburo/log.html

*6:http://www.radagast.se/othello/zebra.html

*7:実際のLogistelloには石の並び以外にも着手可能数などの評価項目もあります。

*8:https://github.com/LeelaChessZero/lc0

*9:設定で変更は可能です