2025 年 6 巻 2 号 p. 60-68
生命科学分野で取得されるデータ集合は、雑多(ヘテロ)な構造になり、ヘテロなデータ構造を扱える理論的な枠組みがもとめられている。本連載では、汎用的なヘテロバイオデータの解析手法である行列・テンソル分解を紹介していく。第7回では、行・列方向でともに和が1となる行列を生成する行列バランス化問題や、近年注目されている最適輸送問題と行列・テンソル分解の関係性について議論する。
これまでの連載[1, 2, 3, 4, 5, 6]で扱った行列分解手法のうち、第4回[4]の分割表の対応分析では(周辺または同時)確率分布、第5回[5]のランダムウォーク手法では遷移確率行列、さらに第5回[5]のラプラシアン固有マップやスペクトラルクラスタリングではランダムウォーク正規化グラフラプラシアンなど、非負値の頻度情報を正規化し、確率値として扱う場面が多く登場した。
遷移確率行列は、さらに行ごとの和が全て1の右確率行列(Right Stochastic Matrix)と、列ごとの和が全て1となる左確率行列(Left Stochastic Matrix)の2種類がある。ここではさらに、行・列の和がともに1となる行列、すなわち二重確率行列(Doubly Stochastic Matrix)を導入し、その導出法である行列バランス化(Matrix Balancing)または行列スケーリング(Matrix Scaling)と呼ばれる手法を紹介する[6]。
あるn×mの非負値行列Xの右確率行列Prowをもとめることは、対角要素にXの行ごとの和の逆数を格納したn×nの対角行列
| (1) |
また、同様にXの左確率行列Pcolをもとめることは、対角要素にXの列ごとの和の逆数を格納したm×mの対角行列
| (2) |
これらに対して、Xの二重確率行列Pdoubleをもとめる際には、n×nの対角行列D1とm×mのD2を用いて、以下のように行にも列にも正規化を適用する必要がある。
| (3) |
ただし、これら2つの正規化ステップは互いに依存するため、
この式(3)の行列バランス化を実現するためのアルゴリズムとして知られているのがSinkhorn-Knopp(SK)アルゴリズムである[7]。まず推定したいPdoubleは、行ごとの和も列ごとの和も全て1であるため、以下が成り立つ。
| (4) |
| (5) |
ただし、1nと1mは各々長さn、mの全要素が1のベクトルである。ここでr=diag(D1)、c=diag(D2)とする(diagは行列の対角要素をベクトルとして返す関数)。これにより、ベクトル表記で式(4)、(5)を各々以下のように書き換えることができる。
| (6) |
| (7) |
なお、x=f(x)を満たすxを、適当な初期値x0からxt+1=f(xt), (t=0, 1, 2, …)というように計算することで逐次的に推定するアルゴリズムを不動点反復法(Fixed Point Iteration)というが、SKアルゴリズムではこの手法を式(6)、(7)に適用することで、以下のようにrとcを逐次的に推定する。
| (8) |
なおrの初期値r0としては単位ベクトルが利用される。
行列バランス化は二重確率行列を得る以外のタスクにも適用でき、式(4)、(5)における1nを任意の正のベクトルo、sにおきかえ、
| (9) |
| (10) |
として、これらを以下のように解くこともできる。
| (11) |
この拡張された定式化を、通常の行列バランス化と区別するために、あえて“一般化”行列スケーリングと呼ぶ場合もある。
行列バランス化は、染色体間相互作用を計測するHi-Cのデータであるコンタクトマップの正規化手法として広く利用されている[7, 8, 9]。また、組織内のある遺伝子の3D空間発現パターンを、x, y, z軸方向でのスライス画像から復元するRNAトモグラフィ解析では、復元の際に用いられるIterative Proportional Fitting(IPF)[10, 11, 12]がSKアルゴリズムと等価の計算となっている。ただし、3Dデータであるため、推定する3D空間発現パターンχ(3階テンソル)に対して、x, y, z軸ごとの和ベクトルo, s, qが事前に与えられていることを想定している。このような問題設定は、行列バランス化のテンソル拡張であるため、テンソルバランス化(Tensor Balancing[13])とも呼ばれる。
ここでは長さnのベクトルuの各要素を長さmのベクトルvの要素に振り分ける操作を「輸送」と定義する(図1)。ただしu, vはいずれも非負値の確率ベクトル、すなわち全要素の和が1であるとする。さらに、u-v間の全ての要素間の組み合わせに対して、輸送のコストが与えられているものとする。この時に全体の輸送コストを最小化するような輸送計画を求める問題が最適輸送(Optimal Transport;OT)問題である[14, 15, 16, 17, 18]。

確率分布uの各要素をvの要素に振り分ける最適輸送において、どのuの要素の確率値をどのvの要素に輸送するのか1対1の対応関係(最適輸送マップ)をもとめるモンジュ型輸送と、どのuの要素の確率値をどのvの要素にどのくらいの割合で輸送するのか(最適輸送プラン)をもとめるカントロビッチ型輸送がある。ただし、ノードの大きさはその要素の確率の値の大きさ、エッジの太さはその確率値を輸送する割合を示している。
最初に提案されたOTは、モンジュによるものであり、uのある要素の確率値は全てvのどこか1要素にのみ輸送することを仮定した割り当て問題(assignment problem)として定式化されていた(図1)。しかしながら、輸送方法が一意に定まらず、解が複数ある、または解が存在しないようなケースがあったため、「uのある要素の確率値を分割して、vの別々の要素に輸送しても良い」というように問題設定の連続緩和(第4回[4]グラフラプラシアン参照)が、その後カントロビッチによってなされ、現在ではこの定式化が主流となっている(図1)。本稿でも以降はカントロビッチの定式化を採用する。
OTを解くためのアルゴリズムは連続最適化に基づくものも、離散最適化に基づくものも多数存在する[15]。ここでそれら全てを紹介することはできないが、勾配を陽に計算するため他のモデルと組み合わせやすい、行列計算として定式化されるためGPUで並列処理しやすく大規模データに適用可能といった理由で、上記の行列バランス化の稿で用いたSKアルゴリズムが広く利用されているため、ここでも紹介する。まず、長さnの確率ベクトルaを長さmの確率ベクトルbに変換するカントロビッチ型のOTの目的関数は、以下の通りである。
| (12) |
ただし、〈,〉は2行列の内積(アダマール積の和)である。ここで、n×mの非負値行列Pは、輸送行列または輸送プラン行列と呼ばれるもので、最適化によって推定されるものである。n×mの行列Cはコスト行列と呼ばれるもので、a-b間の全ての要素の組み合わせで、輸送のコスト(例:ユークリッド距離)をデータから事前に求めたものである。aはPの行ごとの和、bはPの列ごとの和となっており、ともに全ての要素の和は1である。これらは質量保存制約と呼ばれ、特に事前情報が無い場合は一様分布を仮定して、各要素をaは1/n、bは1/mとして設定することが多い。
式(12)のPをもとめる方法としては線形計画法がまず挙げられるが、計算コストが高い、数値的に不安定、解が一意に求まる保証がないといったデメリットがある。そのため、ここでは式(12)にPのエントロピー項H(P)=-∑ijPij(log(Pij)-1)を正則化として加える。
| (13) |
εはPのエントロピーの大きさを調整するハイパーパラメーターであり、値が小さいとPはより疎になり、モンジュの割り当て問題のような1対1対応型の解が選ばれやすくなる。一方、εの値が大きくなるにつれてPは密行列となり、aのある要素における確率値は、bの様々な要素に輸送されるようになる。なお、Pのエントロピーが最大となるのは、P=abTが成立する時、すなわちaとbが互いに独立な場合である[14]。これは第4回[4]の対応分析では、行変数と列変数がなるべく非独立になるように定式化することで、互いに関連する行変数スコアと列変数スコアをもとめていたこととは対照的である。
この目的関数をラグランジュ未定乗数法で以下のように最適化する。
| (14) |
ただし、長さnのベクトルfと長さmのベクトルgはラグランジュ未定乗数である。これをPに関して偏微分すると、
| (15) |
となるため、
| (16) |
となる。これは制約条件P 1m=a、PT 1n=bでの、ガウスカーネルで得られたグラム行列Kをuとvで行列バランス化(一般化行列スケーリング)していることに他ならない。そのため、式(11)でのX, s, t, r, cをK, a, b, u, vに置き換えた上でuとvを交互に推定し、収束後にuKvTを計算することで最適なPを得ることができる。
なお、コスト行列Cの計算には通常ユークリッド距離のp乗が利用されるが、この時の式(12)の目的関数の値の1/p乗をa-b間のp-Wasserstein距離という(図2a)。通常pの値は1または2がよく利用される。統計学の分野では、これはEarth Mover's Distance(EMD)とも呼ばれる[19]。これは元々モンジュによってOTが定式化された際に、土(Earth)の山を運ぶ問題として紹介されたことに由来する。余談だが、自然言語処理分野ではWord Mover's Distance[20, 21]、推薦システム分野ではPreference Mover's Distance[22]、シングルセルオミックス分野ではGene Mover's Distance[23]のように、OT分野では、EMDをもじった名前の手法名や論文のタイトルが多い。

上記の通常のOT(Wasserstein OT)を異なる次元にいるデータ間でも行えるように拡張した手法であるGromov-Wasserstein型のOT(GW-OT)も紹介する。ここではdu次元空間にいるn個のデータを、dv(≠du)次元空間にいるm個のデータに輸送することを考える(図2b)。このように互いに行も列も共有しないデータ行列同士をDiagonalという(第3回[3])。または例え次元のサイズが偶然同じであっても、軸の対応関係が無いデータ同士でも本質的に同じ状況である。ここで問題なのは、異なる次元の空間にいるデータ同士では距離を計算できないため、Wasserstein OTが適用できないことである。GW-OTではこういった場合に、同じ次元の空間にいるデータ同士の類似度行列Ss(n×n)、St(m×m)を利用して、以下のような最適化を考える。
| (17) |
ただし、
OTの各種機械学習タスクへの適用はここ数年でかなり活発化している。特にOTを簡単に実行できるPythonパッケージPython Optimal Transport(POT)[27]をはじめとするオープンソースなツールが公開されるようになってからは、多くのデータサイエンス分野での適用事例が報告されている(図3)。

また、OTにはGW-OT以外にも様々な拡張モデルが提案されており、例えばa, bの質量保存制約を緩和し、輸送しないデータが一部含まれていることもよしとするUnbalanced OT(またはPartial OT)、データ間の構造(例:木)をドメイン知識として取り入れるStructured OT(例:Hierarchical OT)、GW-OTをさらに発展させ、行だけでなく列も同時に輸送するCO-Optimal Transport(COOT)や、Wasserstein OTとGW-OTの目的関数を重み付きで足し合わせたFused Gromov-Wasserstein-OT(FGW-OT)などが挙げられる[28]。
多くの機械学習モデルでは、データの値とモデルの推定値間の誤差を定義した損失関数を最小化するように定式化されるが、このような損失関数はそのままOTとして再設定可能である(図2c)。この方向性の拡張として、本連載の第1回から第3回[1, 2, 3]で紹介した主成分分析(Principal Component Analysis;PCA)、非負値行列分解(Non-negative Matrix Factorization;NMF)、Fisherの線型判別モデル(Linear Discriminant Analysis;LDA)、非負値テンソル分解(Non-negative Tensor Factorization;NTF)[29, 30, 31, 32, 33]などが挙げられる。従来の損失関数(例:Kullback-Leiblerダイバージェンス)と比較したOTによる定式化のメリットとして、1)モデルが返す「完全に正解ではないが近しい予測値」も扱える、2)Wasserstein距離は距離の公理を満たすため、あらゆるデータ間で総当たりのWasserstein距離行列を計算し、その後の次元圧縮やクラスタリングといったタスクに応用できる、3)輸送コストだけでなく輸送プラン行列自体も様々な問題に利用可能、といった点が挙げられる[14]。
ドメイン適応(Domain Adaptation;DA)というタスクでは、ラベルありデータセットと、ラベルなしデータセットを用意し、ラベルなしデータのラベル予測を行う[11](図2d)。この際に、両データセットの分布形状が大きく異なると、うまく予測できない場合があるためOTが利用される。敵対的生成ネットワーク(Generative Adversarial Networks;GAN)をDAに利用したDA-GANでは、両者の特徴分布をOTで近づけるようにGANを拡張している[14]。
また、Wasserstein距離による確率分布間の重心(Wasserstein重心)の計算方法が提案され、図形モーフィングと呼ばれるタスクに応用されており、2つの図形を与えた時に、ある図形を他の図形に輸送する過程の中間状態(両者の1-λ:λ(0≤λ≤1)の比率に位置する内分点)を滑らかに内挿できることが報告されている[14](図2e)。
OTは生命科学において非常に強力なツールとなっており、対応関係が無いデータ間でOTをすることにより、輸送プラン行列を擬似的な対応関係として利用できる。例えば、動きがある生命現象の動画データの解析や、異なる対象物をとらえた画像同士の比較、基準となるリファレンスデータ画像に、他の画像を重ね合わせて知識データベース化していく画像レジストレーション(Image Registration)といったタスクでは、画像のピクセル同士の直接的な比較は意味をなさず、対象物の並行移動や拡大・縮小、回転、反転、変形などを考慮した比較が必要となってくるが、このようなデータの解析にOTを適用することで、データの分布形状に着目した輸送や、Wasserstein距離によるデータ解析が行える(図2a)。こういった方向性の研究として、神経の画像[34]、血管の画像[35]、2次元電気泳動画像[36]、タンパク質3次元構造の電気顕微鏡画像[37]、脳のMRI画像[38]、液体クロマトグラフィー質量分析法(Liquid Chromatography-Mass Spectrometry;LC-MS)のスペクトルデータ[39]の解析にOTが利用された事例がある。
シングルセルオミックスは、時系列計測や摂動実験など複数条件で計測した場合、計測のたびに細胞を殺してしまうため(第3回[3])、バッチ間の細胞の対応関係は失われてしまうが、バッチ間でOTを適用することで対応関係を後から再構築できるため、それにより、時系列に従う遺伝子発現を滑らかに推定したり[40, 41]、摂動により細胞がどの方向に変化したのかを次元圧縮図上のベクトル場で描画する解析が行われている[42, 43]。
また、対応関係が得られることで、Guilt-by-Associationの原理(第5回[5])に基づき、一方のデータに付随する情報を他方に伝播させやすくなる。それにより、細胞型がアノテーションされたscRNA-seqデータを、別のscRNA-seqデータのアノテーションに利用したり[44]、あるいはアノテーション済みのscRNA-seqデータと空間トランスクリプトームをOTで対応付けることで、空間側の各スポットにおける細胞型比率の推定(cell-type deconvolution)や、未計測遺伝子の空間発現の推定(空間再構築解析)といった応用が可能となっている[45]。
Wasserstein重心のような輸送の途中の状態を推定するアプローチは、時系列データ間の未観測データの値を予測(内挿)するのに有用であり、例えばシングルセルオミックス分野での、分化過程の細胞の途中の状態を調べるのに利用されている[40, 41]。また、反応物から生成物になる化学反応過程での遷移状態を推定する際にOTが使われた例がある[46]。Wasserstein重心を発現変動遺伝子の検出に利用した研究もあり[47]、各群のレプリケイトのWasserstein重心を求めて群間比較に利用することで、ノンパラメトリックに群間で変動した遺伝子を検出できている。また、シングルセルデータに対して複数回クラスタリングした結果をマージして、1つのクラスタリング結果にするアンサンブル学習としてもWasserstein重心が利用されている[48]。
OTの各種拡張モデルを生命科学データに応用した事例としては、Partial OTやUnbalance/Structured OTで1細胞RNA-Seqのデータ統合を行ったもの[49, 50]、COOTで1細胞マルチオミックスデータ統合を行ったものが挙げられる[51]。またそれ以外にも、従来のNMFをOT化したものを1細胞マルチオミックスデータの統合解析に利用したもの[52]、フローサイトメトリーデータにおけるDA[53, 54]などがある。
著者が注目しているのはDiagonalなデータにGW-OTやCOOT、FGW-OT(図2b)を適用する方向性であり、これは今後様々な応用が考えられる。既に1細胞RNA-Seqデータと1細胞ATAC-Seqデータの統合[55]、タンパク質相互作用ネットワーク(Protein-Protein Interaction;PPI)間のグラフアライメント[56]といったタスクに利用されている。また、まだ生命科学における応用例はないが、GW-OTやCOOTのテンソル拡張としてOptimal Tensor Transport(OTT)という手法も存在する[57]。
今回は、SKアルゴリズムが行列バランス化や、最適輸送問題に重要な役割を果たしていることを紹介した。最適輸送問題は、生命科学分野においてこれまで比較することが難しかった形状データや、Diagonalデータの統合解析で成果を上げており、今後様々な問題への適用が期待できる。
次回は、イベント発生の時系列性に関する解析である生存時間解析と行列・テンソル分解の関係性について議論する。
Sinkhorn-Knoppアルゴリズム
・IPFIterative Proportional Fitting
・OTOptimal Transport
・EMDEarth Mover's Distance
・GW-OTGromov-Wasserstein OT
・POTPython Optimal Transport
・COOTCO-Optimal Transport
・FGW-OTFused Gromov-Wasserstein OT
・PCAPrincipal Component Analysis
・NMFNon-negative Matrix Factorization
・LDALinear Discriminant Analysis
・NTFNon-negative Tensor Factorization
・DADomain Adaptation
・GANGenerative Adversarial Networks
・LC-MSLiquid Chromatography-Mass Spectrometry
・PPIProtein-Protein Interaction
・OTTOptimal Tensor Transport
![]() |
露崎 弘毅 2015年、東京理科大学生命創薬科学科博士後期課程終了。博士(薬科学)。同年より、理化学研究所(所属3)に在籍し、シングルセルオミックスのデータ解析や解析ツール開発に従事。現在は千葉大学(所属1, 2)で医療データの解析を行なっている。パッケージングに特化したハッカソンBio “Pack” athonを主催。 ホームページ:https://researchmap.jp/kokitsuyuzaki |