我尝试在 java 中使用 matlab fminunc 功能,并在 commons-math 中找到了优化函数。我完全不知道如何正确使用它们,因为该示例没有说明渐变的使用或如何实现它们。任何人都可以给我一些例子吗?
DifferentiableMultivariateFunction 和 NonLinearConjugateGradientOptimizer 的使用
计算科学
优化
matlab
共轭梯度
2021-12-10 05:35:09
2个回答
在这种情况下,梯度通常根据函数的导数进行评估,并与函数最快速减小或增加的方向相关联。如果您正在评估一个定义合理的数学函数,您可能能够自己导出梯度,或者使用 Mathematica 或 SymPy 等工具。一本关于向量微积分的教科书将在这里帮助你理解你在做什么。
如果您知道如何计算雅可比矩阵,请遵循Commons Math 中的综合说明和示例,以派生实现jacobian
和value
方法的类。另一方面,如果您不能对函数进行太多分析,您将需要一个不依赖于梯度的算法,例如 Levenberg-Marquardt,如果您尝试做一些看起来像非线性的事情,这是一个很好的起点最小二乘。
这里有几个使用 Java LevenbergMarquardtOptimizer 类的例子,希望能对你有所帮助。
在一些帮助和提供的答案之后,我为我的目标函数提出了以下代码:
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];
}
private void initCovarianceAndGradients(double[] parameter)
{
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));
//RealMatrix choleskyInput = K.add(k_eye.scalarMultiply(Math.pow(parameter[2], 2)));
RealMatrix choleskyInput = K.add(k_eye);
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 {
this.initCovarianceAndGradients(parameter);
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
public MultivariateVectorialFunction gradient() {
return new MultivariateVectorialFunction() {
public double[] value( double[] parameter) {
initCovarianceAndGradients(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) {
initCovarianceAndGradients(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():初始化边际似然计算和梯度计算所需的矩阵和计算:
在这个函数中,我计算了一些全局的东西,这些东西被 value() 和 gradient() 函数强烈重用。我不太明白的是将 double[] 参数传递给 value() 函数和 gradient() 方法的 value() 函数。优化器是否使用更新的参数调用了这些方法?如果是这种情况,我必须在每次调用 value() 和 gradient() 方法时重新计算全局计算。
感谢您的澄清
更新:该类现在可以了,我在以下过程中运行它:
parameter[2] = 1.0;
parameter[1] = 1.0;
parameter[0] = 1.0;
NonLinearConjugateGradientOptimizer underlying =
new NonLinearConjugateGradientOptimizer(ConjugateGradientFormula.FLETCHER_REEVES);
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();
}
这个过程实际上正在运行,但它收敛于一些负参数(通过方法不可能)。如果实施没问题,我想可能还有其他一些问题,但可以肯定的是,就解决方案发表声明会非常有帮助。
其它你可能感兴趣的问题