- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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 constraintsvoid 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 constraintsvoid 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 valuesdouble 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 = 1271347576Unconstrained Solver:Initial Parameters: a = 3.69707 b = -1.66525Final Parameters: a = 5.33677 b = -2.29319Target Parameters: a = 5.33677 b = -2.29319Results: Iterations = 14 Stop Condition = 3Initial Residuals = 4740.99 Final Residuals = 5.13598e-006Constrained Solver:Initial Parameters: a = 3.69707 b = -1.66525Final Parameters: a = 3.73313 b = -2.03975Target Parameters: a = 5.33677 b = -2.29319Results: Iterations = 11 Stop Condition = 2Initial Residuals = 4740.99 Final Residuals = 3561.57Needless 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 KiazykEdit: my apologies, I cannot seem to get angle brackets to work at all, hopefully the code is still clear.*

Link Copied

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

Reinstating the missing angle brackets in Stephen's in-line source requires more guesswork than I am inclined to risk.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page