ReClaM: QuadraticProgram.h Source File

00001 //===========================================================================
00054 //===========================================================================
00055 
00056 
00057 #ifndef _QuadraticProgram_H_
00058 #define _QuadraticProgram_H_
00059 
00060 
00061 #include <SharkDefs.h>
00062 #include <ReClaM/Model.h>
00063 #include <ReClaM/KernelFunction.h>
00064 #include <Array/Array2D.h>
00065 #include <vector>
00066 
00067 
00069 class QPSolver
00070 {
00071 public:
00072     QPSolver();
00073     virtual ~QPSolver();
00074 };
00075 
00076 
00085 class QPMatrix
00086 {
00087 public:
00089     QPMatrix(unsigned int size);
00090 
00092     virtual ~QPMatrix();
00093 
00102     virtual float Entry(unsigned int i, unsigned int j) = 0;
00103 
00115     virtual void FlipColumnsAndRows(unsigned int i, unsigned int j) = 0;
00116 
00118     inline unsigned int getMatrixSize() const
00119     {
00120         return matrixsize;
00121     }
00122 
00123 protected:
00125     unsigned int matrixsize;
00126 };
00127 
00128 
00135 class KernelMatrix : public QPMatrix
00136 {
00137 public:
00141     KernelMatrix(KernelFunction* kernelfunction,
00142                  const Array<double>& data);
00143 
00145     ~KernelMatrix();
00146 
00148     float Entry(unsigned int i, unsigned int j);
00149 
00151     void FlipColumnsAndRows(unsigned int i, unsigned int j);
00152 
00153 protected:
00155     KernelFunction* kernel;
00156 
00158     Array<ArrayReference<double>* > x;
00159 };
00160 
00161 
00169 class RegularizedKernelMatrix : public KernelMatrix
00170 {
00171 public:
00176     RegularizedKernelMatrix(KernelFunction* kernel,
00177                             const Array<double>& data,
00178                             const Array<double>& diagModification);
00179 
00181     ~RegularizedKernelMatrix();
00182 
00184     float Entry(unsigned int i, unsigned int j);
00185 
00187     void FlipColumnsAndRows(unsigned int i, unsigned int j);
00188 
00189 protected:
00191     Array<double> diagMod;
00192 };
00193 
00194 
00204 class QPMatrix2 : public QPMatrix
00205 {
00206 public:
00212     QPMatrix2(QPMatrix* base);
00213 
00215     ~QPMatrix2();
00216 
00218     float Entry(unsigned int i, unsigned int j);
00219 
00221     void FlipColumnsAndRows(unsigned int i, unsigned int j);
00222 
00223 protected:
00225     QPMatrix* baseMatrix;
00226 
00228     Array<unsigned int> mapping;
00229 };
00230 
00231 
00239 class InputLabelMatrix : public QPMatrix
00240 {
00241 public:
00247     InputLabelMatrix(KernelFunction* kernelfunction,
00248                             const Array<double>& input,
00249                             const Array<double>& target,
00250                             const Array<double>& prototypes);
00251 
00253     ~InputLabelMatrix();
00254 
00256     float Entry(unsigned int i, unsigned int j);
00257 
00259     void FlipColumnsAndRows(unsigned int i, unsigned int j);
00260 
00261 protected:
00263     KernelFunction* kernel;
00264 
00266     Array<ArrayReference<double>* > x;
00267 
00269     std::vector<unsigned int> y;
00270 
00272     Array<double> prototypeProduct;
00273 };
00274 
00275 
00289 class CachedMatrix : public QPMatrix
00290 {
00291 public:
00295     CachedMatrix(QPMatrix* base, unsigned int cachesize = 0x4000000);
00296 
00298     ~CachedMatrix();
00299 
00312     float* Row(unsigned int k, unsigned int begin, unsigned int end, bool temp = false);
00313 
00315     float Entry(unsigned int i, unsigned int j);
00316 
00328     void FlipColumnsAndRows(unsigned int i, unsigned int j);
00329 
00330     inline unsigned int getMaxCacheSize() const
00331     {
00332         return cacheMaxSize;
00333     }
00334 
00335     inline unsigned int getCacheSize() const
00336     {
00337         return cacheSize;
00338     }
00339 
00340     inline unsigned int getCacheRowSize(unsigned int k)
00341     {
00342         return cacheEntry[k].length;
00343     }
00344 
00345     inline void Clear()
00346     { cacheClear(); }
00347 
00348     void CacheRowResize(unsigned int k, unsigned int newsize);
00349     void CacheRowRelease(unsigned int k);
00350 
00351 protected:
00353     QPMatrix* baseMatrix;
00354 
00355     // matrix cache
00356     unsigned int cacheSize;                     // current cache size in floats
00357     unsigned int cacheMaxSize;                  // maximum cache size in floats
00358     struct tCacheEntry                          // cache data held for every example
00359     {
00360         float* data;                            // float array containing a matrix row
00361         int length;                             // length of this matrix row
00362         int older;                              // next older entry
00363         int newer;                              // next newer entry
00364     };
00365     std::vector<tCacheEntry> cacheEntry;        // cache entry description
00366     std::vector<float> cacheTemp;               // single kernel row
00367     int cacheNewest;                            // index of the newest entry
00368     int cacheOldest;                            // index of the oldest entry
00369 
00371     void cacheAppend(int var);
00372 
00374     void cacheRemove(int var);
00375 
00378     void cacheAdd(int var, unsigned int length);
00379 
00381     void cacheDelete(int var);
00382 
00384     void cacheResize(int var, unsigned int newlength);
00385 
00387     void cacheClear();
00388 };
00389 
00390 
00447 class QpSvmCG : public QPSolver
00448 {
00449 public:
00456     static void Solve(const Array<double>& quadratic,
00457             const Array<double>& linear,
00458             const Array<double>& lower,
00459             const Array<double>& upper,
00460             Array<double>& point,
00461             double accuracy = 1e-10);
00462 };
00463 
00464 
00526 class QpSvmDecomp : public QPSolver
00527 {
00528 public:
00531     QpSvmDecomp(CachedMatrix& quadraticPart);
00532 
00534     virtual ~QpSvmDecomp();
00535 
00555     double Solve(const Array<double>& linearPart,
00556                  const Array<double>& boxLower,
00557                  const Array<double>& boxUpper,
00558                  Array<double>& solutionVector,
00559                  double eps = 0.001,
00560                  double threshold = 1e100);
00561 
00576     double ComputeInnerProduct(unsigned int index, const Array<double>& coeff);
00577 
00582     void getGradient(Array<double>& grad);
00583 
00585     inline void setVerbose(bool verbose = false)
00586     {
00587         printInfo = verbose;
00588     }
00589 
00591     inline void setStrategy(const char* strategy = NULL)
00592     {
00593         WSS_Strategy = strategy;
00594     }
00595 
00597     inline void setShrinking(bool shrinking = true)
00598     {
00599         useShrinking = shrinking;
00600     }
00601 
00603     inline SharkInt64 iterations()
00604     {
00605         return iter;
00606     }
00607 
00609     inline bool isOptimal()
00610     {
00611         return optimal;
00612     }
00613 
00617     inline void setMaxIterations(SharkInt64 maxiter = -1)
00618     {
00619         this->maxIter = maxiter;
00620     }
00621 
00625     inline void setMaxSeconds(int seconds = -1)
00626     {
00627         maxSeconds = seconds;
00628     }
00629 
00630 protected:
00632     unsigned int dimension;
00633 
00635     CachedMatrix& quadratic;
00636 
00638     Array<double> linear;
00639 
00641     Array<double> boxMin;
00642 
00644     Array<double> boxMax;
00645 
00647     SharkInt64 iter;
00648 
00652     bool optimal;
00653 
00655     SharkInt64 maxIter;
00656 
00658     int maxSeconds;
00659 
00661     bool printInfo;
00662 
00664     const char* WSS_Strategy;
00665 
00667     bool useShrinking;
00668 
00670     bool(QpSvmDecomp::*currentWSS)(unsigned int&, unsigned int&);
00671 
00673     unsigned int active;
00674 
00676     Array<unsigned int> permutation;
00677 
00680     Array<double> diagonal;
00681 
00684     Array<double> gradient;
00685 
00687     Array<double> alpha;
00688 
00690     bool bFirst;
00691 
00693     unsigned int old_i;
00694 
00696     unsigned int old_j;
00697 
00699     double epsilon;
00700 
00706     bool MVP(unsigned int& i, unsigned int& j);
00707 
00713     bool HMG(unsigned int& i, unsigned int& j);
00714 
00720     bool Libsvm28(unsigned int& i, unsigned int& j);
00721 
00729     virtual bool SelectWorkingSet(unsigned int& i, unsigned int& j);
00730 
00732     void SelectWSS();
00733 
00735     void Shrink();
00736 
00738     void Unshrink(bool complete = false);
00739 
00741     bool bUnshrinked;
00742 
00744     void FlipCoordinates(unsigned int i, unsigned int j);
00745 };
00746 
00747 
00757 class QpBoxDecomp : public QPSolver
00758 {
00759 public:
00762     QpBoxDecomp(CachedMatrix& quadraticPart);
00763 
00765     virtual ~QpBoxDecomp();
00766 
00776     void Solve(const Array<double>& linearPart,
00777                  const Array<double>& boxLower,
00778                  const Array<double>& boxUpper,
00779                  Array<double>& solutionVector,
00780                  double eps = 0.001);
00781 
00795     void Continue(const Array<double>& boxLower,
00796                  const Array<double>& boxUpper,
00797                  Array<double>& solutionVector);
00798 
00813     void Continue(const Array<double>& gradientUpdate,
00814                  Array<double>& solutionVector);
00815 
00817     inline SharkInt64 iterations()
00818     {
00819         return iter;
00820     }
00821 
00823     inline bool isOptimal()
00824     {
00825         return optimal;
00826     }
00827 
00831     inline void setMaxIterations(SharkInt64 maxiter = -1)
00832     {
00833         this->maxIter = maxiter;
00834     }
00835 
00836 protected:
00838     void SMO();
00839 
00841     unsigned int dimension;
00842 
00844     CachedMatrix& quadratic;
00845 
00847     Array<double> linear;
00848 
00850     Array<double> boxMin;
00851 
00853     Array<double> boxMax;
00854 
00856     unsigned int active;
00857 
00859     Array<unsigned int> permutation;
00860 
00863     Array<double> diagonal;
00864 
00867     Array<double> gradient;
00868 
00870     Array<double> alpha;
00871 
00873     double epsilon;
00874 
00875     virtual bool SelectWorkingSet(unsigned int& i);
00876 
00878     void Shrink();
00879 
00881     void Unshrink(bool complete = false);
00882 
00884     bool bUnshrinked;
00885 
00887     SharkInt64 iter;
00888 
00892     bool optimal;
00893 
00895     SharkInt64 maxIter;
00896 
00898     void FlipCoordinates(unsigned int i, unsigned int j);
00899 };
00900 
00901 
00912 class QpBoxAllInOneDecomp : public QPSolver
00913 {
00914 public:
00918     QpBoxAllInOneDecomp(CachedMatrix& kernel);
00919 
00921     virtual ~QpBoxAllInOneDecomp();
00922 
00933     void Solve(unsigned int classes,
00934                 const Array<double>& target,
00935                 const Array<double>& prototypes,
00936                 double C,
00937                 Array<double>& solutionVector,
00938                 double eps = 0.001);
00939 
00957     void Continue(double largerC, Array<double>& solutionVector);
00958 
00960     inline SharkInt64 iterations()
00961     {
00962         return iter;
00963     }
00964 
00966     inline bool isOptimal()
00967     {
00968         return optimal;
00969     }
00970 
00974     inline void setMaxIterations(SharkInt64 maxiter = -1)
00975     {
00976         this->maxIter = maxiter;
00977     }
00978 
00979 protected:
00981     struct tExample
00982     {
00983         unsigned int index;         // example index
00984         unsigned int label;         // label of this example
00985         unsigned int active;        // number of active variables
00986         unsigned int* variables;    // list of variables, active ones first
00987     };
00988 
00990     struct tVariable
00991     {
00992         unsigned int example;       // index into the example list
00993         unsigned int index;         // index into example->variables
00994         unsigned int label;         // label corresponding to this variable
00995     };
00996 
00998     std::vector<tExample> example;
00999 
01001     std::vector<tVariable> variable;
01002 
01004     std::vector<unsigned int> storage;
01005 
01006     unsigned int examples;
01007     unsigned int classes;
01008     unsigned int variables;
01009 
01011     CachedMatrix& kernelMatrix;
01012 
01014     Array<double> linear;
01015 
01017     Array<double> boxMin;
01018 
01020     Array<double> boxMax;
01021 
01023     unsigned int activeEx;
01024 
01026     unsigned int activeVar;
01027 
01030     Array<double> diagonal;
01031 
01034     Array<double> gradient;
01035 
01037     Array<double> alpha;
01038 
01040     double epsilon;
01041 
01043     double Modifier(unsigned int v_i, unsigned int v_j, unsigned int y_i, unsigned int y_j);
01044 
01046     bool SelectWorkingSet(unsigned int& i);
01047 
01049     void Shrink();
01050 
01052     void Unshrink(bool complete = false);
01053 
01055     bool bUnshrinked;
01056 
01058     void DeactivateVariable(unsigned int v);
01059 
01061     void DeactivateExample(unsigned int e);
01062 
01064     Array<double> m_prototypeProduct;
01065 
01067     Array<double> m_prototypeLength;
01068 
01070     SharkInt64 iter;
01071 
01075     bool optimal;
01076 
01078     SharkInt64 maxIter;
01079 };
01080 
01081 
01095 class QpCrammerSingerDecomp : public QPSolver
01096 {
01097 public:
01102     QpCrammerSingerDecomp(CachedMatrix& kernel, const Array<double>& y, unsigned int classes);
01103 
01105     virtual ~QpCrammerSingerDecomp();
01106 
01114     void Solve(Array<double>& solutionVector,
01115                 double beta,
01116                 double eps = 0.001);
01117 
01119     inline SharkInt64 iterations()
01120     {
01121         return iter;
01122     }
01123 
01125     inline bool isOptimal()
01126     {
01127         return optimal;
01128     }
01129 
01133     inline void setMaxIterations(SharkInt64 maxiter = -1)
01134     {
01135         this->maxIter = maxiter;
01136     }
01137 
01138 protected:
01140     struct tExample
01141     {
01142         unsigned int index;         // example index
01143         unsigned int label;         // label of this example
01144         unsigned int active;        // number of active variables
01145         unsigned int* variables;    // list of variables, active ones first
01146     };
01147 
01149     struct tVariable
01150     {
01151         unsigned int example;       // index into the example list
01152         unsigned int index;         // index into example->variables
01153         unsigned int label;         // label corresponding to this variable
01154     };
01155 
01157     std::vector<tExample> example;
01158 
01160     std::vector<tVariable> variable;
01161 
01163     std::vector<unsigned int> storage;
01164 
01165     unsigned int examples;
01166     unsigned int classes;
01167     unsigned int variables;
01168 
01170     CachedMatrix& kernelMatrix;
01171 
01173     Array<double> linear;
01174 
01176     Array<double> boxMin;
01177 
01179     Array<double> boxMax;
01180 
01182     unsigned int activeEx;
01183 
01185     unsigned int activeVar;
01186 
01189     Array<double> diagonal;
01190 
01193     Array<double> gradient;
01194 
01196     Array<double> alpha;
01197 
01199     double epsilon;
01200 
01201     double SelectWorkingSet(unsigned int& i, unsigned int& j);
01202 
01204     void Shrink();
01205 
01207     void Unshrink(bool complete = false);
01208 
01210     bool bUnshrinked;
01211 
01213     void DeactivateVariable(unsigned int v);
01214 
01216     void DeactivateExample(unsigned int e);
01217 
01219     SharkInt64 iter;
01220 
01224     bool optimal;
01225 
01227     SharkInt64 maxIter;
01228 };
01229 
01230 
01232 class QpBoxAndEqCG : public QPSolver
01233 {
01234 public:
01241     static void Solve(const Array<double>& quadratic,
01242             const Array<double>& linear,
01243             const Array<double>& boxMin,
01244             const Array<double>& boxMax,
01245             const Array<double>& eqMat,
01246             Array<double>& point,
01247             double accuracy = 1e-10);
01248 
01249 protected:
01250     static void Orthogonalize(int oc, const Array<double>& ortho, const Array<bool>& active, ArrayReference<double> vec);
01251     static void Orthogonalize(Array<double>& eq, const Array<bool>& active);
01252     static void Project(const Array<bool>& active, Array<double>& eq, const Array<double>& gradient, Array<double>& direction);
01253 };
01254 
01255 
01257 class QpBoxAndEqDecomp : public QPSolver
01258 {
01259 public:
01266     static void Solve(const Array<double>& quadratic,
01267             const Array<double>& linear,
01268             const Array<double>& boxMin,
01269             const Array<double>& boxMax,
01270             const Array<double>& eqMat,
01271             Array<double>& point,
01272             double accuracy = 1e-3);
01273 
01274 protected:
01275     static void Orthogonalize(Array<double>& eq);
01276     static void Project(const Array<double>& eq, const Array<double>& gradient, Array<double>& direction);
01277 };
01278 
01279 
01280 // //!
01281 // //! \brief Quadratic program solver with affine linear constraints
01282 // //!
01283 // //! \par
01284 // //! The purpose of the QpAffine class is two-fold:
01285 // //! First, a QpAffine object is a description of a
01286 // //! quadratic program, that is, an optimization problem
01287 // //! composed of a quadratic objective function and a set
01288 // //! of affine linear constraints. Second, this class is
01289 // //! a solver for this kind of problem.
01290 // //!
01291 // //! \par
01292 // //! This class is designed for the representation and
01293 // //! solution of a quite general family of quadratic
01294 // //! programs. These problems have the following
01295 // //! structure (for \f$ \alpha \in R^{\ell} \f$):
01296 // //!
01297 // //! \par
01298 // //! minimize \f$ W(\alpha) = \frac{1}{2} \alpha^T M \alpha + v^T \alpha\f$<br>
01299 // //! s.t. \f$ E \alpha + e = 0 \f$ (equality constraints)<br>
01300 // //! and \f$ I \leq \alpha + i \leq 0 \f$ component wise (inequality constraints).<br>
01301 // //! Here, E is an \f$ n \times \ell \f$ matrix,
01302 // //! e is an n-dimensional vector, I is an
01303 // //! \f$ m \times \ell \f$ matrix and i is an
01304 // //! m-dimensional vector.
01305 // //! The \f$ \ell \times \ell \f$ matrix M is required
01306 // //! to be symmetric and positive definite, while there
01307 // //! are no constraints for the vector v.
01308 // //! In practive the vector e is not needed as long as
01309 // //! the optimization starts at a feasible point.
01310 // //!
01311 // //! \par
01312 // //! The solution strategy followed by the QpAffine
01313 // //! solver is (probably, I don't have a proof) cubic
01314 // //! in the number of variables, as long as the number
01315 // //! of inequality constraints is linear in the number
01316 // //! of variables. This excludes too large programs.
01317 // //! The algorithm itself imposes an even stronger
01318 // //! constraint, because the Gram-Schmidt
01319 // //! orthogonalization and singular value decomposition
01320 // //! algorithms used are numerically sensitive.
01321 // //! Therefore the solver is suited only for small
01322 // //! problems or, in case of well conditioned problems,
01323 // //! for moderate size problems.
01324 // //!
01325 // //! \par
01326 // //! In theory, the solution computed by QpAffine is
01327 // //! exact. The solver internally solves a transformed
01328 // //! problem and afterwards backtransforms the solution
01329 // //! found. These transformations can decrease the
01330 // //! accuracy. This is usually only visible at the
01331 // //! contraints which are not fulfilled exactly.
01332 // //!
01333 // class QpAffine : public QPSolver
01334 // {
01335 // public:
01336 //  //! Constructor defining the quadratic program.
01337 //  //! Important preprocessing steps are precomputed
01338 //  //! to facilitate storage of the problem and the
01339 //  //! computations necessary to solve the problem.
01340 //  //!
01341 //  //! \param  quadratic  Matrix M, see class description
01342 //  //! \param  linear     vector v, see class description
01343 //  //! \param  eqMat      Matrix E, see class description
01344 //  //! \param  ineqMat    Matrix I, see class description
01345 //  //! \param  ineqVec    Vector i, see class description
01346 //  QpAffine(const Array<double>& quadratic,
01347 //           const Array<double>& linear,
01348 //           const Array<double>& eqMat,
01349 //           const Array<double>& ineqMat,
01350 //           const Array<double>& ineqVec);
01351 // 
01352 //  //! Destructor
01353 //  virtual ~QpAffine();
01354 // 
01355 //  //!
01356 //  //! \brief solve the quadratic program
01357 //  //!
01358 //  //! \par
01359 //  //! This is the core method of the #QpAffine class.
01360 //  //! It computes the solution of the problem defined
01361 //  //! by the QpAffine instance.
01362 //  //!
01363 //  //! \param point   input: initial feasible vector \f$ \alpha \f$; output: solution \f$ \alpha^* \f$
01364 //  //!
01365 //  void Solve(Array<double>& point);
01366 // 
01367 // protected:
01368 //  //! dimension of the original problem
01369 //  int dimension;
01370 // 
01371 //  //! dimension of the transformed problem
01372 //  int freedim;
01373 // 
01374 //  //! problem transformation
01375 //  Array<double> transform;
01376 // 
01377 //  //! inverse problem transformation
01378 //  Array<double> inverse;
01379 // 
01380 //  //! quardatic part of the original problem
01381 //  const Array<double>& quadraticOrig;
01382 // 
01383 //  //! linear part of the transformed problem
01384 //  Array<double> linear;
01385 // 
01386 //  //! inequality constraint matrix of the original problem
01387 //  Array<double> constraintMatOrig;
01388 // 
01389 //  //! inequality constraint matrix of the transformed problem
01390 //  Array<double> constraintMatTrans;
01391 // 
01392 //  //! inequality constraint vector of the original problem
01393 //  const Array<double>& constraintVecOrig;
01394 // 
01395 //  //! Gram-Schmidt algorithm
01396 //  //! \param  ortho  set of orthonormal vectors
01397 //  //! \param  vec    input: vector to process, output: orthogonal vector
01398 //  static void Orthogonalize(const std::vector<double*>& ortho, Array<double>& vec);
01399 // 
01400 //  //! Gram-Schmidt algorithm
01401 //  //! \param  ortho  set of orthonormal vectors
01402 //  //! \param  vec    input: vector to process, output: orthogonal vector
01403 //  static void Orthogonalize(const std::vector<double*>& ortho, ArrayReference<double> vec);
01404 // 
01405 //  //! Gram-Schmidt algorithm
01406 //  //! \param  vec    input: vector to process, output: unit vector
01407 //  //! \return  true on success, false if vec has zero length
01408 //  static bool Normalize(Array<double>& vec);
01409 // 
01410 //  //! Vector Normalization
01411 //  //! \param  vec    input: vector to process, output: unit vector
01412 //  //! \return  true on success, false if vec has zero length
01413 //  static bool Normalize(ArrayReference<double> vec);
01414 // 
01415 //  //! result := M1 * M2
01416 //  static void MatMul(const Array<double>& M1, const Array<double>& M2, Array<double>& result);
01417 // 
01418 //  //! result := Mat * vec
01419 //  static void MatVec(const Array<double>& Mat, const Array<double>& vec, Array<double>& result);
01420 // 
01421 //  //! result := Mat * vec
01422 //  static void MatVec(const Array<double>& Mat, const ArrayReference<double> vec, ArrayReference<double> result);
01423 // 
01424 //  //! result := <vec1, vec2>
01425 //  static double Scp(const Array<double>& vec1, ArrayReference<double> vec2);
01426 // };
01427 
01428 
01429 #endif
01430