optimizer.h

Go to the documentation of this file.

// -*- lsst-c++ -*-
/*
 * LSST Data Management System
 * Copyright 2008-2013 LSST Corporation.
 *
 * This product includes software developed by the
 * LSST Project (http://www.lsst.org/).
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the LSST License Statement and
 * the GNU General Public License along with this program.  If not,
 * see <http://www.lsstcorp.org/LegalNotices/>.
 */

#ifndef LSST_MEAS_MODELFIT_optimizer_h_INCLUDED
#define LSST_MEAS_MODELFIT_optimizer_h_INCLUDED

#include "ndarray.h"

#include "lsst/base.h"
#include "lsst/pex/config.h"
#include "lsst/afw/table/Schema.h"
#include "lsst/meas/modelfit/common.h"
#include "lsst/meas/modelfit/Model.h"

namespace lsst { namespace meas { namespace modelfit {

class Likelihood;
class Prior;
class Optimizer;

class OptimizerObjective {
public:

    int const dataSize;
    int const parameterSize;

    static PTR(OptimizerObjective) makeFromLikelihood(
        PTR(Likelihood) likelihood,
        PTR(Prior) prior = PTR(Prior)()
    );

    OptimizerObjective(int dataSize_, int parameterSize_) :
        dataSize(dataSize_), parameterSize(parameterSize_)
    {}

    void fillObjectiveValueGrid(
        ndarray::Array<Scalar const,2,1> const & parameters,
        ndarray::Array<Scalar,1,1> const & output
    ) const;

    virtual void computeResiduals(
        ndarray::Array<Scalar const,1,1> const & parameters,
        ndarray::Array<Scalar,1,1> const & residuals
    ) const = 0;

    virtual bool differentiateResiduals(
        ndarray::Array<Scalar const,1,1> const & parameters,
        ndarray::Array<Scalar,2,-2> const & derivatives
    ) const {
        return false;
    }


    virtual bool hasPrior() const { return false; }

    virtual Scalar computePrior(ndarray::Array<Scalar const,1,1> const & parameters) const { return 1.0; }

    virtual void differentiatePrior(
        ndarray::Array<Scalar const,1,1> const & parameters,
        ndarray::Array<Scalar,1,1> const & gradient,
        ndarray::Array<Scalar,2,1> const & hessian
    ) const {
        gradient.deep() = 0.0;
        hessian.deep() = 0.0;
    }

    virtual ~OptimizerObjective() {}
};

class OptimizerControl {
public:
    LSST_CONTROL_FIELD(
        noSR1Term, bool,
        "If true, ignore the SR1 update term in the Hessian, resulting in a Levenberg-Marquardt-like method"
    );

    LSST_CONTROL_FIELD(
        skipSR1UpdateThreshold, double,
        "Skip the SR1 update if |v||s| / (|v||s|) is less than this threshold"
    );

    LSST_CONTROL_FIELD(
        minTrustRadiusThreshold, double,
        "If the trust radius falls below this threshold, consider the algorithm converged"
    );

    LSST_CONTROL_FIELD(
        gradientThreshold, double,
        "If the maximum of the gradient falls below this threshold, consider the algorithm converged"
    );

    LSST_CONTROL_FIELD(
        numDiffRelStep, double,
        "relative step size used for numerical derivatives (added to other steps)"
    );

    LSST_CONTROL_FIELD(
        numDiffAbsStep, double,
        "absolute step size used for numerical derivatives (added to other steps)"
    );

    LSST_CONTROL_FIELD(
        numDiffTrustRadiusStep, double,
        "step size (in units of trust radius) used for numerical derivatives (added to relative step)"
    );

    LSST_CONTROL_FIELD(
        stepAcceptThreshold, double,
        "steps with reduction ratio greater than this are accepted"
    );

    LSST_CONTROL_FIELD(
        trustRegionInitialSize, double,
        "the initial trust region will be set to this value"
    );

    LSST_CONTROL_FIELD(
        trustRegionGrowReductionRatio, double,
        "steps with reduction radio greater than this may increase the trust radius"
    );

    LSST_CONTROL_FIELD(
        trustRegionGrowStepFraction, double,
        "steps with length this fraction of the trust radius may increase the trust radius"
    );

    LSST_CONTROL_FIELD(
        trustRegionGrowFactor, double,
        "when increase the trust region size, multiply the radius by this factor"
    );

    LSST_CONTROL_FIELD(
        trustRegionShrinkReductionRatio, double,
        "steps with reduction radio less than this will decrease the trust radius"
    );

    LSST_CONTROL_FIELD(
        trustRegionShrinkFactor, double,
        "when reducing the trust region size, multiply the radius by this factor"
    );

    LSST_CONTROL_FIELD(
        trustRegionSolverTolerance, double,
        "value passed as the tolerance to solveTrustRegion"
    );

    LSST_CONTROL_FIELD(
        maxInnerIterations, int,
        "maximum number of iterations (i.e. function evaluations and trust region subproblems) per step"
    );

    LSST_CONTROL_FIELD(
        maxOuterIterations, int,
        "maximum number of steps"
    );

    LSST_CONTROL_FIELD(
        doSaveIterations, bool,
        "whether to save all iterations for debugging purposes"
    );

    OptimizerControl() :
        noSR1Term(false), skipSR1UpdateThreshold(1E-8),
        minTrustRadiusThreshold(1E-5),
        gradientThreshold(1E-5),
        numDiffRelStep(0.0), numDiffAbsStep(0.0), numDiffTrustRadiusStep(0.1),
        stepAcceptThreshold(0.0),
        trustRegionInitialSize(1.0),
        trustRegionGrowReductionRatio(0.75),
        trustRegionGrowStepFraction(0.8),
        trustRegionGrowFactor(2.0),
        trustRegionShrinkReductionRatio(0.25),
        trustRegionShrinkFactor(1.0/3.0),
        trustRegionSolverTolerance(1E-8),
        maxInnerIterations(20),
        maxOuterIterations(500),
        doSaveIterations(false)
    {}
};

class OptimizerHistoryRecorder {
public:

    OptimizerHistoryRecorder(
        afw::table::Schema & schema,
        PTR(Model) model,
        bool doRecordDerivatives
    );

    explicit OptimizerHistoryRecorder(afw::table::Schema const & schema);

    void apply(
        int outerIterCount,
        int innerIterCount,
        afw::table::BaseCatalog & history,
        Optimizer const & optimizer
    ) const;

    void unpackDerivatives(
        ndarray::Array<Scalar const,1,1> const & nested,
        Vector & gradient,
        Matrix & hessian
    ) const;

    void unpackDerivatives(
        afw::table::BaseRecord const & record,
        Vector & gradient,
        Matrix & hessian
    ) const;

    void unpackDerivatives(
        ndarray::Array<Scalar const,1,1> const & nested,
        ndarray::Array<Scalar,1,1> const & gradient,
        ndarray::Array<Scalar,2,2> const & hessian
    ) const;

    void unpackDerivatives(
        afw::table::BaseRecord const & record,
        ndarray::Array<Scalar,1,1> const & gradient,
        ndarray::Array<Scalar,2,2> const & hessian
    ) const;

    void fillObjectiveModelGrid(
        afw::table::BaseRecord const & record,
        ndarray::Array<Scalar const,2,1> const & parameters,
        ndarray::Array<Scalar,1,1> const & output
    ) const;

    afw::table::Key<int> outer;
    afw::table::Key<int> inner;
    afw::table::Key<int> state;
    ScalarKey objective;
    ScalarKey prior;
    ScalarKey trust;
    ArrayKey parameters;
    ArrayKey derivatives;
};

class Optimizer {
public:

    typedef OptimizerObjective Objective;
    typedef OptimizerControl Control;
    typedef OptimizerHistoryRecorder HistoryRecorder;

    enum StateFlags {
        CONVERGED_GRADZERO = 0x0001,
        CONVERGED_TR_SMALL = 0x0002,
        CONVERGED = CONVERGED_GRADZERO | CONVERGED_TR_SMALL,
        FAILED_MAX_INNER_ITERATIONS = 0x0010,
        FAILED_MAX_OUTER_ITERATIONS = 0x0020,
        FAILED_MAX_ITERATIONS = 0x0030,
        FAILED_EXCEPTION = 0x0040,
        FAILED_NAN = 0x0080,
        FAILED = FAILED_MAX_INNER_ITERATIONS | FAILED_MAX_OUTER_ITERATIONS | FAILED_EXCEPTION | FAILED_NAN,
        STATUS_STEP_REJECTED = 0x0100,
        STATUS_STEP_ACCEPTED = 0x0200,
        STATUS_STEP = STATUS_STEP_REJECTED | STATUS_STEP_ACCEPTED,
        STATUS_TR_UNCHANGED = 0x1000,
        STATUS_TR_DECREASED = 0x2000,
        STATUS_TR_INCREASED = 0x4000,
        STATUS_TR = STATUS_TR_UNCHANGED | STATUS_TR_DECREASED | STATUS_TR_INCREASED,
        STATUS = STATUS_STEP | STATUS_TR,
    };

    Optimizer(
        PTR(Objective const) objective,
        ndarray::Array<Scalar const,1,1> const & parameters,
        Control const & ctrl
    );

    PTR(Objective const) getObjective() const { return _objective; }

    Control const & getControl() const { return _ctrl; }

    bool step() { return _stepImpl(0); }

    bool step(HistoryRecorder const & recorder, afw::table::BaseCatalog & history) {
        return _stepImpl(0, &recorder, &history);
    }

    int run() { return _runImpl(); }

    int run(HistoryRecorder const & recorder, afw::table::BaseCatalog & history) {
        return _runImpl(&recorder, &history);
    }

    int getState() const { return _state; }

    Scalar getObjectiveValue() const { return _current.objectiveValue; }

    ndarray::Array<Scalar const,1,1> getParameters() const { return _current.parameters; }

    ndarray::Array<Scalar const,1,1> getResiduals() const { return _current.residuals; }

    ndarray::Array<Scalar const,1,1> getGradient() const { return _gradient; }

    ndarray::Array<Scalar const,2,2> getHessian() const { return _hessian; }

    void removeSR1Term();

private:

    struct IterationData {
        Scalar objectiveValue;
        Scalar priorValue;
        ndarray::Array<Scalar,1,1> parameters;
        ndarray::Array<Scalar,1,1> residuals;

        IterationData(int dataSize, int parameterSize);

        void swap(IterationData & other);
    };

    friend class OptimizerHistoryRecorder;

    bool _stepImpl(
        int outerIterCount,
        HistoryRecorder const * recorder=NULL,
        afw::table::BaseCatalog * history=NULL
    );

    int _runImpl(HistoryRecorder const * recorder=NULL, afw::table::BaseCatalog * history=NULL);

    void _computeDerivatives();

    int _state;
    PTR(Objective const) _objective;
    Control _ctrl;
    double _trustRadius;
    IterationData _current;
    IterationData _next;
    ndarray::Array<Scalar,1,1> _step;
    ndarray::Array<Scalar,1,1> _gradient;
    ndarray::Array<Scalar,2,2> _hessian;
    ndarray::Array<Scalar,2,-2> _residualDerivative;
    Matrix _sr1b;
    Vector _sr1v;
    Vector _sr1jtr;
};

void solveTrustRegion(
    ndarray::Array<Scalar,1,1> const & x,
    ndarray::Array<Scalar const,2,1> const & F, ndarray::Array<Scalar const,1,1> const & g,
    double r, double tolerance
);

}}} // namespace lsst::meas::modelfit

#endif // !LSST_MEAS_MODELFIT_optimizer_h_INCLUDED