スパースモデリング (3) MATLABで簡単スパースモデリング

# カーブフィッティング
 
手軽にスパースモデリングを体験できる題材として，今回はカーブフィッティングを取り上げます．書籍「<a href="http://amzn.asia/f8DcJhd">スパースモデリング</a>」に掲載されている例題を詳しく解説します．

# 多項式補間の基本
カーブフィッティングで使われる最も単純なテクニックは<font color="red">**多項式補間**</font>です．点と点の間を結ぶ曲線を求めるときに使います．

多項式補間では以下の性質が重要です．
**2点が与えられれば，それらの点を通る直線（1次式）がただ一つ決まる．**
<img width=50%, src="https://qiita-image-store.s3.amazonaws.com/0/202074/833e55c9-bd2d-ba6d-dcf4-891cc06ffa5e.png">
**3点が与えられれば，それらの点を通る曲線（2次曲線）がただ一つ決まる．**
<img width=50%, src="https://qiita-image-store.s3.amazonaws.com/0/202074/78a11a48-bbc2-1452-9873-4f5c559a3d44.png">
**4点が与えられれば，それらの点を通る曲線（3次曲線）がただ一つ決まる．**
<img width=50%, src="https://qiita-image-store.s3.amazonaws.com/0/202074/e7c1c7d6-2d54-9fa6-5cb9-973903dd4a68.png">

以上から，次の一般化が可能です．

>**N個の点が与えられれば，それらの点を通る曲線（N-1次曲線）がただ一つ決まる．**

ここで，N-1次曲線というのは，N-1次の多項式で表される曲線のことです．

**問題：**80次曲線（80次の多項式で表される曲線）を決めるのに，何個の点が必要でしょうか？

上の原則から言えば，81個の点が必要ですよね．これを<font color="red">**たった11点から求められますよ**</font>，というのがスパースモデリングです．ビッグデータならぬ，<font color="red">**スモールデータ（足りないデータ）**</font>から予測しようというわけです．


# スパースモデリングの問題
次の多項式
$$
 y(t) = -t^{80}+t
$$
を考えます．この多項式に
$$
 t = 0,~ 0.1,~ 0.2,\ldots,~ 0.9,~ 1
$$
の11点を代入して11組のデータを求めます．これは，以下のような値になります．

|0|0.1|0.2|0.3|0.4|0.5|0.6|0.7|0.8|0.9|1|
|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|
|0|0.1|0.2|0.3|0.4|0.5|0.6|0.7|0.8|0.9|0|

<details>
$t$の値が1未満のとき，$t^{80}$は極めて小さい値になりますので，
$-t^{80}+t\approx t$となり，表のような結果になります．
</details>

図示すると，以下のようなデータとなります．

<img width=50%, src="https://qiita-image-store.s3.amazonaws.com/0/202074/554c41a4-8c43-c34a-1d15-ed191936549b.png">

なお点線はもとの80次多項式曲線を表しています．

**この11点のデータだけから，もとの80次多項式（点線）を推定しなさい**，というのがここでの問題です．

# MATLAB で解いてみよう
上の問題をMATLABで解いてみましょう．

まずは，多項式の定義から．多項式の係数からなるベクトルを定義します．

```matlab
%% polynomial coefficients
x_orig = [-1,zeros(1,78),1,0]';
```

ここで，`zeros(1,78)`は1行78列（横ベクトル）のゼロだけのベクトルです．$-t^{80}$と$t$の間の係数はすべて0であることに注意してください．

この多項式から，11点のデータを生成します．

```matlab
%% sampling
t = 0:0.1:1;
y = polyval(x_orig,t);
```

`polyval`はMATLAB関数で，多項式の係数ベクトルから，多項式の値を計算する関数です．

### 10次多項式によるカーブフィッティング
まずは，11点を通る10次多項式を求めてみましょう．

```matlab
%% Vandermonde matrix
Phi_v = vander(t);
%% Interpolation polynomial with order 10
x_i = inv(Phi_v)*y';
```
ここで，`vander(t)`というのは，ヴァンデルモンド行列と呼ばれる行列で，これの逆行列を観測データに掛ければ，補間多項式の係数が出てきます．
<details>
詳しくは，書籍「<a href="http://amzn.asia/f8DcJhd">スパースモデリング</a>」の第2.1.2節をご参照ください．
</details>
グラフを表示してみましょう．

```matlab
%% Plot
tt = 0:0.01:1;
figure;
plot(tt,polyval(x_i,tt),'-');
hold on
plot(tt,polyval(x_orig,tt),'--');
stem(t,y);
```
<img width=50%, src="https://qiita-image-store.s3.amazonaws.com/0/202074/19d1d8e1-7c98-7f4d-4514-8d14089115ac.png">

特に，t=0.9〜1.0の間で，もとの80次曲線とずれてしまっています．

### スパースモデリングによる方法
では，スパースモデリングのテクニックを使って，80次曲線を求めてみましょう．

まず補間条件（求める曲線が与えられた点の上を通るための条件）を求めます．

```matlab
% Order of polynomial
M_l = 80;
% Vandermonde matrix
Phi_l=[];
for m=0:M_l
 Phi_l = [t'.^m,Phi_l];
end
```

11点から80次多項式を求めるために，ヴァンデルモンド行列は$11\times 81$の横長の行列になっていることに注意してください．横長の行列なので，当然，逆行列は存在しません．より詳しく言うと，連立方程式

```math
\Phi_l x = y
```

には無数の解が存在します．（$\Phi_l$は$11\times 81$のヴァンデルモンド行列です）．

スパースモデリングは，上の等式を満たす無数の解の中から最もスパースなものを選びます．すなわち，非ゼロの係数が最も少ない多項式を見つけます．

CVXで解いてみましょう．

```matlab
% CVX
cvx_begin
 variable x(M_l+1)
 minimize norm(x,1)
 subject to 
   Phi_l*x == y'
cvx_end
```

簡単に答えが得られて，

```matlab
>> x

x =

   -1.0000
   -0.0000
   -0.0000
   -0.0000
   -0.0000
   -0.0000
   -0.0000
   -0.0000
   -0.0000
   -0.0000
   -0.0000
   -0.0000
   -0.0000
   -0.0000
   -0.0000
   -0.0000
   -0.0000
   -0.0000
   -0.0000
   -0.0000
   -0.0000
   -0.0000
   -0.0000
   -0.0000
   -0.0000
   -0.0000
   -0.0000
   -0.0000
   -0.0000
   -0.0000
   -0.0000
   -0.0000
   -0.0000
   -0.0000
   -0.0000
   -0.0000
   -0.0000
   -0.0000
   -0.0000
   -0.0000
   -0.0000
   -0.0000
   -0.0000
   -0.0000
   -0.0000
   -0.0000
   -0.0000
   -0.0000
   -0.0000
   -0.0000
   -0.0000
    0.0000
    0.0000
    0.0000
    0.0000
    0.0000
    0.0000
    0.0000
    0.0000
    0.0000
    0.0000
    0.0000
    0.0000
    0.0000
    0.0000
    0.0000
    0.0000
    0.0000
    0.0000
    0.0000
    0.0000
   -0.0000
   -0.0000
   -0.0000
   -0.0000
    0.0000
   -0.0000
    0.0000
   -0.0000
    1.0000
    0.0000
```
となります．１番目の係数（$t^{80}$の係数）が$-1$，８０番目の係数（$t$の係数）が$1$となりました．正確にもとの多項式を復元しています．

グラフを描いてみましょう．

```matlab
%% Plot
tt = 0:0.01:1;
figure;
plot(tt,polyval(x,tt),'-');
hold on
plot(tt,polyval(x_orig,tt),'--');
stem(t,y);
```
<img width=50%, src="https://qiita-image-store.s3.amazonaws.com/0/202074/4f8c2138-7ee9-d474-6d5a-107c72397971.png">

完全に一致しています．

----------

スパースモデリングに興味を持った方は，ぜひ，書籍「**スパースモデリング- 基礎から動的システムへの応用 -**」でその仕組みを学んでください．

<img src="https://qiita-image-store.s3.amazonaws.com/0/202074/0c2267dc-c99f-df02-548f-f01a5c26027c.jpeg" width=25%>
[出版社ページ](http://www.coronasha.co.jp/np/isbn/9784339032222/)
[Amazonページ](https://www.amazon.co.jp/gp/product/4339032220/)