【Python】行列を使ってフィボナッチ数を対数時間で解く

April 28th, 2021

フィボナッチ数列とは

フィボナッチ数列は以下のような漸化式で表される数列のこと。


F1=F2=1Fn+2=Fn+1+FnF_1 = F_2 = 1 \\ F_{n + 2} = F_{n + 1} + F_n

再帰的な構造を持つことからプログラミング初学者の間でよく解かされがちなので、見た事がある方も多くいらっしゃるだろう。

Pythonで再帰関数を用いてフィボナッチ数列のNN項目を求めるには、次のように書く。


def fib(n: int) -> int:
    if n < 2:
        return n
    return fib(n-2) + fib(n-1)

再帰とは関数の中で自分自身を呼び出す処理のことで、上の例で言うと受け取ったnnの値が22未満になるまではこのfibfibという関数が呼び出される事になる。

しかし、この方法は指数関数的に計算量が増加する非効率なアルゴリズムであることが知られている。

なぜこのアルゴリズムの効率が悪いかは、次の図のように視覚的に確認すると分かりやすい。

フィボナッチの指数的爆発


fibfib関数に44という引数が与えられた場合、fib(4)fib(4)を計算するためにfib(3)fib(3)fib(2)fib(2)の計算が必要になるといった具合に、この木構造が分岐するたびに計算量が増加することが視覚的に確認できる。

この再帰アルゴリズムの計算量の改善を図る手法として、メモ化というアルゴリズムが存在する。


メモ化再帰

メモ化とは計算結果を保存しておいて、その数値が必要になった際に都度計算し直すのではなくその数値を参照するという手法で、この手法を用いれば計算量はO(N2)O(N^2)からO(N)O(N)まで減らす事ができる。

先程の図の例で言うと、木の最初の分岐点で計算が必要となるfib()fib(2)の値は、fib(3)fib(3)を計算する途中で算出することができるので、同じ計算をするのは二度手間であるから、この計算の値をどこかに保存しておけば、必要な計算はNN回で済むという理屈だ。

メモ化を実際にプログラムに落とす場合は計算結果を連想配列などに記録させる処理が必要なのだが、Pythonでは関数の前に@lru_cacheと書くだけで自動でメモ化を行ってくれる。


from functools import lru_cache

@lru_cache()
def fib(n: int) -> int:
    if n < 2:
        return n
    return fib(n-2) + fib(n-1)

ここら辺までが大体私が大学でアルゴリズムとデータ構造の授業を履修した際に初回の授業で習った事なのだが、 行列を使えば実はもっと高速に(具体的には対数時間で)フィボナッチ数列のNN項目を求める事ができる。


フィボナッチ数列を行列式で求める

フィボナッチ数列は、行列を用いて次のように表す事ができる。


(1110)n=(Fn+1FnFnFn1)\left( \begin{matrix} 1 & 1 \\ 1 & 0 \end{matrix} \right)^n = \left( \begin{matrix} F_{n+1} & F_n \\ F_n & F_{n-1} \end{matrix} \right)

上記の性質は、帰納法を用いれば簡単に証明できる。


(1110)n+1=(1110)n(1110)=(Fn+1FnFnFn1)(1110)=(Fn+1+FnFn+1Fn+Fn1Fn)=(Fn+2Fn+1Fn+1Fn)\left( \begin{matrix} 1 & 1 \\ 1 & 0 \end{matrix} \right)^{n + 1} = \left( \begin{matrix} 1 & 1 \\ 1 & 0 \end{matrix} \right)^{n} \left( \begin{matrix} 1 & 1 \\ 1 & 0 \end{matrix} \right) \\ = \left( \begin{matrix} F_{n+1} & F_n \\ F_n & F_{n-1} \end{matrix} \right) \left( \begin{matrix} 1 & 1 \\ 1 & 0 \end{matrix} \right) \\ = \left( \begin{matrix} F_{n+1} + F_{n} & F_{n+1} \\ F_n + F_{n-1} & F_{n} \end{matrix} \right) \\ = \left( \begin{matrix} F_{n+2} & F_{n+1} \\ F_{n+1} & F_{n} \end{matrix} \right)

このとき(1110)\left( \begin{matrix} 1 & 1 \\ 1 & 0 \end{matrix} \right)QQとおくと、フィボナッチ数列のNN項目はQn1Q^{n-1}の1行1列目と等しい事が分かる。

つまり、FnF^nを求めるにはQn1Q^{n-1}を求めればよいと言うことになる。


再び再帰

ではどうすればQn1Q^{n-1}を対数時間で求めることができるのだろうか。

xのn乗というのは、Exponentiation by squaringという手法を用いれば対数時間で求める事ができる事が知られており、次のように表すことができる。


xn={x(x2)n1/2(nが奇数のとき)(x2)n/2(nが隅数のとき)x^n = \left\{ \begin{matrix}x(x^2)^{n-1 / 2}&(nが奇数のとき)\\ (x^2)^{n / 2}&(nが隅数のとき) \end{matrix} \right.

つまり、QnQ^{n}という行列を求めたければQn/2Q^{n / 2}を求めればよく、Qn/2Q^{n / 2}を求めたければQn/4Q^{n / 4}を求めれば良いことになる。

よって、時間計算量はO(logN)O(log N)となる。

それでは、実際のPythonコードを見てみよう。


Pythonのサンプルコード

まず初めに、行列QQと、行列QQN1N-1乗の1行1列目を返すfib関数を定義する。


Q1 = [1, 1, 
      1, 0]

def fib(n):
    return power(Q1, n-1)[0]

次に、fib関数の中で内部的に使用されているpower関数を定義する。power関数は先ほど紹介したExponentiation by squaringという手法を用いる。

なお、基底部となる行列の00乗にあたる部分だが、行列の00乗は単位行列EE (1001)\left( \begin{matrix} 1 & 0 \\ 0 & 1 \end{matrix} \right) であるので、そちらを用意する。


def power(matrix, m):
    if m == 0:
        return [1, 0, 0, 1]
    elif m == 1:
        return matrix
    else:
        matrix_tmp = matrix
        n = 2
        while n <= m:
            matrix_tmp = multiply(matrix_tmp, matrix_tmp)
            n = n*2
        R = power(matrix, m-n//2)
        return multiply(matrix_tmp, R)

そして最後に、power関数の内部的に使われている行列の掛け算を計算するmultiply関数を定義すればよい。


def multiply(matrix_a, matrix_b):
    a, b, c, d = matrix_a
    x, y, z, w = matrix_b

    return (
        a*x + b*z,
        a*y + b*w,
        c*x + d*z,
        c*y + d*w,
    )

print(fib(50)) # 12586269025

以上が、行列を使用してフィボナッチ数をO(logN)O(log N)で計算する方法である。


参考