non-linear optimization routines

[Petros_Mamales]

Hi,

I am trying to use the trust-region methods w/in mkl 11 and have a couple of questions :

a) what is the layout of the jacobian matrix ?

I expect it to be fortran style matrix with m (mapped -f-  dimension) rows and n columns ( x dimension). Is this correct ?

b) concerning the jacobian calculations :

1) can I completely ignore the provided interface and provide instead with my own evaluation ? Is the jacobian matrix all that I am interested ?

(Unfortunately the examples provided are not very illuminating as far as the part of the API concerning the jacobian valuations. As a result, I still do not understand what it does, besides the obvious/expected ) 

2) if I provide with a NULL pointer, will the jacobian be calculated internally, using numerical differentiation ? -leaving aside for the moment questions of efficiency. I believe not, because a) it is not mentioned in the notes and b) because there is no place to insert the function pointer for the function valuation ( wouldn't that simplify a lot the interface ? The RCI use here seems quite typical, as far as the examples show, at least) - the question still is valid though, since, in this case, the code could request for function values, no ?

3) when calculating the jacobian, one of the inputs is the function values at the point of valuation. As far as I see, from the examples, there is no accompanying evaluation of the function values. Is this always  the case, and I, therefore, can use as fvec input the one held by the pointer passed in the solver routines ?

Finally, is there a more detailed write-up than the one provided in the reference manual ? ( the reference Conn00 comes with a rather heafty toll ;-) ).

Thank you very much in advance, for your help,

Petros

[mecej4]

Petros Mamales wrote:

I am trying to use the trust-region methods w/in mkl 11 and have a couple of questions :

a) what is the layout of the jacobian matrix ?

As the Fortran example shows more clearly, it is an m X n matrix.

I expect it to be fortran style matrix with m (mapped -f-  dimension) rows and n columns ( x dimension). Is this correct ?

Yes.

b) concerning the jacobian calculations :

1) can I completely ignore the provided interface and provide instead with my own evaluation ? Is the jacobian matrix all that I am interested ?

(Unfortunately the examples provided are not very illuminating as far as the part of the API concerning the jacobian valuations. As a result, I still do not understand what it does, besides the obvious/expected )

If you want to provide statements for calculating the Jacobian, you can replace the call to DJACOBI by these statements.

2) if I provide with a NULL pointer, will the jacobian be calculated internally, using numerical differentiation ? -leaving aside for the moment questions of efficiency. I believe not, because a) it is not mentioned in the notes and b) because there is no place to insert the function pointer for the function valuation ( wouldn't that simplify a lot the interface ? The RCI use here seems quite typical, as far as the examples show, at least) - the question still is valid though, since, in this case, the code could request for function values, no ?

Remember that you are using a reverse call interface, except when/if you call DJACOBI. Therefore, NULL pointers make no sense at all.

3) when calculating the jacobian, one of the inputs is the function values at the point of valuation. As far as I see, from the examples, there is no accompanying evaluation of the function values. Is this always  the case, and I, therefore, can use as fvec input the one held by the pointer passed in the solver routines ?

The first argument to DJACOBI is the external subprogram that is called to calculate the objective function. You are not passing the function value, but the address of the function (or function pointer, in C terminology). If DJACOBI needs to use the values in FVEC (as it must, for mathematical reasons) it manages the storage and access of FVEC behind the scenes and you need not be concerned as to the details. I find the Fortran version of the example to be quite clear (as would someone familiar with the hoary RCI).

[Petros_Mamales]

Hi Mecej4,

Thank you very much for the quick reply. Much appreciated!

I am not able to read FORTRAN code any more, which is the reason I never look in the examplesf ;-)

As it turns out, jacobix is more suitable to my needs, namely c++, where I can exploit the user_data of the input.

However, it is not clear from the documentation, the order of the arguments in the function pointer argument (the api that does not contain identifiers and the abscense of  const-qualifiers doesn't help ).

So, is the order :

m, n, fvec, x user_data ?

(the question refers mainly to the floaters).

Thank you again for your help,

Petros

[mecej4]

Petros Mamales wrote:

However, it is not clear from the documentation, the order of the arguments in the function pointer argument (the api that does not contain identifiers and the abscense of  const-qualifiers doesn't help ).

So, is the order :

m, n, fvec, x user_data ?

(the question refers mainly to the floaters).

Ah, I now see the point of doubt from your side, and the explanation is quite simple. The function extendet_powell() in the C example file ex_nlsqp_c.c is an optional function, and its argument list is entirely up to you to select. In fact, instead of defining and calling the function after a return with RCI_Request = 1, you can put statements in line to calculate the objective function. Similarly, you can call another function of your own creation or put statements in line to calculate the jacobian. If, however, you wish to make use of the MKL-provided finite-difference jacobian routine ?jacobi(), you must provide a function that matches the argument list specified for the argument fcn in the documentation for ?jacobi().

In summary, if you wish to write your own code for evaluating the function and the jacobian, you can do all that in-line, with no need for a function such as extendet_powell() and no calls to ?jacobi().

[Petros_Mamales]

Hi Mecej4,

I guess I asked a less than clever question as your kind response indicates ;-)

It seems that, other than conformity with the API, the only thing the jacobix should really use is some function pointer to a function of void * user_data. The rest of the arguments of the function pointer should be superfluous - other than the x and fvec (and the same goes with the arguments of jacobix itself, other than  fjac and x ) since all this could be found in user_data !

Thank you very much for your help,

Petros

PS: In reality, I want to use the mkl numerical jacobian when feasible and something else ( maybe analytical/quasi-analytical ) for other cases.

This can be achieved with the help of a couple of wrappers, that provide the functionality, but do not necessarily adhere to the function calls the api and the examples indicate.

[mecej4]

Here are the lines of code to replace matching lines in the MKL example ex_nlsqp_c.c.

[cpp]        if (RCI_Request == 1) {
            /* recalculate function value */
           fvec[0] = x[0] + 10.0 * x[1];
           fvec[1] = 2.2360679774998 * (x[2] - x[3]);
           fvec[2] = (x[1] - 2.0 * x[2]) * (x[1] - 2.0 * x[2]);
           fvec[3] = 3.1622776601684 * (x[0] - x[3]) * (x[0] -x[3]);
        }
        if (RCI_Request == 2) {
            /* compute jacobi matrix */
            fjac[0]=1.0; fjac[4]=10.0; fjac[8]=0; fjac[12]=0;
            fjac[1]=0; fjac[5]=0; fjac[9]=2.2360679774998;
            fjac[13]=-fjac[9];
            fjac[2]=0; fjac[6]=2*(x[1]-2.0*x[2]); fjac[10]=-fjac[6];
            fjac[14]=0;
            fjac[3]=2*3.1622776601684*((x[0] - x[3])); fjac[7]=fjac[11]=0;
            fjac[15]=-fjac[3];
        }
[/cpp]

The entire function extendet_powell () can be removed. As you can see, when the function and jacobian are calculated explicitly as above, neither that function nor djacobi, jacobix, etc., need to be called.

There are other nonlinear solvers that provide a more flexible interface. The user is asked to provide explicit calculation of as many components of the jacobian as can be programmed with reasonable effort. The solver calculates the remaining components by finite-difference approximation.

[Petros_Mamales]

Hi mecej4,

Yes, I understand. As I mentioned earlier I wanted for myself some of the flexibility and also the possibility to use jacobix, avoiding so to rewrite one by myself ( since there is already fine, fast and fully debugged code that does this for me!).

Thank you very much,

Petros

[Petros_Mamales]

I need one more clarification ;-))

The function pointer that is the input of the jacobix will use the same x and fvec that the rest of the routine will use , correct ?

This I conclude by the assumption that there are no further allocations other the ones the user provides, but still there could be swaps of the

pointers internally. Is the assumption correct ?

The reason I am asking is that, in my environment, I cannot afford raw pointers dangling around. They have to belong to a vector (boost unbounded array, to be exact). My only hope is that I can bundle the x_ and fvec_ structures with the user_data.

For this to work though, the pointers of the function pointer and the function itself have to coincide (in pairs) when used.

Is this correct? This is getting a bit merky..apologies,

TIA,

Petros

[mecej4]

If you consider that the routines such as the solver itself and ?jacobix are written in Fortran, with additional facilities to enable them to be called from C, and you note that Fortran is a higher level language with no such thing as pointers in older versions than F90, you may see that the questions do not make sense. In a Fortran program, there are no pointers as such. You may think of an array name as something similar to a C pointer, but in Fortran you cannot detach that name and reattach it (i.e., pointer-associate with or point at) to another variable. Dangling pointers are next to impossible to create in Fortran 77.

In normal usage from Fortran code, the function pointer would be a constant-valued pointer whose value is filled in at link time. In C, you could have a prototype statement such as "extern void myfcn(int *m, int *n, double *x, double *f,void *udata);" in the caller of ?jacobix, and simply use "myfcn" as the first argument to ?jacobix.

Perhaps you are trying to protect yourself from hazards that do not exist (at least in code written in standard-conforming Fortran 77).

[Petros_Mamales]

Hi mecej,

1)To begin with, I do not know that the library is internally written in FORTRAN, at least in its entirety.

2)You can swap contents (temporarily, for convenience say, w/out swaping pointers (that don't exist) ).

3) dangling pointers cannot be afforded on my side, I am sure mkl has taken care of it.

The real question had to do with wether the function pointer will use the very argument point or a displaced one. his was my worry and I think it is pretty valid ( meaning the buffer).

At any rate, I think I will do something a bit different.

Thank you for your help and have a nice weekend.

Petros