ケィオスの時系列解析メモランダム

時系列解析、生体情報解析などをやわらかく語ります

【統計・時系列解析の基礎としての行列１】行列表現・線形代数が何の役に立つのか？

線形代数 統計 時系列解析 行列 ベクトル

近年、統計解析や時系列解析は、多くの分野で「道具」として広く使われるようになっています。その一方で、そのような手法がどのような構造を仮定し、何を計算しているのかを十分に理解しないまま使っている学生や実務者が増えているとも感じます。統計や時系列解析を、手順どおりに操作するだけのブラックボックスとして扱ってしまうと、結果の妥当性を判断したり、条件を変えたときの振る舞いを正しく理解したりすることが難しくなります。
　今回は、統計・時系列解析の背後にある行列（線形代数）の考え方を理解する方向へ、皆さんを導くことを目指します。そのための共通の言葉として、ベクトルと行列による表現に注目します。
　行列表現は、単に式を簡潔に書くための記法ではありません。多変量データの関係性を一つのまとまりとして扱い、推定、分解、変換といった操作の意味を見通しよく理解するための枠組みです。数式の読み方がわかるようになると、これまで別々の方法だと考えていた手法の間に共通性が見えてきて、全体の見通しが一気に良くなり、理解を効率化することができます。
　例えば、最小二乗法は、「連立一次方程式を解く問題」としてだけでなく、「観測データをある部分空間へ射影する操作」として理解できます。また、共分散行列の固有値分解は、主成分分析（Principal Component Analysis, PCA）の基礎となり、独立成分分析（Independent Component Analysis, ICA）や Second-Order Blind Identification（SOBI）、さらには多変量時系列解析へと自然につながっていきます。 時系列解析においても、自己回帰モデル（Autoregressive Model, AR）は、行列表現によって回帰問題として統一的に扱うことができます。
　以下では、統計および時系列解析において、特に重要で、多くの場面に登場する行列計算を順に取り上げます。今回は概説とし、各節の導出と具体例は別記事で扱います。行列計算を単なる「計算手順」として覚えるのではなく、式の意味を読み取り、手法どうしのつながりを理解するための道具として使えるようになることを、本資料の目標としています。

1. 最小二乗：統計の頻出問題

　統計解析において最も基本的で、かつ頻繁に使用されるのが最小二乗法です。一般の最小二乗法では、次の線形モデルを考えます。

\displaystyle{ \mathbf{y} = \mathbf{X}\, \boldsymbol{\beta} + \boldsymbol{\varepsilon} }

ここで、

  • \mathbf{y} は観測データを縦に並べた N 次元ベクトル
  • \mathbf{X} は説明変数からなる デザイン行列
  • \boldsymbol{\beta} は未知の パラメタベクトル
  • \boldsymbol{\varepsilon} は 誤差ベクトル

を表します。デザイン行列（design matrix）とは各観測点での説明変数の値を並べたものです。
　最小二乗法では、観測データ \mathbf{y} とモデルによる予測値 \mathbf{X}\boldsymbol{\beta} の差、すなわち残差ベクトル

\displaystyle{ \mathbf{y}-\mathbf{X}\, \boldsymbol{\beta} }

の2乗和を最小にするように、パラメタ \boldsymbol{\beta} を決めます。つまり、

\displaystyle{ \lVert\mathbf{y}-\mathbf{X}\boldsymbol{\beta}\rVert^2 = (\mathbf{y}-\mathbf{X}\boldsymbol{\beta})^\top (\mathbf{y}-\mathbf{X}\boldsymbol{\beta}) }

を最小化することが、最小二乗法の基本的な考え方です。

■ 観測 N 点、説明変数 2 個の例

　このような表現に慣れていない人にイメージを持ってもらうために、観測が N 点、説明変数が 2 個の場合を具体的に書いてみます。
　観測が N 点あるとき、観測ベクトルは

\displaystyle{ \mathbf{y} = \begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_N \end{pmatrix} }

という [N] 次元ベクトルです。
　説明変数が 2 個の場合、パラメタベクトルは

\displaystyle{ \boldsymbol{\beta} = \begin{pmatrix} \beta_1 \\ \beta_2 \end{pmatrix} }

という 2 次元ベクトルになります。
　各観測点での説明変数の値を並べると、デザイン行列は

\displaystyle{ \mathbf{X} = \begin{pmatrix} x_{11} & x_{12} \\ x_{21} & x_{22} \\ \vdots & \vdots \\ x_{N1} & x_{N2} \end{pmatrix} }

となります。
　誤差は観測点ごとに 1 つずつあるため、

\displaystyle{ \boldsymbol{\varepsilon} = \begin{pmatrix} \varepsilon_1 \\ \varepsilon_2 \\ \vdots \\ \varepsilon_N \end{pmatrix} }

という N 次元ベクトルになります。
　以上を用いると、線形モデルは

\displaystyle{ \begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_N \end{pmatrix} = \begin{pmatrix} x_{11} & x_{12} \\ x_{21} & x_{22} \\ \vdots & \vdots \\ x_{N1} & x_{N2} \end{pmatrix} \begin{pmatrix} \beta_1 \\ \beta_2 \end{pmatrix} + \begin{pmatrix} \varepsilon_1 \\ \varepsilon_2 \\ \vdots \\ \varepsilon_N \end{pmatrix} }

と書けます。右辺を展開してみると、

\displaystyle{ \begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_N \end{pmatrix} = \begin{pmatrix} x_{11}\beta_1 + x_{12}\beta_2 + \varepsilon_1 \\ x_{21}\beta_1 + x_{22}\beta_2 + \varepsilon_2 \\ \vdots \\ x_{N1}\beta_1 + x_{N2}\beta_2 + \varepsilon_N \end{pmatrix} }

となります。そして、各成分を比較すると、

\displaystyle{ y_i = x_{i1}\beta_1 + x_{i2}\beta_2 + \varepsilon_i \qquad (i=1,\dots,N) }

という、回帰式になっています。

■ パラメタ \boldsymbol{\beta} を求めると

　抽象的なベクトル・行列表現の威力を感じるために、パラメタ \boldsymbol{\beta} を最小二乗法により求めると、一般に、

\displaystyle{ \boldsymbol{\beta} = (\mathbf{X}^\top \mathbf{X})^{-1}\mathbf{X}^\top \mathbf{y} }

となります（ただし、\mathbf{X}^\top\, \mathbf{X} が正則である場合に限る）。この式は、観測が N 点であっても、説明変数の数が増えても、同じ形のまま成り立ちます。この結果を正しく理解するためには、式の展開、ベクトル・行列の微分、逆行列に関する基本的な知識が必要になります。今回は詳細は説明しません。

2. 射影行列：回帰に幾何的視点を導入

　前節では，最小二乗法によってパラメタ \boldsymbol{\beta} が

\displaystyle{ \boldsymbol{\beta} = (\mathbf{X}^\top \mathbf{X})^{-1}\mathbf{X}^\top \mathbf{y} }

と求まることを紹介しました。
　この式を計算公式として使うこともできますが、行列の意味を理解すると、回帰が何をしているのかを幾何学的に理解できるようになります。その鍵となるのが、射影行列（projection matrix）です。

■ 回帰による予測値の行列表現

　最小二乗法で得られたパラメタ \boldsymbol{\beta} を用いると、観測データに対するモデルの予測値は

\displaystyle{ \hat{\mathbf{y}} = \mathbf{X}\boldsymbol{\beta} }

と書けます。ここに、前節で得られた \boldsymbol{\beta} を代入すると、

\displaystyle{ \hat{\mathbf{y}} = \mathbf{X}(\mathbf{X}^\top \mathbf{X})^{-1}\mathbf{X}^\top \mathbf{y} }

となります。ここで、

\displaystyle{ \mathbf{P} = \mathbf{X}(\mathbf{X}^\top \mathbf{X})^{-1}\mathbf{X}^\top }

とおくと、

\displaystyle{ \hat{\mathbf{y}} = \mathbf{P}\mathbf{y} }

と非常に簡潔に書けます。 この行列 \mathbf{P} を 射影行列と呼びます。

■ 射影行列の意味

　射影行列 \mathbf{P} は、観測ベクトル \mathbf{y} を、デザイン行列 \mathbf{X} の列ベクトルが張る部分空間へ写す操作を表しています。つまり、回帰分析とは、「観測データ \mathbf{y} を、説明変数で表現できる部分空間へ直交射影すること」として理解できます。
　この見方に立つと、回帰は「式を当てはめる操作」ではなく、高次元空間における幾何学的な写像であることが分かります。

■ 射影行列の基本的な性質

　射影行列 \mathbf{P} には、次の重要な性質があります。

\displaystyle{ \mathbf{P}^\top = \mathbf{P} \qquad \mathbf{P}^2 = \mathbf{P} }

すなわち、\mathbf{P} は 対称行列であり、同じ射影を 2 回行っても結果は変わらない（冪等性）という特徴を持っています。
　この性質から，予測値 \hat{\mathbf{y}} と残差ベクトル

\displaystyle{ \mathbf{e} = \mathbf{y}-\hat{\mathbf{y}} = (\mathbf{I}-\mathbf{P})\mathbf{y} }

が、互いに直交することも分かります。
　以上をまとめると，観測ベクトル \mathbf{y} は，

\displaystyle{ \mathbf{y} = \hat{\mathbf{y}} + \mathbf{e} = \mathbf{P}\mathbf{y} + (\mathbf{I}-\mathbf{P})\mathbf{y} }

と分解できます。これは、

  • \hat{\mathbf{y}}：モデルで説明できる成分
  • \mathbf{e}：モデルでは説明できない成分

への 直交分解を意味しています。
　このように、射影行列を用いることで、最小二乗回帰の意味を、代数的にも幾何学的にも統一的に理解できるようになります。

3. 共分散行列（対称・正定値）：ばらつき、不確かさの記述

　前節では、回帰分析が射影行列によって記述できることを見ました。次に重要になるのが、誤差やばらつき、不確かさそのものをどう表現するかという問題です。この役割を担うのが、共分散行列（covariance matrix）です。

■ 共分散行列の定義

　平均が 0 の多変量確率変数ベクトル \mathbf{x} に対して、共分散行列は

\displaystyle{ \boldsymbol{\Sigma} = \mathrm{E}\bigl[\mathbf{x}\, \mathbf{x}^\top\bigr] }

と定義されます。成分で書けば、

\displaystyle{ \Sigma_{ij} = \mathrm{E}[x_i\, x_j] }

であり、対角成分は各成分の分散、非対角成分は成分間の共分散を表します。

■ 共分散行列の基本的な性質

　共分散行列 (\boldsymbol{\Sigma}) は、定義から次の重要な性質を持ちます。

\displaystyle{ \boldsymbol{\Sigma}^\top = \boldsymbol{\Sigma} }

すなわち 対称行列です。さらに、任意のベクトル \mathbf{v} に対して，

\displaystyle{ \mathbf{v}^\top \boldsymbol{\Sigma} \mathbf{v} \ge 0 }

が成り立つため、\boldsymbol{\Sigma} は半正定値であり、退化がなければ正定値になります。
　この性質は、「ばらつき」や「不確かさ」が負にならないことを、行列として表現していると理解できます。

■ 回帰と共分散行列の関係

　線形モデル

\displaystyle{ \mathbf{y} = \mathbf{X}\, \boldsymbol{\beta} + \boldsymbol{\varepsilon} }

において、誤差ベクトル (\boldsymbol{\varepsilon}) の共分散行列を

\displaystyle{ \mathrm{E}[\boldsymbol{\varepsilon}\, \boldsymbol{\varepsilon}^\top] = \boldsymbol{\Sigma} }

とすると、\boldsymbol{\Sigma} は「観測誤差の大きさ」と「誤差どうしの相関構造」をまとめて表す量になります。
　誤差が独立で分散が等しい場合には、

\displaystyle{ \boldsymbol{\Sigma} = \sigma^2 \mathbf{I} }

となり、これは通常の最小二乗法の仮定に対応します。一方、誤差に相関がある場合には、共分散行列は一般の対称正定値行列となり、重み付き最小二乗法や一般化最小二乗法へと自然につながっていきます。

■ 共分散行列は「構造の言語」

　このように、共分散行列は単なる分散の集まりではなく、

  • ばらつきの大きさ
  • 成分間の依存関係
  • 不確かさの方向性

を同時に表現するための言語です。

4. 固有値分解と逆平方根：共分散構造から PCA と whitening へ

　前節では、共分散行列が、ばらつきや不確かさ、成分間の依存構造を表す行列であることを見ました。次に重要になるのが，この共分散行列をどのように分解し，どのように利用するかという点です。その中心となる操作が，固有値分解と逆平方根です。

■ 共分散行列の固有値分解

　共分散行列 \boldsymbol{\Sigma} は対称正定値行列であるため，次のように固有値分解できます。

\displaystyle{ \boldsymbol{\Sigma} = \mathbf{V}\, \boldsymbol{\Lambda}\, \mathbf{V}^\top }

ここで、\mathbf{V} は固有ベクトルを列にもつ直交行列、\boldsymbol{\Lambda} は固有値を対角成分にもつ対角行列です。固有ベクトルは，ばらつきが現れる「方向」を表し，固有値はその方向に沿った「大きさ」を表します。

■ 主成分分析との関係

　この固有値分解は、主成分分析（Principal Component Analysis, PCA）の基礎そのものです。データを固有ベクトル方向に回転すると、

  • 成分間の相関が消える
  • 分散が大きい順に成分を並べられる

という性質が得られます。すなわち、PCAとは、共分散行列を固有値分解し、ばらつきの向きと大きさを分離して捉える手法だと理解できます。

■ 共分散行列の逆平方根と whitening

　固有値分解を用いると、共分散行列の逆平方根を次のように定義できます。

\displaystyle{ \boldsymbol{\Sigma}^{-1/2} = \mathbf{V}\, \boldsymbol{\Lambda}^{-1/2}\, \mathbf{V}^\top }

ここで、\boldsymbol{\Lambda}^{-1/2} は、各固有値の平方根の逆数を対角成分にもつ行列です。

　この逆平方根を用いてデータ \mathbf{x} を変換すると，

\displaystyle{ \mathbf{z} = \boldsymbol{\Sigma}^{-1/2}\, \mathbf{x} }

となり、変換後のデータ \mathbf{z} の共分散行列は単位行列になります。この操作を whitening（白色化） と呼びます。

■ PCA と ICA，SOBI へのつながり

白色化（whitening）によって、

  • 分散がすべて 1 にそろう
  • 成分間の相関が除去される

ため、残る自由度は直交変換のみになります。この性質が、独立成分分析（ICA） や SOBI における分離の出発点になります。
　このように、固有値分解と逆平方根は、

  • 共分散行列の理解
  • PCA による次元削減
  • whitening を経由した信号分離

を一つの流れとして結びつける、重要な行列操作です。

5. 特異値分解（SVD）：一般行列と数値安定性の要

　前節では、共分散行列のような対称正定値行列に対して、固有値分解が役立つ道具であることを見ました。しかし、統計や時系列解析で現れる行列は、必ずしも対称行列とは限りません。例えば、デザイン行列 \mathbf{X} や観測データを並べたデータ行列は、一般に非正方行列であり、そのままでは固有値分解を適用できません。
　そのような一般の行列に対して，中心的な役割を果たすのが，特異値分解（Singular Value Decomposition: SVD）です。

■ 特異値分解の定義

　任意の N\times M 行列 \mathbf{X} は、次のように分解できます。

\displaystyle{ \mathbf{X} = \mathbf{U}\, \mathbf{S}\, \mathbf{V}^\top }

ここで、

  • \mathbf{U} は列が互いに直交する N\times N 行列
  • \mathbf{V} は列が互いに直交する M\times M 行列
  • \mathbf{S} は非負の特異値を対角成分にもつ対角行列

です。
　特異値は、行列 \mathbf{X} がどの方向にどれだけ強く伸びているかを表します。

■ 固有値分解との関係

　SVD は、固有値分解の一般化と見ることができます。なぜなら、

\displaystyle{ \mathbf{X}^\top \mathbf{X} = \mathbf{V}\, \mathbf{S}^2\, \mathbf{V}^\top }

となり、\mathbf{X}^\top \mathbf{X} の固有値分解は、SVD から直接得られるからです。つまり、

  • 固有値分解は「対称行列」に対する分解
  • SVD は「任意の行列」に対する分解

だと言えます。

■ 最小二乗法と SVD

　最小二乗の節で見た最小二乗解

\displaystyle{ \boldsymbol{\beta} = (\mathbf{X}^\top \mathbf{X})^{-1}\mathbf{X}^\top \mathbf{y} }

は、\mathbf{X}^\top\, \mathbf{X} が逆行列をもたない場合や、条件が適切でない場合に、数値的に不安定になります。そのような場合も、SVD を用いると、逆行列を明示的に計算することなく、安定に解を求めることができます。
　この意味で、SVD は単なる理論的分解ではなく、実装と数値安定性の要となる行列分解です。

■ ランクと情報量の理解

　SVD を用いると、行列のランク構造が明確になります。小さな特異値に対応する成分は、データの中で情報量が少なく、ノイズに近い成分と解釈できます。
　この考え方は、

  • 次元削減
  • 正則化
  • 数値的不安定性の回避

といった問題で使われます。

■ PCA とのつながり

　共分散行列を経由せず、データ行列そのものに SVD を適用することで、主成分分析を実行することもできます。この方法は、高次元データやサンプル数が少ない場合に特に有効であり、実務では標準的に用いられています。 　このように、SVD は、

  • 一般行列を扱うための基本分解
  • 固有値分解，PCA，最小二乗法を統一する枠組み
  • 数値計算における安定性を保証する道具

として、線形代数の中核をなす概念です。

6. whitening の意味：共分散を単位行列にそろえる

　前節では、SVD によって一般の行列を安定に分解でき、最小二乗や PCA が数値的に安全に実装できることを説明しました。ここでは、それらと深く関係する whitening（白色化） の意味を整理します。

■ whitening とは何か

　多変量データ \mathbf{x} の共分散行列を \boldsymbol{\Sigma} とします。前節で見た固有値分解を用いると、

\displaystyle{ \boldsymbol{\Sigma} = \mathbf{V}\, \boldsymbol{\Lambda}\, \mathbf{V}^\top }

と書けます。このとき、共分散行列 \boldsymbol{\Sigma} が正則であれば、逆平方根を

\displaystyle{ \boldsymbol{\Sigma}^{-1/2} = \mathbf{V}\, \boldsymbol{\Lambda}^{-1/2}\, \mathbf{V}^\top }

と定義できます。これを用いて、データを

\displaystyle{ \mathbf{z} = \boldsymbol{\Sigma}^{-1/2}\, \mathbf{x} }

と変換する操作を whitening と呼びます。なお、\boldsymbol{\Sigma} が特異な場合には、逆平方根の代わりに擬似逆平方根（零に対応する固有値成分を除いたもの）を用いて同様の変換を行います。

■ 共分散が単位行列になる

　この変換後のデータ \mathbf{z} の共分散行列を計算すると、

\displaystyle{ \mathrm{E}[\mathbf{z}\mathbf{z}^\top] = \mathbf{I} }

となります。すなわち、whitening とは、

  • 各成分の分散を 1 にそろえ
  • 成分間の相関をすべて取り除く

操作だと理解できます。この意味で、whitening は、二次統計の意味で共分散構造を取り除いた標準化だと言えます。

■ 幾何学的な意味

　whitening は，次の二つの操作をまとめて行っています。

  1. 固有ベクトル方向への回転
  2. 固有値の平方根でのスケーリング

その結果、もとのデータ分布は、ばらつきの向きや大きさをもたない、等方的な分布へと変換されます。これは，データを「歪みのない座標系」に写し直している操作と考えることができます。

■ whitening 後に残る自由度

　whitening によって、共分散に関する情報はすべて取り除かれます。その結果、whitening 後のデータに対して許される変換は、

\displaystyle{ \mathbf{z} \mapsto \mathbf{Q}\mathbf{z}, \qquad \mathbf{Q}^\top\mathbf{Q}=\mathbf{I} }

のような 直交変換に限られます。
　この「直交変換の不定性」が、独立成分分析や SOBI において、次に解決すべき問題として現れます。

■ なぜ whitening が重要なのか

　whitening を行うことで、

  • スケールや相関の影響を除去できる
  • 問題の自由度を直交変換のみに縮約できる
  • 後続の解析が単純で見通しよくなる

という利点が得られます。このため whitening は，

  • PCA の前処理
  • ICA や SOBI の出発点
  • 多変量解析における標準化操作

として、非常に重要な役割を果たします。

7. まとめ

　研究者は、自身の研究を進める過程で必要となる統計や時系列解析の手法を、独学で身につけてきた方が多いと思います。その理解は、それぞれの試行錯誤や経験に基づくものであり、研究を前に進めるうえでは十分に機能しています。しかし、その知識を他人に伝えようとしたとき、特に基礎的な理解がまだ十分でない学生に教える場面では、思いのほか難しさを感じるのではないでしょうか。
　研究を山登りにたとえるなら、研究者は未踏の山に挑み、独自のルートで登頂を果たす存在だと言えるでしょう。そのルートは、険しく、回り道も多く、必ずしも人に勧められるものではないかもしれません。一方で、後に続く学生は、すでに登頂された山であっても、自らの力を養うために登る必要があります。そのときに重要になるのが、どの道を通れば全体を見渡しながら着実に登れるのかを示すことです。
　このような「適切な登山ルート」を示すことこそが、教育者の役割だと私は考えています。必ずしも最短ルートである必要はありませんが、見通しの悪い谷や、引き返すしかない割れ目に行き詰まることなく、今どこを登っているのか、そしてこの先に何が見えてくるのかを理解できる道筋を、学生に示していきたいと私は考えています。そのために、私自身もいろいろなルートを試行錯誤で学んでいきたいと思います。

登頂ルート

※ もし記事の中で「ここ違うよ」という点や気になるところがあれば、気軽に指摘していただけると助かります。質問や「このテーマも取り上げてほしい」といったリクエストも大歓迎です。必ず対応するとは約束できませんが、できるだけ今後の記事で扱いたいと思います。それと、下のはてなブログランキングはあまり信用できる指標ではなさそうですが (私のブログを読んでいる人は、実際とても少ないです)、押してもらえるとシンプルに励みになります。気が向いたときにポチッとしていただけたら嬉しいです。

ランキング参加中
知識
ランキング参加中
テクノロジー
ランキング参加中
数学・科学・工学