Introduction to Numerical Solutions

///////////////////////////////////////////////////////////////////////////////
//
//  Program.cs
//
//  By Philip R. Braica (HoshiKata@aol.com, VeryMadSci@gmail.com)
//
//  Distributed under the The Code Project Open License (CPOL)
//  http://www.codeproject.com/info/cpol10.aspx
///////////////////////////////////////////////////////////////////////////////

// Using.
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

// Namespace.
namespace SolverDemo
{
    /// <summary>
    /// Linear least squares for 2 variables, 1st and second order,
    /// beyond that requires a decent matrix math package with good pseudo-inverse.
    /// Note that non-linear is more accurate when there is noise but overkill for
    /// this application, this is just a tiny quick fitter for use by the Solver class.
    /// 
    /// See the following for more information:
    /// http://en.wikipedia.org/wiki/Simple_linear_regression
    /// http://en.wikipedia.org/wiki/Invertible_matrix
    /// </summary>
    public class LinearLeastSquares
    {
        /// <summary>
        /// The solution, either:
        ///   y = Solution[0]*x + Solution[1];
        /// or if Solution.Length == 3: 
        ///   y = Solution[0] + Solution[1]*x + Solution[2]*x*x;
        /// </summary>
        public double[] Solution { get; protected set; }

        /// <summary>
        /// Solve for y = mx + b -> 
        /// y = Solution[0]*x + Solution[1];
        /// </summary>
        /// <param name="x"></param>
        /// <param name="y"></param>
        public bool Solve_FirstOrder(double[] x, double[] y)
        {
            // This is the classical software approach to doing
            // this fast, form the sums, do as little math as needed.
            double sx = 0;
            double sxx = 0;
            double sy = 0;
            double sxy = 0;
            double syy = 0;

            // compute the length.
            int n = x.Length < y.Length ? x.Length : y.Length;
            
            // Sum.
            for (int i = 0; i < n; i++)
            {
                double xi = x[i];
                double yi = y[i];
                sx += xi;
                sxx += xi * xi;
                sy += yi;
                sxy += xi * yi;
                syy += yi * yi;
            }

            // y = bx + a
            double b = ((n * sxy) - (sx * sy)) / ((n * sxx) - (sx * sx));
            double a = (sy - (b*sx))/n;

            // Store.
            Solution = new double[2];
            Solution[0] = a;
            Solution[1] = b;

            // Done.
            return true;
        }

        /// <summary>
        /// Solve for y = a0 + a1*x + a2*x^2 -> 
        /// y = Solution[0] + Solution[1]*x + Solution[2]*x*x;
        /// </summary>
        /// <param name="x"></param>
        /// <param name="y"></param>
        public bool Solve_SecondOrder(double[] x, double[] y)
        {
            // This is not done as a matrix operation
            // because it is much faster to calculate this way in part
            // because we know of the inherent matrix symetry and 3x3
            // matrixes aren't soo bad to invert in code this way.
            // 
            // Each is a sum: 
            //  sx = sum of the x values.
            //  sx2 = sum of each x value after squaring that value.
            double sx = 0;
            double sx2 = 0;
            double sx3 = 0;
            double sx4 = 0;
            double sy = 0;
            double sxy = 0;
            double sx2y = 0;
            double sy2 = 0;
            
            // n is the data set length.
            int n = x.Length < y.Length ? x.Length : y.Length;
            
            // Sum everything.
            for (int i = 0; i < n; i++)
            {
                // Avoid repeated de-indexing, and repeated multiplies.
                double xi = x[i];
                double yi = y[i];
                double xi2 = xi * xi;

                // Sum.
                sx += xi;
                sx2 += xi2;
                sy += yi;
                sxy += xi * yi;
                sy2 += yi * yi;
                sx3 += xi2 * xi;
                sx4 += xi2 * xi2;
                sx2y += xi2 * yi;
            }

            // In matrix form, solving for [a0, a1, a2].
            //      n,  sx,  sx2   a0    sy     
            //     sx,  sx2, sx3 * a1  = sxy
            //     sx2, sx3, sx4   a2    sx2y
            //
            // Has a solution if the determinant isn't zero of:
            //     a, b, c   sy             a0
            //     b, d, e * sxy  * 1/det = a1
            //     c, e, f   sx2y           a2
            // det = n * a + sx * b + sx2 * c
            double a = (sx2 * sx4) - (sx3 * sx3);
            double b = (sx3 * sx2) - (sx * sx4);
            double c = (sx * sx3) - (sx2 * sx2);
            double d = (n * sx4) - (sx2 * sx2);
            double e = (sx * sx2) - (n * sx3);
            double f = (n * sx2) - (sx * sx);
            double det = (n * a) + (sx * b) + (sx2 * c);

            // Don't bother trying
            if (det == 0) return false;
            det = 1 / det;
            double a0 = det * ((sy * a) + (sxy * b) + (sx2y * c));
            double a1 = det * ((sy * b) + (sxy * d) + (sx2y * e));
            double a2 = det * ((sy * c) + (sxy * e) + (sx2y * f));
            
            // Store
            Solution = new double[3];
            Solution[0] = a0;
            Solution[1] = a1;
            Solution[2] = a2;

            // Done.
            return true;
        }

        /// <summary>
        /// Compute the estimated Y values for this x.
        /// </summary>
        /// <param name="x"></param>
        /// <returns></returns>
        public double[] ComputeYValues(double[] x)
        {
            double[] y = new double[x.Length];
            if (Solution.Length == 2)
            {
                double b = Solution[0];
                double a = Solution[1];
                for (int i = 0; i < x.Length; i++)
                {
                    y[i] = (a * x[i]) + b;
                }
                return y;
            }
            if (Solution.Length == 3)
            {
                double a0 = Solution[0];
                double a1 = Solution[1];
                double a2 = Solution[2];
                for (int i = 0; i < x.Length; i++)
                {
                    double xi = x[i];
                    y[i] = a0 + (a1 * xi) + (a2 * xi * xi);
                }
                return y;
            }
            return null;
        }

        /// <summary>
        /// Compute the error of this fit.
        /// </summary>
        /// <param name="x"></param>
        /// <returns></returns>
        public double ComputeError(double[] x, double[] y)
        {
            double[] ye = ComputeYValues(x);
            double err = 0;
            for (int i = 0; i < x.Length; i++)
            {
                double de = y[i] - ye[i];
                err += (de * de);
            }
            err = System.Math.Sqrt(err) / x.Length;
            return err;
        }
    }
}