Least Squares in C# QuickStart Sample

Illustrates how to solve least squares problems using classes in the Extreme.Mathematics.LinearAlgebra namespace in C#.

View this sample in: Visual Basic F# IronPython

using System;
// The DenseMatrix and LeastSquaresSolver classes reside in the 
// Extreme.Mathematics.LinearAlgebra namespace.
using Extreme.Mathematics;
using Extreme.Mathematics.LinearAlgebra;

namespace Extreme.Numerics.QuickStart.CSharp
{
    /// <summary>
    /// Illustrates solving least squares problems using the 
    /// LeastSquaresSolver class in the Extreme.Mathematics.LinearAlgebra 
    /// namespace of Extreme Numerics.NET.
    /// </summary>
    class LeastSquares
    {
        static void Main(string[] args)
        {
            // The license is verified at runtime. We're using
            // a demo license here. For more information, see
            // https://numerics.net/trial-key
            Extreme.License.Verify("Demo license");
            // A least squares problem consists in finding
            // the solution to an overdetermined system of
            // simultaneous linear equations so that the
            // sum of the squares of the error is minimal.
            //
            // A common application is fitting data to a
            // curve. See the CurveFitting sample application
            // for a complete example.

            // Let's start with a general matrix. This will be
            // the matrix a in the left hand side ax=b:
            var a = Matrix.Create(6, 4, new double[]
                {
                    1, 1, 1, 1, 1, 1,
                    1, 2, 3, 4, 5, 6,
                    1, 4, 9, 16, 25, 36,
                    1, 2, 1, 2, 1, 2
                }, MatrixElementOrder.ColumnMajor);
            // Here is the right hand side:
            var b = Vector.Create(new double[] { 1, 3, 6, 11, 15, 21 });
            var b2 = Matrix.Create(6, 2, new double[]
                {
                    1, 3, 6, 11, 15, 21,
                    1, 2, 3, 4, 5, 7
                }, MatrixElementOrder.ColumnMajor);
            Console.WriteLine("a = {0:F4}", a);
            Console.WriteLine("b = {0:F4}", b);

            //
            // The LeastSquaresSolver class
            //

            var x = a.LeastSquaresSolve(b);
            var qr = a.GetQRDecomposition();
            qr.LeastSquaresSolve(b);
            // The following creates an instance of the
            // LeastSquaresSolver class for our problem:
            var solver = new LeastSquaresSolver<double>(a, b);
            // We can specify the solution method: normal
            // equations or QR decomposition. In most cases,
            // a QR decomposition is the most desirable:
            solver.SolutionMethod = LeastSquaresSolutionMethod.QRDecomposition;
            // The Solve method calculates the solution:
            x = solver.Solve();
            Console.WriteLine("x = {0:F4}", x);
            // The Solution property also returns the solution:
            Console.WriteLine("x = {0:F4}", solver.Solution);
            // More detailed information is available from
            // additional methods.
            // The values of the right hand side predicted 
            // by the solution:
            Console.WriteLine("Predictions = {0:F4}", solver.GetPredictions());
            // The residuals (errors) of the solution:
            Console.WriteLine("Residuals = {0:F4}", solver.GetResiduals());
            // The total sum of squares of the residues:
            Console.WriteLine("Residual square error = {0}",
                solver.GetResidualSumOfSquares());

            //
            // Direct normal equations
            //

            // Alternatively, you can create a least squares
            // solution by providing the normal equations
            // directly. This may be useful when it is easy
            // to calculate the normal equations directly.
            // 
            // Here, we'll just calculate the normal equation:
            var aTa = SymmetricMatrix<double>.FromOuterProduct(a);
            var aTb = b * a; // a.Transpose() * b;
            // We find the solution by solving the normal equations
            // directly:
            x = aTa.Solve(aTb);
            Console.WriteLine("x = {0:F4}", x);
            // However, properties of the least squares solution, such as
            // error estimates and residuals are not available.

            Console.Write("Press Enter key to exit...");
            Console.ReadLine();
        }
    }
}