系列データにおける正弦波効果・馬蹄効果への省察 ―直鎖グラフにおけるラプラシアン固有マップの視点から―

松本 拡高; 露崎 弘毅

doi:10.11234/jsbibr.2025.2

Abstract

系列データに基づく様々な行列に対し行列分解を行うと、固有ベクトルなどに三角関数のような形状が現れることが経験的に知られている。この様な現象は馬蹄効果・ファントム振動・正弦波効果などと呼ばれ、多様なシナリオで現れることが報告されてきた。本稿では、そのような三角関数様のパターンが現れるメカニズムを、パスグラフに対するラプラシアン固有マップの理論的な結果を軸とし、多様な問題が類似の問題に帰着できることを利用し、直感的に説明することを目指す。

1．イントロダクション

時系列データやゲノム座標と紐づいたデータなど、同質のデータを並べた順序が存在するデータを系列データとよぶ。このような系列データを解析する上で、行列・テンソル分解は基盤的な技術である［1］。また、1細胞RNA-seqにおける分化進行度（擬時間）の推定［2］のように、背後に時間軸が潜んでいるデータに対し次元圧縮とその空間上で細胞系譜を再構築するといった行列分解を用いた発展的な手法も登場してきた。

一方で、このような系列データに対し行列・テンソル分解を行うと、計測データ中には存在しないような振動パターンが自ずと出てくることが多数報告されている。このような傾向は場合によっては馬蹄効果やファントム振動などと呼ばれ、誤った解釈をもたらす危険性があることが指摘されている［3］。また同様の問題として、系列データからパターンを抽出するアプローチとして利用されていた、スライディングウィンドウなどの操作で全体の系列データから抜き出した部分的なデータである部分系列をクラスタリングして得られるパターンとは、入力とした系列データに関わらず正弦波が自ずと出てきて、無意味なものであることが報告・研究されている［4, 5］。

本記事は、本誌同号に同時掲載される「行列・テンソル分解によるヘテロバイオデータ統合解析の数理　―第6回系列データ―」［1］を補足するものとして、系列データの行列分解において三角関数と類似する振動パターンが現れる背景を、時に理論的に、時に直感的に説明を試みるものである。本記事ではまず、n個のノードが直鎖状に連結したパスグラフ（図1a）に対するラプラシアン固有マップ（Laplacian Eigenmap; LE）［6］の結果が三角関数になることを示す。次に、系列データに対するハンケル行列の特異値分解（Singular Value Decomposition; SVD）の結果として、左特異ベクトルにおいて三角関数が現れるメカニズムを、パスグラフのLEの類似性の観点から直感的に説明する。最後に、背後に時系列が存在するようなデータ行列に対するPCAにおいて、振動パターンが容易に現れるメカニズムについても、同様の観点から直感的に説明する。

図１：パスグラフとリンググラフに対するLE

a．パスグラフとそのグラフラプラシアン。b．n＝100としたパスグラフの固有ベクトル。c．パスグラフのLEの第1、2軸。d．2n＝200としたリンググラフとそのグラフラプラシアン。e．パスグラフとリンググラフの対応　f．リンググラフのLEの第1、2軸　g、h、i. パスグラフを模したグラフとそのグラフラプラシアンの固有ベクトルおよびLEの第1、2軸。

2．LEによる系列データの埋め込み

LEとは、ノードi, j間の類似性が非負な重みW_{i, j}で定義される無向グラフが与えられた時、ノードを低次元空間に埋め込むアプローチの1つである。LEでは、D_{i, i}＝∑_jW_{i, j}で定義される対角な次数行列Dを用い、L＝D－Wとして定義されるグラフラプラシアンLの固有値分解を行い、固有値が小さいものから順番に対応する固有ベクトルをデータの要約とした座標として利用するものである¹。

2.1　パスグラフにおけるLE

ここでは、系列データを模したグラフ構造として、n個のノードが直鎖状に連結したパスグラフPを考える（図1a）。つまり、W_{i, i}_＋1＝W_i_＋1,_i＝1 (i∈［1, n－1］)と、それ以外の要素はW_{i, j}＝0となるグラフである。このとき、グラフラプラシアンL_Pは図1a下に示されるものとなる。

詳細は補足説明6.1にて示すが、このグラフラプラシアンのk番目の固有値（k∈［0, n－1］）は、λ_k＝2－2cos(πk/n)であり、その固有ベクトルv_kのi番ノードの値はv_k(i)＝cos(πk(i－1/2)/n)となることが知られている（図1b）²［7］。よって、パスグラフのLEの第1、2軸は図1cとなり、ノードがU字に配置されるという馬蹄効果が現れることがわかる。

2.2　リンググラフにおけるLE

パスグラフにおいて周期的な固有ベクトルが現れるメカニズムを理解するため、2n個のノードが環状に連結したリンググラフRを考える（図1d）。つまり、W_{i, i}_＋1＝W_i_＋1_{, i}＝1 (i∈［1, 2n－1］)とW_{1, 2}_n＝W₂_n_{, 1}＝1で、それ以外の要素はW_{i, j}＝0となるグラフであり、そのグラフラプラシアンL_Rは図1d下で示されるものである。

このグラフラプラシアンのk番目の固有値は、L_Pと同様にλ_k＝2－2cos(πk/n)であり、固有値が重複する固有ベクトルxとyを持つことが知られている［7］（補足説明6.2）。そして、それら固有ベクトルのi番ノードの値はそれぞれx_k(i)＝sin(πki/n)とy_k(i)＝cos(πki/n)となる。よって、リンググラフのLEの第1、2軸は図1fとなり、ノードが円周上に配置されることになる。

ここで、以下の関係式が成り立つことから、パスグラフの固有ベクトルとはリンググラフの固有ベクトルを合成したものであることがわかる。

v k ( i ) = cos ( πk ( i − 1 / 2 ) / n ) = cos ( πk / n ) y k ( i ) + sin ( πk / n ) x k ( i )

これはつまり、パスグラフの固有ベクトルとは、ノード数が2倍のリンググラフにおける固有ベクトルのノード1からnまでの値から合成されたものと見なせるわけである。詳しくは補足説明6.3にて示すが、仮想的にノードiとノード2n－i＋1を対応させると、パスグラフとリンググラフを結びつけることができ、パスグラフとはノード数が2nのリンググラフにおいて1からn番目のノードを抜き出したものであることが証明されている（図1e）。したがって、パスグラフとはリンググラフの部分構造であり、その固有ベクトルにはリンググラフと同様の周期的な構造が現れることになる。そしてパスグラフのk番目の固有ベクトルに現れる周期は、リンググラフの同一固有値の固有ベクトルの周期と一致し、2n/kとなることがわかる。

2.3　パスグラフを模したグラフにおけるLE

本章の最後に、パスグラフを模したグラフにおいても、同様に三角関数様の固有ベクトルが現れることを実験的に示す。ここでは図1gのように、ノードのインデックスが近いもの同士の類似性が高いグラフを考える。ただし、パスグラフとは異なり、非ゼロな重みの辺は一つ隣のノードに限定されず存在するとする。例えば、ノード間の類似度がW_{i, j}＝exp(－|i－j|/0.05n)として与えられるグラフを考えた場合、n＝100におけるグラフラプラシアンの上位の固有ベクトルは図1hとなり、三角関数様のパターンを示すことがわかり、このグラフのLEの第1、2軸は図1iでノードがU字に配置されるという馬蹄効果が現れることがわかる。このように、完全なパスグラフのみでなく、それに類するグラフでもグラフラプラシアンの固有ベクトルに三角関数様のパターンが現れる。

3．ハンケル行列を用いた系列データ解析

3.1　イントロダクション

x＝（x₁, x₂, . . . , x_n_＋_m）という系列データが与えられたとき、長さnの部分系列を列ベクトルとし、以下のように列ベクトルを並べて構築したn×mの行列をハンケル行列と呼ぶ。

H = ( x 1 x 2 . . . x m x 2 x 3 . . . x m + 1 x 3 x 4 . . . x m + 2 . . . . . . x n x n + 1 . . . x n + m )

このようなハンケル行列に対しSVDを行うことは、系列データをパターンに分解し解析する際や、系列データをクラスタリングする際に利用されている。しかし、このような行列に対するSVDで得られる左特異ベクトルにおいて、三角関数が自ずと出てくることが知られている。なお、左特異ベクトルはHH^Tの固有ベクトルと一致することから、以降では基本的にHH^Tのような行列を中心に議論する。

たとえば、n＝100, m＝10000として、平均0、分散1の独立なガウスノイズϵ_iに従うランダムウォークx_i＝x_i_－1＋ϵ_iから生成されたxに対し、HH^Tを計算した一例が図2aである。この行列に対する上位の固有ベクトルは図2bとなり、第2固有ベクトル以降、k番目の固有ベクトルの周期が2n/(k－1)となる三角関数に近い形状が出てくることが見て取れる。

図２：ハンケル行列およびそれを模した行列の形状と固有ベクトル

a、c、e、g. HH^T、JJ^T、H′H′^T、KK^Tの可視化（スケーリング済み）。b、d、f、h. それらの固有ベクトル。

以降ではまず、ハンケル行列を模した行列として独立な系列データを列ベクトルとして並べた行列を考え、列正規化を行なった時に左特異ベクトルに三角関数様のパターンが自ずと出てくることを説明する。次に、それを踏まえて連続系列でも同様の現象が起きることを解説する。また、独立系列において正規化しなかったとしても、三角関数様のパターンが現れることも解説する。

3.2　正規化した独立系列データの場合

まず初めに、計算しやすさのため、独立な系列データを列ベクトルとして並べた行列を考えることにする。ここでは、x＝（x₁, x₂, . . . , x_n）がランダムウォークに従うとし、仮想的な初期値としてx₀＝0を与えたとする。このような長さnの系列データを独立にm個生成し、それらを列ベクトルとて並べた行列を扱うとする。

この行列に対し、各列ベクトルの平均値が0となるように正規化した行列をJとする。実際に、n＝100, m＝10000として生成したJJ^Tを図2cに示す。この行列に対する固有ベクトルは、第k固有ベクトルの周期が2n/kとなる三角関数様のパターンが現れる（図2d）。

詳しくは補足説明6.4にて示すが、mが十分に大きい場合は大数の法則より、JJ^Tのi行j列成分（i≤j）はおおよそ(i－1/2)²＋(j－n－1/2)²で与えられるような構造になっている。以降では、i行j成分（i≤j）が(i－1/2)²＋(j－n－1/2)²で定義される対称行列Aに基づいて議論を進める。この行列Aでは、全ての列方向あるいは行方向に対しその和が(2n³＋n)/6で、最大固有値も(2n³＋n)/6であり、その固有ベクトルは定数ベクトルである。また、行列Aに対し何かしらの係数αを用いてαI－Aとした行列の固有値はAの固有値λ_kに対してα－λ_kが対応し、固有ベクトルは共通であることが一般に成り立つ。したがって、α＝(2n³＋n)/6とした行列αI－Aは、最小の固有値が0となり、その固有ベクトルは定数ベクトルである。ここでAをグラフの重み行列と見なすと、D_{i, i}＝∑_jA_{i, j}＝αであるから、L＝αI－Aとするとグラフラプラシアンの定義とも合っている。

したがって、Aの上位の固有ベクトルを求める操作とは、ノード間の辺の重みをW_ij＝A_ijとした類似性グラフに対しLEを行う問題へと帰着できる。ここで、A_ij＝(i－1/2)²＋(j－n－1/2)²であることを踏まえると、これは図3aのようにしてノード間の類似度を定義していると見なすことができる。すこし奇妙な類似度ではあるが、これは「直線上に並んだノード」にして、先の定義で類似度を計算した類似度グラフがAであると見なすことができる。そう考えると、これは直感的にはパスグラフのLEと類似した状況であり、三角関数様の固有ベクトルが現れることにも納得がいくだろう。実際に、k番目の固有ベクトルの周期は2n/kであり、これはパスグラフにおける周期と一致していることもわかる（図2d）。

図３：類似度の概念図

a．Aに対応する類似度の概念図。b．Bに対応する類似度の概念図。

3.3　連続系列データの場合

次に、1つの系列データに対し長さnの部分系列をスライドさせながら列ベクトルとして並べたハンケル行列の場合を考える。ここでランダムウォークに従う系列データをサンプリングしHH^Tを計算した一例を図2aに示す。この行列の固有ベクトルが図2bであるが、第2固有ベクトル以降は先ほどの正規化を行なった独立系列データの固有ベクトルとおおよそ一致していることがわかる（図2d）。

ここで、HH^Tの第1固有ベクトルはおおよそ定数の値のベクトルとなっていることがわかる（図2b）。というのも、SVDを低ランク近似と捉え、長さnのベクトルuと長さmのベクトルvによる階数1の近似H≃uσv^Tを解くことを考えると、特にトレンドのないランダムウォークを近似する場合には、v_iがHのi番目の列ベクトルのおおよその値を捉えるものになり、それに対するuとは部分系列内の変動を捉えるものではないため、おおよそ定数のベクトルとなる。

ここで、HをSVDしたのち、次の様に第1成分を除去するデフレーションを行った行列 H ' = U , − 1 Σ − 1 , − 1 V , − 1 T を求めたとする³。この行列に対し、H′H′^Tを計算した結果と固有ベクトルは、図2e、fのようになる。この結果は、正規化した独立系列のものとおおよそ一致することがわかるだろう（図2c、d）。というのも、左特異ベクトルがおおよそ定数のベクトルである第1成分をデフレーションする操作とは、Hに対して各列ベクトルを平均0となるように正規化することとおおよそ一致するからである。

以上のことから、正規化のような操作を行わなかったとしても、ウインドウをシフトさせて部分系列を並べたハンケル行列の左特異ベクトルに関して、2番目以降に同様の三角関数様のパターンが生じることがわかる。

3.4　正規化しない独立系列データの場合

最後に余談であるが、独立系列において正規化しなかったとしても、図2hのようにk番目の固有ベクトルにおいて周期が4n/(2k－1)となる三角関数様の形状が出てくることを紹介する。

行列Jにおいて列の正規化を行う前の行列をKとすると、KK^Tは図2gのような形状になる。詳しくは補足説明6.5で示すが、これは以下の行列Bと同様の構造である。

B = ( 1 1 1 . . . 1 1 2 2 . . . 2 1 2 3 . . . 3 . . . . . . . 1 2 3 . . . n )

これはつまり、B_ij＝min(i, j)の行列である。

以降では、このような行列に対しても、固有ベクトルに三角関数様の形状が現れることを直感的に説明する。ここで、以下のようにBを結合した行列B₂を考える。

B 2 = ( B + − B − − B | B )

ただし、B^|＝B_{:, n:1}のように行列の列順を逆にしたもので、B^－は行順を逆にしたもの、B^＋は両方を逆順にしたものを表す。このB₂に関して列の和は0であり、L＝－B₂とするとB₂の対角成分が次数を表し、B₂のi, j成分がノードの類似度を表すグラフのグラフラプラシアンと見なすことができる⁴。詳しくは補足説明6.5にて示すが、B₂の固有値固有ベクトルについては、元のBの固有値λ_kに対し、2λ_kの固有値が存在し、かつB₂の固有ベクトルのn＋1から2nのインデックスの成分がBの固有ベクトルと一致することになる。

いま、ノードi＝1, 2, . . . , 2nに対してノードiの座標をz_i＝i－n－1/2としてz軸上に配置したとする。このとき、ノードi, j間の類似性を以下のように定義する。

W i , j = sign ( z i z j ) * min ( | z i | , | z j | )

これは図3bのような状況である。これも奇妙な類似度の定義ではあるが、これまでと同様に直線上にノードを配置して類似度を測ったものであり、構造的にはパスグラフのLEと同様の問題であることがわかる。したがって、三角関数様のパターンが自ずと現れることが直感的に理解できるだろう。ただし、ここでは元の行列Bの2倍のサイズの行列B₂を考えていたことから、1番目の固有ベクトルの周期はこれまでの2倍の4nとなっている。また、k番目以降の周期は4n/(2k－1)となっている。

このように、系列データに対するSVDの左特異ベクトルには、いろいろな想定や操作によって容易に三角関数様の構造が現れる。そしてそのメカニズムとしては、パスグラフのLEとの類似性を踏まえると、ある程度は直感的に納得できるだろう。

4．主成分分析（PCA）

4.1　イントロダクション

サンプルサイズをnとし、m個の特徴量に関して標準化済みのn行m列のデータ行列Xが与えられた時、XX^Tの固有ベクトルを求めるものがPCAである。ここでは、背後に時系列が存在するような場合に、馬蹄効果がどのように現れるかを説明する。

4.2　背後に系列が存在する場合のPCA

以降では、i番目のサンプルは時刻t_i＝i/nのタイミングでサンプリングされたデータとする。そして、全ての特徴量は時系列に依存した関数によって与えられ、0＜t≤1の範囲において単調増加するものを想定する。すなわち、周期的な変動や振動が個々の特徴量のダイナミクス中には存在しない状況を想定する。以降では、このようなデータ行列のPCAにおいても容易に振動パターンが現れることを実験的に示したのち、その原因をパスグラフのLEと関連づけて説明する。

まず、時間に対し単調増加する単純な特徴量として、f₁(t)＝t, f₂(t)＝t², f₃(t)＝t³の3つの関数に基づくものを想定する（図4a）。最も単純なデータとして、ただ1つの特徴量f₁のみから構成されるn行1列のデータ行列X₁を考える。すると、 X 1 X 1 T は図4bのようになり、その固有ベクトルは図4cとなる。このような場合は、意味のある固有ベクトルはただ一つであり、それは時間軸を表す。次に、2つの特徴量f₁、f₂から構成されたn行2列のデータ行列X₂を考える。このとき、 X 2 X 2 T は図4dのようになり、その固有ベクトルは図4eとなり、第2固有ベクトルに一過的なパターンが現れる。さらに特徴量にf₃を加えたn行3列のデータ行列X₃に対しては、第3固有ベクトルにさらに高周波な振動パターンが現れる（図4f、g）。このように、背後に系列が想定されるデータのPCAにおいては、上位の固有ベクトルに時刻（サンプルのインデックス）に対し振動するようなパターンの固有ベクトルが自然と生じることがわかる。

図４：背後に単純な系列データが存在する場合のPCA

a．正規化後のf₁(t), f₂(t), f₃(t)。b、d、f． X 1 X 1 T , X 2 X 2 T , X 3 X 3 T の可視化。c、e、g． X 1 X 1 T , X 2 X 2 T , X 3 X 3 T の固有ベクトル。

次に、これまでの固有ベクトルに振動パターンが現れる理由を、パスグラフのLEとの類似性から議論する。ここで、XX^Tというのはデータの内積を意味しており、つまりその行列の要素は⟨x_i, x_j⟩＝|x_i||x_j|cosθ_{i, j}として計算される。ここで|x_i|などを無視すると、これはcos類似度を表している。ここで、PCAはXX^Tの固有値が大きいものから固有ベクトルを探索し、LEとは固有値が小さいものから固有ベクトルを探索する問題であることから、負の符号をつけた－XX^Tを考えるとする。このとき、各特徴量は標準化されていることから、－XX^Tの行および列の和は0である。ここで－XX^Tから対角成分を抜き出した行列をD、対角成分以外を抜き出した行列を－Wとし、L＝－XX^T＝D－Wを定義すると、Wがノード間の類似度を表し、Dは次数行列というグラフラプラシアンの定義に合っている。したがって、XX^Tの、i、j成分がノードの類似度を表すグラフ構造のLEと同じ問題へと帰着される⁵。

ここで、図5のように10個のサンプルに対し、3つの特徴量が計測されており、それぞれの特徴量が図5aの表の各列のようにサンプルのインデックスに対し変動するデータを考えてみよう。ただし、簡単のためにそれぞれの特徴量の絶対値は無視し、符号のみを考えるとする。これは、PCAにおいては各特徴量が標準化されているため、サンプルのインデックスに対し単調増加あるいは単調減少する特徴量がある場合、平均値を境にあるサンプルのインデックスの前後でその特徴量の符号が変わることになり、その符号のみを抜き出したトイデータといえる。そして、図5aのデータに対しサンプル間のcos類似度をもとにサンプル間のグラフを構築すると、それは図5bのようなグラフに該当する。そして図5bのグラフとパスグラフの類似性から、そのグラフラプラシアンの固有ベクトルには自ずと三角関数に近い形状が現れることが予想できる（図5c）。また、図5bのグラフを見て分かる通り、図5aではサンプルサイズは10であるが、このデータから構築されるサンプル間の類似性グラフとは、実質的にはノードが4つのパスグラフ様の構造に該当する。したがって、結果として3つの三角関数様の固有ベクトルが現れると予想される。ここで、図4cでは第2固有ベクトル以降に周期的な構造は現れていないが、 X 1 X 1 T のデータ構造はおおよそノードが2つのパスグラフ様になっており、 X 2 X 2 T ではノードが3つのパスグラフ様になっているからだと予想できる。

図５：PCAにおける馬蹄効果のメカニズムの概念図

a．10行3列のトイデータの概念図。b．トイデータに対応するcos類似度に基づくグラフ。黒と赤の辺はそれぞれ正と負の類似度を表し、辺の太さは類似度の大きさを表す。c．そのグラフのLEの結果の概念図。

最後に、パスグラフ様の構造のノード数が増えるような状況で、高周波の振動パターンが現れるかを実験的に確認してみる。まず、シグモイド関数の変曲点を調整し、図6a、cのようにサンプルのインデックスに対し変動する4変量の100行4列のデータ行列X₄を作成すると、その X 4 X 4 T は図6eのようになり、固有ベクトルは図6gとなり、第4固有ベクトルまで三角関数様の形状が現れる。さらに、シグモイド関数の変曲点を一つずつずらして図6b、dのような99変量の100行99列のデータ行列X₅を作成すると、その X 5 X 5 T は図6fのようになり、固有ベクトルは図6hとなる。実際には、6番目以降の固有ベクトルにはさらに高周波な振動パターンが出現する。以上の結果より、データ構造がいくつのノードのパスグラフ様の構造に該当するかによって、どの程度まで高周波な振動パターンが固有ベクトルに現れるかが変わることが予想される。

図６：パスグラフ様構造のノード数の違いの検証

X₄とX₅の可視化（a、b．標準化前、c、d．標準化後）。 X 4 X 4 T , X 5 X 5 T の可視化（e、f）とその固有ベクトル（g、h）。

以上のように、背後に系列情報を含むデータ行列に対してPCAを行った場合、各々の特徴量が時系列などに対し周期的や一過的な変動を含まなかったとしても、その結果には三角関数と類似したパターンが出現することがわかる。これはこれまでの議論と同様に、パスグラフのLEとの類似性を考えると直感的に理解できるだろう。また、PCAにおいては必ずしも高周波な三角関数が発生しないが、これはL＝－XX^Tが想定する直線状のグラフが、いくつのノードのパスグラフに類似するかによって決まると考えると、ある程度は納得がいくだろう。

5．まとめ

本稿では、系列データに関わる行列分解に関して、固有ベクトルなどで三角関数あるいはそれに類似した振動パターンが現れる原因を、パスグラフのLEとの類似性から議論した。特にハンケル行列の場合、周期は部分系列の長さnに依存しており、パラメータの決め方で変わるため注意が必要である。もちろん、系列データ中に周期的なパターンなどが存在すれば、ハンケル行列のSVDによってそれを捉えることができるなど、その有効性は多様な場面で示されているが、このような特徴があることには注意が必要であろう。また、データ解析や機械学習においては行列分解を信号分離の目的で用いることも多く、パターンを分離しそれぞれのパターンおよびそれに寄与する因子などを解釈することも多い。本稿では、背後に系列が存在し、かつ単調増加あるいは減少する特徴量しか存在しないようなデータに対するPCAにおいても、振動パターンが容易に発生することも紹介した。ここで、例えば1細胞RNA-seqデータに対しPCAを行い低次元空間上で可視化し、第2、3主成分に現れる振動パターンを見て、例えば時系列に従って発現が振動するような遺伝子が存在するといった仮説を考えたくなるが、これは誤りである。このように、これら現象は時としてミスリーディングな結果をもたらすことに注意が必要である。

いずれにせよ、系列データに対しても行列分解は非常に有効な手段であることに変わりはない。ただし、本稿で紹介するような現象も存在することを頭の片隅に置き、行列分解を用いたアプローチをより適切に活用することが求められる。例えば、馬蹄効果を強制的に引き伸ばして可視化する除歪対応分析（DCA）［8］などの方法を活用することも有効だろう。また、本稿では正弦波効果や馬蹄効果などの現象の一端を解説したのみで、その全貌は未だ解明されたとは言えない。そのため、これら課題を解決できる新たな理論などの開発も待たれる。

6．補足説明

6.1　パスグラフの固有ベクトル

ここでは、2.1節で示したパスグラフの固有ベクトルが、実際に固有ベクトルの定義を満たすことを確認する。パスグラフのグラフラプラシアンのk番目の固有ベクトルについて、ノードiにおける値をv_k(i)と書くとする。まず、1＜i＜nにおいては、以下の式変形によって固有値・固有ベクトルの関係式が成り立つことが示される。

( L P v k ) ( i ) = 2 v k ( i ) − v k ( i − 1 ) − v k ( i + 1 ) ( L P v k ) ( i ) = 2 cos ( πk ( i − 1 / 2 ) / n ) − cos ( πk ( i − 1 / 2 − 1 ) / n ) − cos ( πk ( i − 1 / 2 + 1 ) / n ) ( L P v k ) ( i ) = ( 2 − 2 cos ( πk / n ) ) cos ( πk ( i − 1 / 2 ) / n ) ( L P v k ) ( i ) = ( 2 − 2 cos ( πk / n ) ) v k ( i )

また、i＝1においては図1eのように0番目のノードを仮想的に考えると、v_k(0)＝cos(πk(－1/2)/n)＝cos(πk(1/2)/n)＝v_k(1)が成り立つことから、以下の関係式を導くことができる。

( L P v k ) ( 1 ) = v k ( 1 ) − v k ( 2 ) = 2 v k ( 1 ) − v k ( 0 ) − v k ( 2 )

以上のことから、先ほどの場合と同様にして固有値・固有ベクトルの関係式が成り立つことが示される。i＝nにおいても同様に（n＋1）番目のノードを仮想的に考えると同様に証明できる。よって、固有値固有ベクトルの関係が全てのノードで満たされており、パスグラフのグラフラプラシアンの固有ベクトルは周期が2n/kの三角関数であることが証明された。

なお、k∈[0, n－1]についてλ_k＝2－2cos(πk/n)の固有値と固有ベクトルv_kが存在するため、L_Pはn個の重複しない固有値を持つことから、これら固有ベクトルは一意に定まる。

6.2　リンググラフの固有ベクトル

パスグラフでは天下り的に固有値・固有ベクトルの証明を行った。以降では、リンググラフのLEにおける固有値・固有ベクトルを導出してみる。L_Rの固有ベクトルをuとすると、L_Ru＝λuが成り立つことから、以下の関係式が成立する。

u ( i − 1 ) − ( 2 − λ ) u ( i ) + u ( i + 1 ) = 0

ここで、u(i)＝u(2n＋i)という周期境界条件を導入すると、i＝1, 2nの場合も含め上記式で説明できる。したがって、以降では周期境界条件を課した三項間漸化式を解くことで、固有値・固有ベクトルを導出する。

この三項間漸化式の特性方程式r²－(2－λ)r＋1＝0の判別式は、Δ＝λ²－4λである。ここで周期境界条件を満たすには、特性方程式の解が|r|＝1となる必要がある。なぜなら、三項間漸化式の解は一般的に2つの解r₁, r₂を用いて u ( i ) = Ar 1 i + Br 2 i として計算されるが、|r|≠1では指数関数的に発散するか、0に収束するかであり、周期性が成り立たないからである。

したがって、特性方程式の2つの解について|r|＝1である必要があるが、これは解が複素共役な場合、つまり判別式がΔ＜0となる場合に相当する。つまり、0＜λ＜4となる必要がある。ここで、－2＜(2－λ)＜2であることから、2cosθ＝(2－λ)として新たな変数θで書き換えることにする。すると、特性方程式はr²－2(cosθ)r＋1＝0である。

ここで、オイラーの公式に基づくとr²－2(cosθ)r＋1＝(r－e^jθ)(r－e^－^jθ)という因数分解が可能であり、特性方程式の解はr₁＝e^jθとr₂＝e^－^jθとなることがわかる（混乱を避けるため、虚数単位をjとして書く）。

よって、この三項間漸化式の実数解はu(i)＝Ae^jθi＋Be^－^jθiであり、整理するとu(i)＝Ccos(θi)＋Dsin(θi)という形で記述できる。ここで周期境界条件から、cos(θi)＝cos(θ(2n＋1)i)かつsin(θi)＝sin(θ(2n＋1)i)を満たす必要がある。これを解くと、θ＝πk/nが得られる（ただし、－2＜(2－λ)＜2であり、θ＝(0, π)の範囲とするためk∈[1, n－1]である）。元の関係式に代入するとλ＝2－2cos(πk/n)となり、この固有値に対応する固有ベクトルとしてはx(i)＝sin(πki/n)とy(i)＝cos(πki/n)があげられる。

なお、判別式がΔ＝0となる場合もある。λ＝0においてu(i－1)－2u(i)＋u(i＋1)＝0を満たすu(i)は定数であり、固有ベクトルが定数ベクトルであることがわかる。これは先ほどの固有値・固有ベクトルにおいてk＝0を代入した場合に該当し、yのみが固有ベクトルとして意味があることに相当している。一方でλ＝4においてu(i－1)＋2u(i)＋u(i＋1)＝0を満たすのはu(i)＝－u(i＋1)と交互に符号が変わるベクトルであることがわかる。これは先ほどの固有値・固有ベクトルにおいてk＝nを代入した場合に該当し、yのみが固有ベクトルとして意味があることに相当している。

したがって、Δ＝0となる場合も含め、固有値はλ＝2－2cos(πk/n)となり、固有ベクトルはx(i)＝sin(πki/n)とy(i)＝cos(πki/n)となる。

また、これまでの議論を力学系の観点から直感的に考えることもできる。先の関係式u(i－1)－(2－λ)u(i)＋u(i＋1)＝0は、以下の2変数の線形力学系で書き換えられる。

( u ( i ) u ( i + 1 ) ) = ( 0 1 − 1 2 − λ ) ( u ( i − 1 ) u ( i ) )

ここで、写像を表す右辺の行列をAとして、その固有値をrとすると、その特性方程式|A－rI|＝0とはr²－(2－λ)r＋1＝0となり、三項間漸化式の特性方程式と同じものであることがわかる。ここで、この線形力学系の解が発散や0に収束せずにu(1)＝u(2n＋1)などの周期境界条件を満たすのは、この力学系が回転運動をしている場合に限られる。そして線形力学系において回転運動を示すのは、特性方程式が複素解を持つ場合であることが一般的に成り立ち、これまでの三項間漸化式を用いた説明と同様の結論に至る。なお、回転運動をする場合のAとは、位相平面上で楕円の軌道を描くような写像を意味している。

6.3　リンググラフとパスグラフの対応

サイズnの単位行列をIとし、単位行列の行（あるいは列）の順序を反転させた行列をI′＝I_1:_n_,_n_:1とする。すると、リンググラフとパスグラフのグラフラプラシアンに関して、以下の関係式が成り立つ。

( I I ' ) L R ( I I ' ) = 2 L P

ここで、L_Rの固有値λ_kに対応する固有ベクトルx, yから、i番目のノードに対してu_k(i)＝cos(πk/n)y_k(i)＋sin(πk/n)x_k(i)＝cos(πk(i－1/2)/n)として合成してベクトルu_kを考えると、これは固有ベクトルから合成されているためL_Ru_k＝λ_ku_kが成り立つ。また、0＜i≤nについてu_k(i)＝u_k(2n－i＋1)が満たされている。

ここで、u_kの1番目からn番目までの値を抜き出したベクトルをv_kとすると、これまでの性質より以下の関係式が成り立つ。

( I I ' ) L R ( I I ' ) v k = ( I I ' ) L R u k = ( I I ' ) λ k u k = 2 λ k v k

ここまでの結果からL_Pv_k＝λ_kv_kが成り立つことが示せ、これがパスグラフの固有ベクトルであることがわかる。これはつまり、リンググラフの固有ベクトルから合成されたu_kから前半部分を抜き出したものがv_kであり、u_k(i)＝u_k(2n－i＋1)の対応から、パスグラフとは頂点iと2n－i＋1を対応させたリンググラフと関連することがわかる。

また、リンググラフはk∈[1, n－1]について固有値λ_k＝2－2cos(πk/n)と固有値が重複する2つの固有ベクトルをもち、そのi番ノードの値はそれぞれx_k(i)＝sin(πki/n)とy_k(i)＝cos(πki/n)となる。k＝0およびk＝nにおいても同様に固有値λ_k＝2－2cos(πk/n)となるが、x_k(i)＝0となり、意味のある固有ベクトルはy_kのみである。したがって、n－1個の重複度2の固有値と2つの重複度1の固有値の、計n個の異なる固有値を持つ。そしてL_Rのn個の異なる固有値とは、L_Pのn個の重複しない固有値と一対一対応する。

6.4　JJ^Tに関して

Jの列ベクトルは独立であり、1つの列ベクトルを扱っても一般性は失わないため、以降では1つの部分系列xを考えるとする。ここで、x_iは以下で書き直すことができる。

x i = ∑ k = 1 i ϵ k

したがって、この部分系列の平均値は以下で与えられる。

μ = 1 n ∑ i = 1 n x i = 1 n ∑ k = 1 n ( n + 1 − k ) ϵ k

よって、xから平均値を引いたものをx′とすると、以下のようになる。

x ' i = x i − μ x ' i = ∑ k = 1 i ϵ k − 1 n ∑ k = 1 n ( n + 1 − k ) ϵ k x ' i = ∑ k = 1 i k − 1 n ϵ k − ∑ k = i + 1 n n + 1 − k n ϵ k

ここで、Jのk番目の列ベクトルに対し同様の操作を行ったものをx′^（^k^）としたとき、 JJ T i , j = ∑ k = 1 m x ' i ( k ) x ' j ( k ) となることから、以降では期待値E(x′_ix′_j) (i≤j)を求めることにする。なお、E(ϵ_iϵ_j)＝δ_ijより、i＝jの項のみが残り、以下の様に計算できる。

E ( x ' i x ' j ) = ∑ k = 1 i ( k − 1 ) 2 n 2 − ∑ k = ( i + 1 ) j ( k − 1 ) ( n + 1 − k ) n 2 + ∑ k = ( j + 1 ) n ( n + 1 − k ) 2 n 2 E ( x ' i x ' j ) = 3 ( i − 1 / 2 ) 2 + 3 ( j − n − 1 / 2 ) 2 − 1 / 2 − n 2 6 n

以上より、厳密な値ではないが、スケールや定数を無視すると、mが十分に大きい場合は大数の法則よりJJ^Tはi行j列成分(i≤j)がおおよそ(i－1/2)²＋(j－n－1/2)²となる行列と同様の構造を持つことがわかる。

6.5　KK^Tに関して

x i = ∑ k = 1 i ϵ k に対し、先と同様に期待値E(x_i x_j) (i≤j)を求めると、k≠lにおいてE[ϵ_k]E[ϵ_l]＝0であるから、残るのは共通するk≤iまでの E [ ϵ k 2 ] である。したがって、期待値は以下の通りである。

E ( x i x j ) = ∑ k = 1 i 1 = i

よって、この場合のKK^Tとは、本文中に示す行列Bと同様の構造を持つことがわかる。

次に、BとB₂の関係に関して議論する。Bにおいて固有値λ_kに対し固有ベクトルv_kを持つとすると、それに対しB₂はu_k＝(－v_k′, v_k)を固有ベクトルに持つ。ここで、v_k′とはベクトルを逆順に並べたものであり、v_k′＝I′v_kで与えられる。補足説明6.2で用いたものと類似の形で、以下の関係式が成り立つことが示せる。

( − I ' I ) B ( − I ' I ) = B 2

そして、以下の式変形が成り立つ。

( − I ' I ) B ( − I ' I ) u k = ( − I ' I ) 2 Bv k = 2 λ k ( − I ' I ) v k = 2 λ k u k

以上のことから、B₂u_k＝2λ_ku_kが成り立つことを意味しており、u_kはB₂の固有ベクトルであり、その固有値は2λ_kである。したがって、本文のようにBを基に構築したB₂の固有ベクトルがわかればBの固有ベクトルがわかることになる。

なお、パスグラフのLEとの類似性からu_k(i)＝cos(πk(i－1/2)/2n)のような形となると期待されるが、このとき、0＜i≤nについて－u_k(i)＝u_k(2n－i＋1)が成り立ち、上述の必要条件が満たされていることがわかる。

謝辞

本記事の執筆にあたり、長崎大学情報データ科学部の植木優夫教授および加葉田雄太郎助教には数学的な内容などに関する議論やアドバイスをいただいた。また、九州大学大学院医学研究院の前原一満准教授にはリンググラフのLEなどに関して提案をいただいた。この場を借りて深謝の意を表する。

脚注

1　ただし、どのような連結グラフについても、ラプラシアンの最小固有値は重複しない0であり、対応する固有ベクトルは常に定数となるので、第1固有ベクトルは意味を持たないので捨てる。

2　補足説明6.1で示す通りn個の固有値は重複がないため、これら固有ベクトルは一意に定まる。

3　U_,－1は左特異ベクトルから第1成分を取り除いたものを意味する。

4　ただし、このような類似度のグラフでは負の重みの辺も存在し、一般的なLEとは少し異なることに注意が必要である。

5　なお、この場合も負の辺を持つことに注意。

References

[1] 露崎弘毅, 松本拡高. (2025). 行列・テンソル分解によるヘテロバイオデータ統合解析の数理―第6回系列データ―. JSBi Bioinformatics Review, 6(1), 11-28.
[2] Trapnell, C., Cacchiarelli, D., Grimsby, J., Pokharel, P., Li, S., Morse, M., ... & Rinn, J.L. (2014). The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells. Nature biotechnology, 32(4), 381-386. https://doi.org/10.1038/nbt.2859
[3] Shinn, M. (2023). Phantom oscillations in principal component analysis. Proceedings of the National Academy of Sciences, 120(48), e2311420120. https://doi.org/10.1073/pnas.2311420120
[4] Keogh, E., & Lin, J. (2005). Clustering of time-series subsequences is meaningless: implications for previous and future research. Knowledge and information systems, 8, 154-177. https://doi.org/10.1007/s10115-004-0172-7
[5] Idé, T. (2006, September). Why does subsequence time-series clustering produce sine waves?. In european conference on principles of data mining and knowledge discovery (pp. 211-222). Berlin, Heidelberg: Springer Berlin Heidelberg. https://doi.org/10.1007/11871637_23
[6] 露崎弘毅. (2022). 行列・テンソル分解によるヘテロバイオデータ統合解析の数理―第4回質的データ，距離，グラフ―．JSBi Bioinformatics Review, 3(2), 33-46. https://doi.org/10.11234/jsbibr.2022.2
[7] Spielman, D. Rings, paths, and Cayley graphs, Spectral Graph Theory, 2014.
[8] Diaconis, P., Goel, S., & Holmes, S. (2008). Horseshoes in multidimensional scaling and local kernel methods.

著者略歴

松本拡高
2016年、東京大学大学院新領域創成科学研究科博士課程終了。博士（科学）。理化学研究所生命機能科学研究センターや革新知能統合研究センターでの研究員を経て、現在は長崎大学情報データ科学部に在籍し、多様な生命現象のデータの解析を行なっている。

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

Corresponding author

Register with J-STAGE for free!