行列・テンソル分解によるヘテロバイオデータ統合解析の数理

―第6回 系列データ―

Abstract

生命科学分野で取得されるデータ集合は、雑多（ヘテロ）な構造になり、ヘテロなデータ構造を扱える理論的な枠組みがもとめられている。本連載では、汎用的なヘテロバイオデータの解析手法である行列・テンソル分解を紹介していく。第6回では、データに順序がある系列データで利用される行列・テンソル分解手法を紹介する。

系列データとは

本稿では、これまでの連載［1, 2, 3, 4, 5］では扱ってこなかった時系列データや空間データといった、データに順序がある系列データに対する行列・テンソル分解の手法を紹介する。生命科学における時系列データとしては、例えば時系列計測したオミックスデータ［6, 7］、神経活動データ［8］、脳波データ［9, 10］などが挙げられる。また、最近では、フローサイトメトリーやシングルセルオミックスデータを一度次元圧縮した後に、最小全域木などのアルゴリズムで全体をなぞるような分化経路を推定し、それを擬似的な時系列とみなすケースも増えている［11, 12］。また、空間データとしては、ゲノム座標（locus）に紐づくデータ［13, 14］や、空間オミックス［15］、細胞形状の画像のピクセルデータ［16］などがある。また、計測年月や計測エリアごとに生物種の個体数をカウントした生態学データ［17］や、時系列計測した脳波データ［18］などは、時間にも空間にも紐づいたデータであるため時空間（Spatio-Temporal）データと呼ばれる。

時間や空間といった変数は、順序を考慮せずにただのカテゴリカルな変数と捉えた場合、1つのモードとして扱うことで素直に行列やテンソルとして表現でき、それにより各種行列・テンソル分解が利用できる（図1）。実際に、そのようにして、時間や空間に関する情報も、他のモードと区別なく扱った事例は多い［6, 7, 8, 9, 10, 13, 14, 15, 16, 17, 18］。一方で、系列データを解析するにあたっては、系列データならではの特徴や問題点があり、それにより通常の行列・テンソル分解がうまく機能しない場合もある。また、系列データが持つデータの並び順や、データ同士の近さといった特徴をうまく組み込んだ、系列データに特化したモデルの方がより良い性能を示すといったこともしばしば起きる。

図1：様々な系列データ（時系列データ、空間データ、時空間データ）

（左上）変数ごとに何かの値を経時計測、または異なる空間で計測した場合、変数×時間、または変数×空間の2階テンソル（行列）になる。（右上）変数ごとに、時間と空間の情報が共に紐づいていた場合、変数×時間×空間の3階テンソルになる。（下）さらに空間が2次元（例：x, y座標）である場合、変数×時間×x座標×y座標の4階テンソルになる。

系列データの問題点の一例としては、系列データは時間や空間座標の変化に対して、観測値が少しずつ連続的に変化するため、隣接または近傍にあるデータは似た値になりやすいという「自己相関性」が挙げられる。これは系列データ間に依存関係があることを示唆しており、i.i.d.^†の仮定から逸脱する性質であることから、i.i.d.を仮定したモデルの性能が落ちる原因となり得る。例えば、自己相関が強い系列データにPrincipal Component Analysis（PCA、第1回［1］）を適用した場合、自己相関が無いデータと比べて、同程度に分散が説明できるまでにより多くの主成分を必要とし、それによりノイズまで拾ってしまうため、説明能力が下がることが報告されている［19］。また、得られたスコアも時間や空間に対して滑らかに変化する保証が無く、値が飛び飛びに変化したパターンとして得られた場合、研究対象の現象によっては、非現実的な解析結果となってしまう。

系列データの特徴や問題点に関しては他にも様々なものが議論されているため、以降はややオムニバス形式とはなってしまうが、系列データに特化した行列・テンソル分解のモデリングの工夫を幾つかのトピックに分けて紹介する。トピック間にも互いに関連するものが多いため、適宜参照することにする。また、今回に関しては、まだ生命科学分野での利用事例が無いモデルが多いことにも注意してほしい。ただし、シングルセルオミックス［11, 12］や空間オミックス［15］など、昨今のハイスループットな実験技術が発展している現状を踏まえると、本稿で紹介した手法が今後利用されるであろうことを著者は期待している。

†：independent and identically distributed（i.i.d.、独立同分布）：これまでの連載［1, 2, 3, 4, 5］で紹介した行列・テンソル分解手法の多くは、確率モデルとして解釈した場合、全てのデータが互いに独立に、かつ同じ分布から生成されるというi.i.d.を仮定している。i.i.d.ではデータ1つ1つに対して同じ目的関数を適用するため、全データ分の目的関数の式の形が簡略され、最適化しやすいメリットがあるため、多くのモデルが暗にこれを仮定している［20］。

時間・空間正則化

行列・テンソル分解における系列データに対する工夫としては、まず時間に関するモードを含んだ行列やテンソルを構築しておき、これを分解する際に、時間に関する因子行列に対しては、隣りあう時間での値の変化が滑らかになるように制約を与える正則化（平滑化）がある（図2）。例えば、空間トランスクリプトームデータにPCAを適用し、スポットのクラスタリングを行う際に、スポット間の物理的な距離に関する正則化を加えることで、より2次元空間上で近傍のスポット同士が同じクラスタになりやすくしたSpatialPCAという手法がある［21］。また、アルツハイマー病の研究を行うAlzheimer’s Disease Neuroimaging Initiativeのチームでは、疾患リスクのある一塩基多型（Single Nucleotide Polymorphism；SNP）と、脳画像データを経時的に計測したコホート研究を行っており、この2種類のデータをCanonical Correlation Analysis（CCA、第2回［2］）で低次元に射影する際に、時間方向では値が急激に変化せず、なるべく滑らかに変化するように正則化項を加えたモデルが提案されている［22］。正則化のための関数としては、時間軸上で隣り合うパラメーターのL2正則化以外にも、グラフ正則化（第2回［2］）が用いられる場合もある［23］。同様に、空間座標（図2）や時空間（時間＋空間）（図2）的に隣接する要素同士に関しても平滑化を行った事例がある［24, 25, 26］。

図2：時間・空間正則化

（上）変数×時間（または空間）の行列を分解することで、変数に関する因子行列Uと、時間（または空間）に関する因子行列Vが得られる。時間（または空間）に関する正則化では、このVの隣り合う列同士の値が急激に変化しないように拘束条件を設定する。（下）変数×時間×空間の3階テンソルでも同様に、変数に関する因子行列Uと、時間に関する因子行列Vと、空間に関する因子行列Wが得られ、VやWの隣り合う列同士の値が急激に変化しないように正則化を行う。

ARモデル

時刻tにおけるテンソルデータX_t（ここではテンソルのオーダーは問わない）が、それよりP個手前のX_t－P（Pは1以上の整数）までのデータX_t－1, X_t－2, . . ., X_t－Pのみに依存するとした仮定に基づくモデルを一般にP次マルコフ過程という。例えば、第5回［5］で紹介したべき乗法やPageRankなどのグラフ上のランダムウォークモデルは1次マルコフなモデルである。また、X_tを目的変数、X_t－1, X_t－2, . . ., X_t－Pを説明変数に設定した回帰モデルはP次自己回帰（Auto-regressive；AR）モデルという。例えば、上記の時空間上で隣り合うデータ同士の影響を考慮した時間・空間正則化は1次のARモデルの一種とも言える。時系列データでは、過去のデータからの影響の蓄積により現在時刻tでの観測値が決定されることが多いため、i.i.d.のモデルよりもこのようなマルコフ性を考慮したモデルの方がより性能が向上しやすい。

ここで、時刻tにおけるテンソルデータX_tを、長さNのベクトルとした多次元時系列データとして考えてみる。ARモデルの入出力をスカラ値からベクトルに拡張したVector Autoregression Model（VAR）では、ある時刻tにおけるN個の変数が、一つ手前の時刻t－1におけるN個の変数に依存すると考える。一般的な重回帰モデルではN個の説明変数から1つの目的変数を予測することになるため、推定する係数は長さNのベクトルであるのに対して、VARではN個の説明変数からN個の目的変数を予測することになるため、N×Nの係数行列を推定しなければならず、Nの数が増えるに従い計算コストが問題となる。

一般にテンソル回帰［27, 28］というアプローチでは、このようにサイズが大きいパラメーターオブジェクトに対して低ランク性を仮定することで、計算コストの問題を解消する。例えば、VARの係数行列に特異値分解（Singular Value Decomposition；SVD、第1回［1］）など行列分解を適用したReduced-rank VAR［29, 30］が挙げられる（図3）。さらにより過去からの影響を捉えるべくP次のVARへモデルを拡張すると、P個のN×Nの係数行列が得られるため、これらを奥行き方向にまとめることでN×N×Pの3階テンソルとなる。Multilinear Low-rank VAR［31］という手法では、このテンソルに低ランク性を仮定し、VARの最適化にテンソル分解の一種であるHigher-order Singular Value Decomposition（HOSVD、第3回［3］）を組み合わせることにより、予測精度を落とさずに学習するパラメーターの数を大幅に削減することに成功している（図3）。

図3：VARと行列・テンソル分解

（上）Reduced-rank VARでは、VARのN×Nの係数行列Aに低ランク性を仮定し、行列分解を適用することにより因子行列UとVの行列積として近似する。（下）Multilinear Low-rank VARでは、P次VARモデルにおける複数のN×Nの係数行列をまとめて3階テンソルとし、低ランク性を仮定してテンソル分解を適用することにより、因子行列U, V, Wとコアテンソル$g$の積として近似する。なお、X_t−1:t−Pは時間tに対して1からP手前のデータを連結した行列である。

Reduced-rank VARと似た考え方ではあるが、状態空間モデルと呼ばれるアプローチでは、係数行列ではなくデータに対して低ランク性を仮定し、N個の変数を一度k（≪N）個のスコア（状態変数）に次元圧縮し、そのスコアにVARモデルを適用する（図4）。このスコアの計算に様々な行列分解が利用されており、例えば、PCA、Independent Component Analysis（ICA、第1回［1］）、Partial Least Squares（PLS、第2回［2］）、CCAを利用したものとして、DiPCA［32］、DWPCA［33］、AR-ICA［34］、DiPLS［35］、DWPLS［36］、DiCCA［37］、DWCCA［38］などがある。なお、テンソル分解においては、このスコアに相当するものはコアテンソルであり、コアテンソルのP次ARモデルとしてMultilinear Orthogonal Autoregressive（MOAR）がある［39］（図4）。

図4：状態空間モデルとしての行列・テンソル分解

（上）状態空間モデルでは行列分解で次元圧縮されたデータのスコアのうち、系列上で隣り合うもの同士でARモデルを仮定する。（下）MOARではテンソル分解で次元圧縮されたテンソルデータのスコア（コアテンソル）のうち、系列上で隣り合うもの同士でARモデルを仮定する。

状態空間モデルの離散値版、すなわち、データが連続値ではなく離散値で、かつデータが生成される際の状態にも離散値の状態変数を設定し、状態ごとにデータの生成確率が変わるとした、隠れマルコフモデル（Hidden Markov Model；HMM）がある［40, 41, 42, 43］。HMMのパラメーター推定法としてはEMアルゴリズムであるBaum-Welchアルゴリズムが利用されるが、面白いことに近年HMMは対称行列版のNon-negative Matrix Factorization（NMF、第1回［1］）であるTriNMFとの等価性が示されている［44］（ただし、実データにおける解の一致性、計算コストなど詳細な違いは不明）。HMMは、バイオインフォマティクス分野においては、遺伝子領域予測、膜貫通領域予測、スプライスサイト予測、配列アライメント、タンパク質ドメイン予測、RNA2次構造予測、など塩基配列やアミノ酸配列に関連したあらゆる解析に広く活用されている。

MAモデル

時系列データの統計的性質（平均、分散、自己相関など）が時間とともに大きく変動する場合、そのデータは「非定常」だと言われる。例えば、トレンドや周期性のあるデータは、時間が経つにつれて平均が変わり、データの変動の仕方も変化するため、このようなデータは非定常とみなされる。非定常なデータをそのまま扱うと、モデルの予測や分析結果に不安定さが生じやすく、データの傾向を掴むことが難しくなる。そのため、時系列データを扱う際には、非定常性を取り除き、データを「定常」に近づけることが望ましい。

時系列データに含まれる非定常性を緩和するためのシンプルな手法として、移動窓を設定し、窓内で平均値をとっていく移動平均（Moving Average；MA［45, 46］）がある（図5）。MAを行うことで、短期的な変動やホワイトノイズが平滑化され、より長期的なパターンが抽出されやすくなる。また類似した手法として過去のデータを辿るほど指数関数的に減衰する重みをかけた指数重み付き移動平均（Exponentially Weighted Moving Average；EWMA）がある。事前に前処理としてこれらの変換を行なっておくことをMA-フィルタ、EWMA-フィルタという。MAやEWMAを利用した行列・テンソル分解として、MA-PCA［47, 48］、MAICA［49］、MA-CP［50］、EWMA-PCA［48］、EWMA-CP［50］などがある。また、上述のARとMAを組み合わせたDynamic PCA［47］、上述のVARとMAを組み合わせて状態空間モデルの状態変数をICAで推定するVARMA-ICAといったモデルも提案されている［51］。

図5：移動窓（MA）と行列・テンソル分解

1次元の時系列データ上で、固定長の移動窓を走査させ、切り出した窓内で平均値をとる。この操作を多変数分行うことで、＋の矢印の方向にデータのサイズが増え、変数×時間の行列が得られる。また、複数個体で同様の行列が得られる場合、テンソルのオーダーが上がり、変数×時間×個体の3階テンソルが得られる。

時間遅れ共分散行列

2つの1次元時系列データa, bを考えてみる。bはaに影響を受けて、値を変化させるものとする。もしこれらの観測値が変化していく間隔よりもデータのサンプリングの間隔の方が小さいと、aの情報がbに伝達され、bの値が変化するまでには遅延（ラグ）がデータに現れる。このラグをどのように扱うかで、様々な行列・テンソル分解のモデルが提案されている。

時間1からTまで計測した多変数時系列データXにおいて、時間1からT－τまでのX₀と1＋τからTまでのX_τとの共分散行列$X_0{X}_{\tau }^T$ を時間遅れ（または時間差）共分散行列という（図6）。τ＝0の時は通常の共分散行列（第1回［1］）と一致する。この行列を利用した次元圧縮のことを遅延埋め込みといい、τ次の自己回帰性を考慮した次元圧縮の一種といえる。時間遅れ共分散行列に対して、各種行列分解を適用することにより、τステップ手前の別のデータとの関係性や、その関係性を踏まえたデータ間の類似関係の把握が可能となる。例えば、時間遅れ共分散ベースのPCAとしては、Dynamic PCA、Moving Dynamic PCA（MDPCA）、Generalized Dynamic PCA（GDPCA）、PCA for time series（TS-PCA）、Generalized PCA for moderately non-stationary vector time series（GTS-PCA）など、実に様々な手法が提案されている［47, 52］。

図6：時間遅れ共分散行列と行列・テンソル分解

時間1からTまで計測した多変数時系列データXにおいて、時間1からT－τまでのX₀と1＋τからTまでのX_τとの共分散行列X₀X_τ^Tを時間遅れ（または時間差）共分散行列という。

ICAベースの手法としては、共分散行列と時間遅れ共分散行列とで共に成り立つ固有ベクトルを一般化固有値問題で解いたTime-structure based Independent Component Analysis（tICA）があり、タンパク質構造に関するMDシミュレーションにおける特徴的な運動パターンの抽出に利用されている［53］。Dynamic Component Analysis（DCA）では、事前にPCAで次元圧縮した時系列データに対して、時間遅れ共分散行列の代わりに時間遅れを考慮した相関行列を作成し、固有値分解（Eigen Value Decomposition；EVD、第1回［1］）を適用しており、MDシミュレーション［54］や経時計測したfMRI［55］のデータ解析に適用されている。また、τ＝1での時間遅れ共分散行列に対してSVDを実行したものは、動的モード分解（Dynamic Mode Decomposition；DMD［56］）といい、近年流体力学分野で発展した行列分解手法である。線虫の体の動きを経時計測した研究［57］では、τの値を複数設定し、τごとに切り出された多変数時系列データを連結し、1つにまとめた行列に対してSVD及びICAで次元圧縮することで、複数時刻間の依存関係を同時に扱っている。

ハンケル行列

ある1次元の系列データに対して移動窓（Sliding window）を適用し、固定長の窓幅で窓をスライドさせて得られたデータセットを束ねて1つの行列としたものを、ハンケル行列（Hankel Matrix）、履歴行列、または軌道行列（Trajectory matrix）などという［27, 58, 59］（図7）。ハンケル行列に対して適用される行列分解の手法としては、まず1次元系列データに含まれる変化点を検出するための手法である、特異スペクトル変換（Singular Spectrum Analysis；SSA［60］）が挙げられる。SSAではまず時系列データを過去の学習データと、新しく生成されたテストデータとに分割する。学習データに対しては、SVDを適用することで、ハンケル行列に頻出するパターンを学習する。テストデータに対しては、学習されたパターンで近似できる度合いを示すスコアによって、これまでにないパターンが含まれていた場合に検知できるようにしている。SSAの生命科学における利用事例としては、時系列計測や空間座標を伴う計測を行なったDNAマイクロアレイ［61］や、著者の松本が線虫の体の動きを解析した事例［62］などがある。

図7：ハンケル行列と行列・テンソル分解

1次元の系列データ上で、固定長の移動窓を走査させ、切り出した窓を縦に連結させる。この操作により得られるフレーム数×窓幅の行列をハンケル行列と呼ぶ。この操作を多変数分行うことで、＋の矢印の方向にデータのサイズが増え、テンソルのオーダーが上がり、フレーム数×窓幅×変数の3階テンソルが得られる。

なお、ハンケル行列は多変数に拡張することで、さらにモードが1つ増え3階テンソルのハンケルテンソルとなる（図7）。これに対してテンソル分解を行ったものとして、Tensor-based SSA［63, 64］がある。また、画像のように2次元座標上にあるデータに関しては、固定サイズの2次元の窓を用意し、これを縦横方向にスキャンする。この操作により、最終的には時系列データよりもさらに1つオーダーが大きい、4階のハンケルテンソルが得られる。テンソル分解の一種であるTucker分解を行うことで欠損値補完に活用したものとしてMissing Slice Recovery（MSR）がある［65］。また動画（＝時系列画像）データに対して、MSRのハンケルテンソルの分解と、MOARの系列コアテンソルの自己回帰性を組み合わせ、さらにAuto Regressive Integrated Moving Average（ARIMA）モデル（AR、MAに加え、隣接データ間の差分をとったもの）へと拡張させたものとしてBlock Hankel Tensor ARIMA（BHT-ARIMA）がある［66］。

シフト、畳み込み

1次元時系列データ同士に特化した距離尺度として、Shape-based distance（SBD［67］）と、動的時間伸長法（Dynamic Time Warping；DTW［68］）がある。ここでは、データ間に依存関係のネットワークがあるものとし、情報伝達にはラグがあるものとする。SBDでは、片方のデータを時間軸に対して並行移動（シフト）させ、ユークリッド距離が最小、もしくはコサイン類似度が最大となるようなシフトを行うことでラグの影響を補正した距離を求める（図8）。一方、DTWは両データでの総当たりの時間ペアでの距離を算出した後に、動的計画法によるペアワイズアライメントよって最も全体のコストが少ない経路を両データの距離とする（図8）。SBDはデータの時間的なずれが一定と考えられる場合に適しており、一方DTWは時間のずれが動的に変化するデータを比較する際に適している。

図8：SBD/DTWと行列・テンソル分解

（左上）SBDでは2つの時系列変数の片方を並行移動させ、最も波形が一致する様にシフト値を探索する。（右上）DTWは両データでの総当たりの時間ペアでの距離を算出した後に、動的計画法によるペアワイズアライメントよって最も全体のコストが少ない経路を求める。（下）各々の距離尺度により、変数×変数の距離行列が取得される。また、複数個体で同様の距離行列が得られる場合、テンソルのオーダーが上がり、変数×変数×個体の3階テンソルが得られる。

DTWはシングルセルオミックスでデファクトスタンダードとなっているツールであるSeuratのVersion 2以降で、複数のバッチを統合する際に、CCAで次元圧縮後の細胞をバッチ間で対応付けさせるのに採用された［69］。著者の露崎が開発した、線虫の神経細胞間の神経活動的な時系列相互作用を検出する手法WormTensor［70］では、SBDまたはDTWによって細胞×細胞の距離行列を求め、これを多個体分揃えることで、細胞×細胞×個体の3階テンソルを算出し、さらにテンソル分解の一種であるMC-MI-HOOI［71］を適用し、細胞のクラスタリングを行うことで、多個体で再現性のある細胞の機能モジュールを同定している。Williamsらは、マウスの神経細胞で類似した解析を行っているが、DTWと同様に時間の対応付けを行うWarping行列を導入し、それ自体もデータから最適化できるようにモデリングしている［72］。さらに同研究チームは、同実験データを複数試行分まとめたデータとして扱えるようにするために、Warping行列とテンソル分解を組み合わせたモデルを提案している［73］。

なお、データにラグが含まれていたとしても、変わらずに同じ解を求めることができるようなモデルは、シフト不変的（Shift-invariant）と呼ばれる［74］。SBDのように、シフト不変的な行列分解では、PCA、ICA、NMFのような通常の行列分解と同様に、変数×時間の行列データを因子行列の積として分解しつつも、抽出する時系列パターンが変数ごとにシフトすることを想定し、シフト変数τも同時に推定する。シフト不変的な行列・テンソル分解としては、SOBI［75］、coroICA［76］、Shifted ICA［77］、Shifted NMF［78］、Shifted CP［79］などが挙げられる。

シフト不変的なモデルと類似した手法として、行列分解を行う際に、抽出する時系列パターンが、異なる時間で繰り返し出現することを想定し、1つずつ時間をずらした時系列パターンの線形和として、時系列データを近似する「畳み込み」型の行列・テンソル分解［80］がある。音源分離のデータ解析においては、室内で音が発された時に、壁や天井などによる反射で、その音がその後も繰り返し聞こえる「残響」を捉えるために利用されている。畳み込み型の行列・テンソル分解としては、Convolutive ICA［81］、Convolutive NMF（またはNon-negative Matrix Factor Deconvolution；NMFD［80, 82］）、Convolutive NTF（またはConvolutive CP［80, 83］）、またConvolutiveかつARなICAとしてCICAAR［84］などが挙げられる。

オンラインアルゴリズム

データが経時的に発生するストリーム型の時系列データを考える。今時刻tでのデータX_tが得られた時に、計測開始時から前回計測時までのデータX₁, X₂, ..., X_t－1を、X_tとマージした上でモデルにフィッティングさせる方針では、新たにデータが発生した時に、その都度過去のデータとマージし直さねばならない。また、データをマージし続けることで、データサイズが増え、いずれメモリに載り切らない状況も考えられる。そのため、X_t－1までで得られたパラメーター（例：ローディング行列U_t－1）に、X_tで得られるパラメーターの変化量を継ぎ足していくことで、パラメーターの値を少しづつ更新し、U_tを得るような定式化が望ましい（図9）。前者の方針はオフラインアルゴリズム、後者はオンラインアルゴリズムという。行列・テンソル分解においては、Online PCA［85］、Online PLS［85］、Dynamic NMF［86］、など様々なアルゴリズムがオンライン化されている。テンソル分解では、Dynamic Tensor Analysis（DTA［87］）、SALS［88］などのモデルが提案されている。また、パラメーターの最適化手法の1つである確率的勾配法では、データの一部分だけを切り出してパラメーター推定に利用するが、時系列データを順番にモデルに入力することで、オンライン学習も可能である。

図9：オンライン型行列・テンソル分解

（上）全てのデータを一度にメモリに載せて計算するオフラインの行列分解とは異なり、オンライン型の行列分解では、系列に従い逐次発生するベクトルデータχ_tにおいて、χ_tと1つ手前のパラメーターU_t−1とを利用して、パラメーターを更新しU_tを得る。（下）行列分解の場合と同様に、オンライン型のテンソル分解では、系列に従い逐次発生するテンソルデータχ_tにおいて、χ_tと1つ手前のパラメーター集合（因子行列：U_t−1、V_t−1、W_t−1、コアテンソル：$g$_t−1）を利用して、パラメーター集合を更新する。

オミックス解析では基本的に既に計測し終えたデータに対して、一度に解析を実施するため、経時的にデータが生成される状況にはなかなかならないが、オンライン型の実装は、全データを一度にメモリ上に展開せずに、メモリに載る分のデータだけを最初に最適化し、その他のデータは逐次的にモデルにフィッティングさせていくことから計算コストが低く、時系列データではないものも含むあらゆる大規模データに対して適用可能である。例えば、大規模なscRNA-SeqデータにOnline PCAを適用した事例としては、著者の露崎が過去に勾配法やべき乗法ベースなど14ものPCAアルゴリズムを実行したもの［89］（JuliaパッケージOnlinePCA.jlとして誰でも利用可能 https://github.com/rikenbit/OnlinePCA.jl）や、UK Biobank上の大規模なSNPデータにオンライン実装のPCA（Randomized SVD、IRAM、第5回［5］）を適用したPCAone［90］などが挙げられる。また、時系列計測したシングルセルマルチオミックスデータに対して、変分オートエンコーダーを用いて逐次的に次元圧縮を適用するMultimodal Integration with Continual Learning（MIRACLE）という手法もある［91］。

時間―周波数解析

1次元の時系列データの信号源に周期性がある場合に、事前に時系列データ行列やハンケル行列に離散フーリエ変換（Discrete Fourier transform；DFT）を適用することで、時間・フレームごとの周波数成分が得られる（図10）。さらに多変数で同様の操作を行うことで、フレーム×周波数×変数の3階テンソルが得られるため、それに対してテンソル分解を適用した事例がある［92, 93］。また、DFTは入力の1次元時系列データをフレーム×周波数のフーリエ基底行列と、長さが周波数分だけある周波数成分ベクトルの積として表現できるため、DFT自体を行列分解、多変数でのDFTをテンソル分解と捉えることもできる［27］。

図10：時間―周波数解析と行列・テンソル分解

1次元の時系列データに対して離散フーリエ変換を適用することにより、時間ごとに周波数成分が取り出され、周波数×時間の行列が得られる。この操作をさらに複数の変数で行うことにより、テンソルのオーダーが上がり、周波数×時間×変数の3階テンソルが得られる。

音声データ解析分野では、複数のマイクを用意して、音声を時系列サンプリングし、DFTを適用することで時間×周波数×マイクの3階テンソルの形でデータが得られる。このデータから元の信号源を分離するブラインド音源分離（Blind Source Separation；BSS）分野では、どれか周波数を固定した上でICAを適用するFrequency-domain ICA（FDICA［94, 95］）や、どれかマイクを固定した上でNMFを行うアプローチ［94, 95］が利用されていたが、前者は異なる周波数間で推定した信号源の並び順を対応づけることが困難であること、後者は複数マイクへの拡張の仕方が不明、抽出したパターンを信号源と結びつけるのが難しいといった課題があった。

現在では、前者の問題に対しては独立ベクトル分析（Independent Vector Analysis；IVA［94, 95］）という手法が提案されており、信号源を別々の1次元ベクトルではなく、複数の信号源を格納した多次元ベクトルとして推定できるように信号源に多変量確率分布を設定することで、異なる周波数間でも並び順が同じになるような工夫がなされている。また、後者の問題に対しては、上記のNMFを複数マイクへ拡張したMultichannel NMF（MNMF［94, 95］）が提案されており、マイク間の位相のずれを考慮しながら、マイク間の依存関係を考慮できるように、因子行列をエルミート行列に拡張した定式化がなされている。また、IVAとMNMFのコンセプトを融合した独立低ランク行列分析（Independent Low-Rank Matrix Analysis；ILRMA［94, 95］）というモデルも提案されている。これらの手法は、さらにARやMAと組み合わせることで、AR-MNMF［96］、ARMA-MNMF［97］、AR-ILRMA［98］といった手法へと発展している。

馬蹄効果

最後に、系列データ特有の問題として知られる馬蹄効果というものを紹介する（図11）。ここでは多変量時系列データを考え、直鎖グラフのように、隣り合うデータ間には依存関係があり、それらは似た値を取りやすいものとする。データに潜在的にこのような構造がある場合、これらをEVDで次元圧縮して得られたスコア（固有ベクトル）は三角関数のように周期性があるものとして推定される。そのため、このスコアを利用して2次元にプロットした際には、データがU字に曲がった状態で配置されることになる。この周期的なパターンは、その形から馬の蹄鉄にちなんで馬蹄（形）効果（Horseshoe effect［99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113］）、またはガットマン効果（Guttman effect［101, 113］）などと呼ばれる。また、Correspondence Analysis（CA、第4回［4］）では、比較的この馬蹄効果が出づらく、U字がやや引き延ばされ、V字になることから、アーチ効果（Arch effect［99, 113］）と呼ばれたり、神経科学や異常検知などの分野における、より長期に計測した系列データにおいては、何度も振動する様なスコアが取り出されることから、正弦波効果［108］、ファントム振動（Phantom oscillations［106］）と呼ばれたりもする。

図11：馬蹄効果

（上）直鎖グラフの様に系列上で隣り合うデータが似通う場合、この行列のグラム行列は、テプリッツ構造を示す。この様な行列に対して固有値分解を行った場合、次元圧縮図上で曲がったデータの並びを示すが、これはデータに周期性があるかには無関係の見かけの周期性であることに注意が必要である。（下）なぜ曲がったプロットになるのかについては、系列上で一定以上離れたデータでは距離が変化しなくなるという飽和特性で説明される。

このようなパターンは、しかしながら、データに実際に周期性があるかどうかには無関係に現れる「見かけの周期性」であり、誤った解釈に繋がる危険性がある。例えば、このデータが「多能性幹細胞→前駆細胞→分化細胞」といった分化経路を示すシングルセルデータだとして、この周期性を根拠として「分化細胞が再び多能性を獲得したのでは？」と勘違いする状況を想像していただければ、これは瑣末な問題ではないことを理解していただけるかと思う（その他「病人が健康になった」「老人が若返った」などで置き換えも可）。そのため、解析者は次元圧縮のスコアの分布だけでなく、変数1つ1つが実際に周期性を持つのかを確認するなど、注意深く結果を解釈する必要がある。

馬蹄効果はこれまでに、PCA［99, 100, 101, 102, 103, 104, 105, 106, 113, 114］やMulti-Dimensional Scaling（MDS、第4回［4］）［101, 102］、Isomap［101］、Local Linear Embedding［101］、カーネルPCA［101］、Laplacian Eigenmap（LE、第4回［4］）［107］、拡散マップ（第4回［4］）［107］、SSA［106, 108］、時系列k-meansクラスタリング［108］など、様々なEVDベースの行列分解手法で観測されている。また、シングルセルオミックス［99, 100, 106, 107, 109, 110］、マイクロバイオーム［102］、神経科学［103］、集団遺伝学［104, 105］、生態学［111］、骨形態計測学［114］など様々な生命科学分野で観測され問題視されている。対数変換などの前処理［99］によってこの馬蹄効果が緩和される場合もあれば、生態学分野では、仮に馬蹄効果が出てしまった場合であっても、強制的に直線上に引き伸ばして可視化するDetrended CA（DCA［111］）という手法も存在する。

馬蹄効果の最大の謎は、本来は直鎖状に配置されるべきデータが、なぜ次元圧縮をするとU字状に曲がってしまうのかである。未だ不明点も多いが、現在までにわかっていることをここでは記す。まず、直感的な説明としては、系列データにおいては、系列上で離れるにつれてそれ以上は距離の値が変化しないという飽和特性（Saturation Property）［102］が挙げられる（図11）。この飽和特性により、あるデータ（例：x₀、図11）を始点として考えた場合、飽和後のデータ点（例：x₈, x₉、図11）は、始点を中心とした同心円状に配置されやすくなり、この特性がx₀以外のどのデータ点にも当てはまると考えると、全体的に辻褄が合う構造として、曲がった構造が出てくるというものである。

より理論的な説明としては、MDSや各種カーネル次元圧縮を統一的に定義した後に、行列摂動理論を適用することにより、最初のスコアは単調増加、2番目のスコアは一峰性となることを示したものがある［101］。ただし、一般に馬蹄効果は3番目以降のスコアにも現れうるが、それに関する説明は特に無い。また、直鎖グラフ様データにおけるグラム行列（第1回［1］）は、図11における、右上から左下へ斜め方向に同じ数値が並ぶ、テプリッツ行列^††（Toeplitz Matrix）と呼ばれる行列となるが、このテプリッツ行列の固有ベクトルは、フーリエ基底すなわち三角関数として近似できることを示したものがある［103］。ただし、この近似が成り立つためには、データ間の類似関係は直線的ではなく、始点と終点が繋がったリング状である必要があり、実際のデータにはそぐわないといった批判もある［106］。

馬蹄効果については、未だ不明点やデータ解析における課題も多い。例えば、馬蹄効果をよりクリアに説明できる理論的な解析や、馬蹄効果とCAのアーチ効果の違い、ファントム振動との関係性、見かけの周期性データの中に真に周期性があるデータが含まれていた場合の両者の識別方法、馬蹄効果にロバストな次元圧縮手法、馬蹄効果を受けた次元圧縮の後のクラスタリングやt-SNE/UMAPへの影響、またDCAのように後から馬蹄効果をデータから取り除く手法…といった研究は、今後も行われると思われる。

なお、図11の帯行列［100, 102, 113］以外にも、ガウス反応型帯行列［113］、リッチネス行列［112］、リニア行列［112］、シフト行列［103］、ハンケル行列［108］など多数の直鎖グラフ様データで馬蹄効果が確認されており、著者の露崎が実際にこれらからU字のプロットが現れることを確認した実験はGoogle Colaboratoryで公開されている［115］。また、LEを軸として、グラフラプラシアン（第4回［4］）、ハンケル行列、系列データ行列における馬蹄効果を考察したものとして、本誌同号に同時掲載されている著者の松本の稿［116］も併せて参照いただきたい。

††：ノイズを含む実データでは、厳密に同じ数値が並ぶわけではないため、「テプリッツ構造を持つ」として区別する場合もある。

おわりに

今回は、系列データのデータ解析をする際には、隣接または近傍にあるデータ同士の自己相関性や、値の連続性、非定常性、ラグ、逐次的なデータの発生、（見かけの）周期性など、系列データならではの特徴や問題点があり、それにより様々なモデリングの工夫が行われていることを紹介した。データ構造としては、これまで通り行列・テンソル分解のモデルをそのまま適用できるものの、こういった工夫を加えることで、より良い予測結果が得られることは、昨今のTransformer型大規模言語モデルの発展からも伺える。

次回は、行・列方向でともに和が1となるように行列を正規化する行列バランス化問題や、近年注目されている、一つの確率分布から別の確率分布へ最小のコストで輸送する方法を求める「最適輸送問題」と、行列・テンソル分解との関係について議論する。

References

• [1]  Tsuyuzaki, K. (2021/4/23). Mathematics for Heterogeneous Biological Data Fusion Analysis with Matrix-Tensor Factorization - Part I. Matrix Factorization -. JSBi Bioinformatics Review, 1(2), 18-25. 10.11234/jsbibr.2021.1
• [2]  Tsuyuzaki, K. (2021/10/5). Mathematics for Heterogeneous Biological Data Fusion Analysis with Matrix-Tensor Factorization - Part II. Simultaneous Matrix Factorization -. JSBi Bioinformatics Review, 2(1), 15-29. 10.11234/jsbibr.2021.6
• [3]  Tsuyuzaki, K. (2022/6/2). Mathematics for Heterogeneous Biological Data Fusion Analysis with Matrix-Tensor Factorization - Part III. Tensor Factorization -. JSBi Bioinformatics Review, 3(1), 20-33. 10.11234/jsbibr.2022.1
• [4]  Tsuyuzaki, K. (2022/11/1). Mathematics for Heterogeneous Biological Data Fusion Analysis with Matrix-Tensor Factorization - Part IV. Categorical Data, Distance, and Graph -. JSBi Bioinformatics Review, 3(2), 33-46. 10.11234/jsbibr.2022.2
• [5]  Tsuyuzaki, K. (2024/6/3). Mathematics for Heterogeneous Biological Data Fusion Analysis with Matrix-Tensor Factorization - Part V. Random Walk -. JSBi Bioinformatics Review, 5(1), 1-15. 10.11234/jsbibr.2024.1
• [6]  Chitforoushzadeh, Z. & Ye, Z. & Sheng, Z. & LaRue, S. & Fry, R.C. & Lauffenburger, D.A. & Janes, K.A. (2016). TNF-insulin crosstalk at the transcription factor GATA6 is revealed by a model that links signaling and transcriptomic data tensors. Sci Signal. 9(431), ra59. 10.1126/scisignal.aad3373
• [7]  Martino, C. & Shenhav, L. & Marotz, C.A. & Armstrong, G. & McDonald, D. & Vázquez-Baeza, Y. & Morton, J.T. & Jiang, L. & Dominguez-Bello, M.G. & Swafford, A.D. & Halperin, E. & Knight, R. (2020). Context-aware dimensionality reduction deconvolutes gut microbial community dynamics. Nat Biotechnol. 39(2), 165-168. 10.1038/s41587-020-0660-7
• [8]  Williams, A.H. & Kim, T.H. & Wang, F. & Vyas, S. & Ryu, S.I. & Shenoy, K.V. & Schnitzer, M. & Kolda, T.G. & Ganguli, S. (2018). Unsupervised Discovery of Demixed, Low-Dimensional Neural Dynamics across Multiple Timescales through Tensor Component Analysis. Neuron, 98(6), 1099-1115.e8. 10.1016/j.neuron.2018.05.015
• [9]  Chatzichristos, C. & Kofidis, E. & Van Paesschen, W. & De Lathauwer, L. & Theodoridis, S. & Van Huffel, S. (2022). Early soft and flexible fusion of electroencephalography and functional magnetic resonance imaging via double coupled matrix tensor factorization for multisubject group analysis. Hum Brain Mapp. 43(4), 1231-1255. doi: 10.1002/hbm.25717
• [10]  Cong, F. & Lin, Q.H. & Kuang, L.D. & Gong, X.F. & Astikainen, P. & Ristaniemi, T. (2015). Tensor decomposition of EEG signals: a brief review. J Neurosci Methods. 2015 Jun 15; 248: 59-69. 10.1016/j.jneumeth.2015.03.018
• [11]  Qiu, P. & Simonds, E.F. & Bendall, S.C. & Gibbs, K.D. Jr. & Bruggner, R.V. & Linderman, M.D. & Sachs, K. & Nolan, G.P. & Plevritis, S.K. (2011). Extracting a cellular hierarchy from high-dimensional cytometry data with SPADE. Nat Biotechnol. 29(10), 886-91. 10.1038/nbt.1991
• [12]  Trapnell, C. & Cacchiarelli, D. & Grimsby, J. & Pokharel, P. & Li, S. & Morse, M. & Lennon, N.J. & Livak, K.J. & Mikkelsen, T.S. & Rinn, J.L. (2014). The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells. Nat Biotechnol. 32(4), 381-386. 10.1038/nbt.2859
• [13]  Zhu, Y. & Chen, Z. & Zhang, K. & Wang, M. & Medovoy, D. & Whitaker, J.W. & Ding, B. & Li, N. & Zheng, L. & Wang, W. (2016). Constructing 3D interaction maps from 1D epigenomes. Nat Commun. 10.1038/ncomms10812.
• [14]  Schreiber, J. & Durham, T. & Bilmes, J. & Noble, W.S. (2021). Avocado: a multi-scale deep tensor factorization method learns a latent representation of the human epigenome. Genome Biol. 21(1), 81. 10.1186/s13059-020-01977-6
• [15]  Townes, F.W. & Engelhardt, B.E. (2023). Nonnegative spatial factorization applied to spatial genomics. Nat Methods. 20(2), 229-238. 10.1038/s41592-022-01687-w
• [16]  Aonishi, T. & Maruyama, R. & Ito, T. & Miyakawa, H. & Murayama, M. & Ota, K. (2022). Imaging data analysis using non-negative matrix factorization. Neurosci Res. 179, 51-56. 10.1016/j.neures.2021.12.001
• [17]  Frelat, R. & Lindegren, M. & Dencker, T.S. & Floeter, J. & Fock, H.O. & Sguotti, C. & Stäbler, M. & Otto, S.A. & Möllmann, C. (2018). Community ecology in 3D: Tensor decomposition reveals spatio-temporal dynamics of large ecological communities. PLoS One. 12(11), e0188205. 10.1371/journal.pone.0188205
• [18]  Onken, A. & Liu, J.K. & Karunasekara, P.P. & Delis, I. & Gollisch, T. & Panzeri, S. (2016). Using Matrix and Tensor Factorizations for the Single-Trial Analysis of Population Spike Trains. PLoS Comput Biol. 12(11), e1005189. 10.1371/journal.pcbi.1005189
• [19]  Vanhatalo, E. & Kulahci, M. (2016) Impact of Autocorrelation on Principal Components and Their Use in Statistical Process Control. Qual. Reliab. Engng. Int., 32, 1483-1500. 10.1002/qre.1858
• [20]  ビショップ，C.M. & 元田浩，栗田多喜夫，樋口知之，松本裕治，村田昇．（2012）．パターン認識と機械学習 下（ベイズ理論による統計的予測），丸善出版
• [21]  Shang, L. & Zhou, X. (2022). Spatially aware dimension reduction for spatial transcriptomics. Nat Commun. 13(1):7203. 10.1038/s41467-022-34879-1.
• [22]  Hao, X. & Li, C. & Yan, J. & Yao, X. & Risacher, S.L. & Saykin, A.J. & Shen, L. & Zhang, D (2017). Alzheimer’s Disease Neuroimaging Initiative. Identification of associations between genotypes and longitudinal phenotypes via temporally-constrained group sparse canonical correlation analysis. Bioinformatics. 33(14), i341-i349. 10.1093/bioinformatics/btx245
• [23]  Yu, H.-F. & Rao, N. & Dhillon, I.S. (2016). Temporal regularized matrix factorization for high-dimensional time series prediction. In Proceedings of the 30th International Conference on Neural Information Processing Systems (NIPS’16), 847-855.
• [24]  Park, J.Y. & Carr, K. & Zheng, S. & Yue, Y. & Yu, R. (2020). Multiresolution Tensor Learning for Efficient and Interpretable Spatial Analysis. Proceedings of the 37th International Conference on Machine Learning, in Proceedings of Machine Learning Research
• [25]  Yokota, T. & Zhao, Q. & Cichocki, A. (2016). Smooth PARAFAC Decomposition for Tensor Completion, in IEEE Transactions on Signal Processing, 64(20), 5423-5436. 10.1109/TSP.2016.2586759
• [26]  竹内孝，鹿島久嗣，上田修功，（2018）．自己回帰テンソル分解による時空間データ予測，2018年度人工知能学会全国大会（第32回），10.11517/pjsai.JSAI2018.0_3L102
• [27]  田中聡久，横田達也，（2024）．テンソルデータ解析の基礎と応用：テンソル表現，縮約計算，テンソル分解と低ランク近似（次世代信号情報処理シリーズ7），コロナ社
• [28]  Guo, W. & Kotsia, I. & Patras, I. (2012). Tensor learning for regression. IEEE Trans Image Process. 21(2), 816-27. 10.1109/TIP.2011.2165291
• [29]  Velu, R.P. & Reinsel, G.C. & Wichern, D.W. (1986). Reduced Rank Models for Multiple Time Series. Biometrika, 73(1), 105-11
• [30]  Reinsel, G.C. & Velu, R.P. & Chen, K. (2022). Multivariate Reduced-Rank Regression: Theory, Methods and Applications (Lecture Notes in Statistics, 225), Springer
• [31]  Wang, D. & Zheng, Y. & Lian, H. & Li, G. (2021). High-Dimensional Vector Autoregressive Time Series Modeling via Tensor Decomposition. Journal of the American Statistical Association, 117(539), 1338-1356. 10.1080/01621459.2020.1855183
• [32]  Dong, Y. & Qin, S.J. (2018). A novel dynamic PCA algorithm for dynamic data modeling and process monitoring, Journal of Process Control, 67, 1-11, 10.1016/j.jprocont.2017.05.002
• [33]  Fei, Z. & Liu, K. (2016). Online process monitoring for complex systems with dynamic weighted principal component analysis, Chinese Journal of Chemical Engineering, 24(6), 775-786, 10.1016/j.cjche.2016.05.038
• [34]  T. Yoshioka, T. Nakatani, M. Miyoshi and H.G. Okuno, Blind Separation and Dereverberation of Speech Mixtures by Joint Optimization, in IEEE Transactions on Audio, Speech, and Language Processing, 19(1), 69-84, 10.1109/TASL.2010.2045183.
• [35]  Dong, Y. & Qin, S.J. (2018). Regression on dynamic PLS structures for supervised learning of dynamic data, Journal of Process Control, 68, 64-72, 10.1016/j.jprocont.2018.04.006
• [36]  Bo, X. & Qinqin, Z. (2020). Online Quality-Relevant Monitoring with Dynamic Weighted Partial Least, Industrial & Engineering Chemistry Research, 59(48), 10.1021/acs.iecr.0c03901
• [37]  Yining Dong, S. & Qin, J. (2018). Dynamic-Inner Canonical Correlation and Causality Analysis for High Dimensional Time Series Data, IFAC-PapersOnLine, 51(18), 476-481, 10.1016/j.ifacol.2018.09.379
• [38]  Zhu, Q. (2020). Auto-regressive modeling with dynamic weighted canonical correlation analysis, Journal of Process Control, 95, 32-44, 10.1016/j.jprocont.2020.09.007
• [39]  Jing, P. & Su, Y. & Jin, X. & Zhang, C. (2019). High-Order Temporal Correlation Model Learning for Time-Series Prediction, in IEEE Transactions on Cybernetics, 49(6), 2385-2397, 10.1109/TCYB.2018.2832085
• [40]  石井健一郎，上田修功．（2014）．続・わかりやすいパターン認識―教師なし学習入門．オーム社
• [41]  阿久津達也，（2024）．バイオインフォマティクスのための人工知能入門：基礎から行列・テンソル分解/深層学習まで．朝倉書店
• [42]  阿久津達也．（2007）．バイオインフォマティクスの数理とアルゴリズム．立出版
• [43]  Pevzner, P.A. & Jones, N.C. & 渋谷哲朗．（2007）．バイオインフォマティクスのためのアルゴリズム入門共立出版
• [44]  Fu, X. & Huang, K. & Sidiropoulos, N.D. & Ma, W. -K. (2019). Nonnegative Matrix Factorization for Signal and Data Analytics: Identifiability, Algorithms, and Applications, in IEEE Signal Processing Magazine, 36(2), 59-80, 10.1109/MSP.2018.2877582.
• [45]  北川源四郎，（2020）．Rによる時系列モデリング入門，岩波書店
• [46]  島田直希，（2021）．時系列解析 自己回帰型モデル・状態空間モデル・異常検知　Advanced Python共立出版
• [47]  Ketelaere, D. & Hubert, B.M. & Schmitt, E. (2015). Overview of PCA-Based Statistical Process-Monitoring Methods for Time-Dependent, High-Dimensional Data. Journal of Quality Technology, 47(4), 318-335. 10.1080/00224065.2015.11918137
• [48]  Jeyanthi, R. & Sahithi, M. & Sireesha, N.V.L. & Srinivasan, MS. & Devanathan, S. & Thampi, S.M. & El-Alfy, E.-S. M. & Trajkovic, L. (2021). Data reconciliation using MA-PCA and EWMA-PCA for large dimensional data. J. Intell. Fuzzy Syst. 41(5), 5731-5736. 10.3233/JIFS-189892
• [49]  Rejer, I. & Górski, P. (2019). MAICA: an ICA-based method for source separation in a low-channel EEG recording. J Neural Eng. 16(5), 056025. 10.1088/1741-2552/ab36db
• [50]  Wang, H. & Zhang, Q. & So, H.Y. & Kwok, A. & Wong, Z.S. (2018). Temporal prediction of in-hospital falls using tensor factorisation. BMJ Innov. 4(2), 75-83. 10.1136/bmjinnov-2017-000221
• [51]  新田益大，杉本謙二. （2007）. 独立成分分析に基づく多入力多出力システムのブラインド同定. 電子情報通信学会論文誌A，J90-A (1)，27-34.
• [52]  Alshammri, F. & Pan, J. (2021). Generalized principal component analysis for moderately non-stationary vector time series, Journal of Statistical Planning and Inference, 212, 201-225, 10.1016/j.jspi.2020.08.007
• [53]  渕上壮太郎，（2014）. 独立成分分析tICAでタンパク質の複雑な運動を解きほぐす（特集 生体高分子の揺らぎとダイナミクス : シミュレーションと実験の統計解析）統計数理，62(2)，243-255
• [54]  Mori T. & Saito, S. (2015). Dynamic heterogeneity in the folding/unfolding transitions of FiP35. J Chem Phys. 142(13), 135101. 10.1063/1.4916641.
• [55]  Thirion, B. & Faugeras, O.D. (2002). Dynamical components analysis of fMRI data. Proceedings IEEE International Symposium on Biomedical Imaging, 915-918.
• [56]  Kutz, J.N. & Brunton, S.L. & Brunton, B.W. & Proctor, J.L. (2017). Dynamic Mode Decomposition: Data-Driven Modeling of Complex Systems, Society for Industrial and Applied Mathematics
• [57]  Wang, Z. & Wang, D. & Li, C. & Xu, Y. & Li, H. & Bao, Z. (2018). Deep reinforcement learning of cell movement in the early stage of C. elegans embryogenesis. Bioinformatics. 34(18), 3169-3177. 10.1093/bioinformatics/bty323
• [58]  Vervliet, N. & Debals, O. & Sorber, L. & Barel, M.V. & Lathauwer, L.D. (2017/5/9). Tensorization techniques - Tensorlab 3.0 documentation, https://www.tensorlab.net/doc/tensorization.html
• [59]  Debals, O. (2017/8). Tensorization and applications in blind source separation. A Doctrial thesis in KU KEUVEN
• [60]  井手剛，井上恵介，（2004）．非線形変換を利用した時系列データからの知識発見，第4回データマイニングワークショップ，日本ソフトウェア科学会データマイニング研究会
• [61]  Hassani, H. & Ghodsi, Z. (2015). A glance at the applications of Singular Spectrum Analysis in gene expression data. Biomol Detect Quantif. 4, 17-21. 10.1016/j.bdq.2015.04.001
• [62]  Ikeda, M. & Matsumoto, H. & Izquierdo, E.J. (2021). Persistent thermal input controls steering behavior in Caenorhabditis elegans. PLoS Comput Biol. 17(1), e1007916. 10.1371/journal.pcbi.1007916
• [63]  Kouchaki, S. & Sanei, S. & Arbon, E.L. & Dijk, D.J. (2015). Tensor based singular spectrum analysis for automatic scoring of sleep EEG. IEEE Trans Neural Syst Rehabil Eng. 23(1), 1-9. 10.1109/TNSRE.2014.2329557
• [64]  Usevich, K., & Golyandina, N. (2010). 2D-extension of Singular Spectrum Analysis: algorithm and elements of theory.
• [65]  Yokota, T. & Erem, B. & Guler, S. & Warfield, S.K. & Hontani, H. (2018). Missing Slice Recovery for Tensors Using a Low-Rank Model in Embedded Space, IEEE/CVF Conference on Computer Vision and Pattern Recognition, 8251-8259, 10.1109/CVPR.2018.00861
• [66]  Shi, Q. & Yin, J. & Cai, J. & Cichocki, A. & Yokota, T. & Chen, L. & Yuan, M. & Zeng, J. (2020). Block Hankel Tensor ARIMA for Multiple Short Time Series Forecasting. AAAI Conference on Artificial Intelligence.
• [67]  Paparrizos, J. & Gravano, L. (2016). K-Shape: Efficient and Accurate Clustering of Time Series. SIGMOD Rec. 45(1), 69-76. 10.1145/2949741.2949758
• [68]  Wang, X. & Mueen, A. & Ding, H. (2013). Experimental comparison of representation methods and distance measures for time series data. Data Min Knowl Disc 26, 275-309. 10.1007/s10618-012-0250-5
• [69]  Butler, A. & Hoffman, P. & Smibert, P. (2018). Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nat Biotechnol 36, 411-420. 10.1038/nbt.4096
• [70]  Tsuyuzaki, K. & Yamamoto, K. & Toyoshima, Y. & Sato, H. & Kanamori, M. & Teramoto, T. & Ishihara, T. & Iino Y. & Nikaido, I. (2023). WormTensor: a clustering method for time-series whole-brain activity data from C. elegans. BMC Bioinformatics. 24(1), 254. 10.1186/s12859-023-05230-2
• [71]  Liu, X. & Ji, S. & Glänzel, W. & De Moor, B. (2013). Multiview Partitioning via Tensor Methods, in IEEE Transactions on Knowledge and Data Engineering, 25(5), 1056-1069, 10.1109/TKDE.2012.95
• [72]  Williams, A.H. & Poole, B. & Maheswaranathan, N. & Dhawale, A.K. & Fisher, T. & Wilson, C.D. & Brann, D.H. & Trautmann, E.M. & Ryu. S. & Shusterman, R. & Rinberg, D. & Ölveczky, B.P. & Shenoy, K.V. & Ganguli, S. (2020). Discovering Precise Temporal Patterns in Large-Scale Neural Recordings through Robust and Interpretable Time Warping. Neuron. 105(2), 246-259.e8, 10.1016/j.neuron.2019.10.020.
• [73]  Williams, A.H. (2020). Combining tensor decomposition and time warping models for multi-neuronal spike train analysis, bioRxiv 10.1101/2020.03.02.974014
• [74]  Mørup, M. (2008). Decomposition methods for unsupervised learning, SeriesDTU Compute PHD
• [75]  Belouchrani, A. & Abed-Meraim, K. & Cardoso, J. & Moulines, E. (1997). A blind source separation technique using second-order statistics, in IEEE Transactions on Signal Processing, 45(2), 434-444, 10.1109/78.554307
• [76]  Pfister, N. & Weichwald, S. & Bühlmann, P. & Scholkopf, B. (2018). Robustifying Independent Component Analysis by Adjusting for Group-Wise Stationary Noise. J. Mach. Learn. Res., 20, 147:1-147:50
• [77]  Mørup, M. & Madsen, K.H. & Hansen, L.K. (2007). Shifted Independent Component Analysis. In: Davies, M.E., James, C.J., Abdallah, S.A., Plumbley, M.D. (eds) Independent Component Analysis and Signal Separation. ICA 2007. Lecture Notes in Computer Science, vol 4666. Springer, 10.1007/978-3-540-74494-8_12
• [78]  Morup, M. & Madsen, K.H. & Hansen, L.K. (2007). Shifted Non-Negative Matrix Factorization, 2007 IEEE Workshop on Machine Learning for Signal Processing, 139-144, 10.1109/MLSP.2007.4414296
• [79]  Mørup, M. & Hansen, L.K. & Arnfred, S.M. & Lim, L.-H. & Madsen, K.H. (2008). Shift-invariant multilinear decomposition of neuroimaging data, NeuroImage, 42(4), 1439-1450, 10.1016/j.neuroimage.2008.05.062
• [80]  Cichocki, A. & Zdunek, R. & Phan, A.H. & Amari, S. (2009/10/12). Nonnegative Matrix and Tensor Factorizations: Applications to Exploratory Multi-way Data Analysis and Blind Source Separation. Wiley
• [81]  Yamazaki, A. & Tajima, T. & Matsuoka, K. Convolutive independent component analysis of EEG data, SICE 2003 Annual Conference (IEEE Cat. No.03TH8734), 1227-1231, 2
• [82]  Smaragdis, P. (2004). Non-negative Matrix Factor Deconvolution; Extraction of Multiple Sound Sources from Monophonic Inputs. In: Puntonet, C.G., Prieto, A. (eds) Independent Component Analysis and Blind Signal Separation. ICA 2004. Lecture Notes in Computer Science, 3195. Springer, 10.1007/978-3-540-30110-3_63
• [83]  Mørup, M. & Schmidt, M.N. (2006). Sparse non-negative tensor factor double deconvolution (SNTF2D) for multi channel time-frequency analysis. Technical report, University of Denmark, DTU
• [84]  Dyrholm, M. & Hansen, L.K. (2004). CICAAR: Convolutive ICA with an Auto-regressive Inverse Model. In: Puntonet, C.G., Prieto, A. (eds) Independent Component Analysis and Blind Signal Separation. ICA 2004. Lecture Notes in Computer Science, 3195. Springer, 10.1007/978-3-540-30110-3_76
• [85]  Arora, R. & Cotter, A. & Livescu, K. & Srebro, N. Stochastic optimization for PCA and PLS, 2012 50th Annual Allerton Conference on Communication, Control, and Computing (Allerton), 861-868, 10.1109/Allerton.2012.6483308.
• [86]  N. Mohammadiha, P. & Smaragdis, G. & Panahandeh Doclo, S. (2015). A State-Space Approach to Dynamic Nonnegative Matrix Factorization,” in IEEE Transactions on Signal Processing, 63(4), 949-959, 10.1109/TSP.2014.2385655
• [87]  Sun, J. & Tao, D. & Faloutsos, C. (2006). Beyond streams and graphs: dynamic tensor analysis. In Proceedings of the 12th ACM SIGKDD international conference on Knowledge discovery and data mining (KDD '06). Association for Computing Machinery, 374-383. 10.1145/1150402.1150445
• [88]  Maehara, T. & Hayashi, K. & Kawarabayashi, K. (2016). Expected Tensor Decomposition with Stochastic Gradient Descent. Proceedings of the AAAI Conference on Artificial Intelligence, 30(1). 10.1609/aaai.v30i1.10292
• [89]  Tsuyuzaki, K. & Sato, H. & Sato, K. & Nikaido, I. (2020). Benchmarking principal component analysis for large-scale single-cell RNA-sequencing. Genome Biology, 21(1):9, 10.1186/s13059-019-1900-3
• [90]  Li, Z. & Meisner, J. & Albrechtsen, A. (2023). Fast and accurate out-of-core PCA framework for large scale biobank data. Genome Res. 33(9), 1599-1608. 10.1101/gr.277525.122
• [91]  Zhou, J. & He, Z. & Wang, J. & Hu, S. & Kan, T. & Dong, G. & Shi, J. & Liu, R. & Ou-Yang, L. & Bo, X. & Ying, X. (2024). Continual integration of single-cell multimodal data with MIRACLE. bioRxiv, 10.1101/2024.09.24.613833
• [92]  Iqbal, N. & Liu, E. & McClellan, J.H. & Al-Shuhail, A. & Kaka, S.I. & Zerguine, A. (2018). “Detection and Denoising of Microseismic Events Using Time-Frequency Representation and Tensor Decomposition,” in IEEE Access, 6, 22993-23006, 10.1109/ACCESS.2018.2830975
• [93]  Hejazi, F. & Joneidi, M. & Rahnavard, N. (2020). A Tensor-based Localization Framework Exploiting Phase Interferometry Measurements, 2020 IEEE International Radar Conference (RADAR), 554-559, 10.1109/RADAR42522.2020.9114560
• [94]  Sawada, H. & Ono, N. & Kameoka, H. & Kitamura, D. & Saruwatari, H. (2019). A review of blind source separation methods: two converging routes to ILRMA originating from ICA and NMF. APSIPA Transactions on Signal and Information Processing. 8, e12. 10.1017/ATSIP.2019.5
• [95]  戸上真人，（2020）．Pythonで学ぶ音源分離 機械学習実践シリーズ，インプレス
• [96]  Sekiguchi, K. & Bando, Y. & Nugraha, A.A. & Fontaine, M. Yoshii, K. (2021). Autoregressive Fast Multichannel Nonnegative Matrix Factorization For Joint Blind Source Separation And Dereverberation, ICASSP 2021-2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 511-515, 10.1109/ICASSP39728.2021.9414857
• [97]  Sekiguchi, K. & Bando, Y. & Nugraha, A.A. & Fontaine, M. & Yoshii, K. & Kawahara, T. (2022). Autoregressive Moving Average Jointly-Diagonalizable Spatial Covariance Analysis for Joint Source Separation and Dereverberation, in IEEE/ACM Transactions on Audio, Speech, and Language Processing, 30, 2368-2382, 10.1109/TASLP.2022.3190734
• [98]  Kagami, H. & Kameoka, H. & Yukawa, M. Joint Separation and Dereverberation of Reverberant Mixtures with Determined Multichannel Non-Negative Matrix Factorization, 2018 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 31-35, 10.1109/ICASSP.2018.8462080
• [99]  Hsu, L.L. & Culhane, A.C. (2020). Impact of Data Preprocessing on Integrative Matrix Factorization of Single Cell Data. Front Oncol. 10, 973. 10.3389/fonc.2020.00973
• [100]  　Shah, N. & Meng, Q. & Zou, Z. & Zhang, X. (2024). Systematic analysis on the horse-shoe-like effect in PCA plots of scRNA-seq data. Bioinform Adv. 4(1), vbae109. 10.1093/bioadv/vbae109
• [101]  　Diaconis, P. & Goel, S. & Holmes, S.P. (2008). HORSESHOES IN MULTIDIMENSIONAL SCALING AND LOCAL KERNEL METHODS. The Annals of Applied Statistics, 2, 777-807.
• [102]  　Morton, J.T. & Toran, L. & Edlund, A. & Metcalf, J.L. & Lauber, C. & Knight, R. (20187). Uncovering the Horseshoe Effect in Microbial Analyses. mSystems. 2(1), e00166-16. 10.1128/mSystems.00166-16
• [103]  　Proix, T. & Perich, M.G. & Milekovic, T. (2022). Misinterpreting the horseshoe effect in neuroscience, bioRxiv, 10.1101/2022.03.04.482986
• [104]  　Novembre, J. & Stephens, M. (2008). Interpreting principal component analyses of spatial population genetic variation. Nat Genet. 40(5), 646-9. 10.1038/ng.139
• [105]  　Frichot, E. & Schoville, S. & Bouchard, G. & François, O. (2012). Correcting principal component maps for effects of spatial autocorrelation in population genetic data. Front Genet. 3, 254. 10.3389/fgene.2012.00254
• [106]  　Shinn, M. (2023). Phantom oscillations in principal component analysis. PNAS. 120(48), e2311420120. 10.1073/pnas.2311420120
• [107]  　小林徹也，杉村薫，舟橋啓，（2020）．実験医学増刊 Vol.38 No.20 機械学習を生命科学に使う! ～シークエンスや画像データをどう解析し、新たな生物学的発見につなげるか?，羊土社
• [108]  　井手剛，（2006）．部分時系列クラスタリングの理論的基礎，2006年度人工知能学会全国大会（第20回），2A1-2, 10.11517/pjsai.JSAI06.0.94.0
• [109]  　Hsu, L.L. & Culhane, A.C. (2023). Correspondence analysis for dimension reduction, batch integration, and visualization of single-cell RNA-seq data. Scientific Report, 13(1197), 10.1038/s41598-022-26434-1
• [110]  　Meng, C. & Zeleznik, O.A. & Thallinger, G.G. & Kuster, B. & Gholami, A.M. & Culhane, A.C. (2016). Dimension reduction techniques for the integrative analysis of multi-omics data. Brief Bioinform. 17(4), 628-41. 10.1093/bib/bbv108
• [111]  　山中武彦，浜崎健児，嶺田拓也，（2005）. 生物・社会調査のための統計解析入門：調査・研究の現場から（その9）―序列化する（対応分析、除歪対応分析、正準対応分析―，農士誌 73(4)，319-324
• [112]  　Podani, J. and Miklós, I. (2002), Resemblance Coefficients and the Horseshoe Effect in Principal Coordinates Analysis. Ecology, 83, 3331-3343. 10.1890/0012-9658(2002)083[3331:RCATHE]2.0.CO;2
• [113]  　Camiz, S. (2004). The Guttman effect: its interpretation and a new redressing method.
• [114]  　Pillay, A.B. & Pathmanathan, D. & Dabo-Niang, S. & Abu, A. & Omar, H. (2024). Functional data geometric morphometrics with machine learning for craniodental shape classification in shrews. Sci Rep. 14(1), 15579. 10.1038/s41598-024-66246-z
• [115]  　系列データにおける馬蹄効果の確認 https://qiita.com/antiplastics/private/72f3ba4f92a1f0f06a3e
• [116]  　松本拡高，露崎弘毅（2025）．系列データにおける正弦波効果・馬蹄効果への省察―直鎖グラフにおけるラプラシアン固有マップの視点から―, JSBi Bioinformatics Review, 6(1), 29-40.

略語リスト

・i.i.d.

independent and identically distributed

・PCA

Principal Component Analysis

・SNP

Single Nucleotide Polymorphism

・CCA

Canonical Correlation Analysis

・AR

Auto-regressive

・VAR

Vector Autoregression Model

・SVD

Singular Value Decomposition

・HOSVD

Higher-order Singular Value Decomposition

・ICA

Independent Component Analysis

・PLS

Partial Least Squares

・MOAR

Multilinear Orthogonal Autoregressive

・HMM

Hidden Markov Model

・NMF

Non-negative Matrix Factorization

・MA

Moving Average

・EWMA

Exponentially Weighted Moving Average

・MDPCA

Moving Dynamic PCA

・GDPCA

Generalized Dynamic PCA

・TS-PCA

PCA for time series

・GTS-PCA

Generalized PCA for moderately non-stationary vector time series

・tICA

Time-structure based Independent Component Analysis

・DCA

Dynamic Component Analysis

・EVD

Eigen Value Decomposition

・DMD

Dynamic Mode Decomposition

・SSA

Singular Spectrum Analysis

・MSR

Missing Slice Recovery

・ARIMA

Auto Regressive Integrated Moving Average

・BHT-ARIMA

Block Hankel Tensor ARIMA

・SBD

Shape-based distance

・DTW

Dynamic Time Warping

・NMFD

Non-negative Matrix Factor Deconvolution

・DTA

Dynamic Tensor Analysis

・MIRACLE

Multimodal Integration with Continual Learning

・DFT

Discrete Fourier transform

・BSS

Blind Source Separation

・FDICA

Frequency-domain Independent Component Analysis

・IVA

Independent Vector Analysis

・MNMF

Multichannel Non-negative Matrix Factorization

・ILRMA

Independent Low-Rank Matrix Analysis

・CA

Correspondence Analysis

・MDS

Multi-Dimensional Scaling

・LE

Laplacian Eigenmap

・DCA

Detrended Correspondence Analys

著者略歴

	露崎 弘毅 2015年、東京理科大学生命創薬科学科博士後期課程終了。博士（薬科学）。同年より、理化学研究所（所属3）に在籍し、シングルセルオミックスのデータ解析や解析ツール開発に従事。現在は千葉大学（所属1, 2）で医療データの解析を行なっている。パッケージングに特化したハッカソンBio“Pack”athonを主催。 ホームページ：https://researchmap.jp/kokitsuyuzaki
松本 拡高 2016年、東京大学大学院新領域創成科学研究科博士課程終了。博士（科学）。理化学研究所生命機能科学研究センターや革新知能統合研究センターでの研究員を経て、現在は長崎大学情報データ科学部に在籍し、多様な生命現象のデータの解析を行なっている。