Fisher Discriminant Orthogonal Decompositionの高速化

Abstract

Fisher Discriminant Orthogonal Decomposition (FDOD)は，Fisher Discriminant Analysis (FDA)に正則化係数と直交分解を適用した判別分析法である．これにより，多変量データの判別分析において過学習を回避すること，グループ数以上の判別軸を求めることが可能になった．しかし，FDODは計算時間やメモリを多く消費する方法である．そこで，分析データを特異値分解して冗長データを取り除くことで，計算時間やメモリを節約する方法，Fast Fisher Discriminant Orthogonal Decomposition (FFDOD)を開発した．7054波数のデータからなる6種のセルロール系繊維の赤外級数スペクトル275個にFFDODを適用した場合，FDODと比較して約1/84の計算時間となった．特異値分解の計算を除外すると，約1/290と顕著な高速化が実現できた．また，FFDODとFDODの結果を比較することとで，計算精度においても同等であることが示された．

Translated Abstract

Fisher Discriminant Orthogonal Decomposition (FDOD) is a discriminant analysis method incorporating regularization coefficient and orthogonal decomposition into ordinary Fisher Discriminant Analysis (FDA). This method makes it possible to avoid overfitting in discriminant analyses of multivariate data and to obtain discriminant axes whose number is greater than that of groups. However, FDOD requires long calculation time and large memory. To solve these problems, a novel technique, Fast Fisher Discriminant Orthogonal Decomposition (FFDOD), has been developed. FFDOD saves calculation time and memory by singular value decomposition of the data to be analyzed to remove redundant data. When FFDOD was applied to 275 infrared spectra of 6 types of cellulosic fibers each of which consists of data at 7054 wavenumbers, the calculation time was reduced to 1/84 of that when using FDOD. If the time required for the singular value decomposition is not considered, a remarkable speedup to about 1/290 was realized. The calculation accuracy of FFDOD has been found equivalent with that of FDOD by comparing the results by FFDOD and FDOD.

1 はじめに

判別分析法は，複数の変数が組になったデータと，そのデータが属するグループインデックスを基にして，変数からグループを推測するモデルを構築する方法である．Fisher Discriminant Analysis (FDA) [1]は，パラメータの射影を求めたときに，グループ内の変動が最も小さくなり，グループ平均間の変動が最も大きくなるような判別軸を求める方法である．有名なアヤメのデータを使った解析 [1]のように変数が4つで観測数がこれを超えている場合，FDAは優れた結果を示す．

しかし，観測数よりも変数の数の方が多い場合，判別分析ではしばしば過学習を引き起こす原因となる．このような場合，正則化係数を適用した判別分析法，Regularized Discriminant Analysis (RDA)によって過学習を防ぐことができる [2, 3]．また，カーネルをFDAに適用したKernel Discriminant Analysis (KDA)でも，過学習を防ぐために正則化係数を用いている [4, 5]．

RDAやKDAでは正則化係数を1つ決定し，グループ数より1つ少ない数の判別軸が得られる．例えば2つのグループを判別する場合，1次元上での判別となる．しかし，より多くの次元の判別軸を用いたほうが，判別性能が高くなる可能性がある．また，3グループ以上の判別では2つ以上の判別軸が得られるが，1つの正則化係数のみなので，繊細な条件調整をすることができない．そこで，我々はRDAを改良し，グループ数以上の判別軸を求めることができ，判別軸ごとに正則化係数を設定できる，Fisher Discriminant Orthogonal Decomposition (FDOD)を開発した．赤外吸収スペクトルをFDODで解析することにより，セルロース系繊維の判別に非常に有効であることが示された [6, 7]，特許を取得した [8]．

FDAやRDAでは一般化固有値問題を1回解けばよいのに対し，FDODでは求めたい判別軸の数だけ一般化固有値問題を解かなくてはならない．正則化係数の最適化のためにはさらに多くの一般化固有値問題を解く必要がある．一般化固有値問題の求解は，通常の固有値問題と比べると計算時間を要する．正則化係数を精緻に決めていくためには，一般化固有値問題の迅速な求解が不可欠である．本研究では，高速化したFDOD，fast FDOD (FFDOD)を開発したので解説する．また，FDODとの性能比較を行う．

2 理論

2.1 Fisher Discriminant Analysis (FDA)

2.1.1 計算方法

m変数からなるデータを，観測数n収めた$n\times m$行列Xを考える．変数がm個あるので1組のデータはm次元空間中の1点となる．つまりXはm次元空間中のn個の点をまとめたものとなっている．m次元の適当な単位ベクトルvを考える．各観測データのv方向成分(スコア)をまとめたベクトルtは，Xのvへの射影であり，   

$\begin{array}{c}t=Xv\end{array}$

(1)のように求まる．

ここで各観測データはG個のグループのどれかに属しているとする．FDAではtのグループ間変動bが大きくなり，かつtのグループ内変動wが小さくなるように，すなわちb / wが最大となるvを求める．Xのグループ間変動行列をB，グループ内変動行列をWとすると，b / wが最大となるvは以下の一般化固有値問題，   

$\begin{array}{c}Bv=\lambda Wv\end{array}$

(2)を解くことで求めることができる．このようにして求まったvを判別軸 (discriminant axis: DA) と称することにする．Bのランクは最大で$G-1$となるので，3グループ以上の判別では固有ベクトルである判別軸が複数見つかる．そのため，3グループ以上での場合は多次元での判別となる．

2.1.2 多変量データにおけるFDAの問題点

Wのランクrank (W)は最大で   

$\begin{array}{c}rank\left(W\right)=\{\begin{array}{c}m\left(m<n-G\right)\\ n-G \left(m\ge n-G\right)\end{array}\end{array}$

(3)となる．スペクトルのように多くの波数や波長からなるデータでは，$m\ge n-G$となることが一般的である．一方，Wのサイズは$m\times m$となるため，$m\ge n-G$となる場合，Wは特異行列となる．この条件下で式(1)を解くと，データがグループごとに1点に集まるようになる．これは一見完全な判別ができたように見える．しかし，得られた判別軸を用いて所属グループ未知のデータの判別を行うと，どのグループからも離れた点となってしまい，判別不能となることが多い．これは，判別軸vをXに特化しすぎたための過学習の状態である．

2.2 Regularized Discriminant Analysis (RDA)

2.2.1 計算方法

FDAにおける過学習は，Wがランク落ちとなる条件であることに起因している．そこで，式(2)の右辺に正則化項を加えた次のような一般化固有値問題を考える．   

$\begin{array}{c}Bv=\lambda \left(W+\zeta I\right)v\end{array}$

(4)

ここで，$\zeta$は非負の正則化係数であり，Iは$m\times m$の単位行列である．正則化項を加えるためW + ζIは正則となり，各グループのデータは拡がりをもって分散することとなる．この拡がりの程度は正則化係数ζで調整することができる．$\zeta \cong 0$では式(2)と等しくなるため，グループごとに1点に集中する．一方，$\zeta \to \infty$の極限ではBの固有値問題を解くことと等しくなる．これはグループ間の変動を最大，すなわち各グループの平均(中心)同士がなるべく離れるようなvが求まることとなる．$\zeta$を調整して過学習を回避しつつ，未知データを適切に判別できるvを見つけることができるのがRDAである．

2.2.2 RDAの問題点

前述したようにBのランクは最大で$G-1$である．つまり2グループの場合は，判別軸vは1つしか見つからない．しかし，正則化項で判別性能を下げているため，1次元だけではなく2次元以上に判別情報が残っていると考えられる．そのため，2グループ判別であっても2次元以上の判別軸も使うことで，より良好な判別が行えると期待できる．

3グループ以上の判別では複数の判別軸が見つかるが，2グループの場合と同様に，より高次元の判別軸を使うことで，より良好な判別が可能であると期待できる．また，$\zeta$は1つしか設定することができないため，判別軸ごとに細かく設定することができない．さらに，複数求まる判別軸は互いに直交ではない$\left(v_i^{\text{T}}{v}_j\ne 0\right)$ため，各判別軸で捉える情報には重複があり無駄が生じることとなる．

2.3 Fisher Discriminant Orthogonal Decomposition (FDOD)

2.3.1 計算方法

FDODではRDAでの課題であった，$G-1$を超えて判別軸が求まらない，判別軸が直交しない，判別軸ごとの正則化係数が設定できない，という3点を解決するために以下のように計算を進める．

準備: 判別モデル校正用試料のスペクトルを収めた行列XをX₁と置く$\left(X_1=X\right)$．

以下の反復計算を$i=1$より始める．

1: X_iのグループ間変動行列B_i, グループ内変動行列W_iを求める．

2: 適当な第i正則化係数${\zeta }_i$を設定し，以下の一般化固有値問題を解く．   

$\begin{array}{c}{B}_i{v}_i={\lambda }_i\left(W_i+{\zeta }_{i}I\right)v_i\end{array}$

(5)解が複数ある場合は，最大固有値とそれに対応する固有ベクトルをそれぞれ，λ_i, v_iとする．

3: スコアt_iを求める．   

$\begin{array}{c}{t}_i=X_i{v}_i\end{array}$

(6)より高次の判別軸が必要ない場合はここで終了する．

4: v_iとt_iを使ってX_iをデフレーションしてX_i+1を求める．   

$\begin{array}{c}{X}_{i+1}=X_i-t_i{v}_i^{\text{T}}\end{array}$

(7)

5: $i\leftarrow i+1$として，手順1に戻る．

このようにして，反復計算を行い必要な数の判別軸を求めることができる．手順2で判別軸ごとに正則化係数を設定している．また手順7で直交分解を行っているため，判別軸は互いに直交するようになる．さらに，手順5のB_iは反復を繰り返すごとに更新されるため，X_iのランク以下にはならない．このようにして，前述の課題が解決されている．

2.3.2 FDODの問題点

FDA, RFDAそしてFDODでは一般化固有値問題を解くことがもっとも重要な過程である．一般化固有値問題は，$Av=\lambda v$の形をとる通常の固有値問題を解くよりも計算に時間を要する．FDAとRFDAでは一般化固有値問題を1回解くだけで終了するが，FDODは必要な判別軸の数だけ解かなくてはならない．さらに，判別軸ごとに正則化係数を調整するために，例えば各${\lambda }_i$を微小量変化させ，評価用試料の判別性能の変化を求めるならば， 確率密度比ると非常に多くの時間を要することになる．

2.4 高速FDOD (fast FDOD: FFDOD)

2.4.1 FFDODの理論

スペクトルを取り扱うときのように$m>n$となる場合，Xのランクは$\text{rank}\left(X\right)\le n$となる．そのためXの各列は互いに線形独立ではなく線形従属となり，Xには重複した情報が含まれていることになる．つまり，BやWは重複した情報から生成した無駄に大きな行列となってしまっている．この冗長性を解消し，WやBのサイズを抑えることができれば，計算速度を改善させることができる．

FFDODではX₁の冗長なデータを取り除くために，まず特異値分解を行う．   

$\begin{array}{c}{X}_1=LS{R}^{\text{T}}\end{array}$

(8)

ここで，X₁のランクをrとすれば，Lはr個の左特異ベクトルからなる$n\times r$行列，Sはr個の特異値を対角にもつ$r\times r$の特異値行列，Rはr個の右特異ベクトルからなる$m\times r$行列となる．特異ベクトルは互いに直交しているので，L^TL = I, R^TR = Iである．ここで，   

$\begin{array}{c}{X}_{s1}:=LS\end{array}$

(9)と定義する．これにより，X₁は   

$\begin{array}{c}{X}_1=X_{\text{s}1}{R}^{\text{T}}\end{array}$

(10)のように表される．これを一般化してiのときについても同様に考え   

$\begin{array}{c}{X}_i=X_{\text{s}i}{R}^{\text{T}}\end{array}$

(11)とする．さらに，判別軸に関して   

$\begin{array}{c}{v}_i=R{v}_{\text{s}i}\end{array}$

(12)となるv_si (|v_si| = 1)を考える．$t_{\text{s}i}:=X_{\text{s}i}{v}_{\text{s}i}$と定義すると，   

$\begin{array}{c}{t}_i=X_i{v}_i=X_{\text{s}i}{R}^{\text{T}}R{v}_{\text{s}i}=X_{\text{s}i}{v}_{\text{s}i}=t_{\text{s}i}\end{array}$

  

$\begin{array}{c}\therefore t_i=t_{\text{s}i}\end{array}$

(13)との関係が得られる．

グループ間変動行列B_iやグループ内変動行列W_iをX_iから求める行列をそれぞれD_B, D_Wとして，式(4)を書き換えると   

$\begin{array}{c}{X}_i^{\text{T}}{D}_{\text{B}}{X}_i{v}_i={\lambda }_i\left(X_i^{\text{T}}{D}_{\text{W}}{X}_i+{\zeta }_i{I}_m\right)v_i\end{array}$

(14)となる．I_mは$m\times m$の単位行列であることを明示するために，mを添字に用いた．式(14)に式(11), (12)を代入し，両辺に左からR^Tをかけると   

$\begin{array}{c}{X}_i^{\text{T}}{D}_{\text{B}}{X}_i{v}_i={\lambda }_i\left(X_i^{\text{T}}{D}_{\text{W}}{X}_i+{\zeta }_i{I}_m\right)v_i\\ R^{\text{T}}R{X}_{\text{s}i}^{\text{T}}{D}_{\text{B}}{X}_{\text{s}i}{v}_{\text{s}i}={\lambda }_i{R}^{\text{T}}\left(R{X}_{\text{s}i}^{\text{T}}{D}_{\text{W}}{X}_{\text{s}i}{R}^{\text{T}}R+{\zeta }_{i}R\right)v_{\text{s}i}\\ X_{\text{s}i}^{\text{T}}{D}_{\text{B}}{X}_{\text{s}i}{v}_{\text{s}i}={\lambda }_i\left(R^{\text{T}}R{X}_{\text{s}i}^{\text{T}}{D}_{\text{W}}{X}_{\text{s}i}+{\zeta }_1{R}^{\text{T}}R\right)v_{\text{s}i}\end{array}$

  

$\begin{array}{c}\therefore X_{\text{s}i}^{\text{T}}{D}_{\text{B}}{X}_{\text{s}i}{v}_{\text{s}i}={\lambda }_i\left(X_{\text{s}i}^{\text{T}}{D}_{\text{W}}{X}_{\text{s}i}+{\zeta }_i{I}_r\right)v_{\text{s}i}\end{array}$

(15)が得られる．I_rは$r\times r$の単位行列である．$X_{\text{s}i}^{\text{T}}{D}_{\text{B}}{X}_{\text{s}i}:=B_{\text{s}i}$, $X_{\text{s}i}^{\text{T}}{D}_{\text{W}}{X}_{\text{s}i}:=W_{\text{s}i}$と定義すると，式(15)は   

$\begin{array}{c}{B}_{\text{s}i}{v}_{\text{s}i}={\lambda }_i\left(W_{\text{s}i}+{\zeta }_i{I}_r\right)v_{\text{s}i}\end{array}$

(16)と書き換えられる．これは，式(11)が成立している条件下では，式(5)を解く代わりに式(16)を解き，得られたv_siを式(12)に適用してv_iを求めることができることを示している．B_iやW_iが$m\times m$であるのに対し，B_siやW_siは$r\times r$となる．$m>n$の条件下ではB_siやW_siの行列サイズが抑えられるため，式(16)を解くための計算時間は式(5)に比べて節約できることになる．

次にデフレーションについて考える．X_iのデフレーションは式(7)で表される．式(7)に式(11), (12)を代入すると，   

$\begin{array}{c}{X}_{i+1}=X_i-t_i{v}_i^{\text{T}}\\ =X_{\text{s}i}{R}^{\text{T}}-t_i{v}_{\text{s}i}^{\text{T}}{R}^{\text{T}}\\ =\left(X_{\text{s}i}-t_i{v}_{\text{s}i}^{\text{T}}\right)R^{\text{T}}\end{array}$

  

$\begin{array}{c}\therefore X_{i+1} =\left(X_{\text{s}i}-t_{\text{s}i}{v}_{\text{s}i}^{\text{T}}\right)R^{\text{T}}\end{array}$

(17)となる．一方，$i+1$のときについて考えると，式(11)は   

$\begin{array}{c}{X}_{i+1}=X_{\text{s}i+1}{R}^{\text{T}}\end{array}$

(18)である．式(17)と式(18)を対比すると，   

$\begin{array}{c}{X}_{\text{s}i+1}=X_{\text{s}i}-t_{\text{s}i}{v}_{\text{s}i}^{\text{T}}\end{array}$

(19)との関係となる．これは，X_siのデフレーションになっている．つまりX_iのデフレーションとX_siのデフレーションは等価であることが分かる．

X₁ではなくX_s1を初期行列としてFDODを適用していくのがFFDODとなる．

2.4.2 FFDODの手順

準備: 判別モデル校正用試料のスペクトルを収めた行列XをX₁と置く$\left(X_1=X\right)$．

1: X₁を特異値分解する．   

$\begin{array}{c}{X}_1=LS{R}^{\text{T}}\end{array}$

(20)

2: X_s1を得る．   

$\begin{array}{c}{X}_{\text{s}1}=LS\end{array}$

(21)以下の反復計算を$i=1$より始める．

3: X_siのグループ間変動行列B_si, グループ内変動行列W_siを求める．

4: 適当な第i正則化係数${\zeta }_i$を設定し，以下の一般化固有値問題を解く．   

$\begin{array}{c}{B}_{\text{s}i}{v}_{\text{s}i}={\lambda }_i\left(W_{\text{s}i}+{\zeta }_{i}I\right)v_{\text{s}i}\end{array}$

(22)が複数ある場合は，最大固有値とそれに対応する固有ベクトルをそれぞれ，λ_i, v_siとする．

5: スコアt_iを求める．   

$\begin{array}{c}{t}_i=t_{\text{s}i}=X_{\text{s}i}{v}_{\text{s}i}\end{array}$

(23)

6: 判別軸v_iを求める．   

$\begin{array}{c}{v}_i=R{v}_{\text{s}i}\end{array}$

(24)より高次の判別軸が必要ない場合はここで終了する．

7: v_siとt_siを使ってX_siをデフレーションしてX_si+1を求める．   

$\begin{array}{c}{X}_{\text{s}i+1}=X_{\text{s}i}-t_{\text{s}i}{v}_{\text{s}i}^{\text{T}}\end{array}$

(25)

8: $i\leftarrow i+1$として，手順1に戻る．

このようにFDODと同様，反復計算を行い必要な数の判別軸を求めていく．FDODとの主な違いは，手順1で特異値分解を行っていることと，手順6で判別軸を求めていることである．いったん求めた判別軸は以後の計算には使われないので，手順6をとばして必要な反復が終わったのちにまとめて求めても構わない．この場合，$V=\left(\begin{array}{ccc}{v}_1& v_2& \cdots \end{array}\right)$, $V_{\text{s}}=\left(\begin{array}{ccc}{v}_{\text{s}1}& v_{\text{s}2}& \cdots \end{array}\right)$とするなら，   

$\begin{array}{c}V=R{V}_{\text{s}}\end{array}$

(26)のように求めることができる．Vを判別軸行列と呼ぶことにする．各判別軸に対応するスコアをまとめたスコア行列Tも以下のように定義する．   

$\begin{array}{c}T:=\left(\begin{array}{ccc}{t}_1& t_2& \cdots \end{array}\right)\end{array}$

(27)

3 スペクトル測定

FFDODによる計算の高速化を検証するために，セルロース繊維の赤外吸収スペクトルを用いた．供試試料として，cotton 40点，linen 39点，ramie 42点，rayon 66点，cupra 47点，lyocell 41点の6種，計275点の繊維布地を用意した．これらの布地の赤外吸収スペクトルは，ATRアタッチメント(JASCO ATRPRO450S)を取り付けた分光光度計(JASCO FT/IR 4700)により測定した．ATRプリズムは ZnSeを用いた．1回反射法によりATR測定を行った．検出器にはTGS検出器を用いた．記録波数領域は4000 − 600 cm⁻¹，波数分解能は2 cm⁻¹，積算回数は16回とした．アポダイゼーション関数はコサイン関数を用いた． ATR 法では波長によるもぐりこみ深さの差異を小さくするためATR 補正(強度のみ)を行った．

4 FDODとFFDODの計算速度比較

得られたATR赤外吸収スペクトル275個(n = 275)に対して，FDODとFFDODによる6グループ判別(G = 6)を行った．本稿では判別アルゴリズムの計算速度比較を目的とするため，スペクトルに対して前処理および正則化係数の最適化は行わなかった．第6判別軸まで求めることとし，正則化係数は各軸において1次から6次判別軸への順に10^−2.5, 10⁻³, 10⁻³, 10⁻³, 10⁻³, 10^−4.5とした．これらはスコアプロットを行った際に，目視によって良好な判別が可能であることより決定した．

解析に供するスペクトルの波数領域は，最低波数600 cm⁻¹は固定し，最高波数を4000 cm⁻¹から50 cm⁻¹程度減らしながら，732 cm⁻¹まで変化させた．4000 − 600 cm⁻¹では波数の数は7054個であり，つまりm = 7054となる．732 − 600 cm⁻¹の場合はm = 276である．

計算にはMATLAB Onlineを用いた．計算時間の計測は変動があることを考慮し，各波数領域条件において10回行った．

5 結果と考察

測定された6種類の繊維種のスペクトルをFigure 1に示した．各繊維種の平均スペクトルが示してある．また，スペクトルの重なりを避けるために各スペクトルに定数を足してずらしてある．どれもセルロース系繊維であるため，スペクトルが似通っていることが分かる．そのため，視覚的にスペクトルから繊維種を判断するのは難しいと言える．

Figure 1.

 Infrared absorption spectra of six kind cellulose fabrics. Spectra are shown with shifted baselines to avoid overlapping each other.

解析波数領域を4000 − 600 cm⁻¹としたとき，FODOおよびFFDODによって得られた，第1−2判別軸2次元スコアプロット(以後"1−2プロット"と記す．他も同様)および5−6プロットをFigure 2に示した．ここでの目的は手法の比較であるため，FDODによる結果を"×"，FFDODによる結果を"+"で表してあり，繊維種ごとの表記は行わなかった．すべてのプロットが"*"のように見えている．これは，2つの手法で各プロットの座標が一致しており，2種類のマーカーが重なったことによる．つまり，2つの手法での判別結果が一致していることを示している．

Figure 2.

 Score plots for (a) 1-2 discriminant axis plane and (b) 5-6 discriminant axis plane.

Each marker looks like [*] because × and + overlap.

1−2プロットでは，天然繊維のうち麻に分類されるlinenとramieが重なってしまい判別不能となった．これは麻同士のスペクトルが似通っているため，第1判別軸と第2判別軸ではlinenとramieのスペクトルの差異を分離できなかったためである．これ以外に関しては，互いに異なる領域にプロットされていた．5−6プロットでは，天然繊維のcottonとramieが重なってしまったものの，linenが異なる領域にプロットされた．再生繊維である，rayon, cupra, lyocellは近い範囲にブロットされ判別は難しい結果となった．1−2プロットでは分離されなかったlinenが，5−6プロットで分離されており，複数の次元でのスコアを用いることの有用性がわかる結果となった．

また，第6判別軸は通常のFDAやRFDAでは得られない次元の軸であり，FDODやFFDODがより高次の判別軸を見つけられることを示している．

前述した1−2プロットや5−6プロットだけでなく，ほか組み合わせのスコアプロットもFDODとFFDODのスコアが一致していた．とくに5−6プロットは反復計算の最後に得られる2つの判別軸である．第5や第6判別軸は，反復初期に求まる第1や第2判別軸に比べ，より多くの計算を経て求まる．そのため，誤差が蓄積しているのでFDODとFFDODの結果に差が生じることが考えられる．

しかしFigure 2を見る限りFDODとFFDODのスコアの間には差が見られなかった．

FDODとFFDODの計算時間の比較のために，計算時間についていくつか定義をしておく．FDOD，FFDODともに，計算を始めてから終了(第6判別軸・DA6を求める)までを総計算時間t_totalとする．特異値分解に要した時間をt_svdとする．FDODには特異値分解はないので，t_svdは含まれない．FDODルーチンの所要時間をt_fdodとする．FDODの場合，FDODルーチンそのものであるので，t_fdod = t_totalである．一般化固有値問題の求解時間をt_gepとする．第6判別軸まで求めているため，t_gepは一般化固有値問題の求解6回の合計となる．

Figure 3-5に変数の数mと計算速度の関係をグラフに示した．Figure 3ではFDODとFFDODのt_totalを比較した．変数の数mが500程度よりも大きいときはFFDODのほうがFDODよりも計算時間が短いことが分かった．mの増加につれて，FDODでは急激にt_totalが増加しているのに対し，FFDODでは変化が緩やかであった．そのため，mが増加するほど計算速度の差が顕著になった．FFDODではm = 1500 − 3000の範囲で計算時間の上昇が見られるが，この現象には再現性があるが，Linux版MATLAB 2019aで計算したときには，このような傾向が別の範囲で見られた．よってこの現象はデータ依存ではなく，MATLAB Onlineの計算システムの特性によるものと考えられる．

Figure 3.

 Comparison of total calculation time (t_total) for FDOD and FFDOD.

Figure 4はFDODおける変数の数の二乗m²とt_totalとt_gepの関係を示した．このようにすると，m²と計算時間の間に比例関係があることが分かる．一般化固有値問題はB_iとW_iを用いるので，これらの行列サイズ$m\times m$によって計算時間が規定されているものと思われる．t_gepはt_totalの5割以上を占めていた．しかし，t_gep以外にも無視できないほどの時間が割かれている．t_totalとt_gepがm²と比例関係にあることから，このその他の部分もm²と比例関係にあることが分かる．この多くはB_iやW_iの準備に費やされている．B_iやW_iはある行列とその転置行列の積から求まる$m\times m$行列である．そのため，生成にはm²に比例した時間を要することになる．

Figure 4.

 Relation between the square of the number of variables and calculation time for FDOD.

Figure 5にはFFDODの計算時間の変化をまとめた．このグラフでは横軸はmである．全体的な傾向では，mの増加に伴いt_totalとt_svdは直線的に上昇していた．前述したm = 1500 − 3000程度の範囲においてt_totalが長くなるのは，この範囲におけるt_svdの影響であることが分かる．全体的な傾向は，特異値分解は元のデータXに対して行われるので，mの変化に影響を受けて直線的に上昇すると説明できる．一方，t_fdodやt_gepは変化に乏しくほぼ一定であった．これはFFDODでは，mに関係なく行列サイズが一定($n\times n$)B_siやW_siに対してFDODを行うので，t_fdodやt_gepの計算時間が一定となるためである．そのため，mがある程度より小さくなると，特異値分解に無駄な時間を費やしてしまい，FDODよりも計算時間が長くなってしまう．正則化係数の検証など計算を多く繰り返すような場合は，数回計算してどちらがより効率的かをあらかじめ求めておくことが必要である．

Figure 5.

 Relation between the number of variables and calculation time for FFDOD.

FDODとFFDODをより詳細に比較するため，最高波数が4000, 2300, 1450, 1025, 813, 732 cm⁻¹の場合について，変数の数，行列サイズ，計算時間，Vの規格直交性についてFDODとFFDODの比較をTable 1, 2にまとめた．各計算時間は10回の試行の平均値と簡略化形式で標準不確かさを示した．また，FDODとFFDODで得られたVやTの同等性も評価した．

Table 1.  Comparison of FDOD and FFDOD calculations part 1.

range (cm−1)	4000 − 600		2300 − 600		1450 − 600
the number of variables	7054		3528		1765
method	FDOD	FFDOD	FDOD	FFDOD	FDOD	FFDOD
	matrix size of B or W
	7054 × 7054	275 × 275	3528 × 3528	275 × 275	1765 × 1765	275 × 275
	calculation time*1(s)
ttotal*2	9.1 (1)	0.108 (4)	2.15 (1)	0.074 (2)	0.498 (7)	0.090 (3)
tsvd*3	–	0.077 (3)	–	0.046 (1)	–	0.060 (2)
tfdod*4	9.1 (1)	0.031 (1)	2.15 (1)	0.028 (1)	0.498 (7)	0.030 (1)
tgep*5	5.8 (1)	0.0168 (7)	1.34 (1)	0.0156 (9)	0.336 (5)	0.0162 (6)
	orthonormality of V
$\max \left\{\left|\left|v_i\right|-1\right|\right\}$ (×10−16)	2.22	2.22	2.22	1.11	1.11	1.11
$\max \left\{\left|v_i^{\text{T}}{v}_j\right|\right\}$ (×10−13)	8.67	1.26	3.50	0.55	3.54	0.59
	maximum difference between each elements of matrices obtained by FDOD and FFDOD*6
T(×10−12)	7.52		2.02		2.16
V(×10−14)	4.84		6.69		4.00

*1: Average of 10 times calculations. The values in parentheses are uncertainty. *2: Total calculation time. *3: Calculation time spent on singular value decomposition routine. *4: Calculation time spent on FDOD routine. *5: Calculation time spent solving general eigenvalue problem. *6: Since the signs corresponding values may be different, the difference was calculated for their absolute values.

Table 2.  Comparison of FDOD and FFDOD calculations part 2.

range (cm−1)	1025 − 600		813 − 600		732 − 600
the number of variables	884		444		276
method	FDOD	FFDOD	FDOD	FFDOD	FDOD	FFDOD
	matrix size of B or W
	884×884	275×275	444×444	275×275	276×276	275×275
	calculation time (s)
ttotal	0.142 (2)	0.065 (2)	0.048 (1)	0.053 (2)	0.030 (1)	0.051 (2)
tsvd	–	0.033 (1)	–	0.0231 (8)	–	0.020 (1)
tfdod	0.142 (2)	0.032 (2)	0.048 (1)	0.030 (2)	0.030 (1)	0.030 (1)
tgep	0.088 (1)	0.018 (1)	0.0267 (8)	0.0167 (9)	0.0167 (9)	0.0172 (7)
	orthonormality of V
$\max \left\{\left|\left|v_i\right|-1\right|\right\}$ (×10−16)	3.33	2.22	2.22	3.33	2.22	3.33
$\max \left\{\left|v_i^{\text{T}}{v}_j\right|\right\}$ (×10−13)	5.53	0.29	4.05	0.69	5.08	1.08
	maximum difference between each elements of matrices obtained by FDOD and FFDOD
T(×10−12)	2.07		0.62		1.05
V(×10−14)	0.66		0.66		0.65

前述したように選択された波数範囲が狭くなる，つまり変数の数が少なくなるにつれてt_totalが短くなる傾向となっている．ただし，波数範囲1450 − 600 cm^–1(m = 1765)の場合，t_total = 0.090 sとなり，2300– 600 cm^–1(m = 3528)の0.074 sよりも長くなっていた(Table 1)．これは，変数の数が前述の計算時間が長くなる1500 − 3000の範囲に入っているため，特異値分解に0.060 sを要していることによる．

m = 7054変数の場合(Table 1)t_totalがFDODでは9.1 sだったのに対し，FFDODでは0.108 sと約1/84と大幅な時間節約となっていた．さらにt_fdodで比較すると，FDODは9.1 sであるが，FFDODでは0.031 sとさらに短くなり，約1/290となった．FDODで最適な判別軸を探索するためには，1つのX₁に対して様々な正則化係数を適用した計算を繰り返す必要がある．この過程には多くの繰り返しが必要なので，変数の数が多くなると多大な時間が必要となる．FFDODを使えば大幅に時間を短縮できるようになる．さらにFFDODでは元のデータ行列X₁に変化がない限り，X_s1も変化しない．つまり，一度特異値分解でX_s1を求めておけば，正則化係数を変化させた計算はX_s1に対して行えばよい．特異値分解を繰り返す必要がないので，さらなる時間の短縮を期待することができる．

FFDODではmに関係なくB_iやW_iの行列サイズが$275\times 275$と一定となるため，t_fdod, t_gepはそれぞれ，0.03 s, 0.017 s程度と一定であった．そのため，mが小さくなると

FDODに対する優位性が低くなっていく．波数範囲が813 − 600 cm^–1の場合，t_totalがFDODでは0.048 s，FFDODでは0.053 sとなっていた．不確かさを考慮するとこの2つの結果の間に有意差があるとは言えないため，FFDODの有意性は認められないこととなった．ただし，FFDODのt_fdodは0.030 sとなり，FDODの0.048 sと比較して37.5%の時間短縮となった．繰り返し計算を行う場合には，FFDODを選択することで効果を得ることができるようになる．

波数範囲が732 − 600 cm^–1(m = 276)場合，t_fdod, t_gepに関しては2つの手法で同等の結果となったが，t_totalではFFDODのほうが劣る結果となった．この場合FFDODを選択する理由はないと言える．

FFDODによって時間が短縮されたとしても，Vの規格直交性が悪かったり，FDODとの結果に差異があったりした場合，FFDODは利用できない．Vの規格直交性を検証するため，$\left|\left|v_i\right|-1\right|$と$\left|v_i^{\text{T}}{v}_j\right| \left(i\ne j\right)$を求めた．$\left|\left|v_i\right|-1\right|$が大きいと規格性が悪く，$\left|v_i^{\text{T}}{v}_j\right|$が大きいと直交性が悪いことが分かる．Table 1, 2には変数の数の条件ごとに，これらの値の最大値$\max \left\{\left|\left|v_i\right|-1\right|\right\}$と$\max \left\{\left|v_i^{\text{T}}{v}_j\right|\right\}$を示した．$\max \left\{\left|\left|v_i\right|-1\right|\right\}$はどの条件においても10⁻¹⁶程度となっており，十分小さく良好な規格性を持っていることが示された．$\max \left\{\left|v_i^T{v}_j\right|\right\}$はどの条件においても${10}^{-13}$程度以下となっていた．この値は$\max \left\{\left|\left|v_i\right|-1\right|\right\}$に比べると大きいが，十分に小さく直交性が保たれていると言える．また，$\max \left\{\left|v_i^T{v}_j\right|\right\}$はFFDODのほうがFDODよりも小さくなっていた．FFDODを用いることで，直交性が保たれるだけでなく，より高い直交性が得られる可能性がある．

たとえ規格直交性が保たれたとしても，得られるTやVがFDODとFFDODの結果の間で異なっていた場合にはFFDODを利用できない．Figure 2に示したように，7054変数の場合の1−2プロットや5−6プロットがFDODとFFDODで一致していることがわかった．視覚的だけでなく，数値的にもこのことを確認するため，スコア行列Tおよび判別軸行列Vについて，FDODとFFDODで対応する

要素同士の差の絶対値を求めた．この差の絶対値のうち，最大値をTable 1, 2に提示した．どの変数条件においても，Tでは10⁻¹²程度，Vでは10⁻¹⁴程度と十分に小さい値になっており，FFDODによってFDODと等しい判別軸とスコアが得られることが数値的にも示された．

6 おわりに

FDODはFDAを基にして，判別軸ごとの正則化係数を適用することで優れた判別性能を持つ方法である．サポートベクターマシン(SVM)やディープラーニング(DL)などの非線形判別法と比較してもいくつかの利点がある．とくに，SVMは多グループ判別に向いていないが，FDAやFDODは多グループ判別が前提となっていることがある．また，SVMやDLは非線形であるがゆえに，判別の要因の解読は難解，もしくは不可能である．これに対しFDAやFDODの判別軸には，どの変量が判別に有効であるかの情報が含まれている．スペクトルのように化学的な情報が含まれるデータの場合，判別軸を調べることで判別に有用な化学的情報を得ることも可能である．

FDODの真価を発揮するためには，適切な正則化係数を見つける必要がある．しかし一般化固有値問題の求解に計算時間の多くが費やされるので，正則化係数をきめ細かく決めるためにFDODを繰り返すと多大な時間がかかることが欠点であった．FDODを行う前に，特異値分解を行うFFDODを用いることで，行列サイズや計算時間を大幅に節約することができた．また，FDODとFFDODOで同じ判別軸やスコアプロットが得られることも示された．

FDODにおける一般化固有値の求解のCPUへの負担は大きく，筆者のPCでも変数が多い条件での連続したFDODの実行は，CPUが高温となりシステム暴走の危険があったため続行を断念した(そのためクラウド上で実行されるMATLAB Onlineを利用した)．これに対し，FFDODでは連続計算であっても，適切な温度が保たれていた．つまり，FFDODの採用は，システムの安定運用という面からも有用であることが分かる．

さらに，近年，小型で安価な近赤外などの分光計が入手可能になってきている．これらの分光計に，シングルボードコンピュータ(SBC)などを組み合わせることで，作業現場や屋外に手軽に持ち運び可能な小型分光システムを構築できる．SBCはCPU性能やメモリ容量などがPCに比べ低く抑えられている．FFDODであれば，このような小型のシステムでであっても採用することが可能である．

FFDODの適用で計算速度を向上させることができるようになったが，最適な正則化係数を得るためには判別モデルのバリデーション法や判別の評価指標も今後の課題である．既報 [7]では，校正用試料と検証用試料を別に用意する外部バリデーションを用いた．また，各グループは正規分布を形成していると仮定し，確率密度を判別指標として用いた．試料が少ない場合などは，クロスバリデーションが有効であると考えられる．評価指標としては，マハラノビス距離や事後確率なども用いることができる．最適な正則化係数の探索法を確立することで，FDODの真価が発揮されると期待される．

謝辞

赤外吸収スペクトルは，東京農工大学農学府修了生の齋藤健悟氏に提供していただいた．ここに感謝の意を表する．

参考文献

• [1]   R. A. Fisher, Ann. Hum. Genet., 7, 179 (1936). doi:10.1111/j.1469-1809.1936.tb02137.x
• [2]   J. H. Friedman, J. Am. Stat. Assoc., 84, 165 (1989). doi:10.1080/01621459.1989.10478752
• [3]   Y. Guo, T. Hastie, R. Tibshirani, Biostatistics, 8, 86 (2007). doi:10.1093/biostatistics/kxj035 PMID:16603682
• [4]   S. Mika, G. Rätsch, Klaus-Robert Müller, Adv. Neural Inf. Process. Syst., 13, 591 (2001).
• [5]   Z. Zhang, G. Dai, C. Xu, M. I. Jordan, J. Mach. Learn. Res., 11, 2199 (2010).
• [6]   K. Sakamoto, M. Kanno, N. Yoshimura, M. Takayanagi, J. Comput. Chem. Jpn., 17, 225 (2018). doi:10.2477/jccj.2018-0058
• [7]   K. Saito, T. Yamagata, M. Kanno, N. Yoshimura, M. Takayanagi, Spectrochim. Acta A Mol. Biomol. Spectrosc., 257, 119772 (2021). doi:10.1016/j.saa.2021.119772 PMID:33887512
• [8]   M. Takayanagi, N. Yoshimura, K. Saito, T. Ando, M. Kanno, International patent publication number, WO/2018/ 074002