C# Iteratively Reweighted Least Sq Example

using System;

using CenterSpace.NMath.Core;


namespace CenterSpace.NMath.Examples.CSharp
{
  /// <summary>
  /// A .NET example in C# demonstrating the features of the classes for solving iteratively reweighted
  /// least squares problems.
  /// </summary>
  class IterativelyReweightedLeastSqExample
  {
    static void Main( string[] args )
    {
      // Set up a least squares problem, Ax = b, with random data.
      var rng = new RandGenUniform( -2, 2, 0x124 );
      int rows = 10;
      int cols = 2;
      var A = new DoubleMatrix( rows, cols, rng );
      var x = new DoubleVector( cols, rng );
      // Fix up the right hand side b so that x
      // is the exact solution, then throw in some outliers.
      DoubleVector b = NMathFunctions.Product( A, x );
      // Throw in a few outliers...
      b[1] = 23;
      b[4] = -10;

      // Create an iteratively reweighted least squares instance
      // and use it to solve the problem using the default settings.
      // The default weighting is DoubleBisquareWeightingFunction which
      // uses the bisquare weighting algorithm.
      var irls = new DoubleIterativelyReweightedLeastSq();

      // Solve. The third parameter below specifies prepending a column of ones to the
      // data in A representing a constant term in the model (which should come out
      // to be zero in the solution from the way we cooked the data). Note that our
      // input matrix A will not actually be changed.
      DoubleVector solution = irls.Solve( A, b, true );

      Console.WriteLine();

      Console.WriteLine( "Solution with bisquare weighting" );
      Console.WriteLine( solution.ToString( "G5" ) );
      Console.WriteLine();
      Console.WriteLine( "||residuals|| = " + irls.Residuals.TwoNorm() );
      Console.WriteLine();
      if ( irls.Iterations >= irls.MaxIterations )
      {
        Console.WriteLine( "The algorithm did not converge in {0} iterations.", irls.MaxIterations );
      }
      else
      {
        Console.WriteLine( "Algorithm converged in {0} iterations.", irls.Iterations );
      }

      // Change some of the settings that control the iteration.
      irls.MaxIterations = 300;
      irls.Tolerance = 1e-7;
      // The convergence function is a delegate that may specified by the user for
      // determining if the algorithm has converged and iteration terminated. The
      // delegate takes as arguments the previous and current solutions and residuals
      // and the tolerance and returns a bool. See the ResidualsChanged function 
      // below.
      var residualsUnchanged =
        new DoubleIterativelyReweightedLeastSq.ToleranceMetFunction( ResidualsUnchanged );
      irls.ConvergenceFunction = residualsUnchanged;

      // Change the weighting function used from the default bisquare weighting to the
      // fair weighting function. See the class DoubleFairWeightingFunction for 
      // particulars.
      irls.WeightsFunction = new DoubleFairWeightingFunction();

      // Solve the problem with the new settings.
      solution = irls.Solve( A, b, true );
      Console.WriteLine();
      Console.WriteLine( "Solution with fair weighting" );
      Console.WriteLine( solution.ToString( "G5" ) );
      Console.WriteLine();
      Console.WriteLine( "||residuals|| = " + irls.Residuals.TwoNorm() );
      irls.Residuals.TwoNorm();
      Console.WriteLine();
      if ( irls.Iterations >= irls.MaxIterations )
      {
        Console.WriteLine( "The algorithm did not converge in {0} iterations.", irls.MaxIterations );
      }
      else
      {
        Console.WriteLine( "Algorithm converged in {0} iterations.", irls.Iterations );
      }

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

    // Convergence function for use in the iteratively reweighted least squares
    // algorithm. This particular convergence function returns true when the residuals
    // from the current iterations are relatively the same as the residuals in the
    // previous iteration.
    static bool ResidualsUnchanged( double tolerance, DoubleVector lastSolution,
      DoubleVector currentSolution, DoubleVector lastResiduals, DoubleVector currentResiduals )
    {
      double maxAbsDiff = NMathFunctions.MaxAbsValue( currentResiduals - lastResiduals );
      return ( maxAbsDiff / NMathFunctions.MaxAbsValue( currentResiduals ) ) < tolerance;
    }
  }
}