qp.hh

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                           */
/*  This file is part of the library KASKADE 7                               */
/*  https://www.zib.de/research/projects/kaskade7-finite-element-toolbox     */
/*                                                                           */
/*  Copyright (C) 2017-2020 Zuse Institute Berlin                            */
/*                                                                           */
/*  KASKADE 7 is distributed under the terms of the ZIB Academic License.    */
/*    see $KASKADE/academic.txt                                              */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#ifndef QP_HH
#define QP_HH
 
#include <memory>
#include <tuple>
 
#include "dune/istl/bvector.hh"
 
#include "linalg/dynamicMatrix.hh"
#include "linalg/threadedMatrix.hh"
 
namespace Kaskade
{
  enum class QPConvexificationStrategy { DONOTHING, INCREMENTALADDITION, EIGENVALUEADDITION };
 
  enum class QPStructure { DENSE, SPARSE };
 
  //-----------------------------------------------------------------------------------------------
  //-----------------------------------------------------------------------------------------------
 
  template <class MatrixA, class MatrixB, class VectorX, class VectorB>
  std::array<typename ScalarTraits<typename MatrixA::field_type>::Real,4>
  checkKKT(MatrixA const& A, MatrixB const& B,VectorX const& c, VectorB const& b,
           VectorX const& x, VectorB const& lambda);
 
 
  //-----------------------------------------------------------------------------------------------
  //-----------------------------------------------------------------------------------------------
  
  template <class Real, class Implementation>
  class QPBoundSolverBase
  {
  public:
    using Vector = Dune::BlockVector<Dune::FieldVector<Real,1>>;
    using Matrix = DynamicMatrix<Dune::FieldMatrix<Real,1,1>>;
 
    virtual ~QPBoundSolverBase()
    {}
    
    std::tuple<Vector,int> solve(Vector const& c, Vector const& b, double tol) const;
 
    std::tuple<Vector,int> solve(Vector const& x, Vector const& c, Vector const& b, double tol) const;
 
    void gradientStep(Vector const& c, Vector const& b, Vector& x) const;
    
    Implementation& setMinSteps(int count);
 
    Implementation& setMaxSteps(int count);
 
  protected:
    QPBoundSolverBase(int dim)
    : dimension(dim) 
    , minSteps(0)
    , maxSteps(200)
    {}
    
    Vector solveSubmatrix(std::vector<int> const& idx, Vector const& b) const;
 
    virtual void updateWithColumn(Real a, int i, Vector& s) const;
    
  private:
    int dimension;      // size of primal variables (i.e. length of x and size of A)
    int minSteps;       // minimal number of steps to perform when solving
    int maxSteps;       // upper bound on iteration count
    
    bool getProjectedGradient(Vector const& x, Vector const& c, Vector const& b, Vector& g, Vector& s) const;
    void getProjectedNewton(Vector const& x, Vector const& c, Vector const& b, Vector& g, Vector& s) const;
    
    bool projectedLineSearch(Vector& x, Vector const& b, Vector& g, Vector& s,
      Vector const* c=nullptr
    ) const;
    void checkNaN(Vector const& v) const;
    
    virtual void printMatrix(std::ostream& out) const = 0;
 
  };
  
  //-----------------------------------------------------------------------------------------------
  //-----------------------------------------------------------------------------------------------
 
  template <class R=double>
  class QPBoundSolver: public QPBoundSolverBase<R,QPBoundSolver<R>>
  {
    using Base = QPBoundSolverBase<R,QPBoundSolver<R>>;
    
  public:
    using Real   = R;
    using Vector = typename Base::Vector;
    using Matrix = DynamicMatrix<Dune::FieldMatrix<Real,1,1>>;
    
    QPBoundSolver(Matrix const& A, QPConvexificationStrategy convex=QPConvexificationStrategy::DONOTHING);
        
  private:
    Matrix A;
 
    friend Base;
    void mv(Vector const& x, Vector& Ax) const { A.mv(x,Ax); }
    Real frobenius_norm() const { return A.frobenius_norm(); }
    Vector solveSubmatrix(std::vector<int> const& idx, Vector const& b) const;
    
    virtual void printMatrix(std::ostream& out) const { out << "A = \n" << A; }
    virtual void updateWithColumn(Real a, int i, Vector& s) const override;
  };
  
  //-----------------------------------------------------------------------------------------------
  //-----------------------------------------------------------------------------------------------
 
  template <class R=double>
  class QPDirectSparse: public QPBoundSolverBase<R,QPDirectSparse<R>>
  {
    using Base = QPBoundSolverBase<R,QPDirectSparse<R>>;
     
    public:
      using Real = R;
      using Vector = typename Base::Vector;       // Dune::BlockVector<Dune::FieldVector<Real,1>>
      using SparseMatrix =  NumaBCRSMatrix<Dune::FieldMatrix<Real,1,1>>;
      
      QPDirectSparse(SparseMatrix const& A, QPConvexificationStrategy convex=QPConvexificationStrategy::DONOTHING);
      
    private:
      SparseMatrix A;
      void mv(Vector const&x, Vector& Ax) const { A.mv(x,Ax); }
      Real frobenius_norm() const;
      friend Base;      
      Vector solveSubmatrix(std::vector<int> const& idx, Vector const& b) const;
      
      virtual void printMatrix(std::ostream& out) const { out << "sparse matrix A in triplet format: \n" << A; }
  };
  
  //--------------------------------------------------------------------------------------------------------------
 
  template <QPStructure sparsity=QPStructure::DENSE, class R=double>
  class QPSlackSolver: public QPBoundSolverBase<R,QPSlackSolver<sparsity,R>>
  {
    using Base = QPBoundSolverBase<R,QPSlackSolver<sparsity,R>>;
    using Entry = Dune::FieldMatrix<R,1,1>;
 
  public:
    using Real = R;
    using Vector = typename Base::Vector;             // Dune::BlockVector<Dune::FieldVector<Real,1>>
    using Matrix = std::conditional_t<sparsity==QPStructure::DENSE,
                                      DynamicMatrix<Entry>,
                                      NumaBCRSMatrix<Entry>>;
    
    QPSlackSolver(Matrix const& A, Matrix const& B, Real gamma, QPConvexificationStrategy convex=QPConvexificationStrategy::DONOTHING);
    
  private:
    Matrix A;
    Matrix B;
    Real gamma;
 
    friend Base;
    void mv(Vector const& x, Vector& Ax) const;
    Real frobenius_norm() const;
    Vector solveSubmatrix(std::vector<int> const& idx, Vector const& b) const;
 
    virtual void printMatrix(std::ostream& out) const { out << "A = \n" << A << "\nB = \n" << B << "\ngamma = " << gamma << "\n"; }
  };
 
  
  //--------------------------------------------------------------------------------------------------------------
  //--------------------------------------------------------------------------------------------------------------
 
   template <int d, class Real=double>
   class QPSolver
   {
     using MatrixM = DynamicMatrix<Dune::FieldMatrix<Real,1,1>>;
 
   public:
     using VectorX = Dune::BlockVector<Dune::FieldVector<Real,d>>;
     using VectorB = Dune::BlockVector<Dune::FieldVector<Real,1>>;
     using MatrixA = DynamicMatrix<Dune::FieldMatrix<Real,d,d>>;
     using MatrixB = DynamicMatrix<Dune::FieldMatrix<Real,1,d>>;
     using Self = QPSolver<d,Real>;
 
     QPSolver(MatrixA const& A, MatrixB const& B, QPConvexificationStrategy convex=QPConvexificationStrategy::DONOTHING);
 
     QPSolver(Self const& qp);
 
     Self& operator=(Self const& qp);
 
     std::tuple<VectorX,VectorB,int> solve(VectorX const& c, VectorB const& b, double tol) const;
 
     MatrixB const& matrixB() const
     {
       return B;
     }
     
      static constexpr QPStructure sparse = QPStructure::DENSE;
 
   private:
     MatrixA Ainv;
     MatrixB B;
     std::unique_ptr<QPBoundSolver<Real>> dual;
   };
 
  //--------------------------------------------------------------------------------------------------------------
  //--------------------------------------------------------------------------------------------------------------
 
   template <int d, QPStructure sparsity = QPStructure::DENSE, class Real=double>
   class QPPenalizedSolver
   {
     using AEntry = Dune::FieldMatrix<Real,d,d>;
     using BEntry = Dune::FieldMatrix<Real,1,d>;
 
   public:
     using VectorX = Dune::BlockVector<Dune::FieldVector<Real,d>>;
     using VectorB = Dune::BlockVector<Dune::FieldVector<Real,1>>;
     using MatrixA = std::conditional_t<sparsity==QPStructure::DENSE,
                                        DynamicMatrix<AEntry>,
                                        NumaBCRSMatrix<AEntry>>;
     using MatrixB = std::conditional_t<sparsity==QPStructure::DENSE,
                                        DynamicMatrix<BEntry>,
                                        NumaBCRSMatrix<BEntry>>;
 
     QPPenalizedSolver(MatrixA const& A, MatrixB const& B);
 
     QPPenalizedSolver(MatrixA const& A, MatrixB const& B, Real gamma);
 
     QPPenalizedSolver<d,sparsity,Real>& setMinSteps(int i);
 
     QPPenalizedSolver<d,sparsity,Real>& setMaxSteps(int i);
 
     QPPenalizedSolver<d,sparsity,Real>& setPenalty(Real gamma);
     
     Real getPenalty() const { return gamma; }
 
     std::tuple<VectorX,VectorB,int> solve(VectorX const& c, VectorB const& b, VectorB const& lambda, double tol) const;
 
 
     std::tuple<VectorX,VectorB,int> solve(VectorX const& c, VectorB const& b,
                                           VectorB lambda, double tol, VectorX const& x) const;
 
     std::tuple<VectorX,int> solve(VectorX const& c, VectorB const& b, double tol, VectorX const& x) const;
     std::tuple<VectorX,int> solve(VectorX const& c, VectorB const& b, double tol) const;
 
     VectorB multiplierEstimate(VectorX const& x, VectorX c) const;
     
      static constexpr QPStructure sparse = sparsity;
 
   private:
     Real gamma = 1e-3;
     MatrixA A;
     MatrixB B;
 
     int minSteps = 0;
     int maxSteps = 10000;
 
     std::vector<int> allCols;
   };
 
   // ---------------------------------------------------------------------------------------------
   // ---------------------------------------------------------------------------------------------
   
   template <int d, QPStructure sparsity = QPStructure::DENSE, class Real=double>
   class QPALSolver
   {
     using AEntry = Dune::FieldMatrix<Real,d,d>;
     using BEntry = Dune::FieldMatrix<Real,1,d>;
 
   public:
     using VectorX = Dune::BlockVector<Dune::FieldVector<Real,d>>;
     using VectorB = Dune::BlockVector<Dune::FieldVector<Real,1>>;
     using MatrixA = std::conditional_t<sparsity==QPStructure::DENSE,
                                        DynamicMatrix<AEntry>,
                                        NumaBCRSMatrix<AEntry>>;
     using MatrixB = std::conditional_t<sparsity==QPStructure::DENSE,
                                        DynamicMatrix<BEntry>,
                                        NumaBCRSMatrix<BEntry>>;
 
     QPALSolver(MatrixA const& A, MatrixB const& B);
 
     QPALSolver<d,sparsity,Real>& setMinSteps(int i);
 
     QPALSolver<d,sparsity,Real>& setMaxSteps(int i);
 
     QPALSolver<d,sparsity,Real>& setPenalty(Real gamma);
 
     Real getPenalty() const { return qp.getPenalty(); }
 
     std::tuple<VectorX,VectorB,int> solve(VectorX const& c, VectorB const& b, double tol) const;
 
 
     std::tuple<VectorX,VectorB,int> solve(VectorX const& c, VectorB const& b, double tol, VectorX x) const;
     
      static constexpr QPStructure sparse = sparsity;
 
   private:
     int minSteps = 0;
     int maxSteps = 10000;
     
     QPPenalizedSolver<d,sparsity,Real> qp;
   };
 
 
}
 
#endif