DifferentiableMultivariateFunction 和 NonLinearConjugateGradientOptimizer 的使用

计算科学 优化 matlab 共轭梯度
2021-12-10 05:35:09

我尝试在 java 中使用 matlab fminunc 功能,并在 commons-math 中找到了优化函数。我完全不知道如何正确使用它们,因为该示例没有说明渐变的使用或如何实现它们。任何人都可以给我一些例子吗?

2个回答

在这种情况下,梯度通常根据函数的导数进行评估,并与函数最快速减小或增加的方向相关联。如果您正在评估一个定义合理的数学函数,您可能能够自己导出梯度,或者使用 Mathematica 或 SymPy 等工具。一本关于向量微积分的教科书将在这里帮助你理解你在做什么。

如果您知道如何计算雅可比矩阵,请遵循Commons Math 中的综合说明和示例,以派生实现jacobianvalue方法的类。另一方面,如果您不能对函数进行太多分析,您将需要一个不依赖于梯度的算法,例如 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();
}

这个过程实际上正在运行,但它收敛于一些负参数(通过方法不可能)。如果实施没问题,我想可能还有其他一些问题,但可以肯定的是,就解决方案发表声明会非常有帮助。