<jupyter_start><jupyter_text>ここで、量子線形システム問題にまつわる少し発展的な話題を紹介する。これらのコラムの内容は本編の理解に必須ではないので、難しくて分からない場合は適宜飛ばして次の章に進んでほしい。 コラム：低ランク行列に対する高速な特異値分解とサンプリング（量子-inspiredアルゴリズム）<jupyter_code>import sys
import time
from collections import Counter

import numpy as np
import sklearn
from sklearn.decomposition import TruncatedSVD
import matplotlib
import matplotlib.pyplot as plt

print("python: %s"%sys.version)
print("numpy: %s"%np.version.version)
print("matplotlib: %s"%matplotlib.__version__)
print("scikit-learn: %s"%sklearn.__version__)<jupyter_output>python: 3.7.5 (v3.7.5:5c02a39a0b, Oct 14 2019, 18:49:57) 
[Clang 6.0 (clang-600.0.57)]
numpy: 1.17.3
matplotlib: 3.1.1
scikit-learn: 0.21.3<jupyter_text>概要低ランクなハミルトニアン $H$ によるダイナミクス $e^{-iHt}$ は量子計算機で効率的にシミュレートすることができる。ハミルトニアンの効率的なシミュレートは量子位相推定アルゴリズムやHHLアルゴリズムを行うための前提条件であるから、量子計算機は低ランクな行列に対して種々の高速なアルゴリズムを自然に実行できる。量子推薦システム (Quantum recommendation system) はそのような応用の一つで、量子計算機は低ランクな行列の任意の行ベクトルに対し、Top-$k$特異値空間へ射影したベクトルからのL2サンプリングを行えるというものである。このサンプリング操作は、例えばWebサービスにおける購買履歴に基づくユーザへのアイテム推薦の操作に対応しており、、量子計算機は高速な推薦システムとして活用することができる。（注：ユーザへの推薦は、行がユーザ・列がアイテム・要素が指定ユーザが指定アイテムを好きかどうか、という嗜好行列（通常は低ランクだと仮定される）に基づいて行われる）しかし、一般に問題に低ランク近似可能性などの強い近似が加わった場合、通常の古典計算機での計算も高速になる点に注意しなければならない。近年、上記の量子推薦システムが行う操作は、種々のデータ構造やアルゴリズムを用いることで通常の計算機でも多項式時間で実現可能だということがわかった。このコラムでは、低ランクな行列に対する打ち切り特異値分解、および特異値空間に射影されたベクトルからのサンプリングが、通常の計算機だけでも行列のサイズに対して従来より指数的に高速になるという点について解説する。この一連の研究は主にEwin Tangによって行われたものである[1][2][3]。この研究は以下の個別の主要要素をつなげたものである。1. セグメント木2. 確率的操作を用いた打ち切り特異値分解3. 棄却サンプリング4. 内積の推定まず、問題設定を提示し、そののち上記の三つについてそれぞれ解説し、これを組みあわせて低ランク近似した特徴量からのサンプリングが可能であることをコードを実行しながら概観する。  問題設定低ランク近似からの特異値空間の射影サンプリングは以下のような問題として定義される。初期状態として、$2^n, 2^m$の大きさの行列$A$を考える。行列$A$は全ての要素が0である。この時、以下の二つの操作がランダムな順序で$q$個要求される。1. 値の更新    $A$の$i$行$j$列の要素を$v$に更新する。2. 射影サンプリング   $A$の$i$行目の成分を$A$のTop-$k$特異値空間に射影したベクトル$v$を得る。ベクトル$v$の$j$番目の要素を$\frac{v_j^2}{\sum_l v_l^2}$の確率で出力する。ただし、$k$は高々$O(n,m)$である。上記に$\epsilon$の誤差を許容した時、高速に処理を行いたい。自明な上記のアルゴリズムは$2^n, 2^m$の配列をスパース形式で確保し、挿入を$O(n,m)$で行い、射影サンプリングは特異値分解を通して${\rm poly}(2^n,2^m)$で行う方法である。これを改善し、${\rm poly}(q,n,m,\epsilon^{-1})$の計算量で行いたい。量子計算機を用いると、上記の計算を${\rm poly}(q,n,m,\epsilon^{-1})$で行うことができる。これは、HHLアルゴリズムの一部を変更し、Quantum projectionというプロトコルを実施することで可能となる。近年のEwin Tangによる論文は、この計算が古典計算機でも同じく${\rm poly}(q,n,m,\epsilon^{-1})$で可能であるというものである。  和のセグメント木通常、我々がプログラムで配列を確保すると、メモリ上の連続的なアドレスに指定した容量の数列が確保される。ランダムアクセスメモリでは読み出しと書き込みはO(1)で完了できるため、通常はこの方法で確保される。一方、配列に対する要求が読み出しと書き込みだけではなく挿入や削除、サンプリングなどがある場合、こうしたデータの保持の仕方は必ずしも高速とは限らない。従って、データに対する要求が単なる読み出しと書き込みに限らない場合、その目的に適した特殊なデータの保持方法が用いられる。これをデータ構造と呼ぶ。和に関するセグメント木はデータ構造の一種であり、データの更新、区間和の計算やサンプリングを高速に行うことができる。対象となるデータの大きさを$N=2^n$とする。２の累乗でない場合は、累乗になるまで末尾を0で埋める。、まず、$2N$の長さの配列を確保する。最初は$2N$個の要素は0で埋まっている。和のセグメント木では以下の4つの操作を行える。1. 値の更新: $i$番目のデータに値$v$を書き込む。2. 値の取得: $i$番目のデータを取り出す。3. 総和の取得: 全データの和を取り出す。4. サンプリング: $i$番目のデータを$\frac{a_i}{\sum_j a_j}$の確率で取り出す。上記を通常の配列で愚直に行うと以下のようなコードになる。<jupyter_code>class VectorSampler():
    def __init__(self,size):
        self.size = size
        self.array = np.zeros(size)
        self.sum = 0
        self.index_list = np.arange(self.size)
        
    def update(self,index,value):
        self.sum += value - self.array[index]
        self.array[index] = value
        
    def get_element(self,index):
        return self.array[index]
    
    def get_sum(self):
        return self.sum
    
    def sampling(self):
        normalized_array = self.array / self.sum
        result = np.random.choice(self.index_list, p = normalized_array)
        return result<jupyter_output><empty_output><jupyter_text>ランダムな読み出し/書き込み等を行って、上記の関数をベンチマークしてみよう。<jupyter_code>n_list_vec = [10**2,10**3,10**4,10**5,10**6,10**7,10**8]
update_time_vec = []
get_element_time_vec = []
get_sum_time_vec = []
sampling_time_vec = []
for n in n_list_vec:
    data_structure = VectorSampler(n)
    
    start_time = time.time()
    rep = 0
    while time.time()-start_time<0.1:
        index = np.random.randint(n)
        value = np.random.rand()
        data_structure.update(index,value)
        rep+=1
    elapsed = (time.time()-start_time)/rep
    update_time_vec.append(elapsed)
    
    start_time = time.time()
    rep = 0
    while time.time()-start_time<0.1:
        index = np.random.randint(n)
        data_structure.get_element(index)
        rep+=1
    elapsed = (time.time()-start_time)/rep
    get_element_time_vec.append(elapsed)
    
    start_time = time.time()
    rep = 0
    while time.time()-start_time<0.1:
        data_structure.get_sum()
        rep+=1
    elapsed = (time.time()-start_time)/rep
    get_sum_time_vec.append(elapsed)
    
    start_time = time.time()
    rep = 0
    while time.time()-start_time<0.1:
        data_structure.sampling()
        rep+=1
    elapsed = (time.time()-start_time)/rep
    sampling_time_vec.append(elapsed)
    
plt.plot(n_list_vec, update_time_vec,label="update")
plt.plot(n_list_vec, get_element_time_vec,label="get_element")
plt.plot(n_list_vec, get_sum_time_vec,label="get_sum")
plt.plot(n_list_vec, sampling_time_vec,label="sampling")
plt.xlabel("size")
plt.ylabel("time (sec)")
plt.xscale("log")
plt.yscale("log")
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>他の関数の時間が定数でスケールする中、samplingだけが線形に時間がかかっていることがわかる。 セグメント木では４つの命令を以下のように実装する。（図でイメージするとわかりやすい；例えばこの[ブログ](https://proc-cpuinfo.fixstars.com/2017/07/optimize-segment-tree/)を参照）1. $i$番目の値の更新   $a_{2^n + i}=v$と更新した上で（$N=2^n$で配列の長さは$2N$）、さらに全ての$j=1,2,\ldots$について $k = \lfloor \frac{2^n + i}{2^j} \rfloor, \: a_k \leftarrow a_k-a_i+v$ という更新を行う。この処理は$O(n)$で終わる。2. $i$番目の値の取得   $a_{2^n + i}$を返す。この処理は$O(1)$で終わる。3. 総和の取得   $a_0$を返す。この処理は$O(1)$で終わる。4. サンプリング   $i=1$とする。$b=0$が$\frac{a_i}{a_i + a_{i+1}}$、$b=1$が$\frac{a_{i+1}}{a_i+a_{i+1}}$の確率で得られるサンプリングを行う。$b=0$となった場合、$i \leftarrow i+2^k$と更新する。$b=1$となった場合、$i \leftarrow i+2^{k+1}$と更新する。$k \leftarrow k+1$とする。$i\geq 2^n$となったら終了し、$i-2^n$を返す。この処理は$O(n)$で終わる。従って、長さ$n$のベクトルに関し、上記四種の命令が$q$個与えられる場合、その実施を高々$O(q \log N)$で行うことができる。実装は以下のようになる。<jupyter_code>class SegmentTree(object):
    def __init__(self, size):
        _size = 1
        while _size<size:
            _size*=2
        self.size = _size*2
        self.array = np.zeros(_size*2)

    def update(self,index, value):
        position = self.size//2 + index
        difference = value - self.array[position]
        while position>0:
            self.array[position]+=difference
            position //= 2
        
    def get_element(self,index):
        return self.array[self.size//2 + index]
    
    def get_sum(self):
        return self.array[1]
                
    def sampling(self):
        current_node = 1
        while current_node<self.size//2:
            left_children = self.array[current_node*2]
            right_children = self.array[current_node*2+1]
            left_weight = left_children / (left_children + right_children)
            if np.random.rand() < left_weight:
                current_node = current_node*2
            else:
                current_node = current_node*2+1
        return current_node - self.size//2<jupyter_output><empty_output><jupyter_text>セグメント木をベンチマークしてみよう。<jupyter_code>n_list_seg = [10**2,10**3,10**4,10**5,10**6,10**7,10**8]
update_time_seg = []
get_element_time_seg = []
get_sum_time_seg = []
sampling_time_seg = []
for n in n_list_seg:
    data_structure = SegmentTree(n)
    
    start_time = time.time()
    rep = 0
    while time.time()-start_time<0.1:
        index = np.random.randint(n)
        value = np.random.rand()
        data_structure.update(index,value)
        rep+=1
    elapsed = (time.time()-start_time)/rep
    update_time_seg.append(elapsed)
    
    start_time = time.time()
    rep = 0
    while time.time()-start_time<0.1:
        index = np.random.randint(n)
        data_structure.get_element(index)
        rep+=1
    elapsed = (time.time()-start_time)/rep
    get_element_time_seg.append(elapsed)
    
    start_time = time.time()
    rep = 0
    while time.time()-start_time<0.1:
        data_structure.get_sum()
        rep+=1
    elapsed = (time.time()-start_time)/rep
    get_sum_time_seg.append(elapsed)
    
    start_time = time.time()
    rep = 0
    while time.time()-start_time<0.1:
        data_structure.sampling()
        rep+=1
    elapsed = (time.time()-start_time)/rep
    sampling_time_seg.append(elapsed)
    
plt.plot(n_list_seg, update_time_seg,label="update")
plt.plot(n_list_seg, get_element_time_seg,label="get_element")
plt.plot(n_list_seg, get_sum_time_seg,label="get_sum")
plt.plot(n_list_seg, sampling_time_seg,label="sampling")
plt.xlabel("size")
plt.ylabel("time (sec)")
plt.xscale("log")
plt.yscale("log")
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>先ほどの素朴なデータ構造の時間と比較してみよう。<jupyter_code>cmap = plt.get_cmap("tab10")
plt.plot(n_list_vec, update_time_vec,label="vector update",color=cmap(0),linestyle="--")
plt.plot(n_list_seg, update_time_seg,label="segtree update",color=cmap(0))
plt.plot(n_list_vec, get_element_time_vec,label="vector get_element",color=cmap(1),linestyle="--")
plt.plot(n_list_seg, get_element_time_seg,label="segtree get_element",color=cmap(1))
plt.plot(n_list_vec, get_sum_time_vec,label="vector get_sum",color=cmap(2),linestyle="--")
plt.plot(n_list_seg, get_sum_time_seg,label="segtree get_sum",color=cmap(2))
plt.plot(n_list_vec, sampling_time_vec,label="vector sampling",color=cmap(3),linestyle="--")
plt.plot(n_list_seg, sampling_time_seg,label="segtree sampling",color=cmap(3))
plt.xlabel("size")
plt.ylabel("time (sec)")
plt.xscale("log")
plt.yscale("log")
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>updateの関数はsegtreeの方が多少遅いものの、samplingのコストが圧倒的に減っていることがわかる。なお、本稿では登場しないが、和のセグメント木は与えられた区間に関する和を求める操作も$O(\log N)$で行うことができる。これは通常のベクトルでは最悪で$O(N)$かかる。また、セグメント木は二つの要素間の演算がモノイド、すなわちゼロ元を持ち結合則が成り立てば様々な対象に対し適用することができる。例えば$f(a,b)=\min(a,b)$とすることで区間最小値を求めるセグメント木を構築できる。セグメント木はベクトルに対するデータ構造だが、多少の拡張をすることで行列をセグメント木の集合として表すことができる。下記のような$m \times n$行列を考えよう。- $i$行目の長さ$n$のベクトル$T_i$をセグメント木$S_i$として保持する。- $i$行目のベクトルの総和を$|T_i|$を並べた長さ$m$のベクトル$C = \{|T_i|\}$をセグメント木$S'$として保持する。上記のようなデータ構造は、以下の操作を全て${\rm poly}(\log n, \log m)$で行う。1. 値の更新: $(i,j)$要素を$v$に更新する。2. 値の取得: $(i,j)$要素を取得する。3. 行ノルムの取得: 与えられた行添え字$i$に対し、$i$行目のノルム$|T_i|$を取得する。4. 行ベクトルのサンプリング: 与えられた行添え字$i$に対し、列添え字$j$を$\frac{T_{ij}}{\sum_j{ T_{ij}}}$の確率で得る。5. 総和の取得: 行列の全要素の和$\sum_{ij} T_{ij}$を取得する。6. 行ノルムのサンプリング: 行添え字$i$を$\frac{|T_i|}{\sum_i |T_i|}$の確率で得る。実装は以下のようになる。<jupyter_code>class SegmentTreeMatrix(object):
    def __init__(self, height, width):
        _height = 1
        while _height<height:
            _height*=2
        _width = 1
        while _width<width:
            _width*=2
        self.width = _width
        self.height = _height
        self.column_norm_segtree = SegmentTree(self.height)
        self.row_segtree_list = []
        for row_index in range(self.height):
            self.row_segtree_list.append(SegmentTree(self.width))
        
    def update(self,i,j, value):
        self.row_segtree_list[i].update(j,value)
        self.column_norm_segtree.update(i, self.row_segtree_list[i].get_sum())
        
    def get_element(self,i,j):
        return self.row_segtree_list[i].get_element(j)
    
    def get_row_norm(self,i):
        return self.row_segtree_list[i].get_sum()
                
    def sampling_column_in_row(self,i):
        return self.row_segtree_list[i].sampling()
    
    def get_sum(self):
        return self.column_norm_segtree.get_sum()
    
    def sampling_row(self):
        return self.column_norm_segtree.sampling()<jupyter_output><empty_output><jupyter_text>確率的な打ち切り特異値分解特異値分解とは、正方行列とは限らない行列$A$について、$A=UDV^{\dagger}$という分解を行う行列操作である。この時、$U,V$ はユニタリ行列であり、$D$は実対角行列である。また、$\dagger$はエルミート共役(転置の複素共役)を表す。$D$の対角要素を$A$の特異値と呼ぶ。行列$A$の幅と高さの小さいほうの値が$n$である場合、特異値は$n$個あり、分解の計算量は$O(n^3)$である。$U$と$V$の$i$番目の列ベクトルを、$A$の$i$番目の特異値の左特異値ベクトル、右特異値ベクトルと呼ぶ。特異値分解はnumpyに標準で搭載されている。<jupyter_code>n = 3
m = 4
A = np.random.rand(n,m)
U,singular_value_list,Vd = np.linalg.svd(A)
D = np.zeros( (n,m) )
np.fill_diagonal(D,singular_value_list)
print("A=\n{}\n".format(A))
print("U=\n{}\n".format(U))
print("D=\n{}\n".format(D))
print("V^d=\n{}\n".format(Vd))
print("UDV^d=\n{}\n".format(U@D@Vd))<jupyter_output>A=
[[0.28518719 0.33357639 0.22698753 0.89912291]
 [0.23855553 0.95448547 0.54962802 0.81592368]
 [0.98739535 0.84929664 0.7374283  0.98272594]]

U=
[[-0.3954486   0.23336956 -0.88834625]
 [-0.55539766  0.70956965  0.43364081]
 [-0.7315421  -0.66486807  0.15098541]]

D=
[[2.40625414 0.         0.         0.        ]
 [0.         0.48352874 0.         0.        ]
 [0.         0.         0.41263792 0.        ]]

V^d=
[[-0.40211518 -0.53333002 -0.38835626 -0.6348558 ]
 [-0.86998319  0.39387582 -0.09786635  0.28002369]
 [-0.00197572  0.59568912  0.35876134 -0.71863821]
 [-0.28535719 -0.45340424  0.84314405  0.04586912]]

UDV^d=
[[0.28518719 0.33357639 0.22698753 0.89912291]
 [0.23855553 0.95448547 0.54962802 0.81592368]
 [0.98739535 0.84929664 0.7374283  0.98272594]]<jupyter_text>計算コストは以下のようになる。<jupyter_code>n_list = np.arange(100,1500,100)
elapsed_times_svd = []
for n in n_list:
    matrix = np.random.rand(n,n)
    start_time = time.time()
    s,d,v = np.linalg.svd(matrix)
    elapsed = time.time()-start_time
    elapsed_times_svd.append(elapsed)

plt.plot(n_list,elapsed_times_svd)
plt.xlabel("size")
plt.ylabel("time (sec)")
plt.show()<jupyter_output><empty_output><jupyter_text>打ち切り特異値分解とは、定数$k$について、特異値の集合のうち大きいものから$k$個の特異値と特異値ベクトルを取り出す操作である。定数$k=n$である場合は通常の特異値分解と一致するため、通常は$k$が$n$よりはるかに小さい場合をしばしば考える。最も素朴な打ち切り特異値分解は、下記の$O(n^2 k)$のアルゴリズムである。1. 各要素が $N(0,1)$ の独立な正規分布に従う$n \times k$の行列$B$を作成する。計算コストは$O(nk)$。2. $C = AB$を計算する。計算コストは$O(n^2 k)$。3. $C$の行ベクトルを直交化し、行列$D$を得る。計算コストは$O(n^2 k)$。4. $E=DA$を計算する。計算コストは$O(n^2 k)$。5. $E$を特異値分解し、$E=U'DV^{\dagger}$を得る。計算コストは$O(n^2 k)$。6. $U = D^{\rm T}U'$を計算する。計算コストは$O(n^2 k)$。結果として得られるU,D,Vは高い確率で$A$の特異値分解と一致する。こうした特異値分解はscikit-learnに実装されている。<jupyter_code>k = 10
n_list_tsvd = np.arange(100,5000,100)
elapsed_times_tsvd = []
for n in n_list_tsvd:
    matrix = np.random.rand(n,n)
    tsvd = TruncatedSVD(n_components=k)
    start_time = time.time()
    tsvd.fit(matrix)
    elapsed = time.time()-start_time
    elapsed_times_tsvd.append(elapsed)

plt.plot(n_list,elapsed_times_svd,label="SVD")
plt.plot(n_list_tsvd,elapsed_times_tsvd,label="truncated SVD")
plt.legend()
plt.xlabel("size")
plt.ylabel("time (sec)")
plt.show()<jupyter_output><empty_output><jupyter_text>全ての特異値を得るFull SVDに比べ、圧倒的に高速であることがわかる。さらに、FKVアルゴリズムというものが存在し、いくつかの条件が満たされると、Truncated SVDからさらなる高速化が可能になる。手順は以下の通りである。1. 行列$A$から行ノルムを重みとして$p$個の行をサンプリングした$B$を作成する。行添え字$i$を行ノルムの重みづけで得るサンプリングのコストを$S_1(A)$として、計算コストは$O(S_1(A) p)$。2. 行列$B$から「ランダムに行を選び、選ばれた行から要素を重みとして列をサンプリングする」操作を$p$回行い$C$を作成する。行ベクトルからの要素のサンプリングコストを$S_2(A)$として、計算コストは$O(S_2(A)p)$。3. $C$を特異値分解して$C=U'DV'^{\dagger}$を得る。計算コストは$O(p^3)$。4. $V = B^T U' D^+$を計算する。計算コストは$O(np^2)$。5. $V = V^{\dagger}D^+A$を計算する。計算コストは$O(n^2 p)$。$p$を$\log n$の多項式としても、結果として得られる$U,D,V$は高い確率で$A$の特異値分解の良い近似となる。手順3. までで得られる、積を取れば$U$が得られる$(B,V',D)$のセットのことを$V$の簡潔表現と呼ぶ。従って、$V$の簡潔表現を得るための計算量は${\rm poly}(S_1(A), S_2(A), \log n)$である。行列の表現としてセグメント木行列を用いると、$S_1(A), S_2(A)$はどちらも${\rm poly}(\log n)$である。従って、FKVアルゴリズムは行列がセグメント木の形式で保持されている場合は効率的に行うことができる。FKVアルゴリズムの実装は以下のようになる。<jupyter_code>def singular_value_decomposition_FKV(segtree_matrix, k, subsample_count):
    # row sampling
    ## sampled row index list
    row_index_list = []
    ## normalization factor for each sampled row index
    row_norm_list = []
    matrix_sum = segtree_matrix.get_sum()
    for _ in range(subsample_count):
        # choose row indices, where its weight is norm of row
        row_index = segtree_matrix.sampling_row()
        row_index_list.append(row_index)
        row_chosen_prob = segtree_matrix.get_row_norm(row_index) / matrix_sum
        row_norm_list.append(subsample_count * row_chosen_prob)

    # column sampling
    ## sampled column index list
    col_index_list = []
    ## normalization factor for each sampled column index
    col_norm_list = []
    ## uniformly choose row index from sampled index for sampling column
    col_sample_row = np.random.choice(row_index_list,size=subsample_count)
    for row_index in col_sample_row:
        # choose column indices, where its weight is element at chosen row/column.
        col_index = segtree_matrix.sampling_column_in_row(row_index)
        col_index_list.append(col_index)
        col_chosen_prob = 0.
        # compute sum of element in chosen column in row-sampled matrix
        for row_index_for_ave in row_index_list:
            col_chosen_prob += segtree_matrix.get_element(row_index_for_ave, col_index) / segtree_matrix.get_row_norm(row_index_for_ave)
        col_norm_list.append(col_chosen_prob)

    # pick normalized element from original matrix
    subsample_matrix = np.zeros( (subsample_count, subsample_count) )
    for sub_row in range(subsample_count):
        for sub_col in range(subsample_count):
            org_row = row_index_list[sub_row]
            org_col = col_index_list[sub_col]
            element = segtree_matrix.get_element(org_row, org_col)
            norm_factor = row_norm_list[sub_row] * col_norm_list[sub_col]
            subsample_matrix[sub_row,sub_col] = element / norm_factor

    # perform SVD for sampled matrix
    ## since segtree matrix contains squared values as elements, take sqrt for svd
    u,d,_ = np.linalg.svd(np.sqrt(subsample_matrix))
    
    # take top-k matrix
    d = d[:k]
    u = u[:,:k]
    return (row_index_list, d, u)<jupyter_output><empty_output><jupyter_text>下のコードで示すように、量子状態のボルンの規則と合わせるため、量子状態の値はL2-normでサンプリングされる。セグメント木は和に対して定義されるため、総和がL2-normの二乗となるよう、セグメント木行列の要素には嗜好行列の要素の事情の値が格納される。従って、サブ行列を作成した後に特異値分解する際には、全要素のsqrtを取る必要がある。このことを踏まえ、早速FKV algorithmを試してみよう。<jupyter_code>n = m = 200
k = 10
p = 50
dense_matrix = np.zeros( (m,n) )
seg_tree = SegmentTreeMatrix(height = m, width = n)
for y in range(m):
    for x in range(n):
        val = np.random.rand()
        seg_tree.update(y,x,val**2)
        dense_matrix[y,x] = val

sampled_row_index_list, d, u = singular_value_decomposition_FKV(seg_tree, k, p)
print("FKV SVD\n",d)
_, d , _ = np.linalg.svd(dense_matrix)
print("Full SVD\n",d[:k])
tsvd = TruncatedSVD(n_components=k, n_iter=1)
tsvd.fit(dense_matrix)
print("Trunc SVD\n",tsvd.singular_values_)<jupyter_output>FKV SVD
 [100.72014493  17.81097382  16.70538197  15.54469501  14.7001673
  14.34469146  13.3012425   12.74757961  12.0585889   11.34122011]
Full SVD
 [100.23396543   7.90274002   7.79754658   7.74914095   7.50160583
   7.48295602   7.40332841   7.34452346   7.29000909   7.23951272]
Trunc SVD
 [100.2339545    7.55592118   7.34886176   7.19372274   7.05660533
   7.01074137   6.93760024   6.77137369   6.6663421    6.4976206 ]<jupyter_text>簡単な検証のため、$200 \times 200$の行列から$50 \times 50$のサブ行列を作成し、Top-10の特異値を表示している。最大特異値はおおむねあっているものの、以降の特異値でずれがひどいことが分かる。次に$200 \times 200$の行列から$1000 \times 1000$のサブ行列(というのか分からないが)を作成し、Top-10の特異値を表示している。当たり前だが、この計算自体は元の行列より大きな行列をサンプリングしているため有用ではない操作である。<jupyter_code>n = m = 200
k = 10
p = 1000
dense_matrix = np.zeros( (m,n) )
seg_tree = SegmentTreeMatrix(height = m, width = n)
for y in range(m):
    for x in range(n):
        val = np.random.rand()
        seg_tree.update(y,x,val**2)
        dense_matrix[y,x] = val

sampled_row_index_list, d, u = singular_value_decomposition_FKV(seg_tree, k, p)
print("FKV SVD\n",d)
_, d , _ = np.linalg.svd(dense_matrix)
print("Full SVD\n",d[:k])
tsvd = TruncatedSVD(n_components=k, n_iter=1)
tsvd.fit(dense_matrix)
print("Trunc SVD\n",tsvd.singular_values_)<jupyter_output>FKV SVD
 [100.96870017   9.44054781   9.11840302   8.99198792   8.63975411
   8.3835851    8.1716375    8.02057196   7.81662004   7.7810329 ]
Full SVD
 [100.94752093   8.02086558   7.80305145   7.77010671   7.60381687
   7.56895073   7.42212447   7.33142897   7.30323491   7.20036441]
Trunc SVD
 [100.94751363   7.57918791   7.35764553   7.22657673   7.13785688
   7.01050535   6.86475243   6.82990868   6.66014474   6.50609105]<jupyter_text>まだ十分とは言えないが、正しい値に漸近的に近づいていることが分かる。実は、FKVで必要となるサンプリング数$p$は$10^7 \times \max \left(\frac{k^{5.5}}{\epsilon^{17}}, \frac{k^4}{\epsilon^{24}} \right)$というとんでもない多項式になっている。従って、$k=10, \epsilon=0.1$に対して理論的な精度保証を得るためには、$p=10^{29.5}$というとんでもない数のサイズの行列が要請される。これが操作として意味を持つには$n,m > p$である必要があり、その行列のサイズは各要素が1bitだとしても$10^{59}$byteとなり、行ですら$10^{14}$PBという到底保持できないサイズとなる。もちろん、上記の多項式や係数は理論的な証明のためのオーバーヘッドを含むため、実際にはそれほど大きくしなくともそこそこの結果は得られると思われるが、いずれにせよ**無条件に多項式だから有用だと言える手法ではない**ことは念頭に入れておきたい。  棄却サンプリング一般に確率分布は高速にサンプリング可能とは限らない。例えば$2^n$個の要素のうち$i$番目の要素が$p_i$の確率で得られる分布で、$i$に対して$p_i$が効率的に計算できるとしても、サンプリングには一般に指数的な時間を要する。棄却サンプリングは上記のような効率的にサンプリングできるとは限らない分布$P$があるとき、$P$と類似している効率的にサンプリング可能な分布$Q$を用いて$P$からのサンプリングを効率的に行う手法である。この時、最終的にサンプリングしたい分布$P$を目的分布、目的分布に近いと期待される効率的にサンプリング可能な分布$Q$を提案分布と呼ぶ。具体的な手続きは以下である。1. $M \geq \max_i \left( \frac{p_i}{q_i} \right)$となる$M$を求める。2. $Q$からサンプリングを行い、$i$を得る。3. $t = \frac{p_i}{Mq_i}$を計算する。4. [0,1]の区間の一様乱数 $r$ を得る。$r<t$なら$i$を出力する。そうでなければ1.からやり直す。上記の手続きで得られるサンプリング操作は$Q$からのサンプリング操作と一致する。上記の操作に必要な計算量は$Q$からのサンプリングの操作が効率的なら$O(M)$の期待値で終了する。効率の最良のケースは確率分布$P$と$Q$が一致している場合であり、$M=1$となる。ただしこの場合はそもそも棄却サンプリングを行う意味がない。最悪のケースは$p_i > 0$でかつ$q_i=0$となるような$i$が存在する場合である。この場合、$M = \infty$となり棄却サンプリングはそもそも実行できない。例として、一様分布に十分近い$n$個のシンボルについての確率分布$\{p_i\}$を考えよう。この分布は全てのindexについて、おおよそ$p_i < 1.1 (1/n)$であることが保障されているとする。$\{p_i\}$からサンプリングするには最悪で$O(n)$が必要になるが、一様分布は$\log n$の速度で高速にサンプリングすることができる。このことを利用して、$\{p_i\}$からのサンプリングを一様分布のサンプリングを用いて棄却サンプリングしてみよう。(以下の実装では規格化のために微妙に$M=1.1$からずれているが、おおむね$M=1.1$となる。)<jupyter_code>def direct_sampling(prob_dist,index_list):
    index = np.random.choice(index_list, p = prob_dist)
    return index

def rejection_sampling(prob_dist, M, n):
    sampling_count = 0
    while True:
        index = np.random.randint(n)
        sampling_count+=1
        thr = prob_dist[index]/(M /n)
        r = np.random.rand()
        if r<thr:
            return index, sampling_count<jupyter_output><empty_output><jupyter_text>まず、棄却サンプリングからサンプリングされる分布が、オリジナルの分布を再現しているかを検証してみよう。<jupyter_code>n = 10
prob_dist = (1./n) * (1.+ (np.random.rand(n)-0.5)*0.1)
prob_dist /= np.sum(prob_dist)
M = np.max(prob_dist / (1./n))
index_list = np.arange(n)
sample_count = 10**6

sample_rejection = []
for _ in range(sample_count):
    sample, _ =rejection_sampling(prob_dist, M, n)
    sample_rejection.append(sample)
counter = Counter(sample_rejection)

rejection_probs = []
for i in index_list:
    rejection_probs.append(counter[i]/sample_count)
plt.plot(index_list, prob_dist, label="target")
plt.plot(index_list, [1/n]*n, label="uniform",linestyle = "--")
plt.errorbar(index_list, rejection_probs, yerr=1./np.sqrt(sample_count), label="rejection")
plt.legend()
plt.xlabel("synbols")
plt.ylabel("probability")
plt.show()<jupyter_output><empty_output><jupyter_text>確かに棄却サンプリングから得られる分布と目的分布が一致していることがわかる。次に、通常のサンプリングとの速度比較のベンチマークを取ってみよう。<jupyter_code>n_list = [10,10**2,10**3,10**4,10**5,10**6,10**7]
direct_sampling_time = []
rejection_sampling_time = []
sampling_count_list = []
M_list = []

for n in n_list:
    prob_dist = (1./n) * (1.+ (np.random.rand(n)-0.5)*0.1)
    prob_dist /= np.sum(prob_dist)
    M = np.max(prob_dist / (1./n))
    index_list = np.arange(n)

    start_time = time.time()
    rep = 0
    while time.time()-start_time < 0.1:
        index = direct_sampling(prob_dist,index_list)
        rep+=1
    elapsed_direct = (time.time()-start_time)/rep

    start_time = time.time()
    rep = 0
    sampling_count_average = 0
    while time.time()-start_time < 0.1:
        index, sampling_count = rejection_sampling(prob_dist, M, n)
        sampling_count_average += sampling_count
        rep+=1
    elapsed_rejection = (time.time()-start_time)/rep
    sampling_count_average /= rep

    direct_sampling_time.append(elapsed_direct)
    rejection_sampling_time.append(elapsed_rejection)
    M_list.append(M)
    sampling_count_list.append(sampling_count_average)
    
plt.plot(n_list, direct_sampling_time,label="direct")
plt.plot(n_list, rejection_sampling_time,label="rejection")
plt.legend()
plt.xlabel("size")
plt.ylabel("time (sec)")
plt.xscale("log")
plt.yscale("log")

print("M={}".format(np.mean(M_list)))
print("average repetition in rejection sampling = {}".format(np.mean(sampling_count_list)))<jupyter_output>M=1.049485039045078
average repetition in rejection sampling = 1.050053426899406<jupyter_text>良いサンプラーを持っている場合、直接的なサンプリングに比べ、棄却サンプリングは圧倒的に高速であることがわかる。また、$M$の値とサンプリング回数がおおむね一致している。棄却サンプリングを行うにあたって注意しなければならない点が二つある。一つ目は、棄却サンプリングを用いる代わりに、目的分布の確率ベクトルに関する累積和やセグメント木などのデータ構造を用いることでも高速にサンプリングが可能になるという可能性が存在する点である。例えばセグメント木を与えられたベクトルで初期化する手続きには$O(n)$が必要となる。サンプリングの回数を$s$とした場合、データ構造を用いて以降のサンプリングを効率化した場合の計算量は$O(n + s \log n)$であり、棄却サンプリングを行った場合は$O(M s \log n)$である。従って、棄却サンプリングを行う利点がある場合は$$\frac{n}{\log n} > (M-1) s$$の場合である。従って、$n$に対して$M$や$s$が非常に大きいような場合は最初にデータ構造を作ってしまうべきである。一方、小さい$M$が保障されておりサンプリングの回数が大きくない場合は棄却サンプリングは有効である。二つ目は$M$の計算コストである。上記では$M$の計算コストは無視した。しかしながら、一般に二つの分布が与えられたとき、その間の$M$がどの程度に収まるかを計算するのは$O(n)$の時間が必要となる。もし$n$の時間をかけて$M$を計算したとすると、棄却サンプリングのコストは$O(n + Ms\log n)$となり常にセグメント木を用いるアプローチより低速である。従って、棄却サンプリングを有効に用いるには、目的分布の既知の性質から近いと確信が持てる提案分布が構成できなければならない。  内積の推定二つの長さ$n$のベクトル$x,y$があるとき、この二つの内積$x\cdot y$を求める計算は$O(n)$を要する。一方、$x,y$についていくつか追加の操作が許されているとき、ある誤差の範囲内でより高速に内積の値を推定することができる。サンプリング可能でノルムが既知のベクトル$x$と、一般のベクトル$y$がある場合、確率$p_i = \frac{x_i^2}{\left|x\right|^2}$の確率で値$z_i = \frac{y_i}{x_i}|x|^2$をとる乱数$z$を考える。この時、$z$の期待値は$E[z] = \sum_i p_i z_i = \sum_i x_i y_i = x \cdot y$となる。また、$z$の分散は$V[z]^2 \leq |x|^2|y|^2$となる。<jupyter_code>n = 100
x = np.random.rand(n)
x_square = x**2
x_square_norm = np.sum(x_square)
x_square_normalized = x_square / x_square_norm
y = np.random.rand(n)
inner_product = np.dot(x,y)
index_list = np.arange(n)

avg = 0
sampling_count_list = np.arange(0,1000000,1000)[1:]
estimation_list = []
for sample_index in range(max(sampling_count_list)):
    index = np.random.choice(index_list, p = x_square_normalized)
    z = y[index]/x[index] * x_square_norm
    avg += z
    if sample_index+1 in sampling_count_list:
        estimation_list.append(avg/(sample_index+1))
plt.plot([min(sampling_count_list), max(sampling_count_list)], [inner_product, inner_product], label = "inner_product")
plt.plot(sampling_count_list, estimation_list,label = "estimation")
plt.legend()
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14327217125382263, 'fraction_numerical': 0.05787461773700306, 'mean_word_length': 5.306894889103182, 'pattern_counts': {'":': 0, '<': 82, '<?xml version=': 0, '>': 68, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '491048', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 15767, 'n_tokens_neox': 14689, 'n_words': 1236, 'path': 'notebooks/7.3c_quantum_inspired_algorithm.ipynb', 'repo_name': 'Hiroya-W/quantum-native-dojo', 'size': '26160'}