C# Iteratively Reweighted Least Sq Example

using System;

using CenterSpace.NMath.Core;
using CenterSpace.NMath.Matrix;

namespace CenterSpace.NMath.Matrix.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.
      RandGenUniform rng = new RandGenUniform(-2, 2, 0x124);
      int rows = 10;
      int cols = 2;
      DoubleMatrix A = new DoubleMatrix(rows, cols, rng);
      DoubleVector 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.
      DoubleIterativelyReweightedLeastSq 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("Solution with bisquare weighting");
      Console.WriteLine(solution);
      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.
      DoubleIterativelyReweightedLeastSq.ToleranceMetFunction 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" + solution);
      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;
    }
  }
}