独立成分分析(ICA)をC++で実装してみる

## 初めに
独立成分分析（ICA: Independent Component Analysis) をC++で実装してみようという記事です．
間違っているところがあれば，コメントお願いします．

(３年くらい前から温めていて，Githubにも載せていなかったのですが，本日高熱出して休養中なので，これを供養します)

## 独立成分分析
### ICAって何？
独立成分分析（ICA: Independent Component Analysis) とは，ざっくり言えば「混ざり合った信号」から「信号源から発せられる信号に別々に分離」することができる手法です．
初めに,以下のような信号があったとします[1]．

$$
s(t) = (s_1(t), s_2(t), ..., s_n(t))^t, t=0,1,2,... 
$$

これらの信号は，干渉はしないかつ，各信号は統計的に独立であるとします．
センサなどから，得られる信号を以下のように定義します．

$$
x(t) = (x_1(t), x_2(t), ..., x_n(t))^t, t=0,1,2,... 
$$

そうすると，これらは，何かしらの行列$A$で混合されているはずです．そうなると以下のように仮定できます．

$$
x(t) = As(t)
$$

ICAでは，この行列$A$を知らずに，分離しようというものです．
ただし, 元の信号$s(t)$について

- 平均は0
- モーメントは適当な次元まで発散しない

である必要があります．

今回は，復元行列$W$というものを仮定して，これを求めていくようなプログラムを書いていきたいと思います．
つまり

$$
s(t) = W^tx(t)
$$

というような$W$を求めていきます．

## FastICA
ICAは一般的に，FastICA([2], [3])と呼ばれるアルゴリズムが使われるらしいです．とりあえずWikiに書いてある方法で実装してみます．

### 高速化？
今回は行列演算を多く行うので，速度の関係からできる限りDSP命令などを使いたいところです．めんどくさいのでライブラリを使うことにします．特に，固有値分解なんかやりたくないですね....

高貴な感じで行くなら，CBLAS/LAPACKだと思います？が，今回は手軽さをもとめてC++の行列計算といえばのEigenを選びました．

## C++での実装
### データのセンタリング
データのセンタリング(中心化)をします．各信号から平均値を引くだけです．
データのレンジがバラバラだと計算がしにくい(深層学習でもデータの標準化はよくやりますよね)のと，平均値は二次統計量です．今回求めたい高次元統計量における独立性には，直接寄与しないはずですので，取り除いてしまっていると理解しています．
（間違っていたら教えてください）

Eigenだとnumpyみたいに簡単に書けます．
```cpp
void __centering() {
    mean_ = x_.rowwise().mean();
    for (int i = 0; i < rows_; i++) {
        x_.row(i).array() -= mean_(i);
    }
}
```

### データのホワイトニング(白色化)
で次に白色化をして，無相関化します．
白色雑音とか思い受けべてもらうといいかも知れません．
ホワイトノイズは全ての周波数成分が等しいパワーで含まれる雑音のことで，その統計的性質として共分散行列が単位行列になること・無相関が特徴だったと思います．

ICAでは非ガウス性を最大化することで，独立成分を見つけます．したがって，白色化して二次統計量を無くしてしまうことで，相関のない状態にし，残された依存性(高次元統計量)から独立成分を見つけるということです．（多分そういうこと）

なんかうだうだ言いましたが，やっていることは単純で，分散共分散行列を計算して，その固有値と固有ベクトルを求めています．

手順ですが，以下の通りです．

```cpp
void __whitening() {
    // 共分散行列の計算: Cov = XX^T / (n-1)
    conv_ = (x_ * x_.transpose()) / (cols_ - 1);

    // 固有値分解: Cov = EDE^T (E=>固有ベクトル行列, D=>固有値の対角行列)
    Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eigen_solver(conv_);
    eigen_ = eigen_solver.eigenvalues();
    e_ = eigen_solver.eigenvectors();

    // 固有値の平方根の逆数を対角行列に変換 => D^(-1/2)
    d_sq_inv_ = Eigen::MatrixXd::Zero(rows_, rows_);
    for (int i = 0; i < rows_; i++) {
        d_sq_inv_(i, i) = 1.0 / std::sqrt(eigen_(i));
    }

    // ホワイトニングの行列=> W = D^(-1/2)E^T
    Eigen::MatrixXd whitening_mat_ = d_sq_inv_ * e_.transpose();

    // で，最後にかけて，ホワイトニングする
    z_ = whitening_mat_ * x_;
}
```

### 独立成分の抽出
で，ここあんまりよくわかっていない部分もあるのですが，不動点アルゴリズムを使用して，非ガウス性を最大化しています．

各成分を反復的に求めて，その前に見つかった成分に対して直交するようにします．そうすると，複数の独立成分が同じ方向に収束することを防ぎ，白色化空間内で互いに無相関な方向が得られるため，互いに独立した成分を順番に抽出することができます．（多分そういうことだと思う）

非線形関数として，`tanh`を使っていますが，ほかにもいろいろあるみたいです．

```cpp

void __ica() {
    w_ = Eigen::MatrixXd(sigs_, rows_);

    for (int c = 0; c < sigs_; c++) {
        // ランダムな初期重みベクトル
        Eigen::VectorXd w = Eigen::VectorXd::Random(rows_);
        w.normalize();

        // 反復計算
        for (int iter = 0; iter < max_iter_; iter++) {
            Eigen::VectorXd w_prv = w;

            // g(x) = tanh(x)
            Eigen::VectorXd wx = w.transpose() * z_;
            Eigen::VectorXd g_wx(cols_);
            Eigen::VectorXd g_p_wx(cols_);

            for (int j = 0; j < cols_; j++) {
                // tanhとその導関数tanh'(x) = 1 - tanh^2(x)を適用する
                g_wx(j) = std::tanh(wx(j));
                g_p_wx(j) = 1.0 - std::pow(std::tanh(wx(j)), 2);
            }

            // 重みベクトルの更新
            w = (z_ * g_wx) / cols_ - (g_p_wx.mean() * w);

            // ベクトルを直交化
            if (c > 0) {
                Eigen::MatrixXd proj = Eigen::MatrixXd::Zero(rows_, rows_);
                for (int i = 0; i < c; i++) {
                    proj += w_.row(i).transpose() * w_.row(i);
                }
                w = w - proj * w;
            }

            w.normalize();

            // 距離見て, toleより小さいなら収束したと判定
            double d = std::abs(1.0 - std::abs(w.dot(w_prv)));
            if (d < tole_) {
                break;
            }

            if (iter == max_iter_ - 1) {
                // TODO: warn log
                std::cerr << "Component " << c + 1 << " did not converge."
                          << std::endl;
            }
        }

        // 収束したやつをとっておく
        w_.row(c) = w.transpose();
    }
}
```

### 計算結果
上記までにも止まったものを掛け合わせればOKです！
```cpp
y_ = w_ * z_;
```
ちなみにも止まった混合行列は
```cpp
// 固有値の平方根を対角行列に変換=> D^(1/2)
Eigen::MatrixXd d_sq = Eigen::MatrixXd::Zero(rows_, rows_);
for (int i = 0; i < rows_; i++) {
    d_sq(i, i) = std::sqrt(eigen_(i));
}

// 混合行列の推定=> A = ED^(1/2)W^(-1)
// だけど，この行列は元の混合行列の推定値であり，スケーリングと順序が異なる場合がある (ICAはこうゆうもん)
transform_ = e_ * d_sq * w_.inverse();
```
とすれば良いです．

### 完成したやつ
殴り書きしたコードなので汚いですが，一応公開します. 
[https://github.com/KikuchiTomo/FastICA](https://github.com/KikuchiTomo/FastICA)

### 開発環境
- OS : MacOS
- コンパイラ : G++ (と思ったけどMacなのでapple clangです...エイリアスはるな!)
- C++ : c++11

## 実行結果
### 元の信号
![original.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/431953/264102bf-9ce8-42b0-add5-840a63a4156b.png)

### 混ぜたやつ
![maxed.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/431953/fa037fb5-75b9-4347-a4fe-abca57e5caa2.png)

### 復元信号
![recomp.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/431953/80384a3f-bea3-4d2f-bf4b-6f34639d0c6d.png)

なんか良さそうですね．

〜 FIN 〜


## 参考文献
[1] 村田昇, 入門　独立成分分析, 東京電機大学出版局，pp28-30, 2005．
[2] Hyvärinen, A.; Oja, E. (2000). "Independent component analysis: Algorithms and applications" (PDF). Neural Networks. 13 (4–5): 411–430. CiteSeerX 10.1.1.79.7003. doi:10.1016/S0893-6080(00)00026-5. PMID 10946390.
[3] Hyvarinen, A. (1999). "Fast and robust fixed-point algorithms for independent component analysis" (PDF). IEEE Transactions on Neural Networks. 10 (3): 626–634. CiteSeerX 10.1.1.297.8229. doi:10.1109/72.761722. PMID 18252563.
[4] FastICA, Wikipedia, 2022/06/28, https://en.wikipedia.org/wiki/FastICA#cite_note-Hyvarinen-1