23/12/30更新

誰もみてないと思いますが、コードにミスがあったのと、効率化の反映で記事を更新しました。
コードのミスの方は、並列化した時の結果が意図せずランダムに削除されるものなので、過去のは参考になりません。

Numbaを使って二次元のQDAを高速に解いた話

訳あって説明変数を2個選んできて、2次元にプロットして、綺麗にまとまっているかを調べる必要ができた。
ただ、その説明変数の種類が20万もあるため、その中から2個選ぶ選び方が
$${}_{200000} C_2 =19999900000$$
今回はこれを解くにあたって、QDA(二次判別)をNumbaで実装し、現実的な時間で計算できるようにしたので、自分のメモ兼ちょっとくらいはおんなじことをしようとしている人の役に立つかもと思い、書き残しておく。
なお環境は、M1macbookair,16Gbメモリ,
python=3.10.13 , numba=0.57.1 , numpy=1.24.4
なお、これを書いているのはpythonに触れて1年くらいの初学者で、大学で使うから頑張って書いたものです。
自分の環境では動いたので、まあそれでいいので他の人が動くかはわかりません。一例としてみてみてください。
また、質問があればぜひ、コメントにて。答えるかわからないし、答えられるのかもわかりません。

QDA(二次判別)

https://qiita.com/m1t0/items/06f2d07e626d1c4733fd
基本的には上のサイトに書いてあるので省く。
ベイズ確率周りを知っているとわかりやすいかも。
概要としては、

import pandas as pd
import numpy as np
from sklearn.preprocessing import MinMaxScaler
from sklearn.svm import SVC
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA
from sklearn.discriminant_analysis import QuadraticDiscriminantAnalysis as QDA
from sklearn.naive_bayes import BernoulliNB,CategoricalNB,ComplementNB,GaussianNB,MultinomialNB
from sklearn.datasets import make_moons,make_classification,make_gaussian_quantiles
#X,Y=make_moons(n_samples=1000,noise=0.1,random_state=4)
#X,Y=make_classification(random_state=12,n_features=2,n_redundant=0,n_informative=1,n_clusters_per_class=1,n_classes=2)
X,Y=make_gaussian_quantiles(random_state=0,n_features=2,n_classes=2)

scaler = MinMaxScaler()
scaler_X = scaler.fit_transform(X)
x1 = np.linspace(X[:,0].min(), X[:,0].max(), 1000)
x2 = np.linspace(X[:,1].min(), X[:,1].max(), 1000)
X1, X2 = np.meshgrid(x1, x2)
plot_X = np.c_[X1.reshape(-1),X2.reshape(-1)]
scaler_plot_X=scaler.fit_transform(plot_X)

#BernoulliNB(),CategoricalNB()は使い方が違う、BernoulliNB()は部分に分けてその中を1色にする?。
model_list=[SVC(kernel="poly"),LDA(),QDA(),ComplementNB(),GaussianNB(),MultinomialNB()]

for i in range(len(model_list)):
    plot_X = np.c_[X1.reshape(-1),X2.reshape(-1)]
    model=model_list[i]
    model.fit(scaler_X,Y)
    plot_y = model.predict(scaler_plot_X)
    plot_X=plot_X.T.reshape(2,1000,-1)
    print(str(model))
    plt.scatter(X[:,0][Y>0],X[:,1][Y>0],c="blue",s=50,alpha=0.8,label="formable")
    plt.scatter(X[:,0][Y<=0],X[:,1][Y<=0],c="red",s=50,alpha=0.8,label="formable")
    plt.contourf(plot_X[0], plot_X[1], plot_y.reshape(1000,1000), alpha=0.4,colors=["red","blue"],antialiased=True,levels=1)
    plt.show()

みたいな感じに分類してくれる。見ると、moonsは限界を感じるが、他はいい感じ。
上のコードでは他にもいろいろな分類器がまとめて比較できる。

Numba

https://qiita.com/gyu-don/items/9d223b007ca620e95abc
こっちは原理は全くよくわからないんですが、インタープリター言語なpythonでコンパイルしてかなり早くしてくれる。
並列化までpythonのコードのまんまできるため、かなーり有難い存在。
ただし、pythonの通常の文法+numpyの一部までしか対応しておらず、エラーがわかりにくく、謎のkernelクラッシュで苦しんだりも、、
もちろん入れていい引数のtypeもnumpyとpythonのみ。
pandasなどのものはnumpyに変換してから入れる。

コードの進み方

解くにあたって、問題点は、
1.かかる時間
2.メモリ
の二つだった。かかる時間は、sklearnの実装を並列化しても、自分が大学を卒業するくらいには時間がかかってしまうみたい。
メモリは工夫次第でなんとかなったが、実装し始めは200Gbくらいかかる計算だった。
最終的にはメモリは1Gb未満、実行時間は全体で時間未満になった(M1macにはファンがないため、扇風機を当てて。)

順序

データ作成
相関係数が高い説明変数を探す
上の相関係数の組でデータ作成を減らす
全ての組み合わせでQDA(二次判別)

1. データ作成

以下のようなデータを想定しています。
説明変数X : numpy.ndarray,shape=(データ数,説明変数の種類数)
目的変数Y : numpy.ndarray,shape=(データ数,)

2. 相関係数が高い説明変数を探す

ここがかなり時間がかかる。
データ数,説明変数の種類数が少ない時にはnp.corrcoefでまとめてやればいいが(圧倒的に最速)、いかんせん20万種もあるため、np.corrcoefの返り値は20万×20万、全部float64、メモリエラーです。
ので、一個一個計算して、返り値を「相関係数がある値(boder)を超えたかどうか」のみに限定してます。
これで、20万×20万(float64)→20万(bool)になるので、メモリ問題を解決しています。

完全に2次元に特化しているため、それ以外には別途拡張が必要です。

import numpy as np
from numba import njit,prange,set_num_threads
@njit(error_model="numpy",fastmath=True)
def sub_corrcoef(X_1,X_2):
    len=X_1.shape[0]
    mean_1=np.sum(X_1)/len
    mean_2=np.sum(X_2)/len
    minus_mean_1=X_1-mean_1
    minus_mean_2=X_2-mean_2
    upper=np.sum(minus_mean_1*minus_mean_2)
    down=np.sqrt(np.sum(minus_mean_1**2)*np.sum(minus_mean_2**2))
    data_corrcoef=np.abs(upper/down)
    return data_corrcoef

@njit(parallel=True,error_model="numpy",fastmath=True)
def sub_corrcoef_drop_njit_nonnan(len_n,index_1,return_list,X,border):
    index_2_list=np.arange(index_1+1,len_n)[~return_list[index_1+1:len_n+1]]
    n_index_2_list=index_2_list.shape[0]
    for i in prange(n_index_2_list):
        index_2=index_2_list[i]
        data_corrcoef=sub_corrcoef(X[:,index_1],X[:,index_2])
        if data_corrcoef>border:
            return_list[index_2]=True
    return return_list

def corrcoef_drop_njit(num_threads,X,border):
    if np.any(np.isinf(X)):
        X=np.where(np.isinf(X), np.nan, X)
    set_num_threads(num_threads)
    len_n=X.shape[1]
    return_list=np.full((len_n),False,dtype="bool")
    for index_1 in range(len_n-1):
        if return_list[index_1]:
            continue
        return_list=sub_corrcoef_drop_njit_nonnan(len_n,index_1,return_list,X,border)
    return np.where(return_list)[0]

num_threads = 8# int(input(f'Enter max number of threads to use: '))
border=0.995#相関係数の絶対値がこれ以上でdrop
drop_list=corrcoef_drop_njit(num_threads,use_X,border)
#drop_list.shape[0]

3. 上の相関係数の組でデータ作成を減らす

2で減らすリストを作ったので、適応するだけ。

use_X=np.delete(use_X, drop_list,axis=1)
index_to_columns=np.delete(index_to_columns,drop_list)

4. 全ての組み合わせでQDA(二次判別)

速さを求めているので、スコア(誤識別率)だけを出すようにしています。
loopする部分QDA_njit_loop関数と1組み合わせでスコアを出すQDA関数で分けていて、QDA関数は上でのQDAの説明で書いたurlの中に出てくる

\begin{gathered}-\frac{1}{2} x^{\prime }Ax+b^{\prime }_{1}x+c_{1}=0\\ A:=\Sigma^{-1}_{k} -\Sigma^{-1}_{l} ,b:=\Sigma^{-1}_{k} \mu_{k} -\Sigma^{-1}_{l} \mu_{l} ,c:=-\frac{1}{2} \mu^{\prime }_{k} \Sigma^{-1}_{k} \mu_{k} +\frac{1}{2} \mu^{\prime }_{l} \Sigma^{-1}_{l} \mu_{l} +\log \frac{\pi_{k} }{\pi_{l} } -\frac{1}{2} \log \frac{|\Sigma_{k} |}{|\Sigma_{l} |} \end{gathered}

の最後を変形して、

x^{T}\Sigma^{-1}_{k} x-2\left( \Sigma^{-1}_{k} \mu_{k} \right)^{T}  x+\mu^{T}_{k} \Sigma^{-1}_{k} \mu_{k} -\left( x^{T}\Sigma^{-1}_{l} x-2\left( \Sigma^{-1}_{l} \mu_{l} \right)^{T}  x+\mu^{T}_{l} \Sigma^{-1}_{l} \mu_{l} \right)  -2\log \frac{\pi_{k} }{\pi_{l} } +\log \frac{|\Sigma_{k} |}{|\Sigma_{l} |} =0

\left( A-B\right)^{T}  =A^{T}-B^{T}\  \  ,\  \  \left( AB\right)^{T}  =B^{T}A^{T}\  \  ,\  \  \left( A^{-1}\right)^{T}  =\left( A^{T}\right)^{-1}  \  \  ,\  \  \Sigma^{T}_{k} =\Sigma_{k}

より、

\left( x-\mu_{k} \right)^{T}  \Sigma^{-1}_{k} \left( x-\mu_{k} \right)  =x^{T}\Sigma^{-1}_{k} x-2\left( \Sigma^{-1}_{k} \mu_{k} \right)^{T}  x+\mu^{T}_{k} \Sigma^{-1}_{k} \mu_{k}

なので、

\left( x-\mu_{k} \right)^{T}  \Sigma^{-1}_{k} \left( x-\mu_{k} \right)  -\left( x-\mu_{l} \right)^{T}  \Sigma^{-1}_{l} \left( x-\mu_{l} \right)  -2\log \frac{\pi_{k} }{\pi_{l} } +\log \frac{|\Sigma_{k} |}{|\Sigma_{l} |} =0

をnumpyで書いたもの。
なお、例によって完全に2次元に特化しているため、それ以外には別途拡張が必要です。

import numpy as np
from numba import njit,prange,set_num_threads
@njit(fastmath=True)#error_model="numpy"
def jit_cov(X):
    n=X.shape[0]
    X_1=X[:,0]-np.sum(X[:,0])/n
    X_2=X[:,1]-np.sum(X[:,1])/n
    var_1=np.sum(X_1*X_1)/n
    var_2=np.sum(X_2*X_2)/n
    cross_cov=np.sum(X_1*X_2)/n
    return var_1,var_2,cross_cov

@njit(error_model="numpy",fastmath=True)
def QDA_nonnan(X,Y):
    class1_data=X[Y]
    class2_data=X[~Y]
    n_1=class1_data.shape[0]
    class1_var_1,class1_var_2,class1_cov=jit_cov(class1_data)
    class2_var_1,class2_var_2,class2_cov=jit_cov(class2_data)
    det_cov1=class1_var_1*class1_var_2-class1_cov**2
    det_cov2=class2_var_1*class2_var_2-class2_cov**2    

    target_x_1=(X-np.array([np.mean(class1_data[:,0]),np.mean(class1_data[:,1])]))
    target_x_2=(X-np.array([np.mean(class2_data[:,0]),np.mean(class2_data[:,1])]))
    value1=class1_var_2*target_x_1[:,0]**2+class1_var_1*target_x_1[:,1]**2-2*class1_cov*target_x_1[:,0]*target_x_1[:,1]
    value1=value1/det_cov1
    value2=class2_var_2*target_x_2[:,0]**2+class2_var_1*target_x_2[:,1]**2-2*class2_cov*target_x_2[:,0]*target_x_2[:,1]
    value2=value2/det_cov2
    value3=-np.log(det_cov1)+np.log(det_cov2)
    value=-value1+value2+value3
    return (value>0)==Y

@njit(fastmath=True)#error_model="numpy",fastmath=True)
def QDA_njit_sub_nonnan(thread_id,loop_array,how_save,len_n,n_data,X,Y):
    border,min_index=0,0
    return_list=np.zeros((how_save+1,3))
    for i in loop_array:
        index_1=int((2*len_n-1-np.sqrt(4*len_n**2-4*len_n+1-8*i))/2)
        index_2=int(i-(index_1*len_n-index_1*(index_1+1)/2)+index_1+1)
        pattern=np.array([index_1,index_2],dtype="int64")
        select_X=X[:,pattern]
        score=np.sum(QDA_nonnan(select_X,Y))/n_data
        if score>border:
            return_list[min_index]=np.array([index_1,index_2,score],dtype="float64")
            min_index=np.argmin(return_list[:,2])
            border=return_list[min_index,2]
    return return_list[np.argsort(return_list[:,2])[::-1]][:-1]

@njit(parallel=True)#,error_model="numpy",fastmath=True)
def QDA_njit_loop(num_threads,X,Y,how_save):
    Y=Y.astype("bool")
    set_num_threads(num_threads)
    len_n=X.shape[1]
    n_data=X.shape[0]
    repeat=int(len_n*(len_n-1)/2)
    return_list_2d=[np.zeros((int(how_save),3),dtype="float64") for _ in range(num_threads)]
    loop_array_list=np.array_split(np.arange(repeat), num_threads)
    for i in prange(num_threads):
        thread_id=i%num_threads
        #thread_id=get_thread_id()
        loop_array=loop_array_list[i]
        return_list_2d[thread_id]=QDA_njit_sub_nonnan(thread_id,loop_array,how_save,len_n,n_data,X,Y)
        
    return_list_2d_=np.zeros((num_threads,int(how_save),3),dtype="float64")
    for i in range(num_threads):
        return_list_2d_[i]=return_list_2d[i]
    return_list=return_list_2d_.reshape(-1,3)
    return_list=return_list[return_list[:,2]!=0]
    return_list=return_list[np.argsort(return_list[:,2])[::-1]]
    return return_list

num_threads = 8#int(input(f'Enter max number of threads to use: '))
how_save=200#int(input(f'Enter max number of threads to use: '))
print(f"data : {use_X.shape[1]}, loop : {use_X.shape[1]*(use_X.shape[1]-1)/2}")
data_ans=QDA_njit_loop(num_threads,use_X,use_Y,how_save)

測定

${}_{3000} C_2$で計測を行いました。
なお環境は、M1macbookair,8Gbメモリ,
python=3.10.13 , numba=0.57.1 , numpy=1.24.4
並列化を行うと基本全てのCPUを使うので、FanのないM1macbookairでは数秒でサーマルスロットリングを起こします。
そこで背中から普通の扇風機を当てて計測していますが、それでもサーマルスロットリングは避けられないので、不安定な結果になっています(ここにきてProにしなかったことを後悔,,,)。
また、メモリ問題はおおよそクリア(使用1Gb未満)になっているのと、測定法がいまいちわからないので、今回なしです。

corrcoef_cutを使用する場合border=0.995にしています。corrcoef_cutで軸種を減らすと一律3000→1487になっています。
QDAではhow_save=40とします。
一応sklearnの結果も載せていますが、これは

QDA_njit_loop()関数のscore=~~の部分をsklearnのQDA()に
sub_corrcoef(9関数をnp.abs(np.corrcoef())に

変えただけのものです。

sklearnはnumba非対応です。実際にはnumba中の並列化はmultiprocessingなどを使用しているそうで、それを使えば並列化は可能ですが、興味が湧かないので、測定はなしで。

		multi	numba	corrcoef_cut	corrcoef_cut_time	QDA_time	total_time
1	自作	T	T	T	333 ms	657 ms	990 ms
2	自作	F	T	T	645 ms	2.84 s	3.49 s
3	自作	F	F	T	24.57 s	80.65 s	105.22 s
4	自作	T	T	F	Nan	2.46 s	2.46 s
5	自作	F	T	F	Nan	11.23 s	11.23 s
6	自作	F	F	F	Nan	331.13 s	331.13 s
7	Sklean	F	F	T	45.94 s	388.10 s	434.04 s
8	Sklean	F	F	F	Nan	1555.75 s	1555.75 s

自作のnumbaなし、sklearnに負けてますね
車輪の再開発的なことになってるんで、まあそうかなとは思っていたんですが、2次元のみの想定で負けるのか、と思ってしまいます。

上の表から、1000000回のloopでの時間を計算してみました。

	データ数	multi	numba	QDA_time	time/1000000 loop
1	1487C2	T	T	657ms	595 ms
2	1487C2	F	T	2.84 s	2.57 s
3	1487C2	F	F	80.65 s	73.00 s
4	3000C2	T	T	2.46 s	547 ms
5	3000C2	F	T	11.23 s	2.50 s
6	3000C2	F	F	331.13 s	73.61 s

データ数1487C2でも3000C2でも1回あたりの実行時間は変わっていませんね。
multiにすると、だいたい4.3倍くらい高速になっているイメージでしょうか。
ここの値はCPUのコア数などによるでしょうし、サーマルスレッドリングにもよるんでしょう
本当にいいPCが欲しい、、、

最後に

実際の用途としては、${}_{200000} C_2$を相関係数が高いもの(border=0.995)で減らしたところ、
corrcoef_cut_timeは4m 58.66sで200000→49181となり、QDA_timeは12m 32.16sで実行できました。
time/1000000 loopは0.622 sなので、まあ想像通りの時間になっているかと思います。
M1macbookでこれなので、サーバーで行えば、1/3ぐらいにはなるのかな？
目的は達成なので、いろいろ実行していきたいと思います。

いろいろと、特にnumbaで躓いていたので、その辺は別ので書き残しておこうかと思います。