Add accessors and mutators for AugLagrangian and L-BFGS. Change int to size_t
where necessary.
This commit is contained in:
@@ -38,16 +38,26 @@ template<typename LagrangianFunction>
|
||||
class AugLagrangian
|
||||
{
|
||||
public:
|
||||
AugLagrangian(LagrangianFunction& function_in, int num_basis);
|
||||
AugLagrangian(LagrangianFunction& function, size_t numBasis);
|
||||
// not sure what to do here yet
|
||||
|
||||
bool Optimize(int num_iterations,
|
||||
bool Optimize(size_t num_iterations,
|
||||
arma::mat& coordinates,
|
||||
double sigma = 0.5);
|
||||
|
||||
//! Get the LagrangianFunction.
|
||||
const LagrangianFunction& Function() const { return function; }
|
||||
//! Modify the LagrangianFunction.
|
||||
LagrangianFunction& Function() { return function; }
|
||||
|
||||
//! Get the number of memory points used by L-BFGS.
|
||||
size_t NumBasis() const { return numBasis; }
|
||||
//! Modify the number of memory points used by L-BFGS.
|
||||
size_t& NumBasis() { return numBasis; }
|
||||
|
||||
private:
|
||||
LagrangianFunction& function_;
|
||||
int num_basis_;
|
||||
LagrangianFunction& function;
|
||||
size_t numBasis;
|
||||
|
||||
/**
|
||||
* This is a utility class, which we will pass to L-BFGS during the
|
||||
@@ -68,11 +78,26 @@ class AugLagrangian
|
||||
|
||||
const arma::mat& GetInitialPoint();
|
||||
|
||||
arma::vec lambda_;
|
||||
double sigma_;
|
||||
//! Get the Lagrangian multipliers.
|
||||
const arma::vec& Lambda() const { return lambda; }
|
||||
//! Modify the Lagrangian multipliers.
|
||||
arma::vec& Lambda() { return lambda; }
|
||||
|
||||
//! Get sigma.
|
||||
double Sigma() const { return sigma; }
|
||||
//! Modify sigma.
|
||||
double& Sigma() { return sigma; }
|
||||
|
||||
//! Get the Lagrangian function.
|
||||
const LagrangianFunction& Function() const { return function; }
|
||||
//! Modify the Lagrangian function.
|
||||
LagrangianFunction& Function() { return function; }
|
||||
|
||||
private:
|
||||
LagrangianFunction& function_;
|
||||
arma::vec lambda;
|
||||
double sigma;
|
||||
|
||||
LagrangianFunction& function;
|
||||
};
|
||||
};
|
||||
|
||||
|
||||
@@ -16,33 +16,33 @@ namespace optimization {
|
||||
|
||||
template<typename LagrangianFunction>
|
||||
AugLagrangian<LagrangianFunction>::AugLagrangian(
|
||||
LagrangianFunction& function_in, int num_basis) :
|
||||
function_(function_in),
|
||||
num_basis_(num_basis)
|
||||
LagrangianFunction& function, size_t numBasis) :
|
||||
function(function),
|
||||
numBasis(numBasis)
|
||||
{
|
||||
// Not sure what to do here (if anything).
|
||||
}
|
||||
|
||||
template<typename LagrangianFunction>
|
||||
bool AugLagrangian<LagrangianFunction>::Optimize(int num_iterations,
|
||||
bool AugLagrangian<LagrangianFunction>::Optimize(size_t numIterations,
|
||||
arma::mat& coordinates,
|
||||
double sigma)
|
||||
{
|
||||
// Choose initial lambda parameters (vector of zeros, for simplicity).
|
||||
arma::vec lambda(function_.NumConstraints());
|
||||
arma::vec lambda(function.NumConstraints());
|
||||
lambda.ones();
|
||||
double penalty_threshold = DBL_MAX; // Ensure we update lambda immediately.
|
||||
|
||||
// Track the last objective to compare for convergence.
|
||||
double last_objective = function_.Evaluate(coordinates);
|
||||
double last_objective = function.Evaluate(coordinates);
|
||||
|
||||
// First, we create an instance of the utility function class.
|
||||
AugLagrangianFunction f(function_, lambda, sigma);
|
||||
AugLagrangianFunction f(function, lambda, sigma);
|
||||
|
||||
// First, calculate the current penalty.
|
||||
double penalty = 0;
|
||||
for (int i = 0; i < function_.NumConstraints(); i++)
|
||||
penalty += std::pow(function_.EvaluateConstraint(i, coordinates), 2);
|
||||
for (size_t i = 0; i < function.NumConstraints(); i++)
|
||||
penalty += std::pow(function.EvaluateConstraint(i, coordinates), 2);
|
||||
|
||||
Log::Debug << "Penalty is " << penalty << " (threshold " << penalty_threshold
|
||||
<< ")." << std::endl;
|
||||
@@ -56,18 +56,18 @@ bool AugLagrangian<LagrangianFunction>::Optimize(int num_iterations,
|
||||
<< ", starting with objective " << last_objective << "." << std::endl;
|
||||
|
||||
// Use L-BFGS to optimize this function for the given lambda and sigma.
|
||||
L_BFGS<AugLagrangianFunction> lbfgs(f, num_basis_);
|
||||
L_BFGS<AugLagrangianFunction> lbfgs(f, numBasis);
|
||||
if (!lbfgs.Optimize(0, coordinates))
|
||||
Log::Debug << "L-BFGS reported an error during optimization."
|
||||
<< std::endl;
|
||||
|
||||
// Check if we are done with the entire optimization (the threshold we are
|
||||
// comparing with is arbitrary).
|
||||
if (std::abs(last_objective - function_.Evaluate(coordinates)) < 1e-10 &&
|
||||
if (std::abs(last_objective - function.Evaluate(coordinates)) < 1e-10 &&
|
||||
sigma > 500000)
|
||||
return true;
|
||||
|
||||
last_objective = function_.Evaluate(coordinates);
|
||||
last_objective = function.Evaluate(coordinates);
|
||||
|
||||
// Assuming that the optimization has converged to a new set of coordinates,
|
||||
// we now update either lambda or sigma. We update sigma if the penalty
|
||||
@@ -75,18 +75,18 @@ bool AugLagrangian<LagrangianFunction>::Optimize(int num_iterations,
|
||||
|
||||
// First, calculate the current penalty.
|
||||
double penalty = 0;
|
||||
for (int i = 0; i < function_.NumConstraints(); i++)
|
||||
penalty += std::pow(function_.EvaluateConstraint(i, coordinates), 2);
|
||||
for (int i = 0; i < function.NumConstraints(); i++)
|
||||
penalty += std::pow(function.EvaluateConstraint(i, coordinates), 2);
|
||||
|
||||
Log::Debug << "Penalty is " << penalty << " (threshold "
|
||||
<< penalty_threshold << ")." << std::endl;
|
||||
if (penalty < penalty_threshold) // We update lambda.
|
||||
{
|
||||
// We use the update: lambda_{k + 1} = lambda_k - sigma * c(coordinates),
|
||||
// We use the update: lambda{k + 1} = lambdak - sigma * c(coordinates),
|
||||
// but we have to write a loop to do this for each constraint.
|
||||
for (int i = 0; i < function_.NumConstraints(); i++)
|
||||
lambda[i] -= sigma * function_.EvaluateConstraint(i, coordinates);
|
||||
f.lambda_ = lambda;
|
||||
for (int i = 0; i < function.NumConstraints(); i++)
|
||||
lambda[i] -= sigma * function.EvaluateConstraint(i, coordinates);
|
||||
f.lambda = lambda;
|
||||
|
||||
// We also update the penalty threshold to be a factor of the current
|
||||
// penalty. TODO: this factor should be a parameter (from CLI). The
|
||||
@@ -100,7 +100,7 @@ bool AugLagrangian<LagrangianFunction>::Optimize(int num_iterations,
|
||||
// parameter (from CLI). The value of 10 is taken from Burer and Monteiro
|
||||
// (2002).
|
||||
sigma *= 10;
|
||||
f.sigma_ = sigma;
|
||||
f.sigma = sigma;
|
||||
Log::Debug << "Updated sigma to " << sigma << "." << std::endl;
|
||||
}
|
||||
}
|
||||
@@ -111,10 +111,10 @@ bool AugLagrangian<LagrangianFunction>::Optimize(int num_iterations,
|
||||
|
||||
template<typename LagrangianFunction>
|
||||
AugLagrangian<LagrangianFunction>::AugLagrangianFunction::AugLagrangianFunction(
|
||||
LagrangianFunction& function_in, arma::vec& lambda_in, double sigma) :
|
||||
lambda_(lambda_in),
|
||||
sigma_(sigma),
|
||||
function_(function_in)
|
||||
LagrangianFunction& functionin, arma::vec& lambdain, double sigma) :
|
||||
lambda(lambdain),
|
||||
sigma(sigma),
|
||||
function(functionin)
|
||||
{
|
||||
// Nothing to do.
|
||||
}
|
||||
@@ -124,16 +124,16 @@ double AugLagrangian<LagrangianFunction>::AugLagrangianFunction::Evaluate(
|
||||
const arma::mat& coordinates)
|
||||
{
|
||||
// The augmented Lagrangian is evaluated as
|
||||
// f(x) + {-lambda_i * c_i(x) + (sigma / 2) c_i(x)^2} for all constraints
|
||||
// f(x) + {-lambdai * c_i(x) + (sigma / 2) c_i(x)^2} for all constraints
|
||||
// Log::Debug << "Evaluating augmented Lagrangian." << std::endl;
|
||||
double objective = function_.Evaluate(coordinates);
|
||||
double objective = function.Evaluate(coordinates);
|
||||
|
||||
// Now loop over constraints.
|
||||
for (int i = 0; i < function_.NumConstraints(); i++)
|
||||
for (int i = 0; i < function.NumConstraints(); i++)
|
||||
{
|
||||
double constraint = function_.EvaluateConstraint(i, coordinates);
|
||||
objective += (-lambda_[i] * constraint) +
|
||||
sigma_ * std::pow(constraint, 2) / 2;
|
||||
double constraint = function.EvaluateConstraint(i, coordinates);
|
||||
objective += (-lambda[i] * constraint) +
|
||||
sigma * std::pow(constraint, 2) / 2;
|
||||
}
|
||||
|
||||
// Log::Warn << "Overall objective is " << objective << "." << std::endl;
|
||||
@@ -146,24 +146,24 @@ void AugLagrangian<LagrangianFunction>::AugLagrangianFunction::Gradient(
|
||||
const arma::mat& coordinates, arma::mat& gradient)
|
||||
{
|
||||
// The augmented Lagrangian's gradient is evaluated as
|
||||
// f'(x) + {(-lambda_i + sigma * c_i(x)) * c'_i(x)} for all constraints
|
||||
// f'(x) + {(-lambdai + sigma * c_i(x)) * c'_i(x)} for all constraints
|
||||
// gradient.zeros();
|
||||
function_.Gradient(coordinates, gradient);
|
||||
function.Gradient(coordinates, gradient);
|
||||
// Log::Debug << "Objective function gradient norm is "
|
||||
// << arma::norm(gradient, 2) << "." << std::endl;
|
||||
// std::cout << gradient << std::endl;
|
||||
|
||||
arma::mat constraint_gradient; // Temporary for constraint gradients.
|
||||
for (int i = 0; i < function_.NumConstraints(); i++)
|
||||
for (int i = 0; i < function.NumConstraints(); i++)
|
||||
{
|
||||
function_.GradientConstraint(i, coordinates, constraint_gradient);
|
||||
function.GradientConstraint(i, coordinates, constraint_gradient);
|
||||
|
||||
// Now calculate scaling factor and add to existing gradient.
|
||||
arma::mat tmp_gradient;
|
||||
tmp_gradient = (-lambda_[i] + sigma_ *
|
||||
function_.EvaluateConstraint(i, coordinates)) * constraint_gradient;
|
||||
tmp_gradient = (-lambda[i] + sigma *
|
||||
function.EvaluateConstraint(i, coordinates)) * constraint_gradient;
|
||||
// Log::Debug << "Gradient for constraint " << i << " (with lambda = "
|
||||
// << lambda_[i] << ") is " << std::endl;
|
||||
// << lambda[i] << ") is " << std::endl;
|
||||
// std::cout << tmp_gradient;
|
||||
gradient += tmp_gradient;
|
||||
}
|
||||
@@ -176,7 +176,7 @@ template<typename LagrangianFunction>
|
||||
const arma::mat& AugLagrangian<LagrangianFunction>::AugLagrangianFunction::
|
||||
GetInitialPoint()
|
||||
{
|
||||
return function_.GetInitialPoint();
|
||||
return function.GetInitialPoint();
|
||||
}
|
||||
|
||||
}; // namespace optimization
|
||||
|
||||
@@ -45,7 +45,7 @@ void AugLagrangianTestFunction::Gradient(const arma::mat& coordinates,
|
||||
gradient[1] = 4 * coordinates[0] + 6 * coordinates[1];
|
||||
}
|
||||
|
||||
double AugLagrangianTestFunction::EvaluateConstraint(int index,
|
||||
double AugLagrangianTestFunction::EvaluateConstraint(size_t index,
|
||||
const arma::mat& coordinates)
|
||||
{
|
||||
// We return 0 if the index is wrong (not 0).
|
||||
@@ -56,7 +56,7 @@ double AugLagrangianTestFunction::EvaluateConstraint(int index,
|
||||
return (coordinates[0] + coordinates[1] - 5);
|
||||
}
|
||||
|
||||
void AugLagrangianTestFunction::GradientConstraint(int index,
|
||||
void AugLagrangianTestFunction::GradientConstraint(size_t index,
|
||||
const arma::mat& coordinates, arma::mat& gradient)
|
||||
{
|
||||
// If the user passed an invalid index (not 0), we will return a zero
|
||||
@@ -107,7 +107,7 @@ void GockenbachFunction::Gradient(const arma::mat& coordinates,
|
||||
gradient[2] = 6 * (coordinates[2] + 3);
|
||||
}
|
||||
|
||||
double GockenbachFunction::EvaluateConstraint(int index,
|
||||
double GockenbachFunction::EvaluateConstraint(size_t index,
|
||||
const arma::mat& coordinates)
|
||||
{
|
||||
double constraint = 0;
|
||||
@@ -130,7 +130,7 @@ double GockenbachFunction::EvaluateConstraint(int index,
|
||||
return constraint;
|
||||
}
|
||||
|
||||
void GockenbachFunction::GradientConstraint(int index,
|
||||
void GockenbachFunction::GradientConstraint(size_t index,
|
||||
const arma::mat& coordinates,
|
||||
arma::mat& gradient)
|
||||
{
|
||||
@@ -203,13 +203,13 @@ void LovaszThetaSDP::Gradient(const arma::mat& coordinates,
|
||||
// std::cout << gradient;
|
||||
}
|
||||
|
||||
int LovaszThetaSDP::NumConstraints()
|
||||
size_t LovaszThetaSDP::NumConstraints() const
|
||||
{
|
||||
// Each edge is a constraint, and we have the constraint Tr(X) = 1.
|
||||
return edges_.n_cols + 1;
|
||||
}
|
||||
|
||||
double LovaszThetaSDP::EvaluateConstraint(int index,
|
||||
double LovaszThetaSDP::EvaluateConstraint(size_t index,
|
||||
const arma::mat& coordinates)
|
||||
{
|
||||
if (index == 0) // This is the constraint Tr(X) = 1.
|
||||
@@ -232,7 +232,7 @@ double LovaszThetaSDP::EvaluateConstraint(int index,
|
||||
return dot(coordinates.col(i), coordinates.col(j));
|
||||
}
|
||||
|
||||
void LovaszThetaSDP::GradientConstraint(int index,
|
||||
void LovaszThetaSDP::GradientConstraint(size_t index,
|
||||
const arma::mat& coordinates,
|
||||
arma::mat& gradient)
|
||||
{
|
||||
|
||||
@@ -29,10 +29,10 @@ class AugLagrangianTestFunction
|
||||
double Evaluate(const arma::mat& coordinates);
|
||||
void Gradient(const arma::mat& coordinates, arma::mat& gradient);
|
||||
|
||||
int NumConstraints() const { return 1; }
|
||||
size_t NumConstraints() const { return 1; }
|
||||
|
||||
double EvaluateConstraint(int index, const arma::mat& coordinates);
|
||||
void GradientConstraint(int index,
|
||||
double EvaluateConstraint(size_t index, const arma::mat& coordinates);
|
||||
void GradientConstraint(size_t index,
|
||||
const arma::mat& coordinates,
|
||||
arma::mat& gradient);
|
||||
|
||||
@@ -62,10 +62,10 @@ class GockenbachFunction
|
||||
double Evaluate(const arma::mat& coordinates);
|
||||
void Gradient(const arma::mat& coordinates, arma::mat& gradient);
|
||||
|
||||
int NumConstraints() const { return 2; };
|
||||
size_t NumConstraints() const { return 2; };
|
||||
|
||||
double EvaluateConstraint(int index, const arma::mat& coordinates);
|
||||
void GradientConstraint(int index,
|
||||
double EvaluateConstraint(size_t index, const arma::mat& coordinates);
|
||||
void GradientConstraint(size_t index,
|
||||
const arma::mat& coordinates,
|
||||
arma::mat& gradient);
|
||||
|
||||
@@ -115,10 +115,10 @@ class LovaszThetaSDP
|
||||
double Evaluate(const arma::mat& coordinates);
|
||||
void Gradient(const arma::mat& coordinates, arma::mat& gradient);
|
||||
|
||||
int NumConstraints();
|
||||
size_t NumConstraints() const;
|
||||
|
||||
double EvaluateConstraint(int index, const arma::mat& coordinates);
|
||||
void GradientConstraint(int index,
|
||||
double EvaluateConstraint(size_t index, const arma::mat& coordinates);
|
||||
void GradientConstraint(size_t index,
|
||||
const arma::mat& coordinates,
|
||||
arma::mat& gradient);
|
||||
|
||||
|
||||
@@ -61,6 +61,41 @@ class L_BFGS
|
||||
*/
|
||||
bool Optimize(const size_t maxIterations, arma::mat& iterate);
|
||||
|
||||
//! Get the memory size.
|
||||
size_t NumBasis() const { return numBasis; }
|
||||
//! Modify the memory size.
|
||||
size_t& NumBasis() { return numBasis; }
|
||||
|
||||
//! Get the Armijo condition constant.
|
||||
double ArmijoConstant() const { return armijoConstant; }
|
||||
//! Modify the Armijo condition constant.
|
||||
double& ArmijoConstant() { return armijoConstant; }
|
||||
|
||||
//! Get the Wolfe parameter.
|
||||
double Wolfe() const { return wolfe; }
|
||||
//! Modify the Wolfe parameter.
|
||||
double& Wolfe() { return wolfe; }
|
||||
|
||||
//! Get the minimum gradient norm.
|
||||
double MinGradientNorm() const { return minGradientNorm; }
|
||||
//! Modify the minimum gradient norm.
|
||||
double& MinGradientNorm() { return minGradientNorm; }
|
||||
|
||||
//! Get the maximum number of line search trials.
|
||||
size_t MaxLineSearchTrials() const { return maxLineSearchTrials; }
|
||||
//! Modify the maximum number of line search trials.
|
||||
size_t& MaxLineSearchTrials() { return maxLineSearchTrials; }
|
||||
|
||||
//! Return the minimum line search step size.
|
||||
double MinStep() const { return minStep; }
|
||||
//! Modify the minimum line search step size.
|
||||
double& MinStep() { return minStep; }
|
||||
|
||||
//! Return the maximum line search step size.
|
||||
double MaxStep() const { return maxStep; }
|
||||
//! Modify the maximum line search step size.
|
||||
double& MaxStep() { return maxStep; }
|
||||
|
||||
private:
|
||||
//! Internal copy of the function we are optimizing.
|
||||
FunctionType function;
|
||||
|
||||
Reference in New Issue
Block a user