セルフィッシュマイニングのインセンティブ分析（マルコフ決定過程による分析）

追記：2018年12月に、より詳細な分析記事を書いています：
https://qiita.com/stnk20/items/297dafff76fd330d8f90

つい先日、Monacoinでセルフィッシュマイニング（Block Withholding Attackといったほうが正確なのかもしれません）が行われていたことがわかりました。事件の詳細は大石さんのこちらの記事が詳しいです。

早速セルフィッシュマイニングに関する分析記事も出ていましたので興味深く読みましたが、分析内容に関しては下記の点が気になりました：

• 記事中の式 $n(p_{n,0}(R+x)−W) $ は、$n(p_{n,0}R−W)+p_{n,0}x$ とすべき
• 撤退オプションが考慮されていない

前者は単なる計算ミスだと思いますが、後者の撤退オプションの取扱いについてはもう少しモデリングに工夫が要ります。

後者について少し説明します。
攻撃中にメインチェーンが先行して伸びてくると攻撃の成功確率は下がりますが、攻撃者がキャンセルしようとしている送金の送金先は、取引所の自分の口座のはずです。ですので、攻撃が失敗しそうなら一旦攻撃をあきらめて、資金を戻して再度チャレンジすることで、より効率的な攻撃を継続することができます。
これはリアルオプションで言うところの撤退オプションです。

2018/05/31追記
Nicehashという、ハッシュレートを一時的に購入できるサービスがあるのですが、Monacoinについてはメインのハッシュレートを圧倒するような量が購入できる状態のようです（参考）。この場合一気にハッシュレートを高めることで攻撃したほうが効果的ですので、実際は本記事のような撤退オプションは考えずに攻撃がなされたと思います。
あくまで読み物としてお楽しみください。

分析は場合分けをコツコツやればできるかもしれませんが、マルコフ決定過程でうまくモデリングできることに気がついたので、これでやってみます。
モデリングに間違いがありそうでしたらコメントお願いします。

MDPがよく使われる機械学習の文脈とは違って、状態遷移確率が既知なので、価値反復法を安直に実行すればよいです。
コードは最下部に置いておきます。

大石さんの記事によると、攻撃対象の送金額が23832モナ、覆ったブロック高が8-24ブロック前後となっています。
この攻撃に十分インセンティブが出てくるようなハッシュレート（シェア）はどの程度でしょうか。
下図は送金額23832、承認数10のときのインセンティブをプロットしたものですが、シェア25%程度からインセンティブが十分に出てくることがわかります。

2018/05/23編集　送金額を1桁多く間違っておりましたので各数値とグラフを修正しております。解析自体の修正ではなくパラメータの修正です。

2018/05/31 100%までの図を追記

ちなみにシェア20%時の最適政策を見てみますと下記のようになっています。右方向が攻撃チェーンの長さ、下方向がメインチェーンの長さ、攻撃しつづけたほうが良い場合は値1、メインチェーンに参加したほうがいい場合は値0になります。
確かにメインチェーンが伸びてきたときに撤退するような戦略になっていることがわかります。

[[1 1 1 1 1 1 1 1 1 1 0]
 [1 1 1 1 1 1 1 1 1 1 0]
 [0 1 1 1 1 1 1 1 1 1 0]
 [0 1 1 1 1 1 1 1 1 1 0]
 [0 0 1 1 1 1 1 1 1 1 0]
 [0 0 1 1 1 1 1 1 1 1 0]
 [0 0 0 1 1 1 1 1 1 1 0]
 [0 0 0 0 1 1 1 1 1 1 0]
 [0 0 0 0 1 1 1 1 1 1 0]
 [0 0 0 0 0 1 1 1 1 1 0]
 [0 0 0 0 0 0 0 0 0 0 0]]

なお、参照先の分析記事と大体同じ条件で分析した結果も貼っておきます。
細かい条件はコードを参照ください。

以下コードです。
読みやすさを優先して計算効率などは考慮してません（承認ブロック数が50以上になると計算に時間がかかります）。

withholding_attack.py

"""
"Block Withholding Attack" をマルコフ決定過程で解析

ノーテーションについて

action: 0->normal mining, 1->selfish mining

T,Rのインデックス:
action, hornest_block(t), attacker_block(t), hornest_block(t+1), attacker_block(t+1)

Vのインデックス:
hornest_block, attacker_block

"""

import numpy as np

def value_function(r,w,q,x,n):
    """
    Args:
        r : 1ブロックあたりの報酬
        w : 1ブロックあたりの全体の計算費用
        q : 攻撃者のハッシュレートの割合
        x : 攻撃者が無効化しようとする送金額
        n : 承認までのブロック数
    Returns:
        V : 状態価値関数 V(honest_block,attacker_block)
        policy : 最適政策
    """
    gamma = 1.0 # MDP割引率
    epsilon = 0.001 # 価値反復の収束判定用
    Niter = 1000

    # transition matrix
    T = np.zeros((2,n+1,n+1,n+1,n+1))
    for i in range(n):
        for j in range(n):
            T[0,i,j,i+1,j  ] = 1.0
            T[1,i,j,i+1,j  ] = 1-q
            T[1,i,j,i  ,j+1] = q
    # terminal state
    for i in range(n):
        T[0,i,n,i,n] = 1.0
        T[1,i,n,i,n] = 1.0
    for j in range(n):
        T[0,n,j,n,j] = 1.0
        T[1,n,j,n,j] = 1.0

    R = np.zeros((2,n+1,n+1,n+1,n+1))
    for i in range(n):
        for j in range(n):
            R[0,i,j,i+1,j  ] = q*r-q*w
            R[1,i,j,i+1,j  ] = -q*w
            if j+1==n:
                R[1,i,j,i  ,j+1] = x+n*r-q*w # withholding attack 成功
            else:
                R[1,i,j,i  ,j+1] = -q*w # withholdingしているあいだは報酬なし

    # value iteration method
    V = np.zeros((n+1,n+1))
    for it in range(Niter):
        t = T*( R+gamma*V.reshape([1,1,1,n+1,n+1]) )
        V1 = np.max( np.sum(t,axis=(-1,-2)), axis=0 )
        if np.sum(np.abs(V1-V)) < epsilon:
            V = V1
            break
        V = V1

    # derive policy
    t = T*( R+gamma*V.reshape([1,1,1,n+1,n+1]) )
    policy = np.argmax( np.sum(t,axis=(-1,-2)), axis=0 )

    return V,policy

def honest_value(r,w,q,n):
    return n*q*(r-w)

import matplotlib.pyplot as plt

## q-n 1d graph
r = 25
w = 20
x = 23832
n = 10

X = []
Y = []
for i in range(1,41):
    q = i/100
    V,P = value_function(r,w,q,x,n)
    print(q)
    print(P)

    incentive = V[0,0]-honest_value(r,w,q,n)

    X.append(q)
    Y.append(incentive)

plt.title("withholding attack incentive")
plt.xlabel("hash share")
plt.ylabel("incentive")
plt.plot(X,Y)
plt.show()


## q-n graph
r = 25
w = 20
x = 100
H = 1750
div = 20
xs,ys = np.meshgrid( np.linspace(1, 1000, div), np.linspace(1, 20, div) )
zs = np.zeros((div,div))
for i in range(div):
    for j in range(div):
        q = xs[i,j]/(xs[i,j]+H)
        n = int(ys[i,j])
        V0 = honest_value(r,w,q,n)
        V,P = value_function(r,w,q,x,n)
        incentive = V[0,0]-V0
        zs[i,j] = incentive
        print(i,j)
fig = plt.figure()
plt.title("withholding attack incentive")
plt.xlabel("hashrate")
plt.ylabel("confirmation")
im = plt.pcolor(xs,ys,zs)
fig.colorbar(im)
plt.show()

## x-n graph
r = 25
w = 20
q = 600/(600+1750)
div = 20
xs,ys = np.meshgrid( np.linspace(1, 1000, div), np.linspace(1, 20, div) )
zs = np.zeros((div,div))
for i in range(div):
    for j in range(div):
        x = xs[i,j]
        n = int(ys[i,j])
        V0 = honest_value(r,w,q,n)
        V,P = value_function(r,w,q,x,n)
        incentive = V[0,0]-V0
        zs[i,j] = incentive
        print(i,j)
fig = plt.figure()
plt.title("withholding attack incentive")
plt.xlabel("transaction amount")
plt.ylabel("confirmation")
im = plt.pcolor(xs,ys,zs)
fig.colorbar(im)
plt.show()