多項式フィッティング - わかばテクノロジ

最適化
多項式フィッティング; Polynomial Fitting

概要

多項式フィッティングとは独立変数\(x_i\)と従属変数\(y_i\)の組で表される\(N\)個のデータ\(i = 1 \cdots N \)があるとき、
\(n\)次多項式\(\tilde{y_i} = \displaystyle \sum_{k=0}^{n} r_k x_i^k\)と\(y_i\)との誤差の総和を最小化する、\(\boldsymbol{r} = (r_0 \cdots r_n)^T \)を求めることである。
誤差には2乗和が用いられる。
\(\quad \displaystyle E = \sum_{i=1}^{N} (y_i - \tilde{y_i})^2 \)

誤差最小値における\(\boldsymbol{r} \)の導出

最小化する対象である誤差の2乗和は以下である。
\(\begin{eqnarray} \quad \displaystyle E &=& \sum_{i=1}^{N} (y_i - \sum_{k=0}^{n} r_k x_i^k)^2 \\ \quad &=& | \boldsymbol{y} - M \boldsymbol{r} |^2 \quad \left( \boldsymbol{y} = (y_1 \cdots y_N)^T, \ \ M_{ki} = x_i^k \right) \\ \quad &=& ( \boldsymbol{y} - M \boldsymbol{r})^T (\boldsymbol{y} - M \boldsymbol{r}) \\ \quad &=& ( \boldsymbol{y}^T - \boldsymbol{r}^T M^T) (\boldsymbol{y} - M \boldsymbol{r}) \\ \quad &=& ( \boldsymbol{y}^T \boldsymbol{y} -\boldsymbol{y}^T M \boldsymbol{r} - \boldsymbol{r}^T M^T \boldsymbol{y} + \boldsymbol{r}^T M^T M \boldsymbol{r}) \\ \quad &=& ( \boldsymbol{y}^T \boldsymbol{y} - (\boldsymbol{y}^T M \boldsymbol{r})^T - \boldsymbol{r}^T M^T \boldsymbol{y} + \boldsymbol{r}^T M^T M \boldsymbol{r}) \quad \cdots \boldsymbol{y}^T M \boldsymbol{r} \in \mathbb{R}\\ \quad &=& ( \boldsymbol{y}^T \boldsymbol{y} -\boldsymbol{r}^T M^T \boldsymbol{y} - \boldsymbol{r}^T M^T \boldsymbol{y} + \boldsymbol{r}^T M^T M \boldsymbol{r}) \\ \quad &=& ( \boldsymbol{y}^T \boldsymbol{y} - 2 \boldsymbol{r}^T M^T \boldsymbol{y} + \boldsymbol{r}^T M^T M \boldsymbol{r}) \end{eqnarray} \)
\(E\)が最小化であるとき、\(\boldsymbol{r}\)についての偏微分は\(\boldsymbol{0}\)になる。
\( \quad \displaystyle \frac{\partial E}{\partial \boldsymbol{r}} = - 2 M^T \boldsymbol{y} + 2 M^T M \boldsymbol{r} = \boldsymbol{0} \\ \quad \rightarrow M^T M \boldsymbol{r} = M^T \boldsymbol{y} \\ \quad \rightarrow \boldsymbol{r} = (M^T M)^{-1} M^T \boldsymbol{y} \\ \)
\((M^T M)^{-1}\)ではパラメータの数(\(n+1\))の正方行列の逆行列を計算することになる。
原点を通る多項式を得たい場合は\(\displaystyle \sum_{k=1}^{n} r_k x_i^k\)としてパラメータの数を1減らし、その後の導出は同様になる。

ソースコード


namespace DataFitting {

    /// <summary>多項式フィッティング</summary>
    public class PolynomialFittingMethod : FittingMethod {
        /// <summary>コンストラクタ</summary>
        public PolynomialFittingMethod(FittingData[] data_list, int degree, bool is_enable_section)
            : base(data_list, degree + (is_enable_section ? 1 : 0)) {

            this.Degree = degree;
            this.IsEnableSection = is_enable_section;
        }

        /// <summary>次数</summary>
        public int Degree {
            get; private set;
        }

        /// <summary>y切片を有効にするか</summary>
        public bool IsEnableSection { get; set; }

        /// <summary>フィッティング値</summary>
        public override double FittingValue(double x, Vector coefficients) {
            if(IsEnableSection) {
                double y = coefficients[0], ploy_x = 1;

                for(int i = 1; i < coefficients.Dim; i++) {
                    ploy_x *= x;
                    y += ploy_x * coefficients[i];
                }

                return y;
            }
            else {
                double y = 0, ploy_x = 1;

                for(int i = 0; i < coefficients.Dim; i++) {
                    ploy_x *= x;
                    y += ploy_x * coefficients[i];
                }

                return y;
            }
        }

        /// <summary>フィッティング</summary>
        public Vector ExecureFitting() {
            Matrix m = new Matrix(data_list.Length, ParametersCount);
            Vector b = Vector.Zero(data_list.Length);

            if(IsEnableSection) {
                for(int i = 0; i < data_list.Length; i++) {
                    double x = data_list[i].X;
                    b[i] = data_list[i].Y;

                    m[i, 0] = 1;

                    for(int j = 1; j <= Degree; j++) {
                        m[i, j] = m[i, j - 1] * x;
                    }
                }
            }
            else {
                for(int i = 0; i < data_list.Length; i++) {
                    double x = data_list[i].X;
                    b[i] = data_list[i].Y;

                    m[i, 0] = x;

                    for(int j = 1; j < Degree; j++) {
                        m[i, j] = m[i, j - 1] * x;
                    }
                }
            }

            return (m.Transpose * m).Inverse * m.Transpose * b;
        }
    }
}

単体テスト


namespace DataFitting.Tests {
    [TestClass()]
    public class PolynomialFittingMethodTests {
        [TestMethod()]
        public void ExecureFittingTest() {
            double[] x_list = { 1, 3, 4, 7, 8, 9, 13, 15, 20 };
            Vector p1 = new Vector(2, -1, 1, 5), p2 = new Vector(4, 3, -1);
 
            FittingData[] data_list1 = new FittingData[x_list.Length], data_list2 = new FittingData[x_list.Length];

            for(int i = 0; i < x_list.Length; i++) {
                double x = x_list[i];

                data_list1[i].X = data_list2[i].X = x;
                data_list1[i].Y = p1[0] + p1[1] * x + p1[2] * x * x + p1[3] * x * x * x;
                data_list2[i].Y = p2[0] * x + p2[1] * x * x + p2[2] * x * x * x;
            }
            
            PolynomialFittingMethod fitting1 = new PolynomialFittingMethod(data_list1, 3, true);
            PolynomialFittingMethod fitting2 = new PolynomialFittingMethod(data_list2, 3, false);

            Assert.AreEqual((fitting1.ExecureFitting() - p1).Norm, 0, 1e-8);
            Assert.AreEqual((fitting2.ExecureFitting() - p2).Norm, 0, 1e-8);
        }
    }
}