Groenen, P.J.F, van de Velden, M. (2016) Multidimensional Scaling by Majorization: A Review. Journal of Statistical Software. 73(8).

　仕事の都合でMDSについて考えていて(滅多にないことである)、Rのsmacofパッケージの実戦投入に先立つ儀式として読んだ論文。
　smacofパッケージについては開発者自身による紹介論文を読んだことがあるのだが、そのときはあまり理解できなかった。この論文はsmacofパッケージの紹介ではあるが、書いているのは第三者のようで、もっとわかりやすいかな、と思って。

　いわく。
　MDSの最初期のアプローチとして、Togerson(1958), Gower(1966)の「古典的尺度化」がある(Gowerは「主座標分析」と呼んだ)。古典的尺度化は固有値分解に基づいており、PCAと密接な関係がある。その後Kruskal(1964)の最小二乗アプローチが登場した[いわゆる非計量的MDSのこと？]。このアプローチではユークリッド距離を(変換した)非類似性に直接フィットさせる。現在最も広く使われているアプローチであり、ここで考えるMDSもこれである。
　最小二乗MDSの基盤となる数値最適化において広く使われているのがJan de LeeuwのSMACOFアルゴリズムである。本論文はde Leeuwの貢献を紹介し、smacofパッケージを使って事例を示す。

　まずは最小二乗MDSについて。
　\(n\)個の対象があるとする。対象\(i\)と\(j\)の対象間非類似性\(\delta_{ij}\)が与えられたとする。対象が\(p\)次元空間において持つ座標を\(n \times p\)行列\(\mathbf{X}\)で表す。この空間におけるユークリッド距離を\(d_{ij}(\mathbf{X}) = \sqrt{\sum_s^p (x_{is}-x{js})^2}\)と書く。
　Kruskalはこう考えた。次式をストレス1と呼ぶ。これが損失関数になる。$$ \sigma_1(\mathbf{X}, \mathbf{\hat{d}}) = \left( \frac{\sum_{i < j} (\hat{d}_{ij} – d_{ij}(\mathbf{X}) )^2}{\sum_{i < j} d^2_{ij}(\mathbf{X})} \right)^{1/2}$$ \(\hat{d}_{ij}\)は所与の\(\delta_{ij}\)の関数で、擬似距離ないしディスパリティと呼ぶ。ここの関数についてKruskalが考えたのは、とにかく順位だけは保たれる、という制約である。この制約を満たす関数はいろいろありうる。線形変換\(\hat{d}_{ij} = a + b \delta_{ij}\)とか、多項式\(\hat{d}_{ij} = b_0 + \sum_q b_q \delta_{ij}^q\)とか、単調スプライン変換とか。

　ここからde Leeuwの考え方に沿って話を進めよう。
　\(\hat{d}_{ij} = \delta_{ij} \geq 0\)とする[おっと… それって計量的MDSではないか。なんだかよくわかんなくなってきた…]。で、ストレス1の分子をちょっと変えた$$ \sigma^2_{raw} (\mathbf{X}, \mathbf{\hat{d}}) = \sum_{i < j} w_{ij} (\hat{d}_{ij} – d_{ij}(\mathbf{X}))^2$$を生ストレスと呼ぶ。ただし\(w_{ij}\)は非負とする。
　ここでのポイントはウェイト\(w_{ij}\)である。1993年、著者のひとりがUCLAの統計学部を訊ねて、De Leeuwに「なんで\(w_{ij}\)を追加したんですか」と訊ねた。先生いわく、このウェイトのおかげで問題が複雑になり、従って、面白くなる。[←ははは]
　ウェイトのスケーリングはなんでもいいので、\(\sum_{i < j} w_{ij} \hat{d}^2_{ij} = n(n-1)/2\)と制約しておく(ウェイトをつけたディスパリティの平均を1にしておく)。こうして、ストレス1の最小化という問題が、生ストレスの制約付き最小化という問題に変わった。
　さて、生ストレス関数を次のように分解しておく: $$ \sigma^2_{raw} (\mathbf{X}, \mathbf{\hat{d}}) = \sum_{i < j} w_{ij} \hat{d}_{ij}^2 + \sum_{i < j} w_{ij} d_{ij}(\mathbf{X})^2 – 2 \sum_{i < j} w_{ij} \hat{d}_{ij} d_{ij}(\mathbf{X})$$ $$= \eta^2_{\hat{d}} + \eta^2(\mathbf{X}) – 2\rho(\mathbf{X})$$ 第1項から順に、全ディスパーション、再現ディスパーション、共ディスパーションと呼ぶ。
　\(\sigma^2_{norm}(\mathbf{X})=\sigma^2_{raw}/\eta^2_{\hat{d}}\)を規準化ストレスと呼ぶ。\(\mathbf{X}\)がlocal minimumのとき、規準化ストレスは0から1の値となる。なぜなら…[説明があったんだけど、学力不足のせいかよく理解できなかった。あとで見直すときのためにメモしておくと、\(\mathbf{X}^{*}\)が生ストレスのlocal mimimumなら、(\(\mathbf{\hat{d}}\)を固定して考えれば) \( \frac{\partial \sigma^2_{raw}(\alpha \mathbf{X}^{*})}{\partial \alpha} = 0\) になるはずだ、というくだりがわからない…]。

　さて、本題。Le Leeuwが提案したSMACOFアルゴリズムは、実のところ、ある一般的な最適化アプローチであった。彼はこれをmajorizationと呼んだ。minimization by majorisation, 略してMMアルゴリズムと呼ばれることもある。別の分野ではgeneralized Weiszfeld’s method とかconvex-concave手続きとも呼ばれている。
　最小化したい関数を\(f(\mathbf{x})\)とする(\(\mathbf{x}\)は\(p\)次元の実数ベクトル)。これを、もっと単純な関数\(g(\mathbf{x}, \mathbf{y})\)で置き換える。この関数をmajorizing functionという。\(\mathbf{y}\)はsupporting pointと呼ばれる既知のベクトルで、たとえば、前のイテレーションで得た点である。
\(g(\mathbf{x}, \mathbf{y})\)は次の3つの性質を持つ必要がある。

• \(f(\mathbf{y}) = g(\mathbf{y}, \mathbf{y}) \)
• \(f(\mathbf{x}) \leq g(\mathbf{x}, \mathbf{y})\)
• 単純であること。ふつうは線形か二次とする。

　いま手元に\(\mathbf{y}\)があるとする。で、\( g(\mathbf{x}^+, \mathbf{y}) \leq g(\mathbf{y}, \mathbf{y})\)である\( \mathbf{x}^+ \)をみつけたとする(\(g\)は簡単な関数なので見つけるのは容易である)。このとき、$$ f(\mathbf{x}^+) \leq g(\mathbf{x}^+, \mathbf{y}) \leq g(\mathbf{y}, \mathbf{y}) = f(\mathbf{y})$$である。De Leeuwはこれをサンドイッチ不等性と呼んだ。ここから、みつけた\(\mathbf{x}^+\)を\(\mathbf{y}\)に代入して繰り返せばよいことがわかる。これは必ず収束することがわかっている。

　話を進める前に記号を導入する。
　平方ユークリッド距離を次のように書き換える。$$ d^2_{ij}(\mathbf{X}) = \sum_s^p (x_{is}-x_{js})^2 $$ $$ = \sum_s^p (\mathbf{x}_s^\top(\mathbf{e}_i – \mathbf{e}_j))^2 = \sum_s^p \mathbf{x}_s^\top (\mathbf{e}_i – \mathbf{e}_j) (\mathbf{e}_i – \mathbf{e}_j)^\top \mathbf{x}_s$$ \(\mathbf{x}_s\)は\(\mathbf{X}\)の\(s\)列目, \(\mathbf{e}_i\)は単位行列\(I\)の\(i\)列目である。一列ずつちまちま処理してないで、どーんと行列演算してあとで必要な要素を足し上げましょう。$$ = tr \mathbf{X}^\top (\mathbf{e}_i – \mathbf{e}_j)(\mathbf{e}_i – \mathbf{e}_j)^\top \mathbf{X} = tr \mathbf{X}^\top \mathbf{A}_{ij} \mathbf{X}$$ \(\mathbf{A}_{ij} \)は、\(a_{ii} = a_{jj} = 1\)、\(a_{ij} = a_{ji} = -1\)、残りが0である。
[頭のいい人には簡単なんだろうけど、私にはちょっと難しい。ええと、\(\mathbf{X}\)は\(3 \times 2\)としよう。たとえば要素1と2の距離を求めるとして、\(s\)列目の差の二乗は $$(x_{1s}-x_{2s})^2 = \left( \left[ \begin{array}{ccc} x_{1s} & x_{2s} & x_{3s} \end{array} \right] \left( \left[ \begin{array}{c} 1 \\ 0 \\ 0\end{array} \right] – \left[ \begin{array}{c} 0 \\ 1 \\ 0\end{array} \right] \right) \right)^2 $$ $$= (\mathbf{x}_s^\top(\mathbf{e}_1 – \mathbf{e}_2))^2 = \mathbf{x}_s^\top (\mathbf{e}_1 – \mathbf{e}_2) (\mathbf{e}_1 – \mathbf{e}_2)^\top \mathbf{x}_s$$と書ける。で、$$ (\mathbf{e}_1 – \mathbf{e}_2) (\mathbf{e}_1 – \mathbf{e}_2)^\top = \left[ \begin{array}{c} 1 \\ -1 \\ 0 \end{array} \right] \left[ \begin{array}{ccc} 1 & -1 & 0 \end{array} \right] = \left[ \begin{array}{ccc} 1 & -1 & 0 \\ -1 & 1 & 0 \\ 0 & 0 & 0 \end{array} \right]$$ を\(\mathbf{A}_{12}\)と呼んでいるわけだ。なるほど…]

　さて、話を戻して、生ストレスの分解$$\sigma^2_{raw}(\mathbf{X}, \mathbf{\hat{d}}) = \eta^2_{\hat{d}} + \eta^2(\mathbf{X}) – 2\rho(\mathbf)$$ をSMACOFで最小化したい。majorizing functionは…
　[ここから約1パージにわたり、3つの要素を行列で書き直して、コーシー・シュワルツの不等式を使って… 云々という議論が続く。力尽きたのでパス…]
　… というわけで、majorizing functionとして下式が手に入った。
$$g(\mathbf{X}, \mathbf{Y}) = \eta^2_{\hat{d}} + tr \mathbf{X}^\top \mathbf{V} \mathbf{X} – 2 tr \mathbf{X}^\top \mathbf{B(Y)Y} $$ $$ \mathbf{V} = \sum_{i < j} w_{ij} \mathbf{A}_{ij} $$ \(\mathbf{B}(\mathbf{Y})\)の要素は $$ b_{ij} = \begin{cases} \frac{ w_{ij} \hat{d}_{ij} }{ d_{ij}(\mathbf{Y}) } & \mathrm{if} \ d_ij(\mathbf{Y}) > 0 \\ 0 & \mathrm{if} \ d_ij(\mathbf{Y}) = 0\end{cases} $$ である。[はいはいわかりました！なんだかしらんが信じますよ！]

　\(\mathbf{X}\)について、たとえば外的変数の\(n \times q\)行列\( \mathbf{H} \)があって\( \mathbf{X} = \mathbf{HC} \)であり \(C\)を最適化したい、というような制約がかかっていたらどうするか。[…スキップするけど、smacofパッケージではそういうこともできるらしい…]

　局所解にどのくらい捕まるか [… 著者のDe Leeuwさんとの個人的思い出話なんかが入っていて面白いんだけど、手に負えないのでスキップ。要するに局所解にがんがん捕まるので、初期値を変えて複数回走らせるべしとのこと。smacofパッケージにはそういう機能もついている由]

　ウェイトについて。De Leeuwさんが導入した\(w_{ij}\)のおかげで、

1. 非類似性行列の欠損を簡単に扱える。その箇所のウェイトを0にするだけでよい。
2. 非類似性の分布が特定の値に固まっているとき、これを一様分布に修正できる。[数値例を使った説明がある。パス]
3. MDSの他の損失関数を模倣できる。
   • Heiser(1988)のSストレス \(\sum_{i < j} (\delta^2_{ij} – d^2_{ij}(\mathbf{X})) \)は[…説明は飛ばして…] \(w_{ij} = 4 \delta^2_{ij}\)とすると近似できる。
   • 多項ストレス \( \sum_{i < j} w_{ij}(\delta^r_{ij} – d^{r}_{ij}(\mathbf{X}))^2 \)は[…説明飛ばして…] \(w_{ij} = r^2 \delta_{ij}^{2r-2} \)とすると近似できる。
   • Ramsay(1977)のmultiscale損失関数 \( \sum_{i < j} ( \log(\delta_{ij}) – \log( d_{ij}(\mathbf{X})))^2 \)は[…スキップ…]\(w_{ij} = \delta_{ij}^{-2}\)で近似できる。
   • Summon(1969)の損失関数は \(w_{ij} = \delta_{ij}^{-1}\)である。[へえ… そうなのか…]
4. ウェイトを非類似性のべき乗にすることで、ウェイトを強調できる。[数値例を挙げて説明している。スキップ]
5. 特定のモデルのために特別なウェイトを与えることができる。
6. 非一意性を避けるためにウェイトを与えることができる。たとえば、フィレンツェの16家族の婚姻ネットワークのデータってあるじゃないですか(networkパッケージに入っている)。パッツィ家はサルヴィアティ家としか繋がっていない[パッツィ家ってメディチ家のロレンツォを暗殺しようとした奴ね…]。つまり、サルヴィアティ家から距離1であればどこに位置してもストレスが変わらないことになる。こういうときは、婚姻関係があったら\(\delta_{ij} = 1\), なかったら\(\delta_{ij} = 1/c\)とコーディングしておいて (\(c\)は0.01というような小さい値)、前者の場合には\(w_{ij} = 1\), 後者の場合には\(w_{ij} = c\)とするとよい。[な・る・ほ・ど… 最後の最後になって、ここはすごく勉強になった]

　… やれやれ、終わった。
　わからないなりに得るところも大きかった。SMACOFというのは要は数値最適化手法なんだから細かいことは気にしなくていいや、と思ってたんだけど、やっぱし勉強してみるもんだな、と。損失関数のなかにウェイトが入っているというのはなにかと便利であることがわかった。
　いっぽう、読んでかえってわかんなくなった面もあって… smacof::mds()でtype = “ordinal”と指定したとき、損失関数のなかの\(\hat{d}_{ij}\)は、入力された非類似性\(\delta_{ij}\)を単に順位に変換したものなの？ つまり、順位を量と見なして最小二乗フィッティングするわけ？ まじで？ KruskalのMDSだと、そういうときにはなんか工夫があったような気がするんだけどな… もう25年くらい前の知識だけど… (自分で書いてて信じられない)