# Usage of DifferentiableMultivariateFunction and NonLinearConjugateGradientOptimizer

I try to use the matlab fminunc functionality in java and found the optimisation functions in commons-math. I have absolutely no clue how to use them right because the example tells nothing about the usage of gradients or how to implement them. Can anybody show me some examples?

In this context, gradients are usually evaluated from derivatives of your function and are associated with the direction that your function is most rapidly decreasing or increasing. If you are evaluating a reasonably well-defined mathematical function, you might be able to derive the gradient yourself or using a tool like Mathematica or SymPy. A textbook on vector calculus will help you a great deal here in understanding what you are doing.

If you know how to compute the Jacobian, follow the comprehensive instructions and example at Commons Math for deriving a class that implements the jacobian and value methods. If on the other hand, you cannot say much analytically about the function, you will need an algorithm that does not rely on the gradient, like Levenberg-Marquardt, which is a good starting point if you are trying to do something that looks like nonlinear least-squares.

Here are several examples of using the Java LevenbergMarquardtOptimizer class which hopefully clarify things a bit for you.

• I already know how to derive the gradients and the what the cost function is. I justz need a proper example how to use the differnt classes to implement them. I don't know in which class the target function goes and where to pass it to the optimizer because the example on the Commons Math homepage is just a plain Least Squares Example which is already provided with the proper funtions. – Andreas May 14 '12 at 7:04
• @Andreas - There's a complete example if you have the ability to compute the Jacobian here. Please try to be more clear in the future when asking your questions in letting us know exactly where you are getting stuck by providing more details. Also, please let us know when you have cross-posted your question to another site, as this is something we try to avoid. – Aron Ahmadia May 14 '12 at 9:09
• Ok thanks. I made my way through the documentation and see some differnces now. Regarding my problem i try to compute the optimizer for a gaussian process where i maximize the marginal liklihood for the partial derivitives of a covariance function. I consider now the DifferentiableMultivariateRealFunction for the marginal likelihood function and its derivitives but i'm not sure. Do you think i should reformulate the question in this context or can this be done also within the context of your answer? – Andreas May 14 '12 at 9:58

After some help and the provided answers i came up with the following code for my target function with gradients :

public class GaussianProcessRegressionMarginalLikelihood implements
DifferentiableMultivariateRealFunction {

RealMatrix K;
RealMatrix invK;

RealMatrix y;
CholeskyDecomposition L;
double L_diag;
RealMatrix alpha;

RealMatrix alpha2;

private double[][] cov;
private double [][][] paramCov;

private double[] targets;
private double[] X;

public GaussianProcessRegressionMarginalLikelihood(double[] targets, double[] X) {

this.targets = targets;
this.X = X;

this.cov = new double[targets.length][targets.length];

}

{
double mean = StatUtils.mean(this.targets);
double[] y_array = new double[this.targets.length];

paramCov = new double[parameter.length][][];

for(int i = 0; i < parameter.length; i++)
{
paramCov[i] = new double[targets.length][targets.length];
}

for(int i = 0; i < this.targets.length;i++)
{
y_array[i] = (targets[i]- mean);
//y_array[i] = targets[i]- mean;
}

for(int i = 0; i < this.X.length;i++)
{

for(int j = 0; j < this.X.length;j++)
{

double covar = (Double)squaredExponential.cov(X[i], X[j],parameter[0],parameter[1],parameter[2]);
this.cov[i][j] = covar;

double dSigmaF = (Double)squaredExponential.dSigmaF(X[i], X[j],parameter[0],parameter[1],parameter[2]);
this.paramCov[0][i][j] = dSigmaF;

double dSigmaN = (Double)squaredExponential.dSigmaN(X[i], X[j],parameter[0],parameter[1],parameter[2]);
this.paramCov[1][i][j] = dSigmaN;

double dL = (Double)squaredExponential.dL(X[i], X[j],parameter[0],parameter[1],parameter[2]);
this.paramCov[2][i][j] = dL;

}
}

y = MatrixUtils.createColumnRealMatrix(y_array);

K =  MatrixUtils.createRealMatrix(cov);

//identity matrix for I
RealMatrix k_eye = MatrixUtils.createRealIdentityMatrix(cov.length);

//k_eye = k_eye.scalarMultiply(Math.pow(Float.MIN_VALUE * 1000, 2));
k_eye = k_eye.scalarMultiply(Math.pow(parameter[2], 2));

CholeskyDecomposition L = null;

try {
L = new CholeskyDecompositionImpl(choleskyInput);
} catch (NonSquareMatrixException e) {
e.printStackTrace();
} catch (NotSymmetricMatrixException e) {
e.printStackTrace();
} catch (NotPositiveDefiniteMatrixException e) {
e.printStackTrace();
}

DecompositionSolver solverLTransponse = new LUDecompositionImpl(L.getLT()).getSolver();
DecompositionSolver solverL = new LUDecompositionImpl(L.getL()).getSolver();

//inverse of Ltranspose for left devision
RealMatrix L_transpose_1 = solverLTransponse.getInverse();
//inverse of Ltranspose for left devision
RealMatrix L_1 = solverL.getInverse();

//L\y
RealMatrix L_y = L_1.multiply(y);

//alpha = L'\(L\y)
alpha = L_transpose_1.multiply(L_y);

double L_diag = 0.0;

for(int i = 0; i < L.getL().getColumnDimension();i++)
{
L_diag += Math.log(L.getL().getEntry(i, i));
}

DecompositionSolver solverK = new LUDecompositionImpl(K).getSolver();
this.invK = solverK.getInverse();

alpha2 = invK.multiply(y);
}

/* log p(y|X,theta) = -(1/2) * y^T*Ky^(-1) * y - 1/2 * log * |Ky| - (n/2) * log(2*pi)
*
*/
@Override
public double value(double[] parameter) throws FunctionEvaluationException,
IllegalArgumentException {

double logpyX = - y.transpose().multiply(alpha).getData()[0][0] * 0.5
- L_diag
- X.length * Math.log(2 * Math.PI) * 0.5;

System.out.println("logPYX: " +logpyX + "param: " + parameter);

return logpyX;
}

@Override
return new MultivariateVectorialFunction() {
public double[] value( double[] parameter) {

RealMatrix innerMatrix = alpha2.multiply(alpha2.transpose()).subtract(invK);

double[] result = new double[paramCov.length];

for(int i = 0; i < paramCov.length; i++)
{
result[i] = parameter[i] * innerMatrix.multiply(MatrixUtils.createRealMatrix(paramCov[i])).getTrace() * 0.5;
}

return result;
}
};
}

@Override
public MultivariateRealFunction partialDerivative(final int k) {
return new MultivariateRealFunction() {
public double value(double[] parameter) {

RealMatrix innerMatrix = alpha2.multiply(alpha2.transpose()).subtract(invK);

double[] result = new double[paramCov.length];

for(int i = 0; i < paramCov.length; i++)
{
result[i] = parameter[i] * innerMatrix.multiply(MatrixUtils.createRealMatrix(paramCov[i])).getTrace() * 0.5;
}

return result[k];
}
};
}


}

initCovarianceAndGradients(): initialisation of matrices and calculations which are needed by both marginal likelihood calculation and gradient calculation:

Within this function i calculate some things globally which are strongly reused by the value() and gradient() functions. What i do not really understand is the passing of the double[] argument to the value() function and the value() function of the gradient() method. Are those methods called by the optimizer with the updated parameters? If this is the case i have to recalculate the global calculations with each call to the value() and gradient() methods.

Thanks for clarification

UPDATE: The Class is now ok and i run it within the following process:

parameter[2] = 1.0;
parameter[1] = 1.0;
parameter[0] = 1.0;

JDKRandomGenerator g = new JDKRandomGenerator();
g.setSeed(753289573253l);

RandomVectorGenerator generator =
new UncorrelatedRandomVectorGenerator(3, new GaussianRandomGenerator(g));

MultiStartDifferentiableMultivariateRealOptimizer optimizer =
new MultiStartDifferentiableMultivariateRealOptimizer(underlying, 10, generator);

GaussianProcessRegressionMarginalLikelihood gprml = new GaussianProcessRegressionMarginalLikelihood(input, X);

RealPointValuePair pair = null;

try {
pair = optimizer.optimize(gprml, GoalType.MAXIMIZE, parameter);
} catch (OptimizationException e) {
// TODO Auto-generated catch block
e.printStackTrace();
} catch (FunctionEvaluationException e) {
// TODO Auto-generated catch block
e.printStackTrace();
}


This procedure is actually running but it converges at some negative parameters (impossible by method). If the implementation is alright i guess there might be some other issues but to be sure it would be very helpful to have a statement about the solution.

• Cool, glad you got it working. – Aron Ahmadia May 14 '12 at 14:43
• Unfortunately this code produces values that are not possible (negative values). I have an implementation in octave using fminunc and this works with the derived gradients. So i guess the direction is right but the conversion of the matrix calculations and covariances to java is either wrong or my selection of the optimizer class it not right for the problem. – Andreas May 14 '12 at 14:50