【C++】【Eigen】Whittaker Smootherを実装する

2024-04-16C++,Eigen,スペクトル,信号処理,数学

スペクトル等のノイズを含むデータをスムージングするために広く使われる手法の1つにSavitzky-Golay法があります。このSavitzky-Golay法と同様に強力なスムージング手法としてWhittaker Smootherがあります。このWhittaker SmootherをC++とEigenを使って実装してみたので紹介します。

Savitzky-Golay法についてはこちら。

Whittaker Smootherとは

Whittaker Smootherは1922年（100年前！）にWhittakerが発表[1]したスムージング手法です。
その後Eilersが2003年の論文[2]で紹介したのち、AsLS法というベースラインの推定手法として本法を応用した手法を発表すると、分光スペクトルやクロマトグラムのベースライン推定手法として注目され、Whittaker Smootherを応用したさまざまなベースライン推定手法が発表されるようになりました。

今回は、Eilersの論文を元にWhittaker SmootherをC++を使って実装してみました。

特徴

  1. 比較的短いコードで実現できる
  2. 疎行列の計算に対応していれば、多数の観測値に対しても比較的短時間でスムージングできる
  3. 欠落データに対応できる
  4. 滑らかさの制御ができる

原理

観測された測定データの列を$\boldsymbol{y} = (y_{1}, \ldots, y_{m})$とします。この値には通常ノイズなどが乗っていることから、データ列の作る線は滑らかな曲線とはなりません。

次に、滑らかなデータ列$\boldsymbol{z} = (z_{1}, \ldots, z_{m})$を考え、$\boldsymbol{z}$を$\boldsymbol{y}$に適合させることを考えます。$\boldsymbol{z}$を$\boldsymbol{y}$に適合させるためには、2つの要素を考慮する必要があります。1つは、$\boldsymbol{y}$に対する$\boldsymbol{z}$の適合度合い、もう1つは$\boldsymbol{z}$の粗さです。
前者は$\boldsymbol{y}$と$\boldsymbol{z}$の差の二乗和、後者は$\boldsymbol{z}$の$s$回差分の二乗和$\sum_{i}(\Delta^{s}z_{i})^{2}$で表すことができます。$\boldsymbol{z}$が滑らかであるほど、$\boldsymbol{y}$からは乖離するため、両者は相反する要素となります。すなわち、ペナルティ付き最小二乗法の問題となります。

$$
\displaystyle
Q = \sum_{i}w_{i}(y_{i} – z_{i})^{2} + \lambda\sum_{i}(\Delta^{s}z_{i})^{2}
\tag{1}
$$

ここで$w_{i}$は重みを表し、$0 \leqq w_{i} \leqq 1$です。欠落データや外れ値などに対して、$0$または$0$に近い値を与えます。また、$\lambda$はペナルティ係数です。$\lambda$を大きくすれば粗さの項が優勢となり、より滑らかな曲線となります。

(1)を行列とベクトルで表すと(2)になります。

$$
\displaystyle
Q = (\boldsymbol{y}-\boldsymbol{z})^{\mathrm{T}}W(\boldsymbol{y}-\boldsymbol{z}) + \lambda \boldsymbol{z}^{\mathrm{T}}(D_{s}^{\mathrm{T}}D_{s})\boldsymbol{z}
\tag{2}
$$

ここで、$W$は重み$\boldsymbol{w} = (w_{1}, \ldots, w_{m})$を対角に配置した対角行列、$D_{s}$は$D_{s}\boldsymbol{z} = \Delta^{s}\boldsymbol{z}$となるような行列です。たとえば、$m = 5, s = 1$では、

$$
\displaystyle
D_{s} = \left(
\begin{array}{rrrrr}
-1 & 1 & 0 & 0 & 0 \\
0 & -1 & 1 & 0 & 0 \\
0 & 0 & -1 & 1 & 0 \\
0 & 0 & 0 & -1 & 1
\end{array}
\right)
$$

$m = 7, s = 2$では、

$$
\displaystyle
D_{s} = \left(
\begin{array}{rrrrrrr}
1 & -2 & 1 & 0 & 0 & 0 & 0 \\
0 & 1 & -2 & 1 & 0 & 0 & 0 \\
0 & 0 & 1 & -2 & 1 & 0 & 0 \\
0 & 0 & 0 & 1 & -2 & 1 & 0 \\
0 & 0 & 0 & 0 & 1 & -2 & 1
\end{array}
\right)
$$

$m = 9, s = 3$では、

$$
\displaystyle
D_{s} = \left(
\begin{array}{rrrrrrrrr}
-1 & 3 & -3 & 1 & 0 & 0 & 0 & 0 & 0 \\
0 & -1 & 3 & -3 & 1 & 0 & 0 & 0 & 0 \\
0 & 0 & -1 & 3 & -3 & 1 & 0 & 0 & 0 \\
0 & 0 & 0 & -1 & 3 & -3 & 1 & 0 & 0 \\
0 & 0 & 0 & 0 & -1 & 3 & -3 & 1 & 0 \\
0 & 0 & 0 & 0 & 0 & -1 & 3 & -3 & 1
\end{array}
\right)
$$

となります。(2)を$\boldsymbol{z}$で偏微分して、$0$と置くと、

$$
\frac{\partial Q}{\partial \boldsymbol{z}} = -2W(\boldsymbol{y} – \boldsymbol{z}) + 2 \lambda(D_{s}^{\mathrm{T}}D_{s})\boldsymbol{z} = 0
$$

$$
(W + \lambda D_{s}^{\mathrm{T}} D_{s}) \boldsymbol{z} = W \boldsymbol{y}
$$

よって、

$$
\boldsymbol{z} = (W + \lambda D_{s}^{\mathrm{T}} D_{s})^{-1} W \boldsymbol{y}
\tag{3}
$$

で滑らかなデータ列$\boldsymbol{z}$が求められます。

Whittakerのオリジナルのものは$s = 3$の場合ですが、$s = 1, 2$の場合も考えることができます。$s = 1$の時はBohlmann Smoothing、$s = 2$の時はHodrick-Prescott Smoothing、$s = 3$の時はWhittaker-Henderson Smoothingとも呼ばれることもあります[3]。ここではひっくるめてWhittaker Smootherと呼ぶことにします。

Eilersによると、$s = 2$の時はうまく機能することが多いようです。

参考文献

  1. Whittaker, E. T. On A New Method of Graduation. Proc. Edinburgh Math. Soc.1922, 41, 63-75.
  2. Eilers, P. H. C. A Perfect Smoother. Anal. Chem.2003, 75(14), 3631–3636.
  3. Orfanidis S. J. Whittaker-Henderson Smoothing, 2018. Applied Optimum Signal Processing. http://eceweb1.rutgers.edu/~orfanidi/aosp/ (accessed Mar 24, 2022)

Eigenによる実装

Whittaker SmootherをEigenを使ってC++で実装します。その前に、Eigenには行列$D_{s}$を作成する関数がありません。MatlabやPython, Rでは相当するdiff()という関数があるため、これも実装する必要があります。

Matlabのdiff()に相当する関数は、以下のエントリで作成していますので、こちらも合わせてご覧ください。

ここで作成したdiff.hを同じディレクトリに置いているものとして実装します。エラーチェックはしていないのでご注意を。コードをGitHubにも公開しましたので、興味のある方はご覧ください。

実装例

  • whittaker.h

    #ifndef __EIGEN_WHITTAKER_H__
#define __EIGEN_WHITTAKER_H__

#include <cassert>
#include "diff.h"

// Whittaker Smootherを実施する
//
// y      : 実際の観測データ
// w      : 各点に対応する重みベクトル(0 <= w <= 1)
// lambda : 平滑化のためのペナルティ係数
// s      : 差分の階数(s = 1, 2, 3)
const Eigen::VectorXd WhittakerSmoother(
  const Eigen::VectorXd & y, const Eigen::VectorXd & w, const double lambda, const unsigned int s
);

#endif // __EIGEN_WHITTAKER_H__

  • whittaker.cpp

    #include "whittaker.h"

const Eigen::VectorXd WhittakerSmoother(
  const Eigen::VectorXd & y, const Eigen::VectorXd & w, const double lambda, const unsigned int s
)
{
  size_t m = y.size();

  Eigen::VectorXd z;
  Eigen::SparseMatrix<double> I, W, D;
  Eigen::SimplicialCholesky< Eigen::SparseMatrix<double> > solver;

  I.resize(m, m);
  I.setIdentity();
  D = Diff(I, s);
  D = lambda * (D.transpose() * D);
  W = w.asDiagonal();

  // W + lambda * D' * D は対称行列なので、Cholesky分解で連立方程式を解く
  solver.compute(W + D);

  assert(solver.info() == Eigen::Success);

  // (W + lambda * D' * D) * z = W * y を解く
  z = solver.solve(W * y);

  assert(solver.info() == Eigen::Success);

  return z;
}

使用例

上記のコードの使用例を示します。コマンドラインで指定したタブ区切り形式のテキストファイルを読み込んだのち、Whittaker Smootherによるスムージングを行った結果を元データとともに標準出力に出力します。

#include <fstream>
#include <iostream>
#include <string>
#include <vector>

#include "whittaker.h"

int main(int argc, char * argv[])
{
    double tx, ty, lambda;
    std::string header_x, header_y;
    std::vector<double> x, y, whittaker_y;

    if(argc < 2){
        return 1;
    }

    std::fstream fin;

    fin.open(argv[1]);

    if(!fin){
        std::cerr << "file not open." << std::endl;
        return 1;
    }

    fin >> header_x >> header_y;

    while(1){
        fin >> tx >> ty;
        if(fin.eof()){
          break;
        }
        x.push_back(tx);
        y.push_back(ty);
    }

    fin.close();

    Eigen::VectorXd vy = Eigen::Map<Eigen::VectorXd>(y.data(), (int)y.size());
    Eigen::VectorXd w  = Eigen::VectorXd::Ones((int)y.size());
    lambda = 1;

    Eigen::VectorXd vz1 = WhittakerSmoother(vy, w, lambda, 1);
    Eigen::VectorXd vz2 = WhittakerSmoother(vy, w, lambda, 2);
    Eigen::VectorXd vz3 = WhittakerSmoother(vy, w, lambda, 3);

    std::cout << header_x << "\t" << header_y;
    std::cout << "\ts=1\ts=2\ts=3" << std::endl;

    for(size_t i = 0; i < x.size(); i++){
        std::cout << x[i] << "\t" << y[i] << "\t";
        std::cout << vz1(i) << "\t" << vz2(i) << "\t" << vz3(i) << std::endl;
    }

    return 0;
}

動作を確認した環境

  • OS
    1. Windows 10
      1. Microsoft C/C++ Optimizing Compiler 19.28.29337 for x64
      2. MinGW-W64 GCC version 8.1.0
      3. Clang version 11.0.1
    2. Linux
      1. GCC version 9.3.0
  • Eigen version 3.3.9

スムージングの例

上記コードを使ってノイズを含む疑似スペクトルをスムージングし、出力データを図示してみました。グラフはJavaScriptのChart.js+chartjs-plugin-zoomを使っています。$x$軸、$y$軸上でマウスのホイールを回転させることで拡大・縮小を、ドラッグでスライドさせることが可能です。また、凡例のラベルをクリックすることで、当該曲線の表示・非表示を変更できます。

平滑化後の元のデータはこちらにあります。








まとめ

Whittaker Smootherについて、原理を含め紹介しました。Savitzky-Golay法と比べると、あまり知られていませんが、比較的短いコードで同様に強力なスムージングをかけられる手法です。

C++,Eigen,スペクトル

Posted by izadori