<|start_filename|>src/test/java/com/manyangled/gibbous/FeasiblePointTest.java<|end_filename|>
/*
Copyright 2018 <NAME>
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/

package com.manyangled.gibbous;

import org.junit.Test;
import static org.junit.Assert.*;

import org.apache.commons.math3.linear.RealMatrix;
import org.apache.commons.math3.linear.RealVector;
import org.apache.commons.math3.linear.Array2DRowRealMatrix;
import org.apache.commons.math3.linear.DiagonalMatrix;
import org.apache.commons.math3.linear.ArrayRealVector;

import org.apache.commons.math3.optim.PointValuePair;
import org.apache.commons.math3.optim.nonlinear.scalar.ObjectiveFunction;
import org.apache.commons.math3.optim.InitialGuess;

import com.manyangled.gibbous.optim.convex.TwiceDifferentiableFunction;
import com.manyangled.gibbous.optim.convex.QuadraticFunction;
import com.manyangled.gibbous.optim.convex.LinearFunction;
import com.manyangled.gibbous.optim.convex.LinearEqualityConstraint;
import com.manyangled.gibbous.optim.convex.InequalityConstraintSet;

import com.manyangled.gibbous.optim.convex.ConvexOptimizer;

import static com.manyangled.gibbous.COTestingUtils.translatedQF;
import static com.manyangled.gibbous.COTestingUtils.eps;

public class FeasiblePointTest {
    private void testFeasibleConstraints(InequalityConstraintSet constraints) {
        PointValuePair pvp = ConvexOptimizer.feasiblePoint(constraints);
        //System.out.format("fkmax= %s\n", pvp.getSecond());
        assertTrue(pvp.getSecond() < 0.0);
        for (TwiceDifferentiableFunction f: constraints.constraints) {
            double y = f.value(pvp.getFirst());
            //System.out.format("y= %s\n", y);
            assertTrue(f.value(pvp.getFirst()) < 0.0);
        }
    }

    @Test
    public void test2DHalfPlane() {
        InequalityConstraintSet hp = new InequalityConstraintSet(
            new LinearFunction(new double[] { -1.0, 0.0 }, 100.0)
        );
        testFeasibleConstraints(hp);
    }


    @Test
    public void test2DSquareRegion() {
        double offset = 100.0;
        InequalityConstraintSet hp = new InequalityConstraintSet(
            new LinearFunction(new double[] { 1.0, 0.0 }, -(5.0 + offset)),
            new LinearFunction(new double[] { -1.0, 0.0 }, 1.0 + offset),
            new LinearFunction(new double[] { 0.0, 1.0 }, -(5.0 + offset)),
            new LinearFunction(new double[] { 0.0, -1.0 }, 1.0 + offset)
        );
        testFeasibleConstraints(hp);
    }

}


<|start_filename|>src/main/java/com/manyangled/gibbous/optim/convex/LinearFunction.java<|end_filename|>
/*
Copyright 2018 <NAME>
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/

package com.manyangled.gibbous.optim.convex;

import org.apache.commons.math3.linear.RealMatrix;
import org.apache.commons.math3.linear.RealVector;
import org.apache.commons.math3.linear.OpenMapRealMatrix;
import org.apache.commons.math3.linear.ArrayRealVector;

/**
 * A linear function b.x + c; The gradient is b, and the Hessian is |0| (all zeros).
 */
public class LinearFunction extends TwiceDifferentiableFunction {
    private final RealVector b;
    private final double c;
    private final int n;

    /**
     * Construct a linear function b.x + c
     *
     * @param b a weight vector
     * @param c a constant
     */
    public LinearFunction(RealVector b, double c) {
        int d = b.getDimension();
        if (d < 1) throw new IllegalArgumentException("Dimension must be nonzero");
        this.b = b;
        this.c = c;
        this.n = d;
    }

    /**
     * Construct a linear function b.x + c
     *
     * @param b a weight vector
     * @param c a constant
     */
    public LinearFunction(double[] b, double c) {
        this(new ArrayRealVector(b), c);
    }

    @Override
    public int dim() { return n; }

    @Override
    public double value(final RealVector x) {
        return c + b.dotProduct(x);
    }

    @Override
    public RealVector gradient(final RealVector x) {
        return b.copy();
    }

    @Override
    public RealMatrix hessian(final RealVector x) {
        // the Hessian is just zero for a linear function
        return new OpenMapRealMatrix(n, n);
    }

    @Override
    public String toString() {
        return String.format("LinearFunction(%g, %s)", c, b.toString());
    }
}


<|start_filename|>src/main/java/com/manyangled/gibbous/optim/convex/SmoothMaxFunction.java<|end_filename|>
/*
Copyright 2018 <NAME>
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/

package com.manyangled.gibbous.optim.convex;

import org.apache.commons.math3.util.Pair;

import org.apache.commons.math3.linear.RealMatrix;
import org.apache.commons.math3.linear.RealVector;
import org.apache.commons.math3.linear.Array2DRowRealMatrix;
import org.apache.commons.math3.linear.ArrayRealVector;

/**
 * Given a set of convex functions f[k] for k = 1,2,..., implements the smooth-max
 * over f[k]: (1/a)log(sum(e^(a f[k]))).
 * See also: http://erikerlandson.github.io/blog/2018/05/27/the-gradient-and-hessian-of-the-smooth-max-over-functions/
 */
public class SmoothMaxFunction extends TwiceDifferentiableFunction {
    private final double alpha;
    private final TwiceDifferentiableFunction[] f;

    /**
     * Construct a smooth-max function over a collection of convex functions
     *
     * @param alpha A weighting parameter
     * @param f A collection of convex functions
     */
    public SmoothMaxFunction(double alpha, TwiceDifferentiableFunction... f) {
        if (f.length < 1) throw new IllegalArgumentException("list of functions must be nonempty");
        this.alpha = alpha;
        this.f = f;
    }

    @Override
    public int dim() { return f[0].dim(); }

    @Override
    public double value(final RealVector x) {
        Pair<Double, double[]> pre = precompute(x);
        double z = pre.getFirst();
        double[] exp = pre.getSecond();
        double s = 0.0;
        for (double e: exp) s += e;
        return z + (Math.log(s) / alpha);
    }

    @Override
    public RealVector gradient(final RealVector x) {
        Pair<Double, double[]> pre = precompute(x);
        double[] exp = pre.getSecond();
        RealVector g = new ArrayRealVector(dim(), 0.0);
        double d = 0.0;
        for (int k = 0; k < f.length; ++k) {
            RealVector gk = f[k].gradient(x);
            d += exp[k];
            g.combineToSelf(1.0, exp[k], gk);
        }
        g.mapDivideToSelf(d);
        return g;
    }

    @Override
    public RealMatrix hessian(final RealVector x) {
        Pair<Double, double[]> pre = precompute(x);
        double[] exp = pre.getSecond();
        int n = dim();
        RealMatrix h = new Array2DRowRealMatrix(n, n);
        double d = 0.0;
        for (int k = 0; k < f.length; ++k) {
            d += exp[k];
            RealVector gfk = f[k].gradient(x);
            RealMatrix hfk = f[k].hessian(x);
            for (int i = 0; i < n; ++i)
                for (int j = 0; j < n; ++j) {
                    h.addToEntry(i, j, exp[k] * hfk.getEntry(i, j));
                    h.addToEntry(i, j, alpha * exp[k] * gfk.getEntry(i) * gfk.getEntry(j));
                }
        }
        d = 1.0 / d;
        RealVector g = gradient(x);
        for (int i = 0; i < n; ++i)
            for (int j = 0; j < n; ++j) {
                h.multiplyEntry(i, j, d);
                h.addToEntry(i, j, (-alpha) * g.getEntry(i) * g.getEntry(j));
            }
        return h;
    }

    // Pre-computes z, the maximum of f[k](x), and also e^(α(f[k]-z)) for each f[k].
    // See: http://erikerlandson.github.io/blog/2018/05/28/computing-smooth-max-and-its-gradients-without-over-and-underflow/
    private Pair<Double, double[]> precompute(final RealVector x) {
        double[] exp = new double[f.length];
        double z = Double.NEGATIVE_INFINITY;
        for (int k = 0; k < f.length; ++k) {
            double fk = f[k].value(x);
            exp[k] = fk;
            z = Math.max(z, fk);
        }
        for (int k = 0; k < f.length; ++k) {
            double fk = exp[k];
            exp[k] = Math.exp(alpha * (fk - z));
        }
        return new Pair<Double, double[]>(z, exp);
    }
}


<|start_filename|>src/main/java/com/manyangled/gibbous/optim/convex/TwiceDifferentiableFunction.java<|end_filename|>
/*
Copyright 2018 <NAME>
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/

package com.manyangled.gibbous.optim.convex;

import org.apache.commons.math3.analysis.MultivariateFunction;
import org.apache.commons.math3.exception.DimensionMismatchException;
import org.apache.commons.math3.linear.RealMatrix;
import org.apache.commons.math3.linear.RealVector;
import org.apache.commons.math3.linear.Array2DRowRealMatrix;
import org.apache.commons.math3.linear.ArrayRealVector;

/**
 * A MultivariateFunction that also has a defined gradient and Hessian
 */
public abstract class TwiceDifferentiableFunction implements MultivariateFunction {
    /**
     * Returns the dimensionality of the function domain.
     * If dim() returns (n) then this function expects an n-vector as its input.
     * @return the expected dimension of the function's domain
     */
    public abstract int dim();
    
    /**
     * Returns the value of this function at (x)
     *
     * @param x a point to evaluate this function at.
     * @return the value of this function at (x)
     */
    public abstract double value(final RealVector x);

    /**
     * Returns the gradient of this function at (x)
     *
     * @param x a point to evaluate this gradient at
     * @return the gradient of this function at (x)
     */
    public abstract RealVector gradient(final RealVector x);

    /**
     * The Hessian of this function at (x)
     *
     * @param x a point to evaluate this Hessian at
     * @return the Hessian of this function at (x)
     */
    public abstract RealMatrix hessian(final RealVector x);
    
    /**
     * Returns the value of this function at (x)
     *
     * @param x a point to evaluate this function at.
     * @return the value of this function at (x)
     */
    @Override
    public double value(final double[] x) {
        return value(new ArrayRealVector(x, false));
    }

    /**
     * Returns the gradient of this function at (x)
     *
     * @param x a point to evaluate this gradient at
     * @return the gradient of this function at (x)
     */
    public RealVector gradient(final double[] x) {
        return gradient(new ArrayRealVector(x, false));
    }

    /**
     * The Hessian of this function at (x)
     *
     * @param x a point to evaluate this Hessian at
     * @return the Hessian of this function at (x)
     */
    public RealMatrix hessian(final double[] x) {
        return hessian(new ArrayRealVector(x, false));
    }
}


<|start_filename|>src/main/java/com/manyangled/gibbous/optim/convex/LinearTransformFunction.java<|end_filename|>
/*
Copyright 2018 <NAME>
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/

package com.manyangled.gibbous.optim.convex;

import org.apache.commons.math3.linear.RealMatrix;
import org.apache.commons.math3.linear.RealVector;

import static com.manyangled.gibbous.optim.convex.ConvexOptimizer.isDense;

/**
 * Applies a linear transform to a function, from f(x) to t(x) = af(x) + b
 */
public class LinearTransformFunction extends TwiceDifferentiableFunction {
    private final double a;
    private final double b;
    private final TwiceDifferentiableFunction f;

    /**
     * construct a linear transform, af(x) + b, of a function f(x)
     * @param a the linear coefficient
     * @param b a constant
     * @param f the function
     */
    public LinearTransformFunction(
        final double a, final double b,
        final TwiceDifferentiableFunction f) {
        this.a = a;
        this.b = b;
        this.f = f;
    }

    @Override
    public int dim() { return f.dim(); }

    @Override
    public double value(final RealVector x) {
        return b + (a * f.value(x));
    }

    @Override
    public RealVector gradient(final RealVector x) {
        RealVector g = f.gradient(x);
        if (a != 1.0) g.mapMultiplyToSelf(a);
        return g;
    }

    @Override
    public RealMatrix hessian(final RealVector x) {
        RealMatrix h = f.hessian(x);
        if (a == 1.0) return h;
        if (isDense(h)) {
            int n = h.getRowDimension();
            for (int j = 0; j < n; ++j)
                for (int k = 0; k < n; ++k)
                    h.multiplyEntry(j, k, a);
            return h;
        } else {
            return h.scalarMultiply(a);
        }
    }
}


<|start_filename|>src/main/java/com/manyangled/gibbous/optim/convex/QuadraticFunction.java<|end_filename|>
/*
Copyright 2018 <NAME>
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/

package com.manyangled.gibbous.optim.convex;

import java.lang.Math;
import java.util.Collection;
import java.util.Iterator;
import java.util.Arrays;
import org.apache.commons.math3.linear.MatrixUtils;
import org.apache.commons.math3.linear.RealMatrix;
import org.apache.commons.math3.linear.DiagonalMatrix;
import org.apache.commons.math3.linear.RealVector;
import org.apache.commons.math3.linear.Array2DRowRealMatrix;
import org.apache.commons.math3.linear.ArrayRealVector;
import org.apache.commons.math3.exception.DimensionMismatchException;

/**
 * Given A, b, c, implements 0.5*(x^T)A(x) + b.x + c;
 * The gradient is A(x) + b, and the Hessian is A
 */
public class QuadraticFunction extends TwiceDifferentiableFunction {
    private final RealMatrix A;
    private final RealVector b;
    private final double c;
    private final int n;

    /**
     * Construct quadratic function 0.5*(x^T)A(x) + b.x + c
     *
     * @param A square matrix of weights for quadratic terms.
     * Typically expected to be positive definite or positive semi-definite.
     * @param b vector of weights for linear terms.
     * @param c a constant
     */
    public QuadraticFunction(RealMatrix A, RealVector b, double c) {
        int d = b.getDimension();
        if (d < 1) throw new IllegalArgumentException("Dimension must be nonzero");
        if (A.getRowDimension() != d)
            throw new DimensionMismatchException(A.getRowDimension(), d);
        MatrixUtils.checkSymmetric(A, 1e-6);
        this.A = A.copy();
        this.b = b.copy();
        this.c = c;
        this.n = d;
    }

    /**
     * Construct quadratic function 0.5*(x^T)A(x) + b.x + c
     *
     * @param A square matrix of weights for quadratic terms.
     * Typically expected to be positive definite or positive semi-definite.
     * @param b vector of weights for linear terms.
     * @param c a constant
     */
    public QuadraticFunction(double[][] A, double[] b, double c) {
        this(new Array2DRowRealMatrix(A), new ArrayRealVector(b), c);
    }

    @Override
    public int dim() { return n; }

    @Override
    public double value(final RealVector x) {
        double v = 0.5 * (A.operate(x)).dotProduct(x);
        v += b.dotProduct(x);
        v += c;
        return v;
    }

    @Override
    public RealVector gradient(final RealVector x) {
        return (A.operate(x)).add(b);
    }

    @Override
    public RealMatrix hessian(final RealVector x) {
        return A.copy();
    }

    /**
     * Create a quadratic function that corresponds to s((x-c).(x-c) &lt; r^2).
     * That is, constrained to an n-dimensional ball of radius r, with scaling factor s.
     *
     * @param center the center of the n-ball
     * @param r a radius, &gt; 0
     * @param s a scaling constant, &gt; 0
     * @return quadratic constraint function for the n-ball constraint
     */
    public static QuadraticFunction nBallConstraintFunction(RealVector center, double r, double s) {
        int n = center.getDimension();
        if (n < 1) throw new IllegalArgumentException("center vector must have dimension > 0");
        if (s <= 0.0) throw new IllegalArgumentException("scale s must be > 0");
        if (r <= 0.0) throw new IllegalArgumentException("radius r must be > 0");
        double[] alls = new double[n];
        java.util.Arrays.fill(alls, s);
        RealMatrix A = new DiagonalMatrix(alls);
        RealVector b = center.mapMultiply(-s);
        double c = 0.5 * ((s * center.dotProduct(center)) - (r * r));
        return new QuadraticFunction(A, b, c);
    }

    /**
     * Create a quadratic function that corresponds to s((x-c).(x-c) &lt; r^2).
     * That is, constrained to an n-dimensional ball of radius r, with scaling factor s.
     *
     * @param center the center of the n-ball
     * @param r a radius, &gt; 0
     * @param s a scaling constant, &gt; 0
     * @return quadratic constraint function for the n-ball constraint
     */
    public static QuadraticFunction nBallConstraintFunction(double[] center, double r, double s) {
        return nBallConstraintFunction(new ArrayRealVector(center, false), r, s);
    }

    /**
     * Create a quadratic function that corresponds to (x-c).(x-c) &lt; r^2.
     * That is, constrained to an n-dimensional ball of radius r
     *
     * @param center the center of the n-ball
     * @param r a radius, &gt; 0
     * @return quadratic constraint function for the n-ball constraint
     */
    public static QuadraticFunction nBallConstraintFunction(RealVector center, double r) {
        return nBallConstraintFunction(center, r, 1.0);
    }

    /**
     * Create a quadratic function that corresponds to (x-c).(x-c) &lt; r^2.
     * That is, constrained to an n-dimensional ball of radius r
     *
     * @param center the center of the n-ball
     * @param r a radius, &gt; 0
     * @return quadratic constraint function for the n-ball constraint
     */
    public static QuadraticFunction nBallConstraintFunction(double[] center, double r) {
        return nBallConstraintFunction(new ArrayRealVector(center), r, 1.0);
    }    
}


<|start_filename|>src/main/java/com/manyangled/gibbous/optim/convex/NewtonOptimizer.java<|end_filename|>
/*
Copyright 2018 <NAME>
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/

package com.manyangled.gibbous.optim.convex;

import org.apache.commons.math3.util.Pair;
import org.apache.commons.math3.optim.OptimizationData;
import org.apache.commons.math3.optim.PointValuePair;
import org.apache.commons.math3.exception.DimensionMismatchException;
import org.apache.commons.math3.exception.ConvergenceException;
import org.apache.commons.math3.linear.RealMatrix;
import org.apache.commons.math3.linear.RealVector;
import org.apache.commons.math3.linear.Array2DRowRealMatrix;
import org.apache.commons.math3.linear.ArrayRealVector;

/**
 * Implements convex optimization, using Newton's method.
 * Supports linear equality constraints and infeasible starting point<p>
 * (Algorithm 10.2) from Convex Optimization, Boyd and Vandenberghe, Cambridge University Press, 2008.
 * <p>
 * {@link BarrierOptimizer} supports the following {@link OptimizationData} parameters as arguments
 * to {@link #optimize(OptimizationData...)}:
 * <ul>
 *   <li>convex objective function: ObjectiveFunction - mandatory: must contain a {@link TwiceDifferentiableFunction} </li>
 *   <li>initial guess: InitialGuess - mandatory: need not satisfy equality constraints </li>
 *   <li>linear equality constraints: {@link LinearEqualityConstraint} - optional </li>
 *   <li>convergence epsilon: {@link ConvergenceEpsilon} - optional </li>
 *   <li>backtracking alpha: {@link BacktrackAlpha} - optional </li>
 *   <li>backtracking beta: {@link BacktrackBeta} - optional </li>
 *   <li>KKT equations solver: {@link KKTSolver} - optional </li>
 * </ul>
 */
public class NewtonOptimizer extends ConvexOptimizer {
    private LinearEqualityConstraint eqConstraint;
    private KKTSolver kktSolver = new SVDSchurKKTSolver();
    private RealVector xStart;
    private double epsilon = ConvergenceEpsilon.CONVERGENCE_EPSILON_DEFAULT;
    private double alpha = BacktrackAlpha.BACKTRACK_ALPHA_DEFAULT;
    private double beta = BacktrackBeta.BACKTRACK_BETA_DEFAULT;
    private HaltingCondition halting;

    public NewtonOptimizer() {
        super();
    }

    @Override
    public PointValuePair optimize(OptimizationData... optData) {
        return super.optimize(optData);
    }

    @Override
    protected void parseOptimizationData(OptimizationData... optData) {
        super.parseOptimizationData(optData);
        for (OptimizationData data: optData) {
            if (data instanceof LinearEqualityConstraint) {
                eqConstraint = (LinearEqualityConstraint)data;
                continue;
            }
            if (data instanceof KKTSolver) {
                kktSolver = (KKTSolver)data;
                continue;
            }
            if (data instanceof ConvergenceEpsilon) {
                epsilon = ((ConvergenceEpsilon)data).epsilon;
                continue;
            }
            if (data instanceof BacktrackAlpha) {
                alpha = ((BacktrackAlpha)data).alpha;
                continue;
            }
            if (data instanceof BacktrackBeta) {
                beta = ((BacktrackBeta)data).beta;
                continue;
            }
            if (data instanceof HaltingCondition) {
                halting = (HaltingCondition)data;
                continue;
            }
        }
        // if we got here, convexObjective exists
        int n = convexObjective.dim();
        if (eqConstraint != null) {
            int nDual = eqConstraint.b.getDimension();
            if (nDual >= n)
                throw new IllegalArgumentException("Rank of constraints must be < domain dimension");
            int nTest = eqConstraint.A.getColumnDimension();
            if ((nDual > 0) && (nTest != n))
                throw new DimensionMismatchException(nTest, n);
        }
        if (this.getStartPoint() != null) {
            xStart = new ArrayRealVector(this.getStartPoint());
            if (xStart.getDimension() != n)
                throw new DimensionMismatchException(xStart.getDimension(), n);
        } else {
            xStart = new ArrayRealVector(n, 0.0);
        }
    }

    @Override
    public PointValuePair doOptimize() {
        final int n = convexObjective.dim();
        if ((eqConstraint == null) || (eqConstraint.b.getDimension() < 1)) {
            // constraints Ax = b are empty
            // Algorithm 9.5: Newton's method (unconstrained)
            RealVector x = xStart;
            double v = convexObjective.value(x);
            while (true) {
                incrementIterationCount();
                RealVector grad = convexObjective.gradient(x);
                RealMatrix hess = convexObjective.hessian(x);
                KKTSolution sol = kktSolver.solve(hess, grad);
                if (sol.lambdaSquared <= (2.0 * epsilon)) break;
                RealVector xDelta = sol.xDelta;
                double gdd = grad.dotProduct(xDelta);
                RealVector tx = null;
                double tv = 0.0;
                boolean foundStep = false;
                for (double t = 1.0; t >= BT_T_LB; t *= beta) {
                    tx = x.add(xDelta.mapMultiply(t));
                    tv = convexObjective.value(tx);
                    if (Double.isInfinite(tv)) {
                        // this is barrier convention for "outside the feasible domain",
                        // so try a smaller step
                        continue;
                    }
                    double vtt = v + (t * alpha * gdd);
                    if (vtt == v) {
                        // (t)(alpha)(gdd) is less than (v)(machine-epsilon)
                        // Further tests for improvement are going to fail
                        break;
                    }
                    if (tv <= vtt) {
                        // This step resulted in an improvement, so halt with success
                        foundStep = true;
                        break;
                    }
                }
                // If there was no forward step to make, that indicates minimum
                if (!foundStep) break;
                // Update x,v for next iteration
                RealVector xprv = x;
                double vprv = v;
                x = tx;
                v = tv;
                // if improvement becomes very small then we are converged
                if (Math.abs(1.0 - (v / vprv)) < epsilon) break;
                // Check halting condition if configured
                if ((halting != null) && halting.checker.converged(
                        getIterations(),
                        new Pair<RealVector, Double>(xprv, vprv),
                        new Pair<RealVector, Double>(x, v))) {
                    break;
                }
            }
            return new PointValuePair(x.toArray(), v);
        } else {
            // constraints Ax = b are non-empty
            // Algorithm 10.2: Newton's method with equality constraints
            final RealMatrix A = eqConstraint.A;
            final RealVector b = eqConstraint.b;
            final RealMatrix AT = A.transpose();
            final int nDual = b.getDimension();
            RealVector x = xStart;
            RealVector nu = new ArrayRealVector(nDual, 0.0);
            double v = convexObjective.value(x);
            while (true) {
                incrementIterationCount();
                RealVector grad = convexObjective.gradient(x);
                double rNorm = residualNorm(x, nu, grad, A, AT, b);
                if (rNorm <= epsilon) break;
                RealMatrix hess = convexObjective.hessian(x);
                KKTSolution sol = kktSolver.solve(hess, A, AT, grad, A.operate(x).subtract(b));
                RealVector xDelta = sol.xDelta;
                RealVector nuDelta = sol.nuPlus.subtract(nu);
                RealVector tx = null;
                RealVector tnu = null;
                double tv = 0.0;
                boolean foundStep = false;
                for (double t = 1.0; t >= BT_T_LB; t *= beta) {
                    tx = x.add(xDelta.mapMultiply(t));
                    tv = convexObjective.value(tx);
                    if (Double.isInfinite(tv)) {
                        // this is barrier convention for "outside the feasible domain",
                        // so try a smaller step
                        continue;
                    }
                    double ftt = 1.0 - (alpha * t);
                    if (ftt == 1.0) {
                        // (t)(alpha) is less than (machine-epsilon)
                        // Further tests for improvement are going to fail
                        break;
                    }
                    tnu = nu.add(nuDelta.mapMultiply(t));
                    RealVector tgrad = convexObjective.gradient(tx);
                    double tNorm = residualNorm(tx, tnu, tgrad, A, AT, b);
                    if (tNorm <= ftt * rNorm) {
                        // This step resulted in an improvement, so halt with success
                        foundStep = true;
                        break;
                    }
                }
                // If there was no forward step to make, that indicates minimum
                if (!foundStep) break;
                // update for next iteration
                RealVector xprv = x;
                double vprv = v;
                x = tx;
                nu = tnu;
                v = tv;
                // if improvement becomes very small then we are converged
                if (Math.abs(1.0 - (v / vprv)) < epsilon) break;
                // check halting condition, if it was configured
                if ((halting != null) && halting.checker.converged(
                        getIterations(),
                        new Pair<RealVector, Double>(xprv, vprv),
                        new Pair<RealVector, Double>(x, v))) {
                    break;
                }
            }
            return new PointValuePair(x.toArray(), v);
        }
    }

    private double residualNorm(
        RealVector x, RealVector nu, RealVector grad,
        RealMatrix A, RealMatrix AT, RealVector b) {
        RealVector r = A.operate(x).subtract(b);
        RealVector rDual = grad.add(AT.operate(nu));
        double rr = r.dotProduct(r) + rDual.dotProduct(rDual);
        return Math.sqrt(rr);
    }

    // I used to use a "> 0" lower bound for t in backtracking, but I noticed a
    // case where t never underflowed to zero, so I am going to set it a bit higher
    // than the minimum positive value
    private static final double BT_T_LB = 1e-300;
}


<|start_filename|>src/main/java/com/manyangled/gibbous/optim/convex/LogBarrierFunction.java<|end_filename|>
/*
Copyright 2018 <NAME>
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/

package com.manyangled.gibbous.optim.convex;

import java.lang.Math;
import java.util.Collection;
import java.util.Iterator;
import java.util.Arrays;
import org.apache.commons.math3.linear.RealMatrix;
import org.apache.commons.math3.linear.RealVector;
import org.apache.commons.math3.linear.Array2DRowRealMatrix;
import org.apache.commons.math3.linear.ArrayRealVector;
import org.apache.commons.math3.exception.DimensionMismatchException;

/**
 * Given a convex objective function f0 and convex constraints f[k],
 * computes the log barrier function: <p>
 * b(x) = (t)f0(x) - sum(log(f[k](x))) <p>
 * returns +inf if any f[k](x) &gt;= 0
 */
public class LogBarrierFunction extends TwiceDifferentiableFunction {
    private final double t;
    private final TwiceDifferentiableFunction f0;
    private final TwiceDifferentiableFunction[] f;
    private final int n;

    /**
     * construct a log-barrier function b(x) = (t)f0(x) - sum(log(f[k](x)))
     * @param t multiplier constant for f0, must be &gt; 0
     * @param f0 a convex objective function
     * @param f a list of convex constraint functions
     */
    public LogBarrierFunction(double t, TwiceDifferentiableFunction f0, TwiceDifferentiableFunction[] f) {
        if (t <= 0.0) throw new IllegalArgumentException("t must be > 0");
        this.t = t;
        this.f0 = f0;
        this.n = f0.dim();
        this.f = f;
        for (TwiceDifferentiableFunction fi: f) {
            if (fi.dim() != n) throw new DimensionMismatchException(fi.dim(), n);
        }
    }

    /**
     * construct a log-barrier function b(x) = (t)f0(x) - sum(log(f[k](x)))
     * @param t multiplier constant for f0, must be &gt; 0
     * @param f0 a convex objective function
     * @param f a list of convex constraint functions
     */
    public LogBarrierFunction(double t, TwiceDifferentiableFunction f0, Collection<TwiceDifferentiableFunction> f) {
        this(t, f0, f.toArray(new TwiceDifferentiableFunction[0]));
    }

    @Override
    public int dim() {
        return n;
    }

    @Override
    public double value(final RealVector x) {
        double v = t * f0.value(x);
        for (TwiceDifferentiableFunction fi: f) {
            double ti = fi.value(x);
            if (ti >= 0.0) {
                return Double.POSITIVE_INFINITY;
            }
            v -= Math.log(-ti);
        }
        return v;
    }

    @Override
    public RealVector gradient(final RealVector x) {
        // g should be dense, due to contributions of barrier functions
        RealVector g = new ArrayRealVector(f0.gradient(x).toArray(), false);
        g.mapMultiplyToSelf(t);
        for (TwiceDifferentiableFunction fi: f) {
            double zi = -1.0 / fi.value(x);
            g.combineToSelf(1.0, zi, fi.gradient(x));
        }
        return g;
    }

    @Override
    public RealMatrix hessian(final RealVector x) {
        // h should be dense, due to contributions of barrier functions
        RealMatrix h = new Array2DRowRealMatrix(f0.hessian(x).getData(), false);
        for (int j = 0; j < n; ++j)
            for (int k = 0; k < n; ++k)
                h.multiplyEntry(j, k, t);
        for (TwiceDifferentiableFunction fi: f) {
            double vi = fi.value(x);
            RealVector gi = fi.gradient(x);
            RealMatrix hi = fi.hessian(x);
            for (int j = 0; j < n; ++j) {
                for (int k = 0; k < n; ++k) {
                    h.addToEntry(j, k, (gi.getEntry(j)*gi.getEntry(k)/(vi*vi)) - hi.getEntry(j, k)/vi);
                }
            }
        }
        return h;
    }
}


<|start_filename|>src/test/java/com/manyangled/gibbous/COTestingUtils.java<|end_filename|>
/*
Copyright 2018 <NAME>
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/

package com.manyangled.gibbous;

import org.junit.Test;
import static org.junit.Assert.assertEquals;
import static org.junit.Assert.assertArrayEquals;

import org.apache.commons.math3.linear.RealMatrix;
import org.apache.commons.math3.linear.RealVector;
import org.apache.commons.math3.linear.Array2DRowRealMatrix;
import org.apache.commons.math3.linear.DiagonalMatrix;
import org.apache.commons.math3.linear.ArrayRealVector;

import com.manyangled.gibbous.optim.convex.QuadraticFunction;

public class COTestingUtils {
    public static final double eps = 1e-8;

    public static QuadraticFunction translatedQF(double h, double[] center) {
        double[] all1 = new double[center.length];
        java.util.Arrays.fill(all1, 1.0);
        RealVector cv = new ArrayRealVector(center);
        RealMatrix A = new DiagonalMatrix(all1);
        RealVector b = cv.mapMultiply(-1.0);
        double c = h + 0.5*cv.dotProduct(cv);
        return new QuadraticFunction(A, b, c);
    }
}


<|start_filename|>src/main/java/com/manyangled/gibbous/optim/convex/InequalityConstraintSet.java<|end_filename|>
/*
Copyright 2018 <NAME>
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/

package com.manyangled.gibbous.optim.convex;

import java.util.ArrayList;
import java.util.Collection;

import org.apache.commons.math3.optim.OptimizationData;

/**
 * A collection of convex constraint functions. A constraint f is satisfied by x if
 * and only if f(x) &lt; 0. Constraints contained in this object are appended to the
 * list of constraints applied by a convex optimizer. Currently used by
 * {@link BarrierOptimizer} and {@link ConvexOptimizer#feasiblePoint(OptimizationData... optData)}
 */
public class InequalityConstraintSet implements OptimizationData {
    public final ArrayList<TwiceDifferentiableFunction> constraints =
        new ArrayList<TwiceDifferentiableFunction>();

    /**
     * Construct an inequality constraint set from a Collection of convex constraint functions
     * @param constraints the Collection of constraint functions to apply
     */
    public InequalityConstraintSet(Collection<TwiceDifferentiableFunction> constraints) {
        this.constraints.addAll(constraints);
    }

    /**
     * Construct an inequality constraint set from an argument list (or array) of convex
     * constraint functions
     * @param constraints the list of constraint functions to apply
     */
    public InequalityConstraintSet(TwiceDifferentiableFunction... constraints) {
        for (TwiceDifferentiableFunction f: constraints)
            this.constraints.add(f);
    }
}


<|start_filename|>src/main/java/com/manyangled/gibbous/optim/convex/ConvergenceEpsilon.java<|end_filename|>
/*
Copyright 2018 <NAME>
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/

package com.manyangled.gibbous.optim.convex;

import org.apache.commons.math3.optim.OptimizationData;

/**
 * Defines an epsilon as a threshold for measuring a small change that approximates convergences of
 * some iterative optimization.
 */
public class ConvergenceEpsilon implements OptimizationData {
    /**
     * The actual epsilon value. Strictly &gt; 0.
     */
    public final double epsilon;

    /**
     * Construct a new ConvergenceEpsilon parameter with the given epsilon value
     *
     * @param eps The epsilon value to use. Must be strictly &gt; 0.
     */
    public ConvergenceEpsilon(double eps) {
        if (eps <= 0.0)
            throw new IllegalArgumentException("epsilon must be > 0");
        this.epsilon = eps;
    }

    /** Default convergence epsilon value */
    public static final double CONVERGENCE_EPSILON_DEFAULT = 1e-9;
}


<|start_filename|>src/main/java/com/manyangled/gibbous/optim/convex/HaltingCondition.java<|end_filename|>
/*
Copyright 2018 <NAME>
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/

package com.manyangled.gibbous.optim.convex;

import org.apache.commons.math3.util.Pair;
import org.apache.commons.math3.optim.OptimizationData;
import org.apache.commons.math3.optim.ConvergenceChecker;
import org.apache.commons.math3.linear.RealVector;

/**
 * Represents a custom halting condition that will cause an optimization to
 * stop iterating if it returns true.
 * <p>
 * An example application is {@link ConvexOptimizer.NegChecker}, which halts
 * an optimization during search for a feasible point if such a point is 
 * discovered.
 * <p>
 * This class is intended primarily for internal use by library algorithms.
 * Use with caution.
 */
public class HaltingCondition implements OptimizationData {
    public final ConvergenceChecker<Pair<RealVector, Double> > checker;

    /**
     * Construct a halting condition from a {@link ConvergenceChecker}
     * @param checker the convergence checker to use as the test for halting condition
     */
    public HaltingCondition(ConvergenceChecker<Pair<RealVector, Double> > checker) {
        this.checker = checker;
    }
}


<|start_filename|>src/main/java/com/manyangled/gibbous/optim/convex/InnerOptimizationData.java<|end_filename|>
/*
Copyright 2018 <NAME>
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/

package com.manyangled.gibbous.optim.convex;

import java.util.ArrayList;
import java.util.Collection;

import org.apache.commons.math3.optim.OptimizationData;

/**
 * A collection of {@link OptimizationData} arguments that are passed to an inner optimizer.
 * <p>
 * For example, the {@link BarrierOptimizer} operates by iterative calls to the
 * {@link NewtonOptimizer}. A user may pass configuration arguments directly to the
 * {@link NewtonOptimizer} by putting them in a {@link InnerOptimizationData} argument
 * while invoking {@link BarrierOptimizer#optimize(OptimizationData... optData)}
 */
public class InnerOptimizationData implements OptimizationData {
    public final ArrayList<OptimizationData> optData =
        new ArrayList<OptimizationData>();

    /**
     * Create an InnerOptimizationData payload from a Collection of OptimizationData
     * @param optData the Collection of OptimizationData arguments
     */
    public InnerOptimizationData(Collection<OptimizationData> optData) {
        this.optData.addAll(optData);
    }

    /**
     * Create an InnerOptimizationData payload from an argument list (or array)
     * of OptimizationData
     * @param optData the list of OptimizationData arguments
     */
    public InnerOptimizationData(OptimizationData... optData) {
        for (OptimizationData data: optData)
            this.optData.add(data);
    }
}


<|start_filename|>src/main/java/com/manyangled/gibbous/optim/convex/ConvexOptimizer.java<|end_filename|>
/*
Copyright 2018 <NAME>
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/

package com.manyangled.gibbous.optim.convex;

import java.util.ArrayList;

import org.apache.commons.math3.optim.OptimizationData;
import org.apache.commons.math3.optim.PointValuePair;
import org.apache.commons.math3.optim.InitialGuess;
import org.apache.commons.math3.optim.nonlinear.scalar.MultivariateOptimizer;
import org.apache.commons.math3.optim.nonlinear.scalar.ObjectiveFunction;
import org.apache.commons.math3.analysis.MultivariateFunction;
import org.apache.commons.math3.util.Pair;
import org.apache.commons.math3.optim.ConvergenceChecker;
import org.apache.commons.math3.linear.RealMatrix;
import org.apache.commons.math3.linear.RealVector;
import org.apache.commons.math3.linear.Array2DRowRealMatrix;
import org.apache.commons.math3.linear.BlockRealMatrix;
import org.apache.commons.math3.linear.ArrayRealVector;

/**
 * An abstract class that represents a {@link MultivariateOptimizer} capable of
 * optimizing convex functions, possibly under a set of convex constraints.
 */
public abstract class ConvexOptimizer extends MultivariateOptimizer {
    protected TwiceDifferentiableFunction convexObjective;

    protected ConvexOptimizer() {
        super(null);
    }

    @Override
    protected void parseOptimizationData(OptimizationData... optData) {
        super.parseOptimizationData(optData);
        for (OptimizationData data: optData) {
            if (data instanceof ObjectiveFunction) {
                MultivariateFunction f = ((ObjectiveFunction)data).getObjectiveFunction();
                if (f instanceof TwiceDifferentiableFunction) {
                  convexObjective = (TwiceDifferentiableFunction)f;
                } else {
                    throw new IllegalArgumentException("TwiceDifferentiableFunction objective required");
                }
                continue;
            }
        }
        if (convexObjective == null)
            throw new IllegalStateException("Expected a TwiceDifferentiableFunction argument");
    }

    @Override
    public PointValuePair optimize(OptimizationData... optData) {
        return super.optimize(optData);
    }

    /**
     * A {@link ConvergenceChecker} specialized to be used as a halting condition when solving
     * for feasible points. Allows the optimization to halt when the objective function becomes
     * negative.
     * <p>
     * NOTE: This class is primarily intended for internal use by
     * {@link #feasiblePoint(OptimizationData... optData)}
     * <p>
     * A useful property is that objective functions having no finite minimum can
     * be handled via a {@link HaltingCondition} having a {@link NegChecker} payload. An example
     * of such a problem is solving a feasible point for a single planar constraint.
     */
    public static class NegChecker implements ConvergenceChecker<Pair<RealVector, Double> > {
        @Override
        public boolean converged(
            int iter,
            Pair<RealVector, Double> prv, Pair<RealVector, Double> cur) {
            return cur.getSecond() < 0.0;
        }
    }

    /**
     * Solve a feasible point for a set of constraints, or determine if the constraints are not feasible.
     * <p>
     * Implements feasible point algorithm described
     * <a href="http://erikerlandson.github.io/blog/2018/06/03/solving-feasible-points-with-smooth-max/">here.</a>
     * <p>
     * Supports the following {@link OptimizationData} parameters as arguments
     * <ul>
     *   <li>initial guess {@link InitialGuess} - optional: defaults to zero-vector. </li>
     *   <li>convex inequality constraints: {@link InequalityConstraintSet} - optional </li>
     *   <li>linear inequality constraints: {@link LinearInequalityConstraint} - optional </li>
     *   <li>linear equality constraints: {@link LinearEqualityConstraint} - optional </li>
     *   <li>convergence epsilon: {@link ConvergenceEpsilon} - optional </li>
     *   <li>inner optimizer parameters: {@link InnerOptimizationData} - optional: passed down to {@link NewtonOptimizer} inner calls. </li>
     * </ul>
     * <p>
     * NOTE: There must be at least one inequality constraint provided, via {@link LinearInequalityConstraint},
     * {@link InequalityConstraintSet}, etc.
     * <p>
     * NOTE: all parameters are also passed to {@link NewtonOptimizer}, and so
     * for example setting {@link ConvergenceEpsilon} here will also set it for inner calls to {@link NewtonOptimizer}. However, any
     * settings passed via {@link InnerOptimizationData} are applied last for {@link NewtonOptimizer}, and so will have precedence.
     * <p>
     * @param optData list of {@link OptimizationData} arguments
     * @return a {@link PointValuePair} where the first element is a feasible point (x), or the point "nearest to feasible"
     * in the sense of minimizing the maximum distance to a constraint surface. The second value is the maximum
     * value f[k](x) over all given constraint functions f[k]. If this value is negative, then (x) is feasible. If the value
     * is &gt; 0, then the constraints cannot be satisfied simultaneously.
     */
    public static PointValuePair feasiblePoint(OptimizationData... optData) {
        double epsilon = ConvergenceEpsilon.CONVERGENCE_EPSILON_DEFAULT;
        RealVector initialGuess = null;
        ArrayList<OptimizationData> barrierArgs = new ArrayList<OptimizationData>();
        ArrayList<TwiceDifferentiableFunction> ineqConstraints =
            new ArrayList<TwiceDifferentiableFunction>();
        final TwiceDifferentiableFunction[] fType = new TwiceDifferentiableFunction[0];
        final OptimizationData[] odType = new OptimizationData[0];
        final ArrayList<OptimizationData> solverArgs = new ArrayList<OptimizationData>();
        final ArrayList<OptimizationData> innerArgs = new ArrayList<OptimizationData>();
        for (OptimizationData data: optData) {
            if (canPassFromMain(data)) {
                solverArgs.add(data);
            }
            if (data instanceof InitialGuess) {
                initialGuess = new ArrayRealVector(((InitialGuess)data).getInitialGuess());
                continue;
            }
            if (data instanceof LinearInequalityConstraint) {
                for (TwiceDifferentiableFunction f: ((LinearInequalityConstraint)data).lcf)
                    ineqConstraints.add(f);
                continue;
            }
            if (data instanceof InequalityConstraintSet) {
                ineqConstraints.addAll(((InequalityConstraintSet)data).constraints);
                continue;
            }
            if (data instanceof ConvergenceEpsilon) {
                epsilon = ((ConvergenceEpsilon)data).epsilon;
                continue;
            }
            if (data instanceof InnerOptimizationData) {
                for (OptimizationData d: ((InnerOptimizationData)data).optData.toArray(odType))
                    if (canPassFromInner(d)) innerArgs.add(d);
                continue;
            }
        }
        // Inner arguments apply after "main" arguments to the inner solver
        solverArgs.addAll(innerArgs);
        // check to see if any linear equality constraints are in effect
        boolean hasLinearInequalityConstraint = false;
        for (OptimizationData data: solverArgs.toArray(odType))
            if (data instanceof LinearEqualityConstraint) hasLinearInequalityConstraint = true;
        if (ineqConstraints.size() < 1)
            throw new IllegalStateException("set of inequality constraints was empty");
        final int n = ineqConstraints.get(0).dim();
        if (initialGuess == null) initialGuess = new ArrayRealVector(n, 0.0);
        final TwiceDifferentiableFunction[] fk = ineqConstraints.toArray(fType);
        // These are free parameters, and might be exposed to a user, but I'm not currently
        // convinced there's a lot of value to tweaking them.
        final double minNBallFactor = Math.log(1e-3);
        // Initialize our location and get the maximum over the constraint functions
        RealVector x = initialGuess;
        double s = fkMax(x.toArray(), fk);
        // If our point is already feasible we are done, unless we need to satisfy
        // linear equality constraints, in which case just run it through the Newton algorithm to make
        // sure they are satisfied
        if (s < 0.0 && !hasLinearInequalityConstraint) return new PointValuePair(x.toArray(), s);
        double alpha = 1.0;
        while (true) {
            // Add the n-ball constraint, to guarantee a non-singular hessian
            // Set the radius to 2x the maximum constraint function value (s). The value (s)
            // corresponds, roughly, to a distance from the current guess to the nearest constraint
            // surface, and so setting the radius to (2s) makes it likely that this ball will overlap
            // the feasible region, to (heuristically) attempt to solve for feasible point in one
            // iteration. Scale 1/(r^2) is set so that the minimum value (center of the ball)
            // remains constant as the radius changes, to prevent it from dominating the computation.
            double radius = Math.max(1.0, 2.0 * s);
            TwiceDifferentiableFunction nbc =
                 QuadraticFunction.nBallConstraintFunction(x, radius, 1.0 / (radius * radius));
            ArrayList<TwiceDifferentiableFunction> augConstraints =
                (ArrayList<TwiceDifferentiableFunction>)(ineqConstraints.clone());
            augConstraints.add(nbc);
            // if necessary, tune the smooth-max alpha to guarantee that our n-ball hessian
            // isn't totally washed out.
            double v0 = nbc.value(x);
            if (v0 < (s + minNBallFactor)) alpha = minNBallFactor / (v0 - s);
            // return the point (x) that minimizes the maximum value of fk(x), and/or (x)
            // where fk(x) is negative for all constraints fk.
            // Smooth-max over fk is always >= to true max, so if smooth-max becomes negative
            // we know max is negative (feasible).
            ArrayList<OptimizationData> args = (ArrayList<OptimizationData>)solverArgs.clone();
            args.add(new InitialGuess(x.toArray()));
            args.add(new ObjectiveFunction(new SmoothMaxFunction(alpha, augConstraints.toArray(fType))));
            args.add(new HaltingCondition(new NegChecker()));
            PointValuePair spvp = (new NewtonOptimizer()).optimize(args.toArray(odType));
            RealVector xprv = x;
            // update our solution x, and the true maximum of constraint functions
            x = new ArrayRealVector(spvp.getFirst());
            s = fkMax(spvp.getFirst(), fk);
            // if our latest x satisfies all contstraints, we can stop
            if (s < 0.0) break;
            RealVector xdelta = x.subtract(xprv);
            // if we are no longer moving, our augmented n-ball constraint is no longer
            // influencing the result, and we've identified our mini-max point, whether
            // it is feasible or not.
            if (xdelta.dotProduct(xdelta) < epsilon) break;
            // increase alpha as we converge, so that smooth-max more closely approximates true max
            // see: http://erikerlandson.github.io/blog/2019/01/02/the-smooth-max-minimum-incident-of-december-2018/
            alpha *= 10.0;
        }
        return new PointValuePair(x.toArray(), s);
    }

    private static boolean canPassFromMain(OptimizationData data) {
        if (data instanceof InitialGuess) return false;
        if (data instanceof ObjectiveFunction) return false;
        if (data instanceof HaltingCondition) return false;
        return true;
    }

    private static boolean canPassFromInner(OptimizationData data) {
        return canPassFromMain(data);
    }

    private static double fkMax(double[] x, TwiceDifferentiableFunction[] fk) {
        double s = Double.NEGATIVE_INFINITY;
        for (TwiceDifferentiableFunction f: fk) {
            double y = f.value(x);
            if (s < y) s = y;
        }
        return s;
    }

    /**
     * Determine if a {@link RealMatrix} subclass is a dense format or not
     * @param M The matrix to test
     * @return true if matrix is known to be a dense format, false otherwise.
     */
    public static boolean isDense(RealMatrix M) {
        if (M instanceof Array2DRowRealMatrix) return true;
        if (M instanceof BlockRealMatrix) return true;
        return false;
    }
}


<|start_filename|>src/test/java/com/manyangled/gibbous/BarrierOptimizerTest.java<|end_filename|>
/*
Copyright 2018 <NAME>
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/

package com.manyangled.gibbous;

import org.junit.Test;
import static org.junit.Assert.*;

import org.apache.commons.math3.linear.RealMatrix;
import org.apache.commons.math3.linear.RealVector;
import org.apache.commons.math3.linear.Array2DRowRealMatrix;
import org.apache.commons.math3.linear.DiagonalMatrix;
import org.apache.commons.math3.linear.ArrayRealVector;

import org.apache.commons.math3.optim.PointValuePair;
import org.apache.commons.math3.optim.nonlinear.scalar.ObjectiveFunction;
import org.apache.commons.math3.optim.InitialGuess;

import com.manyangled.gibbous.optim.convex.*;

import static com.manyangled.gibbous.COTestingUtils.translatedQF;
import static com.manyangled.gibbous.COTestingUtils.eps;

public class BarrierOptimizerTest {
    @Test
    public void testSimpleConstrained2D() {
        double[] center = { 0.0, 0.0 };
        double h = 0.0;
        QuadraticFunction q = translatedQF(h, center);
        double[] ig = { 10.0, 10.0 };
        double[][] A = { { -1.0, -1.0 } }; // constraint x + y > 1
        double[] b = { -1.0 };
        double[] xminTarget = { 0.5, 0.5 };
        double vminTarget = 0.25;
        BarrierOptimizer barrier = new BarrierOptimizer();
        PointValuePair pvp = barrier.optimize(
            new ObjectiveFunction(q),
            new LinearInequalityConstraint(A, b),
            new InitialGuess(ig));
        double[] xmin = pvp.getFirst();
        double vmin = pvp.getSecond();
        assertArrayEquals(xminTarget, xmin, eps);
        assertEquals(vminTarget, vmin, eps);
    }

    @Test
    public void testTranslatedConstrained2D() {
        double[] center = { 3.0, 3.0 };
        double h = 7.0;
        QuadraticFunction q = translatedQF(h, center);
        double[] ig = { 10.0, 10.0 };
        double[][] A = { { -1.0, -1.0 } }; // constraint x + y > 7
        double[] b = { -7.0 };
        double[] xminTarget = { 3.5, 3.5 };
        double vminTarget = 7.25;
        BarrierOptimizer barrier = new BarrierOptimizer();
        PointValuePair pvp = barrier.optimize(
            new ObjectiveFunction(q),
            new LinearInequalityConstraint(A, b),
            new InitialGuess(ig));
        double[] xmin = pvp.getFirst();
        double vmin = pvp.getSecond();
        assertArrayEquals(xminTarget, xmin, eps);
        assertEquals(vminTarget, vmin, eps);
    }

    @Test
    public void testSimpleConstrained3D() {
        double[] center = { 0.0, 0.0, 0.0 };
        double h = 0.0;
        QuadraticFunction q = translatedQF(h, center);
        double[] ig = { 2.0, 2.0, 2.0 };
        double[][] A = { { -1.0, -1.0, -1.0 } }; // constraint x + y + z > 1
        double[] b = { -1.0 };
        double[] xminTarget = { 1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0 };
        double vminTarget = 1.0 / 6.0;
        BarrierOptimizer barrier = new BarrierOptimizer();
        PointValuePair pvp = barrier.optimize(
            new ObjectiveFunction(q),
            new LinearInequalityConstraint(A, b),
            new InitialGuess(ig));
        double[] xmin = pvp.getFirst();
        double vmin = pvp.getSecond();
        assertArrayEquals(xminTarget, xmin, eps);
        assertEquals(vminTarget, vmin, eps);
    }

    @Test
    public void testTranslatedConstrained3D() {
        double[] center = { 10.0, 10.0, 10.0 };
        double h = 10.0;
        QuadraticFunction q = translatedQF(h, center);
        double[] ig = { 15.0, 15.0, 15.0 };
        double[][] A = { { -1.0, -1.0, -1.0 } }; // constraint x + y + z > 31
        double[] b = { -31.0 };
        double[] xminTarget = { 10.0 + (1.0 / 3.0), 10.0 + (1.0 / 3.0), 10.0 + (1.0 / 3.0) };
        double vminTarget = h + (1.0 / 6.0);
        BarrierOptimizer barrier = new BarrierOptimizer();
        PointValuePair pvp = barrier.optimize(
            new ObjectiveFunction(q),
            new LinearInequalityConstraint(A, b),
            new InitialGuess(ig));
        double[] xmin = pvp.getFirst();
        double vmin = pvp.getSecond();
        assertArrayEquals(xminTarget, xmin, eps);
        assertEquals(vminTarget, vmin, eps);
    }

    @Test
    public void testIneqAndEqConstraints2D() {
        QuadraticFunction q = new QuadraticFunction(
            new double[][] { { 1.0, 0.0 }, { 0.0, 1.0 } },
            new double[] { 0.0, 0.0 },
            0.0);
        double[] xminTarget = { 1.0, 1.0 };
        double vminTarget = 1.0;
        BarrierOptimizer barrier = new BarrierOptimizer();
        PointValuePair pvp = barrier.optimize(
            new ObjectiveFunction(q),
            new LinearInequalityConstraint(
                new double[][] { { -1.0, 0.0 } }, // constraint x > 1,
                new double[] { -1.0 }),
            new LinearEqualityConstraint(
                new double[][] { { 0.0, 1.0 } },  // constraint y = 1,
                new double[] { 1.0 }),
            new InitialGuess(new double[] { 10.0, 10.0 }));
        double[] xmin = pvp.getFirst();
        double vmin = pvp.getSecond();
        assertArrayEquals(xminTarget, xmin, eps);
        assertEquals(vminTarget, vmin, eps);
    }

    @Test
    public void testIntegrationWithFPSolver2D() {
        QuadraticFunction q = new QuadraticFunction(
            new double[][] { { 1.0, 0.0 }, { 0.0, 1.0 } },
            new double[] { 0.0, 0.0 },
            0.0);
        double[] xminTarget = { 1.0, 1.0 };
        double vminTarget = 1.0;
        LinearInequalityConstraint ineqc = new LinearInequalityConstraint(
            new double[][] { { -1.0, 0.0 } }, // constraint x > 1,
            new double[] { -1.0 });
        LinearEqualityConstraint eqc = new LinearEqualityConstraint(
            new double[][] { { 0.0, 1.0 } },  // constraint y = 1,
            new double[] { 1.0 });
        PointValuePair fpvp = ConvexOptimizer.feasiblePoint(ineqc, eqc);
        // if not < 0, there is no feasible point
        assertTrue(fpvp.getSecond() < 0.0);
        double[] ig = fpvp.getFirst();
        BarrierOptimizer barrier = new BarrierOptimizer();
        PointValuePair pvp = barrier.optimize(
            new ObjectiveFunction(q),
            ineqc,
            eqc,
            new InitialGuess(ig));
        double[] xmin = pvp.getFirst();
        double vmin = pvp.getSecond();
        assertArrayEquals(xminTarget, xmin, eps);
        assertEquals(vminTarget, vmin, eps);
    }
}


<|start_filename|>src/main/java/com/manyangled/gibbous/optim/convex/BarrierOptimizer.java<|end_filename|>
/*
Copyright 2018 <NAME>
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/

package com.manyangled.gibbous.optim.convex;

import java.util.ArrayList;

import org.apache.commons.math3.util.Pair;
import org.apache.commons.math3.optim.OptimizationData;
import org.apache.commons.math3.optim.nonlinear.scalar.ObjectiveFunction;
import org.apache.commons.math3.optim.InitialGuess;
import org.apache.commons.math3.optim.PointValuePair;
import org.apache.commons.math3.exception.DimensionMismatchException;
import org.apache.commons.math3.linear.RealMatrix;
import org.apache.commons.math3.linear.RealVector;
import org.apache.commons.math3.linear.Array2DRowRealMatrix;
import org.apache.commons.math3.linear.ArrayRealVector;

/**
 * An optimizer that can minimize a convex function in the presence of a set of
 * convex inequality constraints and linear equality constraints.
 * <p>
 * An implementation of the Barrier Method (Algorithm 11.1) from
 * Convex Optimization, Boyd and Vandenberghe, Cambridge University Press, 2008.
 * <p>
 * {@link BarrierOptimizer} supports the following {@link OptimizationData} parameters as arguments
 * to {@link #optimize(OptimizationData...)}:
 * <ul>
 *   <li>convex objective function: {@link ObjectiveFunction} - mandatory: must contain a {@link TwiceDifferentiableFunction} </li>
 *   <li>initial guess: {@link InitialGuess} - mandatory: must be strictly feasible w.r.t. all inequality constraints. </li>
 *   <li>convex inequality constraints: {@link InequalityConstraintSet} - optional </li>
 *   <li>linear inequality constraints: {@link LinearInequalityConstraint} - optional </li>
 *   <li>linear equality constraints: {@link LinearEqualityConstraint} - optional </li>
 *   <li>convergence epsilon: {@link ConvergenceEpsilon} - optional </li>
 *   <li>objective scale factor: {@link BarrierMu} - optional </li>
 *   <li>inner optimizer parameters: {@link InnerOptimizationData} - optional: passed down to {@link NewtonOptimizer} inner calls. </li>
 * </ul>
 * <p>
 * NOTE: all parameters to {@link #optimize(OptimizationData...)} are also passed to {@link NewtonOptimizer}, and so
 * for example setting {@link ConvergenceEpsilon} here will also set it for inner calls to {@link NewtonOptimizer}. However, any
 * settings passed via {@link InnerOptimizationData} are applied last for {@link NewtonOptimizer}, and so will have precedence.
 */
public class BarrierOptimizer extends ConvexOptimizer {
    private ArrayList<TwiceDifferentiableFunction> constraintFunctions =
        new ArrayList<TwiceDifferentiableFunction>();
    private RealVector xStart;
    private double epsilon = ConvergenceEpsilon.CONVERGENCE_EPSILON_DEFAULT;
    private double mu = BarrierMu.BARRIER_MU_DEFAULT;
    private double t0 = BarrierMu.BARRIER_T0_DEFAULT;
    private OptimizationData[] odType = new OptimizationData[0];
    private HaltingCondition halting;
    private ArrayList<OptimizationData> newtonArgs = new ArrayList<OptimizationData>();
    private ArrayList<OptimizationData> innerArgs = new ArrayList<OptimizationData>();

    public BarrierOptimizer() {
        super();
    }

    @Override
    public PointValuePair optimize(OptimizationData... optData) {
        return super.optimize(optData);
    }

    private boolean canPassFromMain(OptimizationData data) {
        if (data instanceof ObjectiveFunction) return false;
        if (data instanceof InitialGuess) return false;
        if (data instanceof HaltingCondition) return false;
        return true;
    }

    private boolean canPassFromInner(OptimizationData data) {
        if (data instanceof HaltingCondition) return true;
        return canPassFromMain(data);
    }

    @Override
    protected void parseOptimizationData(OptimizationData... optData) {
        super.parseOptimizationData(optData);
        // save these for configuring newton optimizers
        for (OptimizationData data: optData) {
            if (canPassFromMain(data)) {
                newtonArgs.add(data);
            }
            if (data instanceof ConvergenceEpsilon) {
                epsilon = ((ConvergenceEpsilon)data).epsilon;
                continue;
            }
            if (data instanceof BarrierMu) {
                mu = ((BarrierMu)data).mu;
                t0 = ((BarrierMu)data).t0;
            }
            if (data instanceof LinearInequalityConstraint) {
                for (TwiceDifferentiableFunction f: ((LinearInequalityConstraint)data).lcf)
                    constraintFunctions.add(f);
                continue;
            }
            if (data instanceof InequalityConstraintSet) {
                constraintFunctions.addAll(((InequalityConstraintSet)data).constraints);
                continue;
            }
            if (data instanceof HaltingCondition) {
                halting = (HaltingCondition)data;
                continue;
            }
            if (data instanceof InnerOptimizationData) {
                for (OptimizationData d: ((InnerOptimizationData)data).optData.toArray(odType))
                    if (canPassFromInner(d)) innerArgs.add(d);
                continue;
            }
        }
        // if we got here, convexObjective exists
        int n = convexObjective.dim();
        if (this.getStartPoint() != null) {
            xStart = new ArrayRealVector(this.getStartPoint());
            if (xStart.getDimension() != n)
                throw new DimensionMismatchException(xStart.getDimension(), n);
        } else {
            xStart = new ArrayRealVector(n, 0.0);
        }
        // append any "inner" args - this overrides anything currently in newtonArgs
        newtonArgs.addAll(innerArgs);
    }

    @Override
    public PointValuePair doOptimize() {
        double m = (double)((constraintFunctions != null) ? constraintFunctions.size() : 0);
        if (m == 0.0) {
            // if there are no inequality constraints, invoke newton's method directly
            ArrayList<OptimizationData> args = (ArrayList<OptimizationData>)newtonArgs.clone();
            args.add(new ObjectiveFunction(convexObjective));
            args.add(new InitialGuess(xStart.toArray()));
            NewtonOptimizer newton = new NewtonOptimizer();
            return newton.optimize(args.toArray(odType));
        }
        RealVector x = xStart;
        for (double t = t0; (t * epsilon) <= m ; t *= mu) {
            TwiceDifferentiableFunction bf = new LogBarrierFunction(t, convexObjective, constraintFunctions);
            NewtonOptimizer newton = new NewtonOptimizer();
            ArrayList<OptimizationData> args = (ArrayList<OptimizationData>)newtonArgs.clone();
            args.add(new ObjectiveFunction(bf));
            args.add(new InitialGuess(x.toArray()));
            PointValuePair pvp = newton.optimize(args.toArray(odType));
            // update for next iteration
            RealVector xprv = x;
            x = new ArrayRealVector(pvp.getFirst());
            if ((halting != null) && halting.checker.converged(
                    getIterations(),
                    new Pair<RealVector, Double>(xprv, convexObjective.value(xprv)),
                    new Pair<RealVector, Double>(x, convexObjective.value(x)))) {
                break;
            }
        }
        return new PointValuePair(x.toArray(), convexObjective.value(x));
    }
}


<|start_filename|>src/main/java/com/manyangled/gibbous/optim/convex/KKTSolution.java<|end_filename|>
/*
Copyright 2018 <NAME>
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/

package com.manyangled.gibbous.optim.convex;

import org.apache.commons.math3.linear.RealVector;

/**
 * This class holds a solution returned from a {@link KKTSolver}.
 * <p>
 * This class is not exposed to the user, and is important only if you are
 * implementing a custom {@link KKTSolver}.
 */
public class KKTSolution {
    /** The delta-x vector solved from the KKT equations */
    public final RealVector xDelta;
    /** The nu+ vector solved from the KKT equations, or null if in (alg 9.5) mode */
    public final RealVector nuPlus;
    /** The lambda-squared value from the KKT equations, or null if in (alg 10.2) mode */
    public final double lambdaSquared;

    /**
     * Construct a KKTSolution for Newton's method with equality constraints.
     * <p>
     * (alg 10.2) of <p>
     * Convex Optimization, Boyd and Vandenberghe, Cambridge University Press, 2008.
     * @param xd the delta-x vector solved from the KKT equation
     * @param nup the nu+ vector solved from the KKT equation
     */
    public KKTSolution(final RealVector xd, final RealVector nup) {
        this.xDelta = xd;
        this.nuPlus = nup;
        this.lambdaSquared = 0.0;
    }

    /**
     * Construct a KKTSolution for Newton's method without equality constraints.
     * <p>
     * (alg 9.5) of <p>
     * Convex Optimization, <NAME> Vandenberghe, Cambridge University Press, 2008.
     * @param xd the delta-x vector solved from the KKT equation
     * @param lsq the lambda-squared value associated with this solution
     */
    public KKTSolution(final RealVector xd, final double lsq) {
        this.xDelta = xd;
        this.lambdaSquared = lsq;
        this.nuPlus = null;
    }
}


<|start_filename|>src/main/java/com/manyangled/gibbous/optim/convex/CholeskySchurKKTSolver.java<|end_filename|>
/*
Copyright 2018 <NAME>
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/

package com.manyangled.gibbous.optim.convex;

import org.apache.commons.math3.linear.RealMatrix;
import org.apache.commons.math3.linear.RealVector;
import org.apache.commons.math3.linear.ArrayRealVector;
import org.apache.commons.math3.linear.CholeskyDecomposition;
import org.apache.commons.math3.linear.DecompositionSolver;

/**
 * Solves KKT conditions for {@link NewtonOptimizer}, using Schur block factorization and
 * Cholesky decomposition.
 * <p>
 * Implements (Algorithm 10.3) and Step 1 of (Algorithm 9.5) from
 * Convex Optimization, Boyd and Vandenberghe, Cambridge University Press, 2008.
 * <p>
 * This is currently the default {@link KKTSolver} for {@link NewtonOptimizer}.
 * <p>
 * See also {@link SVDSchurKKTSolver}
 */
public class CholeskySchurKKTSolver extends KKTSolver {
    // step 1 of algorithm 9.5
    public KKTSolution solve(final RealMatrix H, final RealVector g) {
        DecompositionSolver dsH = (new CholeskyDecomposition(H)).getSolver();
        RealVector v = dsH.solve(g);
        double lsq = g.dotProduct(v);
        v.mapMultiplyToSelf(-1.0);
        return new KKTSolution(v, lsq);
    }

    // Algorithm 10.3
    public KKTSolution solve(
        final RealMatrix H,
        final RealMatrix A, final RealMatrix AT,
        final RealVector g, final RealVector h) {
        DecompositionSolver dsH = (new CholeskyDecomposition(H)).getSolver();
        RealMatrix m1 = dsH.solve(AT);
        RealVector v1 = dsH.solve(g);
        RealMatrix S = A.multiply(m1); // -S relative to 10.3
        DecompositionSolver dsS = (new CholeskyDecomposition(S)).getSolver();
        RealVector w = dsS.solve(h.subtract(A.operate(v1))); // both sides neg, so w same
        RealVector v = dsH.solve(g.add(AT.operate(w))); // this yields -v
        v.mapMultiplyToSelf(-1.0); // correct -v to +v
        return new KKTSolution(v, w);
    }
}


<|start_filename|>src/main/java/com/manyangled/gibbous/optim/convex/BarrierMu.java<|end_filename|>
/*
Copyright 2018 <NAME>
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/

package com.manyangled.gibbous.optim.convex;

import org.apache.commons.math3.optim.OptimizationData;

/**
 * The mu parameter for the Barrier Method used by {@link BarrierOptimizer}.
 * This is the scaling factor for the objective function multiplier (t),
 * as described in (Algorithm 11.1) from <p>
 * Convex Optimization, Boyd and Vandenberghe, Cambridge University Press, 2008.
 */
public class BarrierMu implements OptimizationData {
    /** The value of mu, the scaling factor for the objective function multiplier */
    public final double mu;
    /** The initial value for t, the objective function multiplier */
    public final double t0;

    /**
     * Construct a Barrier Method mu parameter.
     * @param mu the value of the parameter mu, the scaling factor for the objective function multiplier t. Must be &gt; 1.
     * @param t0 the initial value for t, the objective function multiplier. Must be &gt; 0.
     */
    public BarrierMu(double mu, double t0) {
        if (mu <= 1.0) throw new IllegalArgumentException("mu must be > 1");
        if (t0 <= 0.0) throw new IllegalArgumentException("t0 must be > 0");
        this.mu = mu;
        this.t0 = t0;
    }

    /**
     * Construct a Barrier Method mu parameter, the scaling factor for the objective function multiplier t.
     * @param mu the value of the parameter. Must be &gt; 1.
     */
    public BarrierMu(double mu) {
        this(mu, BARRIER_T0_DEFAULT);
    }

    /** Default value for mu */
    public static final double BARRIER_MU_DEFAULT = 15.0;
    /** Default initial value for (t) */
    public static final double BARRIER_T0_DEFAULT = 1.0;
}


<|start_filename|>src/test/java/com/manyangled/gibbous/KKTSolverTest.java<|end_filename|>
/*
Copyright 2018 <NAME>
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/

package com.manyangled.gibbous;

import org.junit.Test;
import static org.junit.Assert.assertEquals;
import static org.junit.Assert.assertArrayEquals;

import org.apache.commons.math3.linear.RealMatrix;
import org.apache.commons.math3.linear.RealVector;
import org.apache.commons.math3.linear.Array2DRowRealMatrix;
import org.apache.commons.math3.linear.ArrayRealVector;

import org.apache.commons.math3.linear.SingularValueDecomposition;

import com.manyangled.gibbous.optim.convex.KKTSolver;
import com.manyangled.gibbous.optim.convex.KKTSolution;
import com.manyangled.gibbous.optim.convex.CholeskySchurKKTSolver;
import com.manyangled.gibbous.optim.convex.SVDSchurKKTSolver;

import static com.manyangled.gibbous.COTestingUtils.eps;

public class KKTSolverTest {
    private RealMatrix inverse(RealMatrix A) {
        SingularValueDecomposition svd = new SingularValueDecomposition(A);
        return svd.getSolver().getInverse();
    }

    private void testNoConstraints(KKTSolver kkts, double[][] Hdata, double[] gdata) {
        RealMatrix H = new Array2DRowRealMatrix(Hdata);
        RealVector g = new ArrayRealVector(gdata);
        KKTSolution sol = kkts.solve(H, g);
        RealVector Hxv = H.operate(sol.xDelta);
        RealVector target = g.mapMultiply(-1.0);
        assertArrayEquals(target.toArray(), Hxv.toArray(), eps);
        RealMatrix Hi = inverse(H);
        double targetLambdaSq = Hi.operate(g).dotProduct(g);
        assertEquals(targetLambdaSq, sol.lambdaSquared, eps);
    }

    private void testWithConstraints(
        KKTSolver kkts,
        double[][] Hdata, double[][] Adata, double[] gdata, double[] hdata) {
        RealMatrix H = new Array2DRowRealMatrix(Hdata);
        RealMatrix A = new Array2DRowRealMatrix(Adata);
        RealVector g = new ArrayRealVector(gdata);
        RealVector h = new ArrayRealVector(hdata);
        RealMatrix AT = A.transpose();
        KKTSolution sol = kkts.solve(H, A, AT, g, h);
        int n = Hdata.length + Adata.length;
        int Hn = Hdata.length;
        double[][] Tdata = new double[n][n];
        for (int j = 0; j < Hn; ++j)
            for (int k = 0; k < Hn; ++k)
                Tdata[j][k] = Hdata[j][k];
        for (int j = 0; j < Adata.length; ++j)
            for (int k = 0; k < Adata[0].length; ++k) {
                Tdata[Hn + j][k] = Adata[j][k];
                Tdata[k][Hn + j] = Adata[j][k];
            }
        RealMatrix T = new Array2DRowRealMatrix(Tdata);
        RealVector t = sol.xDelta.append(sol.nuPlus);
        RealVector Txt = T.operate(t);
        RealVector target = g.append(h).mapMultiply(-1.0);
        assertArrayEquals(target.toArray(), Txt.toArray(), eps);
    }

    @Test
    public void testSchurConstrained1() {
        double[][] H = { { 2.0, 1.0 }, { 1.0, 2.0 } };
        double[][] A = { { 1.0, 1.0 } };
        double[] g = { 1.0, 2.0 };
        double[] h = { 3.0 };
        testWithConstraints(new CholeskySchurKKTSolver(), H, A, g, h);
        testWithConstraints(new SVDSchurKKTSolver(), H, A, g, h);
    }

    @Test
    public void testSchurConstrained2() {
        double[][] H = { { 3.0, 2.0, 1.0 },
                         { 2.0, 3.0, 2.0 },
                         { 1.0, 2.0, 3.0 } };
        double[][] A = { { 1.0, 2.0, 3.0 },
                         { 6.0, 5.0, 4.0 } };
        double[] g = { 1.0, 4.0, 9.0 };
        double[] h = { 3.0, 7.0 };
        testWithConstraints(new CholeskySchurKKTSolver(), H, A, g, h);
        testWithConstraints(new SVDSchurKKTSolver(), H, A, g, h);
    }

    @Test
    public void testSchurUnconstrained1() {
        double[][] H = { { 5.0, 1.0 }, { 1.0, 5.0 } };
        double[] g = { 3.0, 7.0 };
        testNoConstraints(new CholeskySchurKKTSolver(), H, g);
        testNoConstraints(new SVDSchurKKTSolver(), H, g);
    }

    @Test
    public void testSchurUnconstrained2() {
        double[][] H = { { 7.0, 2.0, 1.0 },
                         { 2.0, 7.0, 2.0 },
                         { 1.0, 2.0, 7.0 } };
        double[] g = { 9.0, 4.0, 1.0 };
        testNoConstraints(new CholeskySchurKKTSolver(), H, g);
        testNoConstraints(new SVDSchurKKTSolver(), H, g);
    }
}


<|start_filename|>src/main/java/com/manyangled/gibbous/optim/convex/LinearInequalityConstraint.java<|end_filename|>
/*
Copyright 2018 <NAME>
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/

package com.manyangled.gibbous.optim.convex;

import org.apache.commons.math3.linear.RealMatrix;
import org.apache.commons.math3.linear.RealVector;
import org.apache.commons.math3.linear.Array2DRowRealMatrix;
import org.apache.commons.math3.linear.ArrayRealVector;
import org.apache.commons.math3.optim.OptimizationData;
import org.apache.commons.math3.exception.DimensionMismatchException;

/**
 * A set of linear inequality constraints expressed as Ax &lt; b
 */
public class LinearInequalityConstraint implements OptimizationData {
    /** The corresponding set of individual linear constraint functions */
    public final LinearFunction[] lcf;

    /**
     * Construct a set of linear inequality constraints from Ax &lt; B
     * @param A A matrix linear coefficient vectors
     * @param b A vector of constants
     */
    public LinearInequalityConstraint(final RealMatrix A, final RealVector b) {
        int k = A.getRowDimension();
        if (b.getDimension() != k)
            throw new DimensionMismatchException(b.getDimension(), k);
        this.lcf = new LinearFunction[k];
        for (int j = 0; j < k; ++j)
            lcf[j] = new LinearFunction(A.getRowVector(j), -b.getEntry(j));
    }

    /**
     * Construct a set of linear inequality constraints from Ax &lt; B
     * @param A A matrix linear coefficient vectors
     * @param b A vector of constants
     */
    public LinearInequalityConstraint(final double[][] A, final double[] b) {
        this(new Array2DRowRealMatrix(A), new ArrayRealVector(b));
    }
}


<|start_filename|>src/main/java/com/manyangled/gibbous/optim/convex/LinearEqualityConstraint.java<|end_filename|>
/*
Copyright 2018 <NAME>
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/

package com.manyangled.gibbous.optim.convex;

import org.apache.commons.math3.linear.RealMatrix;
import org.apache.commons.math3.linear.RealVector;
import org.apache.commons.math3.linear.Array2DRowRealMatrix;
import org.apache.commons.math3.linear.ArrayRealVector;
import org.apache.commons.math3.optim.OptimizationData;
import org.apache.commons.math3.exception.DimensionMismatchException;

/**
 * Represents a set of linear equality constraints given as Ax = b.
 */
public class LinearEqualityConstraint implements OptimizationData {
    public final RealMatrix A;
    public final RealVector b;

    /**
     * Construct a set of linear equality constraints Ax = b.
     * Represents equations A[i].x = b[i], for each row of A.
     * @param A the matrix of linear weights
     * @param b the vector of constants
     */
    public LinearEqualityConstraint(final RealMatrix A, final RealVector b) {
        int k = A.getRowDimension();
        if (b.getDimension() != k)
            throw new DimensionMismatchException(b.getDimension(), k);
        this.A = A;
        this.b = b;
    }

    /**
     * Construct a set of linear equality constraints Ax = b.
     * Represents equations A[i].x = b[i], for each row of A.
     * @param A the matrix of linear weights
     * @param b the vector of constants
     */    
    public LinearEqualityConstraint(final double[][] A, final double[] b) {
        this(new Array2DRowRealMatrix(A), new ArrayRealVector(b));
    }
}


<|start_filename|>src/main/java/com/manyangled/gibbous/optim/convex/BacktrackAlpha.java<|end_filename|>
/*
Copyright 2018 <NAME>
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/

package com.manyangled.gibbous.optim.convex;

import org.apache.commons.math3.optim.OptimizationData;

/**
 * The backtracking alpha parameter, used in {@link NewtonOptimizer}.
 * As described in (Algorithm 9.2) of
 * Convex Optimization, Boyd and Vandenberghe, Cambridge University Press, 2008.
 */
public class BacktrackAlpha implements OptimizationData {
    public final double alpha;

    /**
     * Construct a backtracking alpha parameter
     *
     * @param a the value for alpha. Must be on interval (0,1/2)
     */
    public BacktrackAlpha(double a) {
        if ((a <= 0.0) || (a >= 0.5))
            throw new IllegalArgumentException("alpha must be on (0,1/2)");
        this.alpha = a;
    }

    public static final double BACKTRACK_ALPHA_DEFAULT = 0.4;
}


<|start_filename|>src/main/java/com/manyangled/gibbous/optim/convex/KKTSolver.java<|end_filename|>
/*
Copyright 2018 <NAME>
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/

package com.manyangled.gibbous.optim.convex;

import org.apache.commons.math3.linear.RealMatrix;
import org.apache.commons.math3.linear.RealVector;
import org.apache.commons.math3.optim.OptimizationData;

/**
 * Solves the KKT conditions for {@link NewtonOptimizer}.
 * <p>
 * Default solvers make no assumptions about matrix structure.
 * Custom subclasses can be implemented to take advantage of structural
 * knowledge about matrix sparsity.
 */
public abstract class KKTSolver implements OptimizationData {
    /**
     * solve block factored matrix equation:
     * <pre>
     * | H AT | | v | = -| g |
     * | A  0 | | w |    | h |
     * </pre>
     * <p>
     * where (v, w) are primal/dual delta-x and "nu+" from algorithm 10.2 of <p>
     * Convex Optimization, Boyd and Vandenberghe, Cambridge University Press, 2008.
     * @param H Hessian matrix
     * @param A coefficient matrix of equality constraints
     * @param AT transpose of A
     * @param g gradient, corresponding to H
     * @param h constant vector block corresponding to A
     * @return solution (delta-x, nu+)
     */
    public abstract KKTSolution solve(
        final RealMatrix H,
        final RealMatrix A, final RealMatrix AT,
        final RealVector g, final RealVector h);

    /**
     * Solve constraint-free system Hv = -g <p>
     * returns delta-x (aka v) and lambda-squared from Algorithm 9.5 of <p>
     * Convex Optimization, <NAME> Vandenberghe, Cambridge University Press, 2008. <p>
     * delta-nu (aka w, aka the dual) is returned as null
     * @param H Hessian matrix
     * @param g gradient, corresponding to H
     * @return solution delta-x with lambda-squared
     */
    public abstract KKTSolution solve(final RealMatrix H, final RealVector g);
}


<|start_filename|>src/test/java/com/manyangled/gibbous/NewtonOptimizerTest.java<|end_filename|>
/*
Copyright 2018 <NAME>
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/

package com.manyangled.gibbous;

import org.junit.Test;
import static org.junit.Assert.assertEquals;
import static org.junit.Assert.assertArrayEquals;

import org.apache.commons.math3.linear.RealMatrix;
import org.apache.commons.math3.linear.RealVector;
import org.apache.commons.math3.linear.Array2DRowRealMatrix;
import org.apache.commons.math3.linear.DiagonalMatrix;
import org.apache.commons.math3.linear.ArrayRealVector;

import org.apache.commons.math3.optim.PointValuePair;
import org.apache.commons.math3.optim.nonlinear.scalar.ObjectiveFunction;
import org.apache.commons.math3.optim.InitialGuess;

import com.manyangled.gibbous.optim.convex.NewtonOptimizer;
import com.manyangled.gibbous.optim.convex.QuadraticFunction;
import com.manyangled.gibbous.optim.convex.LinearEqualityConstraint;

import static com.manyangled.gibbous.COTestingUtils.translatedQF;
import static com.manyangled.gibbous.COTestingUtils.eps;

public class NewtonOptimizerTest {
    @Test
    public void testSimpleQuadratic2D() {
        double[] center = { 0.0, 0.0 };
        double h = 0.0;
        QuadraticFunction q = translatedQF(h, center);
        double[] ig = { 10.0, 10.0 };
        NewtonOptimizer nopt = new NewtonOptimizer();
        PointValuePair pvp = nopt.optimize(
            new ObjectiveFunction(q),
            new InitialGuess(ig));
        double[] xmin = pvp.getFirst();
        double vmin = pvp.getSecond();
        assertArrayEquals(center, xmin, eps);
        assertEquals(h, vmin, eps);
    }

    @Test
    public void testTranslatedQuadratic2D() {
        double[] center = { 3.0, 7.0 };
        double h = 11.0;
        QuadraticFunction q = translatedQF(h, center);
        NewtonOptimizer nopt = new NewtonOptimizer();
        PointValuePair pvp = nopt.optimize(
            new ObjectiveFunction(q));
        double[] xmin = pvp.getFirst();
        double vmin = pvp.getSecond();
        assertArrayEquals(center, xmin, eps);
        assertEquals(h, vmin, eps);        
    }

    @Test
    public void testTranslatedQuadratic3D() {
        double[] center = { 3.0, 7.0, 11.0 };
        double h = -13.0;
        QuadraticFunction q = translatedQF(h, center);
        NewtonOptimizer nopt = new NewtonOptimizer();
        PointValuePair pvp = nopt.optimize(
            new ObjectiveFunction(q));
        double[] xmin = pvp.getFirst();
        double vmin = pvp.getSecond();
        assertArrayEquals(center, xmin, eps);
        assertEquals(h, vmin, eps);        
    }

    @Test
    public void testSimpleConstrained2D() {
        double[] center = { 0.0, 0.0 };
        double h = 0.0;
        QuadraticFunction q = translatedQF(h, center);
        double[] ig = { 10.0, 10.0 };
        double[][] A = { { 1.0, 1.0 } }; // equality constraint x + y = 1
        double[] b = { 1.0 };
        double[] xminTarget = { 0.5, 0.5 };
        double vminTarget = 0.25;
        NewtonOptimizer nopt = new NewtonOptimizer();
        PointValuePair pvp = nopt.optimize(
            new ObjectiveFunction(q),
            new LinearEqualityConstraint(A, b),
            new InitialGuess(ig));
        double[] xmin = pvp.getFirst();
        double vmin = pvp.getSecond();
        assertArrayEquals(xminTarget, xmin, eps);
        assertEquals(vminTarget, vmin, eps);
    }

    @Test
    public void testTranslatedConstrained2D() {
        double[] center = { 3.0, 3.0 };
        double h = 7.0;
        QuadraticFunction q = translatedQF(h, center);
        double[] ig = { 10.0, 10.0 };
        double[][] A = { { 1.0, 1.0 } }; // equality constraint x + y = 7
        double[] b = { 7.0 };
        double[] xminTarget = { 3.5, 3.5 };
        double vminTarget = 7.25;
        NewtonOptimizer nopt = new NewtonOptimizer();
        PointValuePair pvp = nopt.optimize(
            new ObjectiveFunction(q),
            new LinearEqualityConstraint(A, b),
            new InitialGuess(ig));
        double[] xmin = pvp.getFirst();
        double vmin = pvp.getSecond();
        assertArrayEquals(xminTarget, xmin, eps);
        assertEquals(vminTarget, vmin, eps);
    }

    @Test
    public void testSimpleConstrained3D() {
        double[] center = { 0.0, 0.0, 0.0 };
        double h = 0.0;
        QuadraticFunction q = translatedQF(h, center);
        double[][] A = { { 1.0, 1.0, 1.0 } }; // equality constraint x + y + z = 1
        double[] b = { 1.0 };
        double[] xminTarget = { 1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0 };
        double vminTarget = 1.0 / 6.0;
        NewtonOptimizer nopt = new NewtonOptimizer();
        PointValuePair pvp = nopt.optimize(
            new ObjectiveFunction(q),
            new LinearEqualityConstraint(A, b));
        double[] xmin = pvp.getFirst();
        double vmin = pvp.getSecond();
        assertArrayEquals(xminTarget, xmin, eps);
        assertEquals(vminTarget, vmin, eps);
    }

    @Test
    public void testTranslatedConstrained3D() {
        double[] center = { 10.0, 10.0, 10.0 };
        double h = -100.0;
        QuadraticFunction q = translatedQF(h, center);
        double[][] A = { { 1.0, 1.0, 1.0 } }; // equality constraint x + y + z = 31
        double[] b = { 31.0 };
        double[] xminTarget = { 10.0 + (1.0 / 3.0), 10.0 + (1.0 / 3.0), 10.0 + (1.0 / 3.0) };
        double vminTarget = -100.0 + (1.0 / 6.0);
        NewtonOptimizer nopt = new NewtonOptimizer();
        PointValuePair pvp = nopt.optimize(
            new ObjectiveFunction(q),
            new LinearEqualityConstraint(A, b));
        double[] xmin = pvp.getFirst();
        double vmin = pvp.getSecond();
        assertArrayEquals(xminTarget, xmin, eps);
        assertEquals(vminTarget, vmin, eps);
    }
}


<|start_filename|>src/main/java/com/manyangled/gibbous/optim/convex/BacktrackBeta.java<|end_filename|>
/*
Copyright 2018 <NAME>
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/

package com.manyangled.gibbous.optim.convex;

import org.apache.commons.math3.optim.OptimizationData;

/**
 * The backtracking beta parameter, used in {@link NewtonOptimizer}.
 * As described in (Algorithm 9.2) of
 * Convex Optimization, Boyd and Vandenberghe, Cambridge University Press, 2008.
 */
public class BacktrackBeta implements OptimizationData {
    public final double beta;

    /**
     * Construct a backtracking beta parameter
     *
     * @param b the value for beta. Must be on interval (0,1)
     */
    public BacktrackBeta(double b) {
        if ((b <= 0.0) || (b >= 1.0))
            throw new IllegalArgumentException("beta must be on (0,1)");
        this.beta = b;
    }

    public static final double BACKTRACK_BETA_DEFAULT = 0.8;
}
