Auto-suggest helps you quickly narrow down your search results by suggesting possible matches as you type.

Showing results for

- Intel Community
- Software
- Software Development SDKs and Libraries
- Intel® oneAPI Math Kernel Library
- Eps for Jacobian calculation

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

Shigeo_Kobayashi

Beginner

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

12-17-2013
09:34 PM

194 Views

Eps 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

= 0, then the trial step meets the required precision (≤ 1.0D-10)." is too hard for me to understand.*eps*(6)

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...

1 Solution

mecej4

Black Belt

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

12-20-2013
04:23 AM

194 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}.

Link Copied

10 Replies

Alexander_K_Intel2

Employee

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

12-18-2013
09:04 PM

194 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

Shigeo_Kobayashi

Beginner

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

12-19-2013
01:10 AM

194 Views

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.

mecej4

Black Belt

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

12-19-2013
05:41 AM

194 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.

Gennady_F_Intel

Moderator

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

12-19-2013
09:36 AM

194 Views

thanks mecej, we need to slighly improve C test.

Gennady_F_Intel

Moderator

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

12-19-2013
09:36 AM

194 Views

thanks mecej, we need to slighly improve C test.

Shigeo_Kobayashi

Beginner

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

12-19-2013
05:54 PM

194 Views

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.

mecej4

Black Belt

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

12-19-2013
07:39 PM

194 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.]

Shigeo_Kobayashi

Beginner

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

12-19-2013
10:44 PM

194 Views

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.

mecej4

Black Belt

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

12-20-2013
04:23 AM

195 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}.

Shigeo_Kobayashi

Beginner

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

12-21-2013
04:49 AM

194 Views

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)

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

For more complete information about compiler optimizations, see our Optimization Notice.