2019-08-25更新：ノイズ強度$q_t$を動的に変更するように変更。
2020-03-21更新：予測更新の式の見直しほか、ノーテーションの整理など。
導入
難易度調整アルゴリズム（Difficulty Adjustment Algorithm: DAA）についての論文"A Lucas Critique to the Difficulty Adjustment Algorithm of the Bitcoin System" (日本語解説)を読みました。PoWではブロックタイムを保つために難易度調整を行いますが、難易度を動かすとマイナーの収益性が変化するためにハッシュレート自体も動いてしまいます（こういう状況を弾力的であると表現します）。このような状態の理論解析と観測データからのシミュレーションが展開されており、結果としてブロックタイムが不安定になる場合があるという主張をしています。
これを読んでからDAAについて少し調べたのですが、様々な手法が提案されているにもかかわらず、その理論的な背景についてはあまり語られていない気がします（ちゃんと調べたわけではないので知っている方は教えてください）。
そこで今回は、ハッシュレートの逐次ベイズ推定（近似あり）により難易度調整アルゴリズムを導出し、簡単なシミュレーションを行ってみました。
仮定
ブロック$t$において、ハッシュレート$h_t$は状態変数$l_t$で表現することにします
\log h_t = l_t
$l_t$は、コイン価格、マイニング報酬、電気代、採掘機材のシェアなどに依存しますが、これらはプロトコル内から見えないので、ランダムに変動するものとして大雑把にとらえて、時変のガウス分布に従うとします。時間変動は次の状態方程式で記述します（$T$ブロックごとの更新を想定します）
l_{t+T} = l_t + v_t
ここで、$v_t$は次の分布から得られるガウシアンホワイトノイズとします
v_t \sim N(0,q_t)
PoWによるブロック発見はおおよそポアソン過程とみなせます。このとき$T$ブロックの合計ブロックタイム$B_t$の分布はガンマ分布となり、（ハッシュ・時間あたりの）当選確率を$w_t$として、次のようになります。
B_t \sim \frac{(w_t h_t)T}{\Gamma(T)} B_t^{T-1} \exp(-w_t h_t B_t)
状態推定の方法
予測更新
状態遷移による状態変数の平均・分散の動きは状態方程式から導くことができます
\begin{align}
\bar{l}_{t+T} &= \hat{l}_t \\
\bar{s}_{t+T} &= \hat{s}_t + q_t
\end{align}
このモデルは実時間では一定間隔の更新にならないため、$q_t$は経過時間に比例する形にしておきます
q_t = q B_t
ハッシュレート急落時はブロック発見が遅く、急騰時はブロック発見が速くなりますが、上記の取扱いにより時間あたりの更新幅を調整するようながあります。
ここはハッシュレート急落時の調整・急騰時の調整どちらを重視するかによって、変更してもよいかもれません。
計測更新
$T$ブロックの合計ブロックタイム$B_t$を計測したときの状態変数の事後分布は、ベイズの定理から、
P(l_t|B_t) = \eta P(B_t|l_t)P(l_t) = \eta' \exp(-J_t)
\\
J_t = -T l_t
+ w_t \exp(l_t) B_t
+ \frac{(l_t-\bar{l}_t)^2}{2\bar{s}_t}
右辺第二項を状態量の予測値周辺で2次まで近似すると、
\exp(l_t) \simeq
\exp(\bar{l}_t) \{1 + (l_t - \bar{l}_t) + \frac{1}{2}(l_t - \bar{l}_t)^2 \}
事後分布の平均値$\hat{l}_t$は、$J_t$が極小となる条件から求まります
\frac{\partial J_t}{\partial{l_t}} = 0
また、事後分布の分散$\hat{S}_t$は$J_t$の二階偏微分から求まります
\frac{\partial^2 J}{{\partial l_t}^2} = \hat{s}_t^{-1}
これらの式をまとめて、計測更新式が得られます
\hat{s}_t = (\bar{s}_t^{-1} + w_t \exp(\bar{l}_t) B_t)^{-1}
\\
\hat{l}_t = \bar{l}_t + (T-w_t \exp(\bar{l}_t) B_t)\hat{s}_t
難易度更新
上記の更新式で得られたハッシュレートの推定値に対する最適な当選確率を次の当選確率とします（$B_0$は基準ブロックタイム）
w_{t+T} = \frac{1}{\exp(\hat{l}_t) B_0}
参考実装
pythonで実装するとこんな感じです
class LevelDAA(object):
def __init__(self, q, w0, s0, T=1, b0=600):
self.T = T
self.count = 0
self.bsum = 0
self.b0 = b0
self.q = q
self.l = np.log(1/w0/b0)
self.s = s0
def winning_rate(self, w, b):
self.bsum += b
self.count = (self.count + 1)%self.T
if self.count == 0:
sp = self.s+self.q*b/self.b0
r = w*np.exp(self.l)
self.s = 1/(1/sp+r*self.bsum)
self.l += self.s*(self.T-r*self.bsum)
w = 1/np.exp(self.l)/self.b0
self.bsum = 0
return w
簡易シミュレーション
シミュレーションに用いたコードはこちらですので、細かい点はコードを参照してください。
均衡する当選確率がある時点で急上昇するような条件設定で、各アルゴリズムの当選確率の推移を図示します。オレンジ色の線が正解値です。
DAA1, DAA2はそれぞれBitcoin, BitcoinCashのアルゴリズムに類するもので、名称は冒頭の論文での記載にならっています。
今回考えたLevelDAAは、Tが小さいと即応性がありますが、Tが大きいと外れ値の影響が小さくなります。ここでは$T=1$と$T=36$の場合を用意しました。
ハッシュレートが弾力的でない場合
この場合はどのアルゴリズムでも問題なく動作しています。
ハッシュレートが弾力的な場合
DAA1は発散してしまいました。DAA2も健闘していますがオーバーシュートがあります。
ハッシュレートが弾力的な場合の仮想ステップ応答
冒頭の論文では理論的な議論がありますが、簡単にシミュレーションしてみます。
まとめ
無理のない仮定からシンプルで実用的と思われる難易度調整アルゴリズムを導出することができました。
実際に利用するにあたっては、様々な変動や攻撃を想定してパラメータの決定や安定性の分析をする必要があるかと思います。
（おまけ）
ちなみに、今回の定式化とほぼ同様の手続きで、ハッシュレートと弾力性を同時に推定するアルゴリズムを作ることができますが、推定が遅かったり精度が悪かったりして、安定して動作するものにはなりませんでした。