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

Eps for Jacobian calculation

Shigeo_Kobayashi
Beginner
852 Views

The MKL reference manual does not say much about eps to be specified for Jacobian calculation.
The examples such as ex_nlsqp_c.c etc actually specify the same eps (an array of  6 elements) for dtrnlsp_init()  to djacobi().
Which element of eps[] is used for the  Jacobian calculation?
Is it eps[5]  (The trial step precision.)?  
 "If eps(6) = 0, then the trial step meets the required precision ( 1.0D-10)." is too hard for me to understand.
I'm sure djacobi/djacobix internally allocate the arrays(2 arrays?) holding function values each time they are called(and free them on exit).  I do not think this is a good idea...

 

0 Kudos
1 Solution
mecej4
Honored Contributor III
852 Views

Kobayashi-san.:

If the elements of x were all about 1.0 in magnitude, the eps_Jac could be added/subtracted directly. However, when different elements of x are of different magnitudes (and, even different units, such as kg, m and joules), scaling is necessary. Therefore, what is usually used is dx_i = eps_jac • |x_i|, with special handling for those elements of x that are close to zero. Secondly, as you may have observed, only one component of x is changed at a time, with the others held fixed.

In the example code of sjacobi_rci_c, eps_jacobi = 10-4, and each x component = 10. Thus, the perturbations used are 10-3, as you can see for yourself by printing out the x values before calling the extended_powell() function.

For normal use, the proper value to use for eps_jacobi with central differencing is proportional to the cube root of machine epsilon, which works out to about 10-5.

View solution in original post

0 Kudos
10 Replies
Alexander_K_Intel2
852 Views

Hi,

In example you mentioned i founded following line:

            if (djacobi (extendet_powell, &n, &m, fjac, x, eps) != TR_SUCCESS)

Based on MKL documentation eps is pointer on double precision variable so djacobi functionality used only first element of array eps - eps[0]

Thanks,

Alexander Kalinkin

0 Kudos
Shigeo_Kobayashi
Beginner
852 Views

Hi Alexander,
Thank you for your answer. Yes,it must be eps[0]!  I was little bit confused.
But, eps[0] (trust region criterion) is suitable for jacobi matrix calculation or not is another problem...  
Thank you again anyway.

0 Kudos
mecej4
Honored Contributor III
852 Views

The C example that you cited may cause needless confusion. The EPS argument to DJACOBI/JACOBI is a scalar, whereas the argument of the same name for the solver is a 6-element array. The two are not related. The Fortran version of the example keeps the two separate.

From an algorithmic viewpoint, if the objective function can be computed to full machine precision, the optimum value for EPS_JAC is proportional to SQRT(machine_epsilon) if forward differences are used, and cube_root(machine_epsilon) if central differences are used.

It would be good to revise the C example source code(s) accordingly.

0 Kudos
Gennady_F_Intel
Moderator
852 Views

thanks mecej, we need to slighly improve C test.

0 Kudos
Gennady_F_Intel
Moderator
852 Views

thanks mecej, we need to slighly improve C test.

0 Kudos
Shigeo_Kobayashi
Beginner
852 Views

Hi mecej,
Thank you for your comments on EPS_JAC which I really wanted to know.
The reference says that arrays f1 and f2 specified to sjacobi_solve() etc contain f(x+eps) and f(x-eps) respectively.
Is this true ?
If this is true then rci cycle "while (successful == 0)" in the example sjacobi_rci_c.c is useless.

0 Kudos
mecej4
Honored Contributor III
852 Views

Shigeo Kobayashi wrote:
If this is true then rci cycle "while (successful == 0)" in the example sjacobi_rci_c.c is useless.

The integer flag "successful" is maintained and updated in the user program. It is zero initially, and is set to 1 when rci_request becomes 0. The variable rci_request itself could have been used as the flag for loop termination.

MKL provides two different interfaces for computing the Jacobian using numerical approximation. In the first, the function pointer is passed as an argument to ?jacobi, and the MKL ?jacobi routine itself calls that function as many times as it wants. However, the interface to the function has to be fixed, to agree with whatever interface was declared when the ?jacobi routine itself was compiled.

The second interface uses RCI instead. You may call another function or do the calculations in-line to find f(x1) and f(x2). Therefore, there is no restriction as to the function interface. In fact, you could as a friend to look up the function value in a printed table and type that value in!

Thus, (at least) three calls to the ?jacobi_solve routine are needed:

   Call-1: The ?jacobi_solve routine returns to you with x1; you compute f(x1)

   Call-2: You send f(x1) to ?jacobi_solve, and it returns to you with  x2, You compute f(x2).

   Call-3: You send f(x1) and f(x2) to ?jacobi_solve, and it computes and stores the Jacobian using f(x1) and f(x2), and returns to you hinting "OK, I have the Jacobian computed and stored, proceed with whatever is the next step in your optimization algorithm."

[Note that you would have made a Call-0 earlier -- to ?jacobi_init, setting it up with the value of x. The MKL routines use that value of x to obtain the slightly perturbed vectors x1 and x2.]

0 Kudos
Shigeo_Kobayashi
Beginner
852 Views

Hi mecej4,
Thank you for detailed explanation.
My question might little be ambigous.,sorry.
As far as I tried the example,I watched  x, and I found the perturbation value is not the same as eps as the reference describes.
That is,when rci_request == 1 (rci_request == 2 the same),  the perturbation value(say dx) is not the same as eps.
sjacobi_solve() actually sets x+dx  (dx !=  eps)  for computing jacobi matrix.
Pls note eps = 0.0001 but dx = 0.001.  The behavior dx != eps is what I expected.
And what I want to know is how dx is obtained from eps,otherwise I can not estimate a better eps value for different optimization cases.

0 Kudos
mecej4
Honored Contributor III
853 Views

Kobayashi-san.:

If the elements of x were all about 1.0 in magnitude, the eps_Jac could be added/subtracted directly. However, when different elements of x are of different magnitudes (and, even different units, such as kg, m and joules), scaling is necessary. Therefore, what is usually used is dx_i = eps_jac • |x_i|, with special handling for those elements of x that are close to zero. Secondly, as you may have observed, only one component of x is changed at a time, with the others held fixed.

In the example code of sjacobi_rci_c, eps_jacobi = 10-4, and each x component = 10. Thus, the perturbations used are 10-3, as you can see for yourself by printing out the x values before calling the extended_powell() function.

For normal use, the proper value to use for eps_jacobi with central differencing is proportional to the cube root of machine epsilon, which works out to about 10-5.

0 Kudos
Shigeo_Kobayashi
Beginner
852 Views

Hi mecej4,
Thank you!  I'm clear now!! . What you commented is what I expected.
Perhaps if the jacobian is singular (or |f1|==|f2|) then perturbation value dx is increased and rci cycle keeps going ?
Thank you again for your valuable comments,and I hope this issue is added to the documentation.
("Kobayashi-san" ... Yes, Merry Christmas from Kanagawa Japan)

0 Kudos
Reply