Non-linear regression in C#

I used the MathNet.Iridium release because it is compatible with .NET 3.5 and VS2008. The method is based on the Vandermonde matrix. Then I created a class to hold my polynomial regression

using MathNet.Numerics.LinearAlgebra;

public class PolynomialRegression
{
    Vector x_data, y_data, coef;
    int order;

    public PolynomialRegression(Vector x_data, Vector y_data, int order)
    {
        if (x_data.Length != y_data.Length)
        {
            throw new IndexOutOfRangeException();
        }
        this.x_data = x_data;
        this.y_data = y_data;
        this.order = order;
        int N = x_data.Length;
        Matrix A = new Matrix(N, order + 1);
        for (int i = 0; i < N; i++)
        {
            A.SetRowVector( VandermondeRow(x_data[i]) , i);
        }

        // Least Squares of |y=A(x)*c| 
        //  tr(A)*y = tr(A)*A*c
        //  inv(tr(A)*A)*tr(A)*y = c
        Matrix At = Matrix.Transpose(A);
        Matrix y2 = new Matrix(y_data, N);
        coef = (At * A).Solve(At * y2).GetColumnVector(0);
    }

    Vector VandermondeRow(double x)
    {
        double[] row = new double[order + 1];
        for (int i = 0; i <= order; i++)
        {
            row[i] = Math.Pow(x, i);
        }
        return new Vector(row);
    }

    public double Fit(double x)
    {
        return Vector.ScalarProduct( VandermondeRow(x) , coef);
    }

    public int Order { get { return order; } }
    public Vector Coefficients { get { return coef; } }
    public Vector XData { get { return x_data; } }
    public Vector YData { get { return y_data; } }
}

which then I use it like this:

using MathNet.Numerics.LinearAlgebra;

class Program
{
    static void Main(string[] args)
    {
        Vector x_data = new Vector(new double[] { 0, 1, 2, 3, 4 });
        Vector y_data = new Vector(new double[] { 1.0, 1.4, 1.6, 1.3, 0.9 });

        var poly = new PolynomialRegression(x_data, y_data, 2);

        Console.WriteLine("{0,6}{1,9}", "x", "y");
        for (int i = 0; i < 10; i++)
        {
            double x = (i * 0.5);
            double y = poly.Fit(x);

            Console.WriteLine("{0,6:F2}{1,9:F4}", x, y);
        }
    }
}

Calculated coefficients of [1,0.57,-0.15] with the output:

    x        y
 0.00   1.0000
 0.50   1.2475
 1.00   1.4200
 1.50   1.5175
 2.00   1.5400
 2.50   1.4875
 3.00   1.3600
 3.50   1.1575
 4.00   0.8800
 4.50   0.5275

Which matches the quadratic results from Wolfram Alpha.

Edit 1 To get to the fit you want try the following initialization for x_data and y_data:

Matrix points = new Matrix( new double[,] {  {  1, 82.96 }, 
               {  2, 86.23 }, {  3, 87.09 }, {  4, 84.28 }, 
               {  5, 83.69 }, {  6, 89.18 }, {  7, 85.71 }, 
               {  8, 85.05 }, {  9, 85.58 }, { 10, 86.95 }, 
               { 11, 87.95 }, { 12, 89.44 }, { 13, 93.47 } } );
Vector x_data = points.GetColumnVector(0);
Vector y_data = points.GetColumnVector(1);

which produces the following coefficients (from lowest power to highest)

Coef=[85.892,-0.5542,0.074990]
     x        y
  0.00  85.8920
  1.00  85.4127
  2.00  85.0835
  3.00  84.9043
  4.00  84.8750
  5.00  84.9957
  6.00  85.2664
  7.00  85.6871
  8.00  86.2577
  9.00  86.9783
 10.00  87.8490
 11.00  88.8695
 12.00  90.0401
 13.00  91.3607
 14.00  92.8312

@ja72 code is pretty good. But I ported it on the present version of Math.NET (MathNet.Iridium is not supported for now as I understand) and optimized code size and performance (For instance, Math.Pow function is not used in my solution because of slow performance).

public class PolynomialRegression
{
    private int _order;
    private Vector<double> _coefs;

    public PolynomialRegression(DenseVector xData, DenseVector yData, int order)
    {
        _order = order;
        int n = xData.Count;

        var vandMatrix = new DenseMatrix(xData.Count, order + 1);
        for (int i = 0; i < n; i++)
            vandMatrix.SetRow(i, VandermondeRow(xData[i]));
        
        // var vandMatrixT = vandMatrix.Transpose();
        // 1 variant:
        //_coefs = (vandMatrixT * vandMatrix).Inverse() * vandMatrixT * yData;
        // 2 variant:
        //_coefs = (vandMatrixT * vandMatrix).LU().Solve(vandMatrixT * yData);
        // 3 variant (most fast I think. Possible LU decomposion also can be replaced with one triangular matrix):
        _coefs = vandMatrix.TransposeThisAndMultiply(vandMatrix).LU().Solve(TransposeAndMult(vandMatrix, yData));
    }

    private Vector<double> VandermondeRow(double x)
    {
        double[] result = new double[_order + 1];
        double mult = 1;
        for (int i = 0; i <= _order; i++)
        {
            result[i] = mult;
            mult *= x;
        }
        return new DenseVector(result);
    }

    private static DenseVector TransposeAndMult(Matrix m, Vector v)
    {
        var result = new DenseVector(m.ColumnCount);
        for (int j = 0; j < m.RowCount; j++)
        {
            double v_j = v[j];
            for (int i = 0; i < m.ColumnCount; i++)
                result[i] += m[j, i] * v_j;
        }
        return result;
    }

    public double Calculate(double x)
    {
        return VandermondeRow(x) * _coefs;
    }
}

It's also available on github:gist.

I don't think you want non linear regression. Even if you are using a quadratic function, it is still called linear regression. What you want is called multivariate regression. If you want a quadratic you just add a x squared term to your dependent variables.