drwnQPSolver.h

/******************************************************************************
** DARWIN: A FRAMEWORK FOR MACHINE LEARNING RESEARCH AND DEVELOPMENT
** Distributed under the terms of the BSD license (see the LICENSE file)
** Copyright (c) 2007-2017, Stephen Gould
** All rights reserved.
**
******************************************************************************
** FILENAME:    drwnQPSolver.h
** AUTHOR(S):   Stephen Gould <stephen.gould@anu.edu.au>
**
*****************************************************************************/

#pragma once

#include <cstdlib>
#include <vector>

#include "Eigen/Core"
#include "Eigen/Cholesky"
#include "Eigen/LU"

#include "drwnBase.h"

using namespace std;
using namespace Eigen;

// drwnQPSolver -------------------------------------------------------------

class drwnQPSolver {
 public:
    // line search parameters
    static double alpha; 
    static double beta;  
    
 protected:
    // objective
    MatrixXd _mP;  
    VectorXd _q;   
    double _r;     

    // constraints and bounds
    MatrixXd _mA;  
    VectorXd _b;   
    MatrixXd _mG;  
    MatrixXd _h;   
    VectorXd _l;   
    VectorXd _u;   

    VectorXd _x; 

 public:
    drwnQPSolver();
    drwnQPSolver(const MatrixXd& P, const VectorXd& q, double r = 0.0);
    virtual ~drwnQPSolver();

    // define problem
    void setObjective(const MatrixXd& P, const VectorXd& q, double r = 0.0);
    void setEqConstraints(const MatrixXd& A, const VectorXd& b);
    void setIneqConstraints(const MatrixXd& G, const VectorXd& h);
    void setBounds(const VectorXd& lb, const VectorXd& ub);    
    void clearEqConstraints();
    void clearIneqConstraints();
    void clearBounds();

    virtual void initialize(const VectorXd& x);

    // solve
    virtual double solve();
    double objective() const;
    double objective(const VectorXd& x) const;
    inline VectorXd solution() const { return _x; }

    // accessors
    inline int size() const { return _x.rows(); }    
    inline double operator[](unsigned i) const { return _x[i]; }

 protected:
    // compute gradient: _mP x + _q
    inline VectorXd gradient() const { return _mP * _x + _q; }
    inline VectorXd gradient(const VectorXd& x) const { return _mP * x + _q; }
    
    // check strictly feasibility
    bool isFeasiblePoint(const VectorXd& x) const;

    // special case solvers
    double solveOnlyBounds();
    double solveSingleEquality();
    double solveSimplex();
    double solveNoBounds(); // only equality constraints
    double solveGeneral();

    // line search
    double lineSearchNoBounds(const VectorXd& x, const VectorXd& dx,
        const VectorXd& nu, const VectorXd& dnu) const;
    double lineSearchGeneral(const VectorXd& x, const VectorXd& dx,
        const VectorXd& nu, const VectorXd& dnu) const;

    // solve the system of equations:
    //    | H_x  A^T | | x | = - | c |
    //    |  A   H_y | | y |     | b |
    void solveKKTSystem(const MatrixXd& Hx, const MatrixXd& Hy,
        const MatrixXd& A, const VectorXd& c, const VectorXd& b,
        VectorXd& x, VectorXd& y) const;

    //    | H_x  A^T | | x | = - | c |
    //    |  A   0   | | y |     | b |
    void solveKKTSystem(const MatrixXd& Hx, const MatrixXd& A, 
        const VectorXd& c, const VectorXd& b,
        VectorXd& x, VectorXd& y) const;
};

// drwnLogBarrierQPSolver ---------------------------------------------------

class drwnLogBarrierQPSolver : public drwnQPSolver {
 public:
    // update parameters
    static double t0; // > 0
    static double mu; // > 1
    
 public:
    drwnLogBarrierQPSolver();
    drwnLogBarrierQPSolver(const MatrixXd& P, const VectorXd& q, double r = 0.0);
    virtual ~drwnLogBarrierQPSolver();

    virtual bool findFeasibleStart();

    virtual void initialize(const VectorXd& x);

    // solve
    virtual double solve();
};