From ffa8c7a757efba5b377a2d1f4f91674b9ff75159 Mon Sep 17 00:00:00 2001 From: Manish Kumar Date: Sun, 28 Jan 2018 16:10:21 +0530 Subject: [PATCH 1/6] Optimized LRSDP function computations. --- .../optimizers/sdp/lrsdp_function_impl.hpp | 27 ++++++++++++++----- 1 file changed, 20 insertions(+), 7 deletions(-) diff --git a/src/mlpack/core/optimizers/sdp/lrsdp_function_impl.hpp b/src/mlpack/core/optimizers/sdp/lrsdp_function_impl.hpp index 30cb237ffe..974180d632 100644 --- a/src/mlpack/core/optimizers/sdp/lrsdp_function_impl.hpp +++ b/src/mlpack/core/optimizers/sdp/lrsdp_function_impl.hpp @@ -47,8 +47,8 @@ LRSDPFunction::LRSDPFunction(const size_t numSparseConstraints, template double LRSDPFunction::Evaluate(const arma::mat& coordinates) const { - const arma::mat rrt = coordinates * trans(coordinates); - return accu(SDP().C() % rrt); + // For computation optimization we will be taking CR first. + return trace((trans(coordinates) * SDP().C()) * coordinates); } template @@ -64,11 +64,13 @@ double LRSDPFunction::EvaluateConstraint( const size_t index, const arma::mat& coordinates) const { - const arma::mat rrt = coordinates * trans(coordinates); + // For computation optimization we will be taking AR first. if (index < SDP().NumSparseConstraints()) - return accu(SDP().SparseA()[index] % rrt) - SDP().SparseB()[index]; + return trace((trans(coordinates) * SDP().SparseA()[index]) * coordinates) + - SDP().SparseB()[index]; const size_t index1 = index - SDP().NumSparseConstraints(); - return accu(SDP().DenseA()[index1] % rrt) - SDP().DenseB()[index1]; + return trace((trans(coordinates) * SDP().DenseA()[index1]) * coordinates) + - SDP().DenseB()[index1]; } template @@ -96,6 +98,9 @@ UpdateObjective(double& objective, for (size_t i = 0; i < ais.size(); ++i) { // Take the trace subtracted by the b_i. + // Here taking AR first is not recommended as we are already + // using pre-computed RR^T. Taking AR first will result in increase + // in number of computations. const double constraint = accu(ais[i] % rrt) - bis[i]; objective -= (lambda[lambdaOffset + i] * constraint); objective += (sigma / 2.) * constraint * constraint; @@ -116,6 +121,9 @@ UpdateGradient(arma::mat& s, { for (size_t i = 0; i < ais.size(); ++i) { + // Here taking AR first is not recommended as we are already + // using pre-computed RR^T. Taking AR first will result in increase + // in number of computations. const double constraint = accu(ais[i] % rrt) - bis[i]; const double y = lambda[lambdaOffset + i] - sigma * constraint; s -= y * ais[i]; @@ -137,13 +145,18 @@ EvaluateImpl(const LRSDPFunction& function, // Let's start with the objective: Tr(C * (R R^T)). // Simple, possibly slow solution-- see below for optimization opportunity // - // TODO: Note that Tr(C^T * (R R^T)) = Tr( (CR)^T * R ), so + // Note that Tr(C^T * (R R^T)) = Tr( (CR)^T * R ), so // multiplying C*R first, and then taking the trace dot should be more memory // efficient // // Similarly for the constraints, taking A*R first should be more efficient + // + // For computation optimization we will be taking CR first. const arma::mat rrt = coordinates * trans(coordinates); - double objective = accu(function.SDP().C() % rrt); + + // Optimized objective function. + double objective = trace((trans(coordinates) * function.SDP().C()) + * coordinates); // Now each constraint. UpdateObjective(objective, rrt, function.SDP().SparseA(), From fd4016f191a93c313506380bf69d27ec711c05ad Mon Sep 17 00:00:00 2001 From: Manish Kumar Date: Tue, 30 Jan 2018 10:18:47 +0530 Subject: [PATCH 2/6] Comments Updated --- .../optimizers/sdp/lrsdp_function_impl.hpp | 22 ++++++++++--------- 1 file changed, 12 insertions(+), 10 deletions(-) diff --git a/src/mlpack/core/optimizers/sdp/lrsdp_function_impl.hpp b/src/mlpack/core/optimizers/sdp/lrsdp_function_impl.hpp index 974180d632..1c65bb8751 100644 --- a/src/mlpack/core/optimizers/sdp/lrsdp_function_impl.hpp +++ b/src/mlpack/core/optimizers/sdp/lrsdp_function_impl.hpp @@ -47,7 +47,7 @@ LRSDPFunction::LRSDPFunction(const size_t numSparseConstraints, template double LRSDPFunction::Evaluate(const arma::mat& coordinates) const { - // For computation optimization we will be taking CR first. + // For computation optimization we will be taking R^T * C first. return trace((trans(coordinates) * SDP().C()) * coordinates); } @@ -64,7 +64,7 @@ double LRSDPFunction::EvaluateConstraint( const size_t index, const arma::mat& coordinates) const { - // For computation optimization we will be taking AR first. + // For computation optimization we will be taking R^T * A first. if (index < SDP().NumSparseConstraints()) return trace((trans(coordinates) * SDP().SparseA()[index]) * coordinates) - SDP().SparseB()[index]; @@ -98,8 +98,8 @@ UpdateObjective(double& objective, for (size_t i = 0; i < ais.size(); ++i) { // Take the trace subtracted by the b_i. - // Here taking AR first is not recommended as we are already - // using pre-computed RR^T. Taking AR first will result in increase + // Here taking R^T * A first is not recommended as we are already + // using pre-computed R * R^T. Taking R^T * A first will result in increase // in number of computations. const double constraint = accu(ais[i] % rrt) - bis[i]; objective -= (lambda[lambdaOffset + i] * constraint); @@ -121,8 +121,8 @@ UpdateGradient(arma::mat& s, { for (size_t i = 0; i < ais.size(); ++i) { - // Here taking AR first is not recommended as we are already - // using pre-computed RR^T. Taking AR first will result in increase + // Here taking R^T * A first is not recommended as we are already + // using pre-computed R * R^T. Taking R^T * A first will result in increase // in number of computations. const double constraint = accu(ais[i] % rrt) - bis[i]; const double y = lambda[lambdaOffset + i] - sigma * constraint; @@ -146,12 +146,14 @@ EvaluateImpl(const LRSDPFunction& function, // Simple, possibly slow solution-- see below for optimization opportunity // // Note that Tr(C^T * (R R^T)) = Tr( (CR)^T * R ), so - // multiplying C*R first, and then taking the trace dot should be more memory - // efficient + // multiplying C * R first, and then taking the trace dot should be more + // memory efficient. // - // Similarly for the constraints, taking A*R first should be more efficient + // Similarly for the constraints, taking R^T * A first should be + // more efficient. // - // For computation optimization we will be taking CR first. + // For computation optimization we will be taking R^T * C first. + // Objective function = Tr((R^T * C) * R) const arma::mat rrt = coordinates * trans(coordinates); // Optimized objective function. From 0b7ba2e4dac9c4b8d529a9cbdcb650a8e1aaac6a Mon Sep 17 00:00:00 2001 From: Manish Date: Sat, 10 Feb 2018 00:54:24 +0530 Subject: [PATCH 3/6] Using Cached RRT() --- .../core/optimizers/sdp/lrsdp_function.hpp | 12 +++++++ .../optimizers/sdp/lrsdp_function_impl.hpp | 35 +++++++++++++++---- 2 files changed, 41 insertions(+), 6 deletions(-) diff --git a/src/mlpack/core/optimizers/sdp/lrsdp_function.hpp b/src/mlpack/core/optimizers/sdp/lrsdp_function.hpp index 00c39f2993..36778ba2aa 100644 --- a/src/mlpack/core/optimizers/sdp/lrsdp_function.hpp +++ b/src/mlpack/core/optimizers/sdp/lrsdp_function.hpp @@ -88,12 +88,24 @@ class LRSDPFunction //! Modify the SDP object representing the problem. SDPType& SDP() { return sdp; } + //! Get R*R^T matrix. + const arma::mat& RRT() const { return rrt; } + + /** + * Update R*R^T matrix. Note caching R*R^T provide computation + * optimization by reducing redundant R*R^T calculations. + */ + void UpdateRRT(const arma::mat& newrrt) const; + private: //! SDP object representing the problem SDPType sdp; //! Initial point. arma::mat initialPoint; + + //! R*R^T matrix. + arma::mat rrt; }; // Declare specializations in lrsdp_function.cpp. diff --git a/src/mlpack/core/optimizers/sdp/lrsdp_function_impl.hpp b/src/mlpack/core/optimizers/sdp/lrsdp_function_impl.hpp index 1c65bb8751..6c27524012 100644 --- a/src/mlpack/core/optimizers/sdp/lrsdp_function_impl.hpp +++ b/src/mlpack/core/optimizers/sdp/lrsdp_function_impl.hpp @@ -29,6 +29,9 @@ LRSDPFunction::LRSDPFunction(const SDPType& sdp, Log::Warn << "LRSDPFunction::LRSDPFunction(): solution matrix will have " << "more columns than rows. It may be more efficient to find the " << "transposed solution." << std::endl; + + // Initialize R*R^T matrix. + rrt = initialPoint * trans(initialPoint); } template @@ -42,13 +45,17 @@ LRSDPFunction::LRSDPFunction(const size_t numSparseConstraints, Log::Warn << "LRSDPFunction::LRSDPFunction(): solution matrix will have " << "more columns than rows. It may be more efficient to find the " << "transposed solution." << std::endl; + + // Initialize R*R^T matrix. + rrt = initialPoint * trans(initialPoint); } template double LRSDPFunction::Evaluate(const arma::mat& coordinates) const { - // For computation optimization we will be taking R^T * C first. - return trace((trans(coordinates) * SDP().C()) * coordinates); + // Note: We don't require to update the R*R^T matrix here as the current + // is used by AugLagrangian, which do not update the coordinate matrix. + return accu(SDP().C() % rrt); } template @@ -64,11 +71,15 @@ double LRSDPFunction::EvaluateConstraint( const size_t index, const arma::mat& coordinates) const { - // For computation optimization we will be taking R^T * A first. + // Note: We don't require to update the R*R^T matrix here as the current + // is used by AugLagrangian, which do not update the coordinate matrix. + + // Using cached R*R^T gives better optimization for sparse matrices. if (index < SDP().NumSparseConstraints()) - return trace((trans(coordinates) * SDP().SparseA()[index]) * coordinates) - - SDP().SparseB()[index]; + return accu(SDP().SparseA()[index] % rrt) - SDP().SparseB()[index]; const size_t index1 = index - SDP().NumSparseConstraints(); + + // For computation optimization we will be taking R^T * A first. return trace((trans(coordinates) * SDP().DenseA()[index1]) * coordinates) - SDP().DenseB()[index1]; } @@ -83,6 +94,12 @@ void LRSDPFunction::GradientConstraint( << "for arbitrary optimizers!" << std::endl; } +template +void LRSDPFunction::UpdateRRT(const arma::mat& newrrt) const +{ + *(arma::mat*)(&rrt) = newrrt; +} + //! Utility function for calculating part of the objective when AugLagrangian is //! used with an LRSDPFunction. template @@ -156,6 +173,12 @@ EvaluateImpl(const LRSDPFunction& function, // Objective function = Tr((R^T * C) * R) const arma::mat rrt = coordinates * trans(coordinates); + // Update R*R^T matrix. + // Note that we can only use this optimization in case of L-BFGS optimizer + // or any other optimizer which calls Evaluate() before Gradient() + // with same coordinates matrix. + function.UpdateRRT(rrt); + // Optimized objective function. double objective = trace((trans(coordinates) * function.SDP().C()) * coordinates); @@ -183,7 +206,7 @@ GradientImpl(const LRSDPFunction& function, // with // S' = C - sum_{i = 1}^{m} y'_i A_i // y'_i = y_i - sigma * (Trace(A_i * (R R^T)) - b_i) - const arma::mat rrt = coordinates * trans(coordinates); + const arma::mat rrt = function.RRT(); arma::mat s(function.SDP().C()); UpdateGradient( From 1153ee59d0f6bb3127cc2363f3af33c7e10e01a1 Mon Sep 17 00:00:00 2001 From: Manish Date: Sat, 10 Feb 2018 15:08:09 +0530 Subject: [PATCH 4/6] Code Optimization --- src/mlpack/core/optimizers/sdp/lrsdp_function.hpp | 7 ++----- .../core/optimizers/sdp/lrsdp_function_impl.hpp | 12 ++++++++---- 2 files changed, 10 insertions(+), 9 deletions(-) diff --git a/src/mlpack/core/optimizers/sdp/lrsdp_function.hpp b/src/mlpack/core/optimizers/sdp/lrsdp_function.hpp index 36778ba2aa..81f419940a 100644 --- a/src/mlpack/core/optimizers/sdp/lrsdp_function.hpp +++ b/src/mlpack/core/optimizers/sdp/lrsdp_function.hpp @@ -91,11 +91,8 @@ class LRSDPFunction //! Get R*R^T matrix. const arma::mat& RRT() const { return rrt; } - /** - * Update R*R^T matrix. Note caching R*R^T provide computation - * optimization by reducing redundant R*R^T calculations. - */ - void UpdateRRT(const arma::mat& newrrt) const; + //! Modify R*R^T matrix. + arma::mat& RRT() { return rrt; } private: //! SDP object representing the problem diff --git a/src/mlpack/core/optimizers/sdp/lrsdp_function_impl.hpp b/src/mlpack/core/optimizers/sdp/lrsdp_function_impl.hpp index 6c27524012..0e27aeac7d 100644 --- a/src/mlpack/core/optimizers/sdp/lrsdp_function_impl.hpp +++ b/src/mlpack/core/optimizers/sdp/lrsdp_function_impl.hpp @@ -94,10 +94,14 @@ void LRSDPFunction::GradientConstraint( << "for arbitrary optimizers!" << std::endl; } +//! Utility function for updating R*R^T matrix. +//! Note: Caching R*R^T provide computation optimization by reducing +//! redundant R*R^T calculations. template -void LRSDPFunction::UpdateRRT(const arma::mat& newrrt) const +void UpdateRRT(LRSDPFunction& function, + const arma::mat& newrrt) { - *(arma::mat*)(&rrt) = newrrt; + function.RRT() = newrrt; } //! Utility function for calculating part of the objective when AugLagrangian is @@ -149,7 +153,7 @@ UpdateGradient(arma::mat& s, template static inline double -EvaluateImpl(const LRSDPFunction& function, +EvaluateImpl(LRSDPFunction& function, const arma::mat& coordinates, const arma::vec& lambda, const double sigma) @@ -177,7 +181,7 @@ EvaluateImpl(const LRSDPFunction& function, // Note that we can only use this optimization in case of L-BFGS optimizer // or any other optimizer which calls Evaluate() before Gradient() // with same coordinates matrix. - function.UpdateRRT(rrt); + UpdateRRT(function, rrt); // Optimized objective function. double objective = trace((trans(coordinates) * function.SDP().C()) From f4a28e07c96febca6bf5b173b8c9e0f61d11fbac Mon Sep 17 00:00:00 2001 From: Manish Date: Sat, 10 Feb 2018 16:12:27 +0530 Subject: [PATCH 5/6] Updated Comments --- .../core/optimizers/sdp/lrsdp_function.hpp | 2 +- .../optimizers/sdp/lrsdp_function_impl.hpp | 24 ++++++++++++++----- 2 files changed, 19 insertions(+), 7 deletions(-) diff --git a/src/mlpack/core/optimizers/sdp/lrsdp_function.hpp b/src/mlpack/core/optimizers/sdp/lrsdp_function.hpp index 81f419940a..084031050e 100644 --- a/src/mlpack/core/optimizers/sdp/lrsdp_function.hpp +++ b/src/mlpack/core/optimizers/sdp/lrsdp_function.hpp @@ -101,7 +101,7 @@ class LRSDPFunction //! Initial point. arma::mat initialPoint; - //! R*R^T matrix. + //! Cache R*R^T matrix. arma::mat rrt; }; diff --git a/src/mlpack/core/optimizers/sdp/lrsdp_function_impl.hpp b/src/mlpack/core/optimizers/sdp/lrsdp_function_impl.hpp index 0e27aeac7d..978c7da33f 100644 --- a/src/mlpack/core/optimizers/sdp/lrsdp_function_impl.hpp +++ b/src/mlpack/core/optimizers/sdp/lrsdp_function_impl.hpp @@ -54,7 +54,8 @@ template double LRSDPFunction::Evaluate(const arma::mat& coordinates) const { // Note: We don't require to update the R*R^T matrix here as the current - // is used by AugLagrangian, which do not update the coordinate matrix. + // function is only used by AugLagrangian, which do not update the coordinates + // matrix. return accu(SDP().C() % rrt); } @@ -72,7 +73,8 @@ double LRSDPFunction::EvaluateConstraint( const arma::mat& coordinates) const { // Note: We don't require to update the R*R^T matrix here as the current - // is used by AugLagrangian, which do not update the coordinate matrix. + // function is only used by AugLagrangian, which do not update the coordinates + // matrix. // Using cached R*R^T gives better optimization for sparse matrices. if (index < SDP().NumSparseConstraints()) @@ -95,8 +97,9 @@ void LRSDPFunction::GradientConstraint( } //! Utility function for updating R*R^T matrix. -//! Note: Caching R*R^T provide computation optimization by reducing -//! redundant R*R^T calculations. +//! Note: Caching R*R^T provide significant computation optimization +//! by reducing redundant R*R^T calculations in case of functions are not used +//! updating coordinates matrix, hence leaving R*R^T unchanged. template void UpdateRRT(LRSDPFunction& function, const arma::mat& newrrt) @@ -175,12 +178,19 @@ EvaluateImpl(LRSDPFunction& function, // // For computation optimization we will be taking R^T * C first. // Objective function = Tr((R^T * C) * R) + + // Calculate R*R^T for updating cache. const arma::mat rrt = coordinates * trans(coordinates); // Update R*R^T matrix. // Note that we can only use this optimization in case of L-BFGS optimizer - // or any other optimizer which calls Evaluate() before Gradient() - // with same coordinates matrix. + // or any other similar optimizer which calls Evaluate() before Gradient() + // with same coordinates matrix and uses only Evaluate() to update + // coordinates matrix. + + // Note: In case optimizer also uses Gradient() for updating coordinates + // matrix than the same line of code can be used to update R*R^T through + // Gradient(). UpdateRRT(function, rrt); // Optimized objective function. @@ -210,6 +220,8 @@ GradientImpl(const LRSDPFunction& function, // with // S' = C - sum_{i = 1}^{m} y'_i A_i // y'_i = y_i - sigma * (Trace(A_i * (R R^T)) - b_i) + + // Directly reterive R*R^T from cache. const arma::mat rrt = function.RRT(); arma::mat s(function.SDP().C()); From 58f096bf7824a620b2d30ac15b1a8cd4322b8729 Mon Sep 17 00:00:00 2001 From: Manish Kumar Date: Tue, 13 Feb 2018 19:17:32 +0530 Subject: [PATCH 6/6] Added comment for LRSDPfunction --- src/mlpack/core/optimizers/sdp/lrsdp_function.hpp | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/src/mlpack/core/optimizers/sdp/lrsdp_function.hpp b/src/mlpack/core/optimizers/sdp/lrsdp_function.hpp index 084031050e..1ebb7369dd 100644 --- a/src/mlpack/core/optimizers/sdp/lrsdp_function.hpp +++ b/src/mlpack/core/optimizers/sdp/lrsdp_function.hpp @@ -22,6 +22,15 @@ namespace optimization { /** * The objective function that LRSDP is trying to optimize. + * + * Note: LRSDPfunction is designed and implemented to specifically work + * with a combination of AugLagrangian and L-BFGS optimizer. Implemenation + * of LRSDPFunction includes caching R * R^T matrix in order to avoid redundant + * computations(Specifically with above two optimizers) of R * R^T matrix + * through AugLagrangian optimizer call, as only L-BFGS takes part in updating + * R(coordinates) matrix. So, be careful while using LRSDP with some other + * optimizer. You may need to modify caching process of R * R^T matrix. + * See EvaluateImpl() in lrsdp_function_impl.hpp for more details. */ template class LRSDPFunction