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):

//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[i] = 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; }

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.