Non-linear optimization fails when bounds constraints are applied

[stephen_kiazyk]

Greetings,

I am currently developing an application which requires the use of a non-linear solver. In our application, it would be beneficial to be able to place loose but definite bounds constraints on the possible solution set, and hence we were intending to use the bound-constrained version of the non-linear solver. However, after running some test cases, we have found that the bound-constrained solver will often fail to reach an appropriate solution, even in cases where the unconstrained solver would have no problem, and both the initial values and the desired solution were well within the bound-constraints placed on the solver. I am not an expert on non-linear solvers nor on MKL, so I am not sure if these undesired results are due to my misunderstanding of the problem formulation, or misapplication of the MKL library to this problem, or some other factor; if possible, could somebody perhaps indicate what is being done incorrectly in this situation.

Here is an example test we have run (it is an intentionally trivial program, this is in no way reflective of our intended usage of the mkl solver in the application proper. I have provided a custom random function and pre-defined seed only so that it could immediately give a failure case for the purpose of illustration):
[bash]//general includes
#include <C:/Program Files/Intel/MKL/10.2.4.032/include/mkl.h>
#include
#include
#include
using namespace std;

#define MKL_DIR "C:/Progra~1/Intel/MKL/10.2.4.032"

//linking to the mkl libraries
#ifdef _WIN64

#pragma comment( lib, MKL_DIR "/em64t/lib/mkl_intel_thread.lib" )
#pragma comment( lib, MKL_DIR "/em64t/lib/mkl_core.lib" )
#pragma comment( lib, MKL_DIR "/em64t/lib/libguide.lib" )
#pragma comment( lib, MKL_DIR "/em64t/lib/mkl_intel_lp64.lib" )

#else

#pragma comment( lib, MKL_DIR "/ia32/lib/mkl_intel_thread.lib" )
#pragma comment( lib, MKL_DIR "/ia32/lib/mkl_intel_c.lib" )
#pragma comment( lib, MKL_DIR "/ia32/lib/mkl_core.lib" )
#pragma comment( lib, MKL_DIR "/ia32/lib/libguide.lib" )

#endif

//custom random function provided to allow illustration of this case
//regardless of target system configuration
//from POSIX.1-2001
static unsigned int nextSeed = 1;

void CustomSRand( unsigned int seed )
{
nextSeed = seed;
}

int CustomRand()
{
nextSeed = nextSeed*1103515245 + 12345;
return int((nextSeed/65536) % RAND_MAX);
}

//configuration parameters to the non-linear solvers (see mkl docs for details)
struct MKLNonLinearConfig
{
double epsilon[6];

MKL_INT maxIterations;
MKL_INT maxTrialIterations;

double initialStepBound;
double jacobianPrecision;
};

//output data from the non-linear solvers
struct MKLNonLinearResults
{
MKL_INT iterationsUsed;
MKL_INT stopCriterion;
double initialResiduals;
double finalResiduals;
};

//data used inside the objective function
struct NonLinearData
{
double targetAValue;
double targetBValue;
double lowerXBounds;
double upperXBounds;
};

//a reasonably simple non-linear function: 'a' and 'b' are the parameters
//I want to solve for, x is the position to evaluate the function at
double NonLinearFunction( double a, double b, double x )
{
return a*std::exp( -b*x );
}

//this objective will sample a number of points between the bounds specified in the user data, and give the residuals for
//the solver's current parameters compared to the target parameters
void MKLObjective( MKL_INT* pNumResiduals, MKL_INT* pNumParams, double* params, double* residuals, void* userData )
{
const NonLinearData* data = static_cast( userData );

const double stepIncrement = (data->upperXBounds - data->lowerXBounds) / ((*pNumResiduals)-1);

assert( *pNumParams == 2 );

double currentX = data->lowerXBounds;

for ( size_t i = 0; i < size_t(*pNumResiduals); ++i )
{
residuals = NonLinearFunction( data->targetAValue, data->targetBValue, currentX ) - NonLinearFunction( params[0], params[1], currentX );
currentX += stepIncrement;
}
}

//solver without bounds constraints
void MKLNonLinearSolve( MKL_INT numParams, MKL_INT numResiduals, MKLNonLinearConfig* config, NonLinearData* userData, double* inoutParams, MKLNonLinearResults* results )
{
double* residuals = new double[numResiduals];
double* jacobian = new double[numParams*numResiduals];

_TRNSP_HANDLE_t mklOperationHandle = 0;

MKL_INT initResult = dtrnlsp_init( &mklOperationHandle, &numParams, &numResiduals, inoutParams, config->epsilon, &config->maxIterations, &config->maxTrialIterations, &config->initialStepBound );
assert( initResult == TR_SUCCESS );

MKL_INT rciRequest = 0;

do
{
MKL_INT solveResult = dtrnlsp_solve( &mklOperationHandle, residuals, jacobian, &rciRequest );
assert( solveResult == TR_SUCCESS );

switch( rciRequest )
{
case 1:
MKLObjective( &numResiduals, &numParams, inoutParams, residuals, static_cast(userData) );
break;
case 2:
MKL_INT jacobianResult = djacobix( &MKLObjective, &numParams, &numResiduals, jacobian, inoutParams, &config->jacobianPrecision, static_cast(userData) );
assert( jacobianResult == TR_SUCCESS );
break;
}
}
while( rciRequest >= 0 );

MKL_INT getResult = dtrnlsp_get( &mklOperationHandle, &results->iterationsUsed, &results->stopCriterion, &results->initialResiduals, &results->finalResiduals );
assert( getResult == TR_SUCCESS );

MKL_INT deleteResult = dtrnlsp_delete( &mklOperationHandle );
assert( deleteResult == TR_SUCCESS );

delete [] residuals;
delete [] jacobian;
}

//solver with bounds constraints
void MKLNonLinearSolveWithConstraints( MKL_INT numParams, MKL_INT numResiduals, MKLNonLinearConfig* config, NonLinearData* userData, double* lowerBounds, double* upperBounds, double* inoutParams, MKLNonLinearResults* results )
{
double* residuals = new double[numResiduals];
double* jacobian = new double[numParams*numResiduals];

_TRNSPBC_HANDLE_t mklOperationHandle = 0;

MKL_INT initResult = dtrnlspbc_init( &mklOperationHandle, &numParams, &numResiduals, inoutParams, lowerBounds, upperBounds, config->epsilon, &config->maxIterations, &config->maxTrialIterations, &config->initialStepBound );
assert( initResult == TR_SUCCESS );

MKL_INT rciRequest = 0;

do
{
MKL_INT solveResult = dtrnlspbc_solve( &mklOperationHandle, residuals, jacobian, &rciRequest );

assert( solveResult == TR_SUCCESS );

switch( rciRequest )
{
case 1:
MKLObjective( &numResiduals, &numParams, inoutParams, residuals, static_cast(userData) );
break;
case 2:
MKL_INT jacobianResult = djacobix( &MKLObjective, &numParams, &numResiduals, jacobian, inoutParams, &config->jacobianPrecision, static_cast(userData) );
assert( jacobianResult == TR_SUCCESS );
break;
}
}
while( rciRequest >= 0 );

MKL_INT getResult = dtrnlspbc_get( &mklOperationHandle, &results->iterationsUsed, &results->stopCriterion, &results->initialResiduals, &results->finalResiduals );
assert( getResult == TR_SUCCESS );

MKL_INT deleteResult = dtrnlspbc_delete( &mklOperationHandle );
assert( deleteResult == TR_SUCCESS );

delete [] residuals;
delete [] jacobian;
}

//a utility function to create random double values
double randomDouble( double lowerBound, double upperBound )
{
return lowerBound + ((upperBound - lowerBound) * (double(CustomRand()) / double(RAND_MAX)));
}

///Entry Point///
int main( int argc, char** argv )
{
//this seed will give a failure condition on the first iteration
unsigned int randomSeed = 1271347576;
CustomSRand( randomSeed );

//specify the configuration parameters
MKLNonLinearConfig config;

config.maxIterations = 1000;
config.maxTrialIterations = 1000;
config.epsilon[0] = 0.00001;
config.epsilon[1] = 0.00001;
config.epsilon[2] = 0.00001;
config.epsilon[3] = 0.00001;
config.epsilon[4] = 0.00001;
config.epsilon[5] = 0.00001;
config.initialStepBound = 100.0;
config.jacobianPrecision = 0.00001;

//specify the data to be used in the objective function
NonLinearData userData;

userData.targetAValue = randomDouble( -10.0, 10.0 );
userData.targetBValue = randomDouble( -3.0, 3.0 );
userData.lowerXBounds = -3.0;
userData.upperXBounds = 3.0;

//set up the starting paramters based on the target data
//(I keep a copy of the starting data for comparison at the end)
double startingParameters[2];

double unconstrainedSolverResults[2];
double constrainedSolverResults[2];

startingParameters[0] = userData.targetAValue * randomDouble( 0.5, 1.5 );
startingParameters[1] = userData.targetBValue * randomDouble( 0.5, 1.5 );

//place loose constraints on the parameters
double lowerBounds[2] = { -100.0, -100.0 };
double upperBounds[2] = { 100.0, 100.0 };

MKLNonLinearResults unconstrainedResults;
MKLNonLinearResults constrainedResults;

std::memcpy( unconstrainedSolverResults, startingParameters, sizeof(double)*2 );
MKLNonLinearSolve( 2, 10, &config, &userData, unconstrainedSolverResults, &unconstrainedResults );

std::memcpy( constrainedSolverResults, startingParameters, sizeof(double)*2 );
MKLNonLinearSolveWithConstraints( 2, 10, &config, &userData, lowerBounds, upperBounds, constrainedSolverResults, &constrainedResults );

std::cout << "Seed = " << randomSeed << std::endl;
std::cout << std::endl;
std::cout << "Unconstrained Solver:" << std::endl;
std::cout << "Initial Parameters: a = " << startingParameters[0] << "\tb = " << startingParameters[1] << std::endl;
std::cout << "Final Parameters: a = " << unconstrainedSolverResults[0] << "\tb = " << unconstrainedSolverResults[1] << std::endl;
std::cout << "Target Parameters: a = " << userData.targetAValue << "\tb = " << userData.targetBValue << std::endl;
std::cout << "Results: Iterations = " << unconstrainedResults.iterationsUsed << " Stop Condition = " << unconstrainedResults.stopCriterion << std::endl;
std::cout << "Initial Residuals = " << unconstrainedResults.initialResiduals << " Final Residuals = " << unconstrainedResults.finalResiduals << std::endl;
std::cout << std::endl;
std::cout << "Constrained Solver:" << std::endl;
std::cout << "Initial Parameters: a = " << startingParameters[0] << "\tb = " << startingParameters[1] << std::endl;
std::cout << "Final Parameters: a = " << constrainedSolverResults[0] << "\tb = " << constrainedSolverResults[1] << std::endl;
std::cout << "Target Parameters: a = " << userData.targetAValue << "\tb = " << userData.targetBValue << std::endl;
std::cout << "Results: Iterations = " << constrainedResults.iterationsUsed << " Stop Condition = " << constrainedResults.stopCriterion << std::endl;
std::cout << "Initial Residuals = " << constrainedResults.initialResiduals << " Final Residuals = " << constrainedResults.finalResiduals << std::endl;

return 0;
}[/bash]
When run, it gives the following results:


Seed = 1271347576

Unconstrained Solver:
Initial Parameters: a = 3.69707 b = -1.66525
Final Parameters: a = 5.33677 b = -2.29319
Target Parameters: a = 5.33677 b = -2.29319
Results: Iterations = 14 Stop Condition = 3
Initial Residuals = 4740.99 Final Residuals = 5.13598e-006

Constrained Solver:
Initial Parameters: a = 3.69707 b = -1.66525
Final Parameters: a = 3.73313 b = -2.03975
Target Parameters: a = 5.33677 b = -2.29319
Results: Iterations = 11 Stop Condition = 2
Initial Residuals = 4740.99 Final Residuals = 3561.57


Needless to say, I'm rather confused about why the addition of bounds constraints changed the result of the solver, and prevented it from finding an appropriate solution, when otherwise all other paramters to the solver were the same. I have also had similar problems when I trying to tightly constrain some of the parameters by setting the upper and lower bounds to be equal, in which case the solver would frequently outright ignore not only those tightly constrained bounds, but all bounds constraints provided and return with an absolutely insane solution. Is there something I am doing incorrectly, or some aspect of the solver that I am not correctly taking into account? And is it incorrect to specify the an equal value for both the upper and lower bounds constraints on a given parameter, or must there always besome non-empty range between them for the solver to work correctly?

Thank you for your time
- Stephen Kiazyk

Edit: my apologies, I cannot seem to get angle brackets to work at all, hopefully the code is still clear.

[Kate_D_]

I ran into problems with dtrnlspbc_solve ignoring bounds, and looking on the forum came across this topic. I have run Stephen's code and can reproduce the odd behaviour with MKL version 11.0. Can anyone from Intel address Stephen's (and now my identical) questions?

[mecej4]

Please attach a file with the problematic source code -- source that can be compiled with no compilation errors.
Reinstating the missing angle brackets in Stephen's in-line source requires more guesswork than I am inclined to risk.

[Kate_D_]

Here is a text file containing Stephen's source code. (.h and .cpp files in one text file)

[mecej4]

Thanks for posting the code.

In general, a nonlinear unconstrained optimization routine may fail to converge within the specified number of iterations for certain initial values of the unknown variables. Likewise, a bound-constrained optimization routine may fail to converge within the specified number of iterations for the same reason and, in addition, if the bounds are too tight.

Many successful optimization algorithms do not guarantee that every trial point will satisfy user-specified bounds, whereas other algorithms exist that maintain feasibility (i.e., constraints are satisfied) at all times. In either case, if the user specifies unreasonable bounds -- in particular, if the bounds are tight and/or do not contain the (or one of the) solutions, the algorithm may fail to converge.

For the specific example in the source code, I found that the bound-constrained optimizer yields the same solution as the unconstrained optimizer if the initial point is given as (5.2, 2.2), the bounds are specified as -100 and +100 for both unknowns, and the "target values" are a = 5.33677 b = -2.29319.

The observed failure(s) in the initial post may be caused partly because of trial values being too far from the solution, incorrect bounds, or in using randomly generated numbers for generating the data and other random numbers as initial parameter values.

The results:

--------------- Unconstrained Solver --------------- Initial Parameters: a = 5.2 b = 2.2 Target Parameters: a = 5.33677 b = 2.29319 Final Parameters: a = 5.33677 b = 2.29319 Stop Condition = 3 Final Residuals = 2.3441e-008 --------------- Constrained Solver ----------------- Initial Parameters: a = 5.2 b = 2.2 Target Parameters: a = 5.33677 b = 2.29319 Final Parameters: a = 5.33677 b = 2.29319 Stop Condition = 3 Final Residuals = 2.34438e-008

[Kate_D_]

Thank you very much mecej4, I too get that answer when I tweaked the userData.lowerXBounds/.upperXBounds, rather than the upper and lower bounds given to the solver. I perhaps don't understand what is happening in the MKLObjectiveFunciton, I'll look at it again. Thanks again for replying to my questions.