2. Another example

The example below is a small constrained quadratic minimization. It is actually a tiny SVM problem that uses only inequality constraints, dropping the equality constraint as described here. Note that the Vector type is loaded with operators, which allow natural mathematic syntax describing a problem.

var g = new Vector(new double[] { -1, -1, -1, -1, -1 });
 
var gramMatrix = 
  new double[,]
  {
    { 3, 1, 3, 0, 1 },
    { 1, 3, 3, 0, 1 },
    { 3, 3, 5, 1, 3 },
    { 0, 0, 1, 2, 3 },
    { 1, 1, 3, 3, 5 }
  };
 
var Q = Vector.GetTensor(gramMatrix);
 
ScalarFunction f = // Goal function.
  x => 0.5 * x * Q(x) + g * x;
 
VectorFunction df = // Gradient of the goal function.
  x => Q(x) + g;
 
int P = g.Length;
 
var range = Enumerable.Range(0, P);
 
VectorFunction d2fd = // Diagonal of the Hessian of the goal.
  x =>
    range.Select(i => gramMatrix[i, i]);
 
int constraintsCount = 2 * P;
 
// Constraint functions: 0 <= fc(i)(x) <= 10
Func<int, ScalarFunction> fc =
  i =>
    x =>
      i < P ? 
      -x[i] : 
      x[i - P] - 10.0;
 
Func<int, VectorFunction> dfc = // Gradients of constraints.
  i =>
    x =>
      i < P ?
      range.Select(j => (j == i) ? -1.0 : 0.0) :
      range.Select(j => (j == i - P) ? 1.0 : 0.0);
 
var zeroVector = new Vector(P);
 
Func<int, VectorFunction> d2fcd = // Diagonals of the constraints' Hessians.
  i =>
    x => zeroVector;
 
// Unility function that returns Jacobi 
// preconditioner Tensor function.
var M = ConjugateGradient.ConstrainedMinimizeOptions.GetJacobiPreconditioner(
  d2fd,
  P,
  fc,
  dfc,
  d2fcd);
 
var certificate = ConjugateGradient.LineSearchConstrainedMinimize(
  f,
  df,
  new Vector(new double[] { 0.5, 0.5, 0.5, 0.5, 0.5 }), // Initial pt.
  constraintsCount,
  fc,
  dfc,
  new ConjugateGradient.LineSearchConstrainedMinimizeOptions()
  {
    DualityGap = 1e-15
  },
  M);
 
// What we expect.
var xe = new Vector(new double[] { 0.25, 0.25, 0.0, 0.5, 0.0 });
 
Assert.IsTrue(
  (certificate.Optimum - xe).Norm2 < ε, 
  String.Format(
    "Should have found {0}, found {1} instead", 
    xe, 
    certificate.Optimum));