C# Constrained Least Squares Example

using System;

using CenterSpace.NMath.Core;


namespace CenterSpace.NMath.Examples.CSharp
{
  /// <summary>
  /// A .NET example in C# showing how to use ConstrainedLeastSquares class to
  /// solve the constrained least squares problem
  /// Cx = d, subject to the constraints
  /// Ax &lt; b
  /// </summary>
  public class ConstrainedLeastSquaresExample
  {
    static void Main( string[] args )
    {
      // Solve, in the least squares sense, Cx = d, subect to Ax <= b 
      // and -0.1 <= x[i] <= 2.0
      var C = new DoubleMatrix( "5x4 [0.9501    0.7620    0.6153    0.4057 " +
                                              "0.2311    0.4564    0.7919    0.9354 " +
                                              "0.6068    0.0185    0.9218    0.9169 " +
                                              "0.4859    0.8214    0.7382    0.4102 " +
                                              "0.8912    0.4447    0.1762    0.8936]" );

      var d = new DoubleVector( 0.0578, 0.3528, 0.8131, 0.0098, 0.1388 );

      // Constraint coefficient matrix
      var A = new DoubleMatrix( "3x4[0.2027    0.2721    0.7467    0.4659 " +
                                             "0.1987    0.1988    0.4450    0.4186 " +
                                             "0.6037    0.0152    0.9318    0.8462]" );

      // Constraints right hand sides.
      var b = new DoubleVector( 0.5251, 0.2026, 0.6721 );

      // Create the constrained least squares problem for minimizing
      // || Cx - d||^2 subject to Ax <= b and -0.1 <= x[i] <= 2.0
      // We first construct the problem object from the matrix C and the
      // vector d, we then add the constraints.
      var problem = new ConstrainedLeastSquaresProblem( C, d );
      // Add the inequality constraints Ax <= b using a constraint tolerance
      // of 0.00001. This allows for small violations of the constraints. 
      // Specifically the constraints will be considered satisfied for a 
      // vector x if
      // Ax <= b + 0.00001
      double constraintTolerance = 0.00001;
      for ( int i = 0; i < A.Rows; i++ )
      {
        problem.AddUpperBoundConstraint( A.Row( i ), b[i], constraintTolerance );
      }
      // All variable values for the solution must satisfy the bounds
      // -0.1 <= x[i] <= 2.0
      var lb = new DoubleVector( problem.NumVariables, -0.10 );
      var ub = new DoubleVector( problem.NumVariables, 2.0 );
      for ( int i = 0; i < problem.NumVariables; i++ )
      {
        problem.AddBounds( i, lb[i], ub[i], .00001 );
      }

      // Create the solver instance.
      var solver = new ConstrainedLeastSquares();
      // The ConstrainedLeastSquares solver uses a QP (Quadratic Programming) solver
      // to solve the constrained least squares problem.
      // The current default QP solver is the NMath active set quadratic programming
      // solver with default options.
      bool success = solver.Solve( problem );
      Console.WriteLine( "Default solver success = " + success );
      Console.WriteLine( "Default solver solution x = " + solver.X.ToString( "0.0000" ) );
      Console.WriteLine( "Default solver residual norm = " + solver.ResidualNorm.ToString( "0.0000" ) );
      Console.WriteLine( "Default solver performed {0} iterations", solver.Iterations );

      // You can pass in an instance of a quadratic programming solver for the 
      // constrained least squares class to use. This allows you to set 
      // option on the QP solver and inspect results of the QP
      // solver.
      var interiorPointQp = new InteriorPointQPSolver();
      var solverParams = new InteriorPointQPSolverParams
      {
        MaxIterations = 10000,
        PresolveLevel = InteriorPointQPSolverParams.PresolveLevelOption.None
      };
      solver.Solve( problem, interiorPointQp, solverParams );
      Console.WriteLine( "\nInterior point solver success = " + success );
      Console.WriteLine( "Interior point QP result = " + interiorPointQp.Result );
      Console.WriteLine( "Interior point QP solver iteration count = " + solver.Iterations );
      Console.WriteLine( "Interior point solver solution x = " + solver.X.ToString( "0.0000" ) );
      Console.WriteLine( "Interior point solver residual norm = " + solver.ResidualNorm.ToString( "0.0000" ) );

      // If you use the active set QP solver you can determine which constraints
      // are active in the solution by accessing the active set QP solvers 
      // Lagrange multiplier property. A constraint is active if its corresponding
      // Lagrange multiplier is non-zero.
      var activeSetQP = new ActiveSetQPSolver();
      success = solver.Solve( problem, activeSetQP );
      Console.WriteLine( "\nActive set solver success = " + success );
      Console.WriteLine( "Active set solver solution x = " + solver.X.ToString( "0.0000" ) );
      Console.WriteLine( "Active set solver residual norm = " + solver.ResidualNorm.ToString( "0.0000" ) );
      // Print out the active constraints.
      for ( int i = 0; i < activeSetQP.LagrangeMultiplier.Length; i++ )
      {
        if ( activeSetQP.LagrangeMultiplier[i] != 0.0 )
        {
          Console.WriteLine( "Constraint {0} = {1} is active", i, problem.Constraints[i].ToString() );
        }
      }

      Console.WriteLine();
      Console.WriteLine( "Press Enter Key" );
      Console.Read();
    }
  }
}