Optimo  0.1.0
A C++ header library for optimization
 All Classes Functions Variables Pages
sparse_primal_dual_lp_api.h
1 // Copyright (C) 2014 Victor Fragoso <vfragoso@cs.ucsb.edu>
2 // All rights reserved.
3 //
4 // Redistribution and use in source and binary forms, with or without
5 // modification, are permitted provided that the following conditions are
6 // met:
7 //
8 // * Redistributions of source code must retain the above copyright
9 // notice, this list of conditions and the following disclaimer.
10 //
11 // * Redistributions in binary form must reproduce the above
12 // copyright notice, this list of conditions and the following
13 // disclaimer in the documentation and/or other materials provided
14 // with the distribution.
15 //
16 // * Neither the name of the University of California, Santa Barbara nor the
17 // names of its contributors may be used to endorse or promote products
18 // derived from this software without specific prior written permission.
19 //
20 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
21 // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
22 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
23 // ARE DISCLAIMED. IN NO EVENT SHALL VICTOR FRAGOSO BE LIABLE FOR ANY DIRECT,
24 // INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
25 // (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
26 // LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
27 // ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
28 // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
29 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
30 //
31 
32 #ifndef OPTIMO_SOLVERS_SPARSE_PRIMAL_DUAL_LP_API_H_
33 #define OPTIMO_SOLVERS_SPARSE_PRIMAL_DUAL_LP_API_H_
34 
35 #include <glog/logging.h>
36 #include <algorithm>
37 #include <vector>
38 #include <Eigen/Core>
39 #include <Eigen/SparseCore>
40 #include <Eigen/SPQRSupport>
41 #include "optimo/solvers/solver.h"
42 
43 namespace optimo {
44 namespace solvers {
45 
46 using Eigen::Matrix;
47 using Eigen::DiagonalMatrix;
48 using Eigen::Dynamic;
49 using Eigen::SparseMatrix;
50 
51 // Sparse Version
53 
80 template <typename Scalar = double>
81 class SparsePrimalDualLP : public Solver<Scalar> {
82  public:
84  class Params {
85  friend SparsePrimalDualLP;
86  public:
88  Params(const uint n ,
89  const uint m ,
90  const uint p ) :
91  n_(n), m_(m), p_(p), l_(n + m + p),
92  beq(p), bin(m), c(n), F_(n + m + p, n + m + p) {
93  bin.setConstant(0.0);
94  beq.setConstant(0.0);
95  c.setConstant(0.0);
96  }
97 
98  // Destructor
99  virtual ~Params(void) { }
100 
102  inline void reserve(void) {
103  const int nnz = 2*(m_ + p_)*n_ + m_*m_; // Estimated nnz elements
104  F_.reserve(nnz); // Reserve memory
105  }
106 
108  inline bool reserve(const std::vector<uint>& col_sizes) {
109  if (col_sizes.size() != l_) return false;
110  F_.reserve(col_sizes);
111  return true;
112  }
113 
114  // Public members
115  Matrix<Scalar, Dynamic, 1> beq;
116  Matrix<Scalar, Dynamic, 1> bin;
117  Matrix<Scalar, Dynamic, 1> c;
118 
119  // Placing matrices Ain, Aeq
121  inline bool setAeqElement(const uint i, const uint j, const Scalar& val) {
122  if (i >= p_ || j >= n_) return false;
123  if (val == static_cast<Scalar>(0.0)) return false;
124  const uint f_col = n_ + m_ + i; // It is transposed in F
125  const uint f_row = j; // It is transposed
126  F_.coeffRef(f_col, f_row) = val; // To first block-row of F
127  F_.coeffRef(f_row, f_col) = val; // To third block-row of F
128  return true;
129  }
130 
132  inline bool setAinElement(const uint i, const uint j, const Scalar& val) {
133  if (i >= this->m_ || j >= this->n_) return false;
134  if (val == static_cast<Scalar>(0.0)) return false;
135  const uint f_col = n_ + i;
136  const uint f_row = j;
137  F_.coeffRef(f_row, f_col) = val;
138  return true;
139  }
140 
141  private:
142  SparseMatrix<Scalar> F_; // KKT system
143  const uint n_; // Num. of variables
144  const uint m_; // Num. of ineq. constraints
145  const uint p_; // Nun. of eq. constraints
146  const uint l_; // l = n + m + p
147  };
148 
149  // Constructor
150  SparsePrimalDualLP(void) { }
151 
152  // Destructor
153  virtual ~SparsePrimalDualLP(void) {
154  F_ = nullptr;
155  beq_ = nullptr;
156  bin_ = nullptr;
157  c_ = nullptr;
158  y_ = nullptr;
159  }
160 
162  TERMINATION_TYPE
163  operator()(Params* params ,
164  Matrix<Scalar, Dynamic, 1>* y ,
165  Scalar* min_value );
166 
167  protected:
168  inline bool feasibleStartingPoint(const Matrix<Scalar, Dynamic, 1>& y) {
169  // Useful blocks
170  const auto& x = y.block(0, 0, n_, 1); // Primal variables
171  const auto& Ain_t = F_->block(0, n_, n_, m_); // Ain transpose
172  const auto& Aeq_t = F_->block(0, n_ + m_, n_, p_); // Aeq transpose
173  Matrix<Scalar, Dynamic, 1> ineq_x = Ain_t.transpose()*x - *bin_;
174  for (int i = 0; i < m_; i++) if (ineq_x(i) > 0.0) return false;
175  // // TODO(vfragoso): Check if we can ignore the eq. constraints
176  // Matrix<Scalar, Dynamic, 1> eq_x = Aeq_t.transpose()*x - *beq_;
177  // for (int i = 0; i < p_; i++) {
178  // if (fabs(eq_x(i)) > this->options.eps_feas_) {
179  // return false;
180  // }
181  // }
182  return true;
183  }
184 
185  bool backTracking(const Scalar t,
186  const Scalar r_norm,
187  const Matrix<Scalar, Dynamic, 1>& ynt,
188  Matrix<Scalar, Dynamic, 1>* y_plus,
189  Matrix<Scalar, Dynamic, 1>* r_plus,
190  Scalar* s);
191 
192  void calculateResiduals(
193  const Matrix<Scalar, Dynamic, 1>& y,
194  const DiagonalMatrix<Scalar, Dynamic, Dynamic>& diag,
195  const Matrix<Scalar, Dynamic, 1> & fx,
196  const Scalar t_inv,
197  Matrix<Scalar, Dynamic, 1>* residuals);
198 
199  void buildKKTSystem(const double mu,
200  double* eta,
201  double* t,
202  Matrix<Scalar, Dynamic, 1>* residuals,
203  double* rd_norm,
204  double* rp_norm);
205 
206  inline TERMINATION_TYPE
207  mapResult(const Eigen::SPQR<Eigen::SparseMatrix<Scalar> >& solver) {
208  switch (solver.info()) {
209  case Eigen::NumericalIssue:
210  return NUMERICAL_ISSUE;
211  case Eigen::NoConvergence:
212  return NO_CONVERGENCE;
213  case Eigen::InvalidInput:
214  return INVALID_ARGUMENTS;
215  case Eigen::Success:
216  return SOLVED;
217  }
218  }
219 
220  inline bool isF_empty(void) {
221  Scalar acc = static_cast<Scalar>(0.0);
222  for (int i = 0; i < F_->outerSize(); i++) {
223  for (typename SparseMatrix<Scalar>::InnerIterator it(*F_, i); it; ++it) {
224  acc += static_cast<Scalar>(fabs(it.value()));
225  }
226  }
227  return acc == static_cast<Scalar>(0.0);
228  }
229 
230  private:
231  uint n_; // Num. of variables
232  uint m_; // Num. of ineq. constraints
233  uint p_; // Nun. of eq. constraints
234  uint l_; // l = n + m + p
235 
236  // Members so that the user can fill out the matrices
237  // J(x) = c'*x
238  SparseMatrix<Scalar>* F_ = nullptr; // KKT system
239  Matrix<Scalar, Dynamic, 1>* beq_ = nullptr; // Vector of equalities
240  Matrix<Scalar, Dynamic, 1>* bin_ = nullptr; // Vector of equalities
241  Matrix<Scalar, Dynamic, 1>* c_ = nullptr; // Cost vector
242  Matrix<Scalar, Dynamic, 1>* y_ = nullptr; // Unknowns [x lambdas nu]'
243 };
244 } // solvers
245 } // optimo
246 #endif // OPTIMO_SOLVERS_SPARSE_PRIMAL_DUAL_LP_API_H_
Abstract class for a Solver algorithm.
Definition: solver.h:53
Matrix< Scalar, Dynamic, 1 > c
Cost vector of the LP.
Definition: sparse_primal_dual_lp_api.h:117
TERMINATION_TYPE operator()(Params *params, Matrix< Scalar, Dynamic, 1 > *y, Scalar *min_value)
Solves the LP problem.
bool setAeqElement(const uint i, const uint j, const Scalar &val)
Method that inserts an element to matrix .
Definition: sparse_primal_dual_lp_api.h:121
bool setAinElement(const uint i, const uint j, const Scalar &val)
Method that inserts an element to matrix .
Definition: sparse_primal_dual_lp_api.h:132
void reserve(void)
This function will reserve memory for the sparse matrix (optional).
Definition: sparse_primal_dual_lp_api.h:102
Sparse LP params.
Definition: sparse_primal_dual_lp_api.h:84
Matrix< Scalar, Dynamic, 1 > bin
Vector for inequalities.
Definition: sparse_primal_dual_lp_api.h:116
Solves a constrained linear program (LP)
Definition: sparse_primal_dual_lp_api.h:81
Matrix< Scalar, Dynamic, 1 > beq
Vector for equalities.
Definition: sparse_primal_dual_lp_api.h:115
bool reserve(const std::vector< uint > &col_sizes)
Reserve memory per column for the sparse matrix (optional).
Definition: sparse_primal_dual_lp_api.h:108
Params(const uint n, const uint m, const uint p)
LP parameters constructor.
Definition: sparse_primal_dual_lp_api.h:88