Myrkgard logo
Series logo MATH001 Best Fitting Curves
Sat, 18 Nov 2017 05:00:00 +0000

Best Fitting Curves

In this episode we'll learn how to find the best fitting curve through a given set of \(x\), \(y\) values. At its core this means finding an approximate solution of an overdetermined system. We'll see how to find the best fitting straight line and, very similar to this, how to find the best fitting parabola. The solution presented relies on the Eigen C++ library's built-in function for computing the singular-value decomposition.

The generic problem behind this is called total least squares problem. The best fitting straight line is a linear regression and the best fitting parabola is a quadratic regression.

The Singular-Value Decomposition

For each real matrix \(\mathbf{A}\) of size \(m \times n\) and rank \(r\) there exists a unique decomposition \[\mathbf{A} = \mathbf{U} \boldsymbol{\Sigma} \mathbf{V}^T \] with \[\boldsymbol{\Sigma} = \left(\begin{array}{ccc|ccc} \sigma_1 & & & & & \\ & \ddots & & & 0 & \\ & & \sigma_r & & & \\ \hline & & & & & \\ & 0 & & & 0 & \\ & & & & & \\ \end{array}\right)\] and \[\sigma_1 \geq \cdots \geq \sigma_r > 0 .\]

\(\mathbf{U} \boldsymbol{\Sigma} \mathbf{V}^T\) is called the singular-value decomposition (SVD) of \(\mathbf{A}\).

Best Fitting Straight Line

Let's say we have a set of \(x\), \(y\) value pairs, e.g. measurement data from some lab experiment: \[\begin{pmatrix} x_1 & y_1 \\ \vdots & \vdots \\ x_n & y_n \\ \end{pmatrix}\] We want to find the best fitting straight line \[f(x) = a_1 x + a_0\] that passes through our \(x\), \(y\) value point cloud.

Step one is to create a matrix \(\mathbf{A}\) from the \(x\) values: \[\mathbf{A} = \begin{pmatrix} x_1 & 1 \\ \vdots & \vdots \\ x_n & 1 \\ \end{pmatrix} .\]

Next we create a vector \(\mathbf{b}\) from the \(y\) values: \[\mathbf{b} = \begin{pmatrix} y_1 \\ \vdots \\ y_n \\ \end{pmatrix} ,\]

Now we need to calculate the SVD of \(\mathbf{A}\). That will give us the matrices \(\mathbf{U}\), \(\boldsymbol{\Sigma}\) and \(\mathbf{V}\).

Using \(\mathbf{U}\) and \(\mathbf{b}\), compute \[\mathbf{d} = \mathbf{U}^T \mathbf{b} .\]

And finally calculate \[\boldsymbol{\alpha} = \begin{pmatrix} \alpha_1 \\ \vdots \\ \alpha_n \\ \end{pmatrix} = \mathbf{V} \begin{pmatrix} d_1 / \sigma_1 \\ \vdots \\ d_r / \sigma_r \\ \hline 0 \\ \end{pmatrix} .\]

This last vector \(\boldsymbol{\alpha}\) contains the coefficients of our best fitting straight line: \[f(x) = \alpha_1 x + \alpha_2 .\]

Best Fitting Parabola

Again, say we have a set of \(x\), \(y\) values: \[\begin{pmatrix} x_1 & y_1 \\ \vdots & \vdots \\ x_n & y_n \\ \end{pmatrix} .\] But this time we want to find the best fitting parabola \[f(x) = a_2 x^2 + a_1 x + a_0 .\] Everything works just like before. The only difference is that our matrix \(\mathbf{A}\) now contains one additional column: \[\mathbf{A} = \begin{pmatrix} x_1^2 & x_1 & 1 \\ \vdots & \vdots & \vdots \\ x_n^2 & x_n & 1 \\ \end{pmatrix} .\] And the function we are looking for evaluates to \[f(x) = \alpha_1 x^2 + \alpha_2 x + \alpha_3 .\]

Eigen's SVD Class

Eigen comes with its
JacobiSVD
template class. It has the name Jacobi in it to indicate which algorithm it uses internally. Following the straight line example above, code using
JacobiSVD
could look like the following snippet.
Eigen::JacobiSVD<Eigen::MatrixXd> svd = A.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV);
Eigen::VectorXd s = svd.singularValues();
Eigen::MatrixXd U = svd.matrixU();
Eigen::MatrixXd V = svd.matrixV();
Eigen::VectorXd d = U.transpose() * b;

Eigen::VectorXd tmp;
tmp.resize(2);
tmp(0) = d(0) / s(0);
tmp(1) = d(1) / s(1);

Eigen::VectorXd alpha = V * tmp;
Fortunately, as this is a very common problem, the Eigen developers provide a shortcut, the
solve()
function. So the code that's actually needed shrinks to a single line:
Eigen::VectorXd alpha = A.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(b);

Code Example, Straight Line

#include <iostream>
#include <Eigen/Dense>
#include <Eigen/SVD>

Eigen::VectorXd bestFitLinear(Eigen::MatrixXd const &data)
{
    auto const samples = data.rows();
    
    Eigen::MatrixXd A;
    A.resize(samples, 2);
    A.col(0) = data.col(0);
    A.col(1) = Eigen::VectorXd::Constant(samples, 1);
    
    Eigen::VectorXd b = data.col(1);

    return A.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(b);
}

int main()
{
    auto const samples = 8;
    Eigen::MatrixXd data;
    data.resize(samples, 2);
    //       x    y
    data << 17., 51.,
            19., 47.,
            22., 48.,
            24., 46.,
            26., 41.,
            27., 42.,
            29., 38.,
            31., 36.;

    std::cout << bestFitLinear(data) << std::endl;
    /* output:
       -1.03051
        68.7437
     */
}
plot

Code Example, Parabola

#include <iostream>
#include <Eigen/Dense>
#include <Eigen/SVD>

Eigen::VectorXd bestFitQuadratic(Eigen::MatrixXd const &data)
{
    auto const samples = data.rows();

    Eigen::MatrixXd A;
    A.resize(samples, 3);
    A.col(0) = data.col(0).array().pow(2);
    A.col(1) = data.col(0);
    A.col(2) = Eigen::VectorXd::Constant(samples, 1);

    Eigen::VectorXd b = data.col(1);  

    return A.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(b);
}

int main()
{
    auto const samples = 8;
    Eigen::MatrixXd data;
    data.resize(samples, 2);
    //        x    y
    data << 2.5, 173.,
            3.1, 170.,
            3.8, 122.,
            5.2, 128.,
            6.4, 151.,
            7.0, 194.,
            7.5, 268.,
            8.3, 297.;

    std::cout << bestFitQuadratic(data) << std::endl;
    /* output:
         13.747
       -126.428
        411.912
     */
}
plot