Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.

Nonlinear Least Squares Problem with Linear (Bound) Constraints

Shokriafra__Mohammad
847 Views

Hi,

I noticed that the nonlinear solver uses central difference to calculate the Jacobi matrix, I was wondering is it possible to pass the analytical Jacobi elements to the code or not? I mean, in my problem, I want to use analytical Jacobi instead of numerical Jacobi. How can I do that?

Thanks
Javad

0 Kudos
3 Replies
mecej4
Honored Contributor III
847 Views

The nonlinear solver uses the reverse call interface. When RCI_REQUEST is equal to 2, you have to update the Jacobian matrix FJAC and call the solver again. The method used for calculating the Jacobian is entirely in your hands.  You can enter formulae for the elements of the Jacobian in line, or in a separate subroutine that you write and call; you can use finite difference approximations for some or all elements, or obtain the derivatives from a finite-field extension using complex variables.

You can, for example, take the furnished testcase ex_nlsqp_f.f and modify it so that FJAC is calculated analytically instead of calling DJACOBI. 

0 Kudos
Shokriafra__Mohammad
831 Views

Thanks for your reply,

I am using this solver in C++, I used the "ex_nlsqp_bc_c.c" example and tried to implement what you mentioned as the following.

I defined the function as

extern void jacobian(MKL_INT *, MKL_INT *, double *, double *);

----

I called the Jacobian function as

if (RCI_Request == 2)
{
jacobian(&m, &n, x, fjac);
}

----

and defined the Jacobian like the following:

void jacobian(MKL_INT * m, MKL_INT * n, double *x, double *fj)
{
fj[0] = 2 * x[0];
fj[1] = 2 * x[1];
fj[2] = 1;
fj[3] = -1;

return;
}

----

The function was defined as:

void extended_powell(MKL_INT * m, MKL_INT * n, double *x, double *f)
{
f[0] = pow(x[0],2) + pow(x[1],2) - 9;
f[1] = x[0] - x[1] + 3;

return;
}

----

When I call the jacobian, unfortunately it does not work, but if I call djacobi it works.

I tracked the iteration of jacobian, it does around 20 iterations and becomes close to the solution, but I do not know why it fails after several iteration while it is getting closer to the correct solution.

The correct solution for this problem is "zero" and "3", and in the following I printed the two last iterations of the solution when I am using jacobian.  

x1= 0.10554969128064530                        x2= 2.9935339376683032

x1= 0.10554965354656229                        x2= 2.9935339407828785

 

Can anyone help me?

 

 

0 Kudos
mecej4
Honored Contributor III
806 Views

I wrote a detailed reply, but the forum software marked it as spam and removed it. Here is a short replacement:

Your Jacobian vector arranges the elements of the 2-dimensional matrix in row-major order. The MKL routine may expect it to be in column-major (as in Fortran) order. Try printing out the values of fjac from the finite-difference calculation and compare those with the analytically calculated values for accuracy and matching order.

0 Kudos
Reply