# ブログ（次元削減について）

## 次元削減

## 概要

### 次元削減とは？

画像、音響などの高次元データの集まりの構造または分布の全体像を人が解釈するのは不可能です。よって次元削減を用いて３次元または、それ以下の次元に落とせば容易にデータの解釈が行なえます。それを含めて他にも、次元削減の主な目的として挙げられるのは、

* データの可視化。
* データ圧縮。
* 低次元特徴量の作成。

例えば128x128ピクセルの白黒画像を考えた時、各画像のデータは128x128=16384次元のデータとなります。つまり、このサイズの画像データには16384次元の自由度があるわけですが、実際のデータはこの広大な空間のごく一部の体積しか占めておらず、データの分布には元の次元より遥かに低い内在的な次元があると仮定されています。この仮定は、「多様体仮説」とも言われ、次元削減における多様体学習の基盤となっております。([1] I. Goodfellow et al.)では高次元から一様的にサンプルされた画像に現実的なデータが現れることがないことや、動画などのデータの連続性などから多様体仮説の信憑性を裏付けています。

次元削減の説明によく用いられる具体的な例は次のスイスロールのデータです。
データは３次元ベクトルですが、データの分布（多様体）はそれより低い２次元の面です。良い次元削減アルゴリズムを用いれば、スイスロール広げて平たい２次元のデータへ圧縮することができます。
しかし、相応しくないアルゴリズムを用いればスイスロールを広げずに２次元空間へデータをそのまま射影するのみのため、データ点どうしの位置関係がごちゃごちゃになります。
データの解釈または、低次元特徴量として用いる際も位置関係が変わっていしまえば、データの誤認識、予測モデルの精度低下などが生じます。ここで重要なのは次元削減後のデータ点どうしの繋がりが元の高次元空間でのデータ点どうしの繋がりをできるだけ忠実に表現することです。それを行うには、アルゴリズムがデータの形を学習する必要性があります。よって多様体「学習」という用語が含まれています。

### 多様体とは？

機械学習においての多様体という用語の言い回しは数学上での正確な定義とは異なり、データ点どうしの繋がりから形成されるデータの形と解釈して良いと思います。しかし、論文・資料によっては正式な数学の定義で扱われることもあるので、機械学習を勉強している方は知っておいても損はないと思います。多様体とは幾何学が由来の用語であり、大まかには、トポロジカル多様体と可微分多様体があります。前者は距離などの概念は考慮せず物の繋がり方のみに焦点を置き、後者は繋がりのみだけでなく距離などを含め、滑らか（可微分）な形の物、例えばボールやドーナッツなど、を扱います。高校で習うような幾何学との大きな違いの一つとしては、座標の考え方で、座標というものは、人による「幻想」に過ぎないということがより強く明記されているように思えます。力学などでも、主にある座標を基準に物事を考えます。しかし、多様体を扱う微分幾何学などでは座標にこだわらず、それに依存しない（任意の滑らかな座標変換に対して不変である。）公式が書けるので物理法則など観測者に依存してはならないものの方程式が書ける便利な言語です。

機械学習の話に戻りましょう。実際の所、機械学習で多様体と呼んでいるものを統計学用語では、ただの「分布」です。ではなぜこの様な曲がりくねったような用語を用いるのでしょうか？ただ単に「分布学習」で良いじゃんと思われる方もいらっしゃると思います。

実際の歴史的な用語の区別の理由は存じませんが、その意図は、次をご覧いただけると何となくお分かりいただけると思えます。

https://www.autodeskresearch.com/publications/samestats

これらの２次元データは、異なるように見えますが、これらのデータの統計的値となる平均と分散、そしてデータの相関係数などを比較してみると「全て」が同じ値を持っています。データの形または意味が異なっても、これらの値において統計上、区別がつかないものになってしまいます。つまり、統計的なデータの表現のみでは多くの情報を見逃してしまう可能性を示唆しています。そこで、分布を考える際に統計学が焦点を置く値とは別に、より形に焦点をおいた多様体という用語が使われているように思えます。

さらに、低次元データに現れるガウス分布や他の解析的に統計学で定義されるような分布も、現実に扱われる高次元データにおいては殆ど見られません。そのうえ音響データなどは非定常なものが多く、形は多種多様であります。よって、とある分布を帰納的バイアスとした古典的モデル、例えば混合ガウス分布モデルなどの精度が低下する現象が見られる要因の一つだと考えられます。それに対して、よりゆるい帰納的バイアスを持った深層学習は、多種多様なデータの形に対応できるため、高次元かつビッグデータ領域では圧倒的な勝者となります。

### 次元削減の種類：

次元削減削減法は大まかに２種類あります。線形と非線形次元削減です。今回の記事で数学的な概要を紹介するのは、前者のほうですが、先程紹介した多様体学習は殆どの場合、後者を意味しています。アルゴリズムとしては、次のようなものが例として挙げられます。

線形|非線形
|:---:|:---:|
Principal Component Analysis (PCA) <br> Muti-Dimensional Scaling (MDS) <br> Non-negative Matric Factorization (NMF) <br> Independent Component Analysis (ICA) <br> Auto-Encoder (AE) | Kernel PCA <br> Isomap <br> Local Linear Embedding (LLE) <br> Laplacian Eigenmap <br> T-SNE <br> Uniform Manifold Assumption Projection (UMAP) <br> PHATE<br> Deep Auto-Encoder (DAE)

一見多い様ですが、これらは次元削減手法の本の一部です。個人的によく使うのは、PCA, [T-SNE](https://distill.pub/2016/misread-tsne/), [UMAP](https://pair-code.github.io/understanding-umap/)となりますが、どの手法を使ったら良いかは、ユーザーの目的、データの大きさや構造などによって異なりますので、どれがベストであるかは一意的に言えないです。線形、非線形の手法を比較した際、大まかなメリットやデメリットが存在します。それは、線形的手法は「大域的」なデータ構造を低次元で表現するのを得意とすることと、逆に非線形的手法は「局所的」な構造の表現を得意とすることです。非線形的手法の多くはデータ点の「近傍」の位置関係の低次元表現に焦点を置いているため大域的な構造が大分変わってしまうことがあります。例としてよく言われるのは、T-SNEは大域的な構造を現すのが不得意だと言われており、UMAPではそれがある程度改善されていると言われています。一方、Isomapは非線形的な手法ですが、大域的な手法でもあります。線形的な手法は、データの大まかな特徴を保てる部分ベクトル空間にデータを射影しているのみなので、細かく入り組んだ構造の情報は失われたり、誤認を生むような形に変形してしまいます。

### 次元削減アルゴリズムにて重要な機能

主な目的である次元削減という機能とは別にも、アルゴリズムによっては、次の幾つかの機能を搭載しています。

1. 学習データを圧縮された低次元から、元の高次元データに復元できること。
2. 学習データにない未知のデータを圧縮できること。
3. 圧縮された未知のデータを高次元データへ復元できること。

１は圧縮されたデータを解凍するのには必須であります。低次元の特徴量を用いて帰納的学習を行うには、２が必須となります。異常検知などでは、３ができないことを上手く利用したモデルもあります。多くの線形の手法は１から３を満たし、逆に多くの非線形の手法は１をも満たしません。アルゴリズムによっては２が可能でありつつも、未知のデータが学習データと独立同分布でなければ妥当な結果が得られないものもあります。他にも、目的によっては、次元削減後の座標に意味づけができるか否かなどが重要になります。多くの非線形アルゴリズムによる低次元表現の座標軸に人が解釈できる意味はなく、解析などの意味づけされた座標軸が必要である場合は、NMFなどが用いられます。

### 次元削減による可視化の注意点

#### 1. アルゴリズムよって生じる問題点

アルゴリズムによっては不適切なパラメータを用いるとデータ点どうしの繋がりをちぎってしまい、元のデータにはない隔離したクラスタが現れたりします。他にも、元のデータから密度やスケール感が大分変わってしまう表現がなされるアルゴリズムが多くありますので、各アルゴリズムの限度や弱点を把握しておくと良いと思います。可視化の信頼性をチェックする方法の一つは、クラスラベルや時刻などのメタデータを付与して可視化を行うことです。音響データなどの物理データには、時間方向に連続的な部分が多いので、不自然にデータが隔離していたりすれば、そこで実際に突発的な音響があるのかをチェックできます。なければ、可視化が上手くいっていない可能性があります。疑わしい時は複数のアルゴリズムを試してみたり、パラメータチューニングを行って見て下さい。

#### 2. データの形による限度
データの形によっては、可視化できる次元へ落とし込むこと（３次元またはそれ以下）が適切でない場合があります。冒頭でのスイスロールの件では２次元の多様体を２次元に落とし込みましたが、決して多様体の次元が２次元でも同じ２次元に落とし込むことができるとは限りません。有名な[向き付け不可能](https://ja.wikipedia.org/wiki/%E5%90%91%E3%81%8D%E4%BB%98%E3%81%91%E5%8F%AF%E8%83%BD%E6%80%A7#%E5%A4%9A%E6%A7%98%E4%BD%93%E3%81%AE%E5%90%91%E3%81%8D%E4%BB%98%E3%81%91%E5%8F%AF%E8%83%BD%E6%80%A7)な２次元多様体の例としては、クラインの壺や実射影平面です。


[Klein bottle - Wikipedia](https://en.wikipedia.org/wiki/Klein_bottle)
上のクラインの壺の図では分かりにくいのですが２次元の面で構成されています。自己接点が見られますが、実際の壺は４次元空間に存在し、接点は有りません。３次元空間では、壺の面の繋がりが正しく表現できず、接点では間違った構造の表現がなされるので、可視化に相応しくありません。よって、多様体の次元数が$m$だとしても決して同じ$m$次元の空間に落とせないことが分かります。ここで、「[ホイットニー埋め込み定理](https://en.wikipedia.org/wiki/Whitney_embedding_theorem)」によれば、どんな多様体の形でも埋め込む空間の次元数が最高でも$2m$あれば繋がりの関係を崩さずに表現できる事が保証されています。埋め込みに関しては他にも面白い証明がなされています。

* 向き付け可能な$m$次元の多様体は$2m-1$次元の空間に埋め込める。
* 多様体の次元数が$m$が２の整数乗でなければ、どれも$2m-1$次元に埋め込める。

幾何学的にいくら何をやっても可視化次元で元のデータの構造を上手く表現できないようなデータの形があることを頭に置いておくと良いと思います。

#### 次元削減のより正式な定義

ここからは、より数学チックな話しとなります。

元のデータの次元数を$d$とします。そして、データの多様体$\mathcal{M}$の次元数を$m$として冒頭で述べた多様体仮説を数学チックに書き換えると次となります。

>多様体仮説：
$d$次元のデータ点$\{\boldsymbol{x}_i\}_{i=1}^{n}$は、$\mathbb{R}^{d}$へ写像 $g: \mathcal{M} \to \mathcal{N} \subset \mathbb{R}^d$によって埋め込まれた$m$次元の多様体$\mathcal{N}$からサンプルされるとした時、$m<<d$であることが多様体仮説となります。

>次元削減：
次元削減のゴールとしては、もとの$\mathcal{M}$を復元することでありそのためには、$g$の逆写像$g^{-1}: \mathcal{N} \to \mathcal{M}$,　$g^{-1}({\boldsymbol{x}_i})\mapsto \boldsymbol{z}_i$を算出する事です。

これは、構造を「ありのまま」に保った次元削減と解釈できますが、３次元またはそれ以下のデータの可視化などでは、殆どこれを満たすはできないので、この記事では次のゆるい目的に書き換える事にしましょう。

>次元削減（その２）：
次元削減された$k$次元のデータ点$\{\boldsymbol{z}_i \}_{i=1}^n$は多様体$\mathcal{L} \subset \mathbb{R}^k$からサンプルされるとし、全てまたは部分的にな構造を保つ埋め込み写像$f:\mathcal{N} \to \mathcal{L}\subset \mathbb{R}^k, f({\boldsymbol{x}_i})\mapsto \boldsymbol{z}_i$を算出する事とします。

ものすごくゆるい定義となりました。どういう構造を保つかはアルゴリズムや目的によって異なる事とします。今回紹介する線形次元削減アルゴリズムたちは、目的が異なるように見えて結局、実は「皆同じ」構造を保っているという話しの流れになっています。

#### 本記事に現れる変数のリスト

これから代表的な線形次元削減の紹介をしていきます。それは、PCA、MDSとSVDとなります。そこで扱われる変数は次のようになります。

変数|定義|サイズ
|:---:|:---:|:---:|
$d$|入力データの次元数|なし|
$n$|入力データのサンプル数|なし|
$k$|次元削減後の次元数$\lt d$|なし|
$X$|$n$個の$d$次元入力ベクトルからなる行列。| $d\times n$
$\tilde{X}$|　$X$の$n$個のデータの平均によって中心化した行列 | $d\times n$
$U$|$d$個の$d$次元正規直交ベクトルからなる行列 |$d\times d$
$V$|$n$個の$n$次元正規直交ベクトルからなる行列 |$n\times n$
$\Sigma$|特異値分解に用いる特異値を対角成分に持つ行列 | $d\times n$
$Z$|次元削減後の$n$個の$k$次元ベクトルからなる行列 | $k\times n$|
$U_k$|$k$次元へ次元削減後の$U$ |$d\times k$|
$V_k$|$k$次元へ次元削減後の$V$ |$n\times k$|
$\Sigma_k$|$k$次元へ次元削減後の$\Sigma$ |$k\times k$|

## PCA
PCAについて説明します．
PCAとは "Principle Component Analysis" の略称であり，主成分分析と呼ばれています．
その名の通り観測したデータの中で重要な要素 "主成分" を探し出すアルゴリズムです．
主成分で構成したデータは観測したデータを少ない次元数で表現することができます．
では，どのようにして， "主成分" を見つけることができるのでしょうか？

データの中で重要な要素をどのように見つけることができるのか考えてみましょう．
これは，データの**情報量をなるべく保存**したまま，新しい点へ射影することができれば良いと考えることができます．射影とは物体の影をイメージしてもらえると分かりやすいと思います．例えば，N太くんを適当な方向（正面・真上・真横）に射影してみましょう．

N太くんで重要な要素はおそらく "ヒト" だと思われます．真上は "ヒト" の要素は全く捉えられていません．真横・正面は "ヒト" の要素を捉えていており，特に正面は良く "ヒト" の要素が捉えられています．
では，データの情報量を保つということはどのような方向なのでしょうか．これは，ばらつき（分散）が大きくなる方向にデータを射影することになります．上のN太くんを3次元のデータの集合と考えるとデータのばらつきは直感的に正面$＞$真横$＞$真上となりそうです．

PCAでは，
データの分散が最も大きい方向 => 第一主成分
データの分散が2番目に大きい方向 => 第二主成分
....
のように主成分を構成していきます．

では，実際に分散を最大化して，主成分を導出してみましょう．

### 1. 導出

$d$ 次元 $n$ 個のデータとして
$$
X= (\boldsymbol{x}_{1}, \boldsymbol{x}_{2}, ..., \boldsymbol{x}_{n}) =
\begin{pmatrix}
x_{11} & \cdots & x_{1j} & \cdots & x_{1n}\\
\vdots & \ddots & & & \vdots \\
x_{i1} & & x_{ij} & & x_{in} \\
\vdots & & & \ddots & \vdots \\
x_{d1} & \cdots & x_{dj} & \cdots & x_{dn}
\end{pmatrix}
$$
とします．

データの平均ベクトルと共分散行列をそれぞれ，
$$
\bar{\boldsymbol{x}} = \frac{1}{N}\sum_{n=1}^{N}(\boldsymbol{x_n}) \qquad \Sigma = \frac{1}{N}\sum_{n=1}^{N}(\boldsymbol{x_n}-\bar{\boldsymbol{x}})(\boldsymbol{x_n}-\bar{\boldsymbol{x}})^T
$$
と表すようにします．

主成分方向に写すため，新しいデータ空間が正規直交基底
$$
U=\boldsymbol{u}_{1}, \boldsymbol{u}_{2}, ... , \boldsymbol{u}_{d}
$$
で張られているとします.

1次元空間 $\boldsymbol{u}_{1}$ へのデータを射影することを考えてみましょう．
データ $\boldsymbol{x}_{n}$ を $\boldsymbol{u}_{1}$ へ写すことは内積によって，$\boldsymbol{u_1}^T\boldsymbol{x_n}$ と表すことができ，その平均は $\boldsymbol{u_1}^T\boldsymbol{\bar{x}}$です．
したがって，$\boldsymbol{u}_{1}$ に移したデータの分散は
$$
\frac{1}{N}\sum_{n=1}^{N}\{ \boldsymbol{u_1}^T\boldsymbol{x_n}-\boldsymbol{u_1}^T\boldsymbol{\bar{x}} \}^2 =
\frac{1}{N}\sum_{n=1}^{N}\{ \boldsymbol{u_1}^T(\boldsymbol{x_n}-\boldsymbol{\bar{x}}) \}^2 =
\boldsymbol{u_1}^T\Sigma\boldsymbol{u_1}
$$
となります．
ここまでで，新しいデータ空間へ射影した分散を導くことができました．PCAでは**分散が最大となるような方向にデータを射影していくこと**でした．したがって，先程導いた分散 $\boldsymbol{u_1}^T\Sigma\boldsymbol{u_1}$ を $\boldsymbol{u_1}$ に対して最大化ししていきましょう．しかし，この最大化は $\boldsymbol{||{u_1}||}\rightarrow\infty$ を防ぐような制約付きの最大化にする必要があります．$\boldsymbol{u_1}$ は正規直交であるため，$\boldsymbol{u_1}^T\boldsymbol{u_1}=1$ の制約を課して最大化を行います．このような制約付き最適化問題はラグランジュ乗数 $\lambda$ を用いて，量
$$
\boldsymbol{u_1}^T\Sigma\boldsymbol{u_1} + \lambda(1-\boldsymbol{u_1}^T\boldsymbol{u_1})
$$
を制約なしに最大化することで解くことができます(詳しい解法省きます)．$\boldsymbol{u_1}$ に関する微分を $0$ とすると以下の結果を得ることができます．
$$
\Sigma\boldsymbol{u_1} = \lambda\boldsymbol{u_1}
$$
左から， $\boldsymbol{u_1^T}$ を掛けて， $\boldsymbol{u_1^Tu_1}=1$ を使うと分散は
$$
\boldsymbol{u_1^T}\Sigma\boldsymbol{u_1} = \lambda
$$
と表すことができます．この式によれば分散と固有値は比例しています．分散を最大のものを選ぶ場合，射影した新しい基底として最も良いものは，固有ベクトルのうち，固有値が最大のものとなります．この固有ベクトルが第一主成分となります．
第二，第三主成分も，すでに得られた主成分と直交するという条件を加えて，同様の方法で解くことができます．

### 2.PCAによる次元削減
今度は先程求めた，固有ベクトルと固有値を使って次元削減を行いたいと思います．
前章で求めた $\lambda$ が大きいものから $k$ 個とってくるとすると，正規直交基底を並べた $u_1, ... , u_k$ で張られた空間は元の $d$ 次元ベクトルが $k$ 次元ベクトルとして表現することができます．正規直交基底を並べた $d\times k$ の行列を $U_{k} = \{u_1, ... , u_k\}$ とすると，
$$
Z = U_{k}^T\tilde{X}
$$
により $d\rightarrow k$ へ次元削減することができます．

#### 引用メモ
- [PRML下巻(p277~286)](https://www.amazon.co.jp/%E3%83%91%E3%82%BF%E3%83%BC%E3%83%B3%E8%AA%8D%E8%AD%98%E3%81%A8%E6%A9%9F%E6%A2%B0%E5%AD%A6%E7%BF%92-%E4%B8%8B-%E3%83%99%E3%82%A4%E3%82%BA%E7%90%86%E8%AB%96%E3%81%AB%E3%82%88%E3%82%8B%E7%B5%B1%E8%A8%88%E7%9A%84%E4%BA%88%E6%B8%AC-C-M-%E3%83%93%E3%82%B7%E3%83%A7%E3%83%83%E3%83%97/dp/4621061240/)
-

## Dual PCA
データ数に対して次元数が多い場合（例えば，遺伝子データなど）に対してPCAを行うのは計算コストの面から良くないです．そこで，PCAを $\boldsymbol{d}>>\boldsymbol{n}$ となるデータに対して計算コストを抑えたのがDual PCA です．では実際にPCAとDualPCAにどのような関係があるのか導出してみます．
以下では中心化したデータ行列として $\tilde{X}$を用います．
### 導出
PCAで導いた固有方程式に対して左から ${\tilde{X}}^T$ を掛けると
$$
\begin{eqnarray}
\Sigma\boldsymbol{u_1} &=& \lambda\boldsymbol{u_1}\\
\Leftrightarrow \tilde{X}^{T}\tilde{X}(\tilde{X}^{T}\boldsymbol{u_1}) &=& \lambda (\tilde{X}^{T}\boldsymbol{u_1})\\
\Leftrightarrow G\boldsymbol{v_1} &=& \lambda \boldsymbol{v_1}
\end{eqnarray}
$$
となります．$G$は $n \times n$ のグラム行列で ，$\boldsymbol{v_1}$ は $\boldsymbol{v_1} = \alpha \tilde{X}^T\boldsymbol{u_1}(\alphaは定数)$ としています．
上式をみるとPCAを求めたときの固有方程式と固有ベクトルが対応していることがわかります．この関係を導いてみます．$\boldsymbol{v_1}$を単位ベクトルとすると，以下のように $\alpha$ を求めることができます．

$$
\begin{eqnarray}
||\boldsymbol{v_1}||^2 = (\alpha \tilde{X}^{T} \boldsymbol{u_1})^2 &=& \alpha^2(\tilde{X}^{T}\boldsymbol{u_1})^T(\tilde{X}^T\boldsymbol{u_1}) \\
&=& \alpha^2 \boldsymbol{u_1^T}\tilde{X}\tilde{X}^T\boldsymbol{u_1} \\
&=& \alpha^2 \boldsymbol{u_1^T}\tilde{X}\tilde{X}^T\boldsymbol{u_1} \\
&=& \alpha^2 \lambda |\boldsymbol{u_1^T}\boldsymbol{u_1}| = 1 \\
\\
\Leftrightarrow \alpha = \frac{1}{\sqrt{\lambda}}
\end{eqnarray}
$$

この $\alpha$ を用いると以下のような関係を導くことができます．
$$
\boldsymbol{u_i} = \frac{1}{\sqrt{\lambda}}\tilde{X}\boldsymbol{v_i}
$$
このように，PCAとDualPCAの固有値問題の固有値は一致しており，この関係を用いることで，固有ベクトル同士も相互変換することができます．

### Dual PCA を用いた次元圧縮
では，Dual PCAの固有値問題で求めた固有値と固有ベクトルを用いて次元削減を行いたいと思います．PCAと同様に$\lambda$ が大きいものから $k$ 個とってきた規格化された固有ベクトル $\boldsymbol{v}_1, ... , \boldsymbol{v}_k$ で張られた空間は，元の $d$ 次元ベクトルを $k$ 次元ベクトルとして表現することができます．
固有ベクトルを並べた行列を $V_{k} = \{\boldsymbol{v}_1, ... , \boldsymbol{v}_k\}$ とし，先程導いた関係 $\boldsymbol{u_i} = \frac{1}{\sqrt{\lambda}}\tilde{X}\boldsymbol{v_i}$ により $U=\{\boldsymbol{u_1}, ..., \boldsymbol{u_k}\}$ へ変換すると，PCA と同様に，

$$
Z = U_{k}^T\tilde{X}
$$

によって $d\rightarrow k$ へ次元削減することができます．

さらに，この式を展開していくと，

$$
\begin{eqnarray}
Z = U_{k}^T\tilde{X} &=& \frac{1}{\sqrt{\lambda}}(\tilde{X}V_k)^T\tilde{X}\\
&=& \frac{1}{\sqrt{\lambda}}V_k^T\tilde{X}^T\tilde{X} \\
&=& \Sigma_kV_k^T
\end{eqnarray}
$$

のように表すこともできます．$ \Sigma_k$ は $\Sigma_k = \Lambda_k^{1/2} = \text{diag}(\lambda_1,...,\lambda_k)^{1/2}$ と表される固有値を並べた対角行列です．

## MDS

Multi-Dimensional Scaling (MDS)では、元のデータ点どうしの距離を考えます。例えば、場所どうしの距離の情報からマップを描くようなイメージです。結論から言うと、この距離がユークリッド距離である場合は、DualPCAと同じ結果となります。距離を変えることで非線形の次元削減も可能です。例としてのIsomapは、距離を測地線上で計る、MDSが元となる非線形次元削減アルゴリズムです。

MDSの目的関数は、元の$d$次元であるデータ点どうしの距離関係と、低次元$k$での距離関係を近似させることで、次を最小化する$\{\boldsymbol{z}_i\}_{i=1}^n \subset \mathbb{R}^k$探すことです。

$$
\displaystyle\min_{\boldsymbol{z}_i \in \mathbb{R}^k}\displaystyle\sum_{i,j}(\underbrace{\tilde{d}(\boldsymbol{z}_i,\boldsymbol{z}_j)^2}_{\small \text{低次元での距離}\\\small\text{（非類似度）}} - \underbrace{d(\boldsymbol{x}_i,\boldsymbol{x}_j)^2}_{\small \text{入力次元での距離}\\\small\text{（非類似度）}} )^2
$$

距離がユークリッドの場合は、これを内積の行列、グラム行列$G$を用いて次の目的関数に書き換えることができます。
$$
\displaystyle\min_{Z \in \mathbb{R}^{k \times n}}\|Z^TZ - G \|_F^2
$$
[Eckart-Young-Mirskyの定理](https://en.wikipedia.org/wiki/Low-rank_approximation#Proof_of_Eckart%E2%80%93Young%E2%80%93Mirsky_theorem_(for_Frobenius_norm))を用いれば、半正定値であるグラム行列$G$をスペクトル分解し、トップ$k$の固有値による$G$のランク$k$近似で解が導き出せることが分かります。
### 導出
#### Step① 距離行列をグラム行列に書き換える。
グラム行列はデータの内積の行列です。つまり、距離と内積の関係が分かれば書き換えることができます。
$d$次元ベクトル$x_i$と$x_j$のユークリッド距離の二乗を次のように書くと距離と内積の関係が式に現れます。
$$
d_{ij}^2= \| \boldsymbol{x}_i - \boldsymbol{x}_j \|^2 = (\boldsymbol{x}_i-\boldsymbol{x}_j )^T(\boldsymbol{x}_i-\boldsymbol{x}_j ) \\ = \boldsymbol{x}_i^T\boldsymbol{x}_i+\boldsymbol{x}_j^T\boldsymbol{x}_j - 2\boldsymbol{x}_i^T\boldsymbol{x}_j
$$

この関係を用いり,$n \times n$のデータ点どうしの内積を要素とした行列、$\boldsymbol{k} = X^TX$と$n \times n$の要素がすべて１である行列、$\boldsymbol{1}$を定義すると次のように書けます。

$$
D: = (d_{ij}^2) = \boldsymbol{k} \cdot \boldsymbol{1}^T + \boldsymbol{1} \cdot \boldsymbol{k}^T - 2\boldsymbol{k}
$$

[エルミート](https://ja.wikipedia.org/wiki/%E3%82%A8%E3%83%AB%E3%83%9F%E3%83%BC%E3%83%88%E8%A1%8C%E5%88%97#:~:text=%E7%B7%9A%E5%9E%8B%E4%BB%A3%E6%95%B0%E5%AD%A6%E3%81%AB%E3%81%8A%E3%81%91%E3%82%8B%E3%82%A8%E3%83%AB%E3%83%9F%E3%83%BC%E3%83%88,%E3%82%88%E3%81%86%E3%81%AA%E3%82%82%E3%81%AE%E3%82%92%E8%A8%80%E3%81%86%E3%80%82)である、中心化行列、$H = I - \frac{1}{n} \boldsymbol{1}\cdot\boldsymbol{1}^T$を用いると、元のデータ$X$の平均を０に移せます。これによって、$X$の中心化、$\tilde{X}:=XH$を行列の演算として定義する事ができます。そして、$n \times n$のグラム行列$G$は中心化された元のデータ$X$の内積となり、上記の距離行列と次のような関係を持ちます。

$$
G:=\tilde{X}^T\tilde{X} = HX^TXH = - \frac{1}{2}HDH
$$

#### Step②　行列の対角化を計算し、トップｋの固有値を選択する。

上記で述べたように次のコスト関数の最小解は、スペクトル分解によって導けます。それは、元の高次元空間での内積行列を近似する次元削減後の内積行列を求めるのに等価しています。
$$
\displaystyle\min_{ \text{rank}(G_k)=k} \| G - G_k\|_F^2 \iff
\displaystyle\min_{{Z \in \mathbb{R}^{k \times n}}} \|\underbrace{\tilde{X}^T\tilde{X}}_{\small\text{入力の次元での内積}} - \underbrace{Z^TZ}_{\small \text{低次元での内積}} \|_F^2
$$
行列$G$の固有値と固有ベクトルを求め、固有値の多きいものから順番に並べていきます。ここで、$G$は半正定値であるので、非負の固有値を持つ事に注意します。その後、トップ$k$個の固有値とそれに対応する固有ベクトル、$V_k=(\boldsymbol{v}_1,...,\boldsymbol{v}_k)$、用いれば解は次となります。

$$
G_k=V_k\Lambda_k V_k^T, \qquad \Lambda^k=\text{diag}(\lambda_1,...,\lambda_k)
$$
$k$次元への削減は、特異値分解に現れる$\Sigma_k$を用いれば次の用にも書けます。
$$
Z=\Lambda_k^{1/2}V_k^T=\Sigma_kV_k^T
$$

データを中心化した、ユークリッド距離によるMDSの解は、DualPCAと全く同じ解となっております。よって、同じように$d>>n$の時にPCAよりも計算効率が良く、次元数が多くデータ数が少ない遺伝子データなどを扱う場合はこちらのほうが速く計算できます。

## SVD

特異値分解(Singular Value Decomposition (SVD) )とは，三角化や，対角化といった行列分解の一つです．幾何的な意味を確認しつつ，導出方法を確認していきましょう．

### 1. 導出

まず，任意の$d \times n$の行列で、最良の一次元部分空間を見つけることを考えてみましょう．
「最良」とは，各点における部分空間への距離が最小である，ということです．(下図でいうと，距離$b$に当たります．)

(簡略のため，$d=2$と仮定)

ここで図を見ると，三平方の定理より，各点の位置($a$)と，部分空間への射影後の位置($c$)がわかれば，$b$を導けることがわかります．これを全点において考えると，
$$
\sum_{i=1}^{n}\left(b_{i 1}^{2}+b_{i 2}^{2}+\cdots+b_{i d}^{2}\right) = \sum_{i=1}^{n}\left(a_{i 1}^{2}+a_{i 2}^{2}+\cdots+a_{i d}^{2}\right) - \sum_{i=1}^{n}\left(c_{i 1}^{2}+c_{i 2}^{2}+\cdots+c_{i d}^{2}\right)
$$
が成り立ちます．右辺の第一項は定数なので，左辺を最小化するには，右辺第二項を最大化すれば良いこととなります．

上記の最良についての定義と説明では行列を$n$個の$d$次元のベクトル（列ベクトル空間）として扱いましたが、次からは、横の方向に見て、行列を$d$個の$n$次元ベクトル（行ベクトル空間）とした観点から話しを進めて行きたいと思います。
任意の$d \times n$行列$X$の行ベクトル空間を考えたとき、各行ベクトル$\boldsymbol{\xi_j}$は任意の$n$個の$n$次元の正規直交行列、$\{\boldsymbol{v}_i\}_{i=1}^n$で次のように展開できます。
$$
X　=\left[\begin{array}{c}
\boldsymbol{\xi_{1}^{T}} \\
\boldsymbol{\xi_{2}^{T}} \\
\vdots \\
\boldsymbol{\xi_{d}^{T}}
\end{array}\right] =\left[\begin{array}{c}
\displaystyle\sum_{i=1}^n\left\langle \boldsymbol{\xi_{1}, v_i}\right\rangle \boldsymbol{v}_i^T\\
\displaystyle\sum_{i=1}^n\left\langle \boldsymbol{\xi_{2}, v_i}\right\rangle \boldsymbol{v}_i^T\\
\vdots \\
\displaystyle\sum_{i=1}^n\left\langle \boldsymbol{\xi_{d}, v_i}\right\rangle \boldsymbol{v}_i^T
\end{array}\right]
$$

$X$はこのように任意の正規直交基底で展開できますが、SVDではどのような正規直交基底を選べば,「最良」の度合いで基底をランク付けできるかを考えます。そのため、$\{\boldsymbol{\xi}_j\}_{j=1}^d$の基底ベクトルへの射影を考えます。

ここで、$\{\boldsymbol{v}_i\}_{i=1}^n$のうちの一つを考え、それを$n \times 1$の単位ベクトル、$\boldsymbol{v}$とします．
すると，行列$X$の$\boldsymbol{v}$への射影は，以下のように表すことができます．
$$
X \boldsymbol{v}=\left[\begin{array}{c}
\boldsymbol{\xi_{1}^{T}} \\
\boldsymbol{\xi_{2}^{T}} \\
\vdots \\
\boldsymbol{\xi_{d}^{T}}
\end{array}\right] \boldsymbol{v}=\left[\begin{array}{c}
\boldsymbol{\xi_{1}^{T} v} \\
\boldsymbol{\xi_{2}^{T} v} \\
\vdots \\
\boldsymbol{\xi_{d}^{T} v}
\end{array}\right]=\left[\begin{array}{c}
\left\langle \boldsymbol{\xi_{1}, v}\right\rangle \\
\left\langle \boldsymbol{\xi_{2}, v}\right\rangle \\
\vdots \\
\left\langle \boldsymbol{\xi_{d}, v}\right\rangle
\end{array}\right]
$$
またこのとき射影の長さの二乗は、
$$
\|X \boldsymbol{v}\|^{2}=\left|\left\langle \boldsymbol{\xi_{1}, v}\right\rangle\right|^{2}+\left|\left\langle \boldsymbol{\xi_{2}, v}\right\rangle\right|^{2}+\cdots+\left|\left\langle \boldsymbol{\xi_{d}, v}\right\rangle\right|^{2} \\
=\boldsymbol{v}^TX^TX\boldsymbol{v}
$$
と表すことができます。そして、最良となるとなる基底は、
$$
\boldsymbol{v_1}: = \underset{\|\boldsymbol{v}\|=1}{\text{argmax }} \boldsymbol{v}^TX^TX\boldsymbol{v}
$$
となります。上記の目的関数は、PCAで記載されているものと形式が同じであるので、解は、最大の固有値を持つ$X^TX$の固有ベクトルとなり、これを第一特異ベクトルと言います。よって、ノルム$\|X\boldsymbol{v_1}\|$は、その固有値のルートをとったものとなり、これを第一特異値、$\sigma_1$と言います．

次に，二次元部分空間$\underset{\boldsymbol{v}_1\perp\boldsymbol{v},\|\boldsymbol{v}\|=1}{\text{span}}\{\boldsymbol{v_1},\boldsymbol{v}\}$ を考え、それを最良とする$\boldsymbol{v}$を第二特異ベクトルと言います。
$$
\boldsymbol{v_2}: = \underset{\boldsymbol{v}_1\perp\boldsymbol{v},\|\boldsymbol{v}\|=1}{\text{argmax }} \boldsymbol{v}^TX^TX\boldsymbol{v}
$$

第二特異ベクトル$\boldsymbol{v_2}$は第一と同様に$X$の射影の長さが第二特異値$\sigma_2$となります。以降の右特異ベクトルや特異値も同様な形で定義し算出することができます．
$$
\boldsymbol{v_l}: = \underset{\boldsymbol{v}_1...\boldsymbol{v}_{l-1}\perp\boldsymbol{v},\|\boldsymbol{v}\|=1}{\text{argmax }} \boldsymbol{v}^TX^TX\boldsymbol{v}, \quad \sigma_l = \|X\boldsymbol{v_l} \|
$$

さて，ここで行列$V=[\boldsymbol{v_1} \, \boldsymbol{v_2} \, ... \, \boldsymbol{v_n}]$を考えます．
すると，
$$
XV=[X\boldsymbol{v_1} \, X\boldsymbol{v_2} \, ... \, X\boldsymbol{v_n}]
$$

ここでの注意点は、行列$X$のサイズ（$d\times n$）によっては特異値がゼロとなる固有ベクトルが必然的に存在することです。今回は、$d<n$であることを想定しましょう。そうなると、$\text{rank}(X^TX)=\text{rank}(X)\le \min\{d,n\}=d$となり[階数・退化次数の定理](https://ja.wikipedia.org/wiki/%E9%9A%8E%E6%95%B0%E3%83%BB%E9%80%80%E5%8C%96%E6%AC%A1%E6%95%B0%E3%81%AE%E5%AE%9A%E7%90%86)から$\text{nullity}(X^TX)\ge n-d$となります。$\text{nullity}(X^TX)$は固有値がゼロとなる固有ベクトルの数に値します。そして、
$$
X^TX\boldsymbol{v}=0 \implies \boldsymbol{v}X^TX\boldsymbol{v}=\|X\boldsymbol{v}\|^2=0 \implies \|X\boldsymbol{v}\|=0 \implies X\boldsymbol{v}=0\\
\boldsymbol{v} \in \text{ker}(X)
$$
となりますので、この場合の$XV$の$d$以上の要素はゼロとなります。
$$
XV=[X\boldsymbol{v_1} \, X\boldsymbol{v_2} \, ... \, X\boldsymbol{v_d},0,...,0]
$$
ここで，正規化された$X\boldsymbol{v_j}(1 \le j \le d)$を$\boldsymbol{u_j}$とおくと，
$$
X\boldsymbol{v_j} = \sigma_j\boldsymbol{u_j}
$$

が成り立ちます．ただし、$X(X^TX\boldsymbol{v_j})=\sigma_j^2X\boldsymbol{v_j}\implies XX^T\boldsymbol{u_j}=\sigma_j^2\boldsymbol{u_j}$、からも$\boldsymbol{u_j}$が$XX^T$の固有ベクトルであることが分かり、これらを左特異ベクトルと言います。

上記をすべての$j$についてまとめると，
$$
XV = [\boldsymbol{u_1} \, \boldsymbol{u_2} \, ... \, \boldsymbol{u_d}]
\left[\begin{array}{ccc}
\sigma_{1} & & & O & 0 & \dots &0\\
& \sigma_{2} & & &0 & &0\\
& & \ddots & &\vdots & &\vdots\\
O & & & \sigma_{d} & 0 & \dots &0
\end{array}\right] = U\Sigma
$$
(ただし，$U=[\boldsymbol{u_1} \, \boldsymbol{u_2} \, ... \, \boldsymbol{u_d}]$, $\Sigma = \left[\begin{array}{ccc}
\sigma_{1} & & & O & 0 & \dots &0\\
& \sigma_{2} & & &0 & &0\\
& & \ddots & &\vdots & &\vdots\\
O & & & \sigma_{d} & 0 & \dots &0
\end{array}\right]$)

以上より，
$$
X=U \Sigma V^{T}
$$
と特異値分解することができました．

### 2. SVDによる行列の近似
ここで，$\left\|X-A\right\|_{F}$を最小にするような行列$A$による$X$の低ランク近似($\text{rankA}=k$)($0<k<\text{rank(X)}$)について考えます．SVDでは、最良の度合いによって要素がランク付けされていますので、ランク削減の際に最下位の要素を省くことで、情報損失を最小限に抑えた、１ランク下の行列の表現が可能となります。よって次の[フロベニウスノルムによるEckart-Young-Mirskyの定理](https://en.wikipedia.org/wiki/Low-rank_approximation#Proof_of_Eckart%E2%80%93Young%E2%80%93Mirsky_theorem_(for_Frobenius_norm))が成り立ちます。

$$
\min _{\mathrm{rank} A = k}\|X-A\|_{F}=\left\|X-X_{k}\right\|_{F}
$$

つまり，$X$に対して特異値分解を用いて$k$次元($0<k<d$)へランク削減を行った行列$X_k$($\text{rank}(X_k)=k$)が，フロベニウスノルムの意味で$X$の最良の近似となっています．実際にこの解は、フロベニウスノルムのみに限らずスペクトルノルムや核型ノルムでも成立します。この定理の証明に関しては、シンプルな証明が[math stackexchange](https://math.stackexchange.com/questions/759032/proof-of-eckart-young-mirsky-theorem)に挙げられていますので良かったらご覧下さい。

## PCA, MDSとSVDの関係

今回紹介した次元削減の手法を次にまとめました。

<br>|PCA|DualPCA or Euclidean MDS
|:---:|:---:|---:|
学習データの次元削減| $Z=U_k^TX$ | $Z=\Sigma_k V_k^T$
学習データの復元| $X_k = U_k Z$ | $X_k = XV_kV_k^T$
未知なテストデータの次元削減| $\boldsymbol{z}=U_k^T \boldsymbol{x}$ |$\boldsymbol{z}=\Sigma_k^{-1} V_k^T X^T \boldsymbol{x}$
未知なテストデータの復元| $\boldsymbol{x_k}=U_k \boldsymbol{z}$ | $\boldsymbol{x_k}= XV_k \Sigma_k^{-2} V_k^T X^T \boldsymbol{z}$
ここで注意するのは、PCAは、左特異値ベクトルのみを扱い、MDSは、右特異値ベクトルのみを扱う事です。そして、SVDによる、ランクｋ行列の近似を用いることで、これらの手法は次のように繋がります。

$$
\Large{X \approx X_k =} \underbrace{\Large{\quad U_k \quad}}_{\small\text{PCAによる復元演算子}}\underbrace{\Large{\quad　\Sigma_kV_k^T　\quad}}_{\small\text{MDSによる次元削減データ}}
$$

これによって、目的が異なったこの2つの手法は統一されます。同じ結果がでるアルゴリズムを、このように2つの見方が可能であるので、これらの手法は互いに「双対」であると言われます。その他にも、PCAによる$X \approx X_k = U_kU_k^TX$の誤差がフロベニウスノルムの最小化解であることを鑑みると、これは、三層のニューラルネットである線形オートエンコーダと等価することとなります。

## まとめ

今回は、次元削減の紹介と概要、そして、幾つかの線形的次元削減の詳細を紹介し、それらが等価することを示しました。キーポイントは次となります。

* 可視化による誤認識が起きる可能性もあるので要注意！
* 一番良い次元削減法と言えるものはなく、それぞれ弱点や長所あり、目的によっても異ってくる。
* PCAは分散の最大化を用いて、データの低次元表現を行う手法。
* MDSは低次元での類似度を低次元で近似することで、データの低次元表現を行う手法。
* SVDは、フロベニウスノルムによる誤差において、低ランク行列による行列の近似の最適解を導ける。（実際には、他の幾つかのノルムで成立する。）
* 中心化されたデータとユークリッド距離を用いた場合は、PCAとMDSは双対の関係にある。
* PCAは特徴量の次元がデータ数より低いさいに計算効率がよく、MDSは、その逆。

次回は非線形次元削減の手法をご紹介できたら良いと思います。

## 参考文献
[1] I. Goodfellow, Y. Bengio and A. Courville, *DEEP LEARNING*, Massachusetts Institute of Technology, 2016.
[2] Y. Ma and Y. Fu, *Manifold Learning Theory and Application*, Taylor & Francis Group, 2012.
[3] C. M. Bishop, *Pattern Recognition and Machine Learning*, Springer-Verlag New York, 2006.
http://seetheworld1992.hatenablog.com/entry/2018/08/16/155013
http://seetheworld1992.hatenablog.com/entry/2017/04/25/114244
https://qiita.com/kidaufo/items/0f3da4ca4e19dc0e987e