Questions about document of Optimization Solver Routines

[simbalee]

I am trying to use Optimization Solver Routines ( OSR ) of MKL 10.0.3.021.
No of the reference manual is: 630813-027US.
On page 2727, there are descriptions about "fvec" parameter of dtrnlsp_solve. It is said that:
fvec( i ) = yi - fi(x)
I understand that here "yi" is the i-th known values to be approximated through optimized values of "x",
and "fi(x)" is the i-th values that are computed with "x" of the current stage and a given calculate-routine.
So the fvec(i) here is a kind of residual error between approximate value "fi(x)" and expected value "yi" for me.
Am I correct?

However, on page 2726, it is also said that:
-------------------------------------------
RCI_Request= 1 that indicates the requirement to recalculate the function at vector X and put the result into fvec.
-------------------------------------------
It seems to me that instead of yi - fi(x), it is the fi(x) that should be put into fvec by following this instruction.

Now you may understand how confused I am.

The example given on page 2737 to page 2743 is consistent with the descriptions on page 2726.
Unfortunately the example itself is not quite self-explaining.
Say, the dimension of function value ( m ) is the same with number of function variables ( n ), which is atypical.
And the fucntion "extendet_powell" is too simple. The coefficents for calculating f(x) are all hard-coded in the "extendet_powell".
What if I need an array of un-redictable values to calculate f(x) ?
How do I pass pointer of the array to the calculating routine?
Is it impossible to do so?
Or do I need to attach the array behind f(x) and use a fixed shift in memory to access the array?

Your help will be deeply appreciated.

[ArturGuzik]

Quoting - simbalee

I am trying to use Optimization Solver Routines ( OSR ) of MKL 10.0.3.021.
No of the reference manual is: 630813-027US.
On page 2727, there are descriptions about "fvec" parameter of dtrnlsp_solve. It is said that:
fvec( i ) = yi - fi(x)
I understand that here "yi" is the i-th known values to be approximated through optimized values of "x",
and "fi(x)" is the i-th values that are computed with "x" of the current stage and a given calculate-routine.
So the fvec(i) here is a kind of residual error between approximate value "fi(x)" and expected value "yi" for me.
Am I correct?

However, on page 2726, it is also said that:
-------------------------------------------
RCI_Request= 1 that indicates the requirement to recalculate the function at vector X and put the result into fvec.
-------------------------------------------
It seems to me that instead of yi - fi(x), it is the fi(x) that should be put into fvec by following this instruction.

Now you may understand how confused I am.

The example given on page 2737 to page 2743 is consistent with the descriptions on page 2726.
Unfortunately the example itself is not quite self-explaining.
Say, the dimension of function value ( m ) is the same with number of function variables ( n ), which is atypical.
And the fucntion "extendet_powell" is too simple. The coefficents for calculating f(x) are all hard-coded in the "extendet_powell".
What if I need an array of un-redictable values to calculate f(x) ?
How do I pass pointer of the array to the calculating routine?
Is it impossible to do so?
Or do I need to attach the array behind f(x) and use a fixed shift in memory to access the array?

Your help will be deeply appreciated.

Hi,

(1) the description is OK, everything depends on how you define your problem. What you do is minimization of the residual, so you minimize yi-fi(x). The doc reads:

RCI_Request= 1 that indicates the requirement to recalculate the function at vector X and put the result into fvec.

but just below there is a description of input argument as:

fvec DOUBLE PRECISION. Array of size m. Contains the function values at X, where fvec(i) = (yi fi (x)).

So it looks consistent. Simply f(x) doesn't mean fvec, as well as function doesn't always mean f(x). Here means residual.

(2) I have to say I don't know what you mean. You just define your function as external. define input and output etc. (containing fvec). and that's all. If the problem is somewhere else please elaborate what you intend to do.


A.

[simbalee]

Quoting - ArturGuzik

Hi,

(1) the description is OK, everything depends how you define your problem. What you do is minimizing the residual, so you minimize yi-fi(x). If you take a look at docs it reads:

RCI_Request= 1 that indicates the requirement to recalculate the function at vector X and put the result into fvec.

but below there is description of input as:

fvec DOUBLE PRECISION. Array of size m. Contains the function values at X, where fvec(i) = (yi fi (x)).

So it looks consistent. Simply f(x) doesn't mean fvec

(2) I have to say I don't know what you mean. You just define your function as external. define input and output etc. (containing fvec). and that's all. If the problem is somewhere else please elaborate what you intend to do.


A.

Thank you Artur.
(1) I have come to understand that the input "fvec" is yi - fi(x) and the meaning of f(x) is changed from an approximation function itself to its resual error within 3 pages. I trust you that your clarification is an official comment on this issue.

(2) The task I am dealing with is a calibration work of camera. That is to say, I have two kinds data,
A. 3D positions (xw, yw, zw ) of some hundreds points in 3d space
B. 2D posititons ( XI, yI ) of projection on a image of those points in 3D space.

What I am going to do is to estabished a way ( f_cam ) to project points in 3D space to 2D plane,

f_cam = ( xw, yw, zw ) -> ( XI, yI )

To use dtrnlsp, I need to define a external function to calculate fvec with a fixed prototype as you and the manual described. The prototype is defined as this way:
fcn( int *m, int *n, double *x, double *f )
here:
m should be totale number of all points times 2 for both XI and YI.
n is number of camera parameters.
x is the array of camera parameters.
f is output, which should be y - f(x). And here y is ( XI, YI ), fx is calculated projected position of ( xw, yw, zw ).

You see, there is no place to input either ( XI, YI ) or ( xw, yw, zw ) in the interface of fcn.
That is the question I asked.

Thank you for your patient.

[Nikita_S_Intel]

Quoting - simbalee

Thank you Artur.
(1) I have come to understand that the input "fvec" is yi - fi(x) and the meaning of f(x) is changed from an approximation function itself to its resual error within 3 pages. I trust you that your clarification is an official comment on this issue.

(2) The task I am dealing with is a calibration work of camera. That is to say, I have two kinds data,
A. 3D positions (xw, yw, zw ) of some hundreds points in 3d space
B. 2D posititons ( XI, yI ) of projection on a image of those points in 3D space.

What I am going to do is to estabished a way ( f_cam ) to project points in 3D space to 2D plane,

f_cam = ( xw, yw, zw ) -> ( XI, yI )

To use dtrnlsp, I need to define a external function to calculate fvec with a fixed prototype as you and the manual described. The prototype is defined as this way:
fcn( int *m, int *n, double *x, double *f )
here:
m should be totale number of all points times 2 for both XI and YI.
n is number of camera parameters.
x is the array of camera parameters.
f is output, which should be y - f(x). And here y is ( XI, YI ), fx is calculated projected position of ( xw, yw, zw ).

You see, there is no place to input either ( XI, YI ) or ( xw, yw, zw ) in the interface of fcn.
That is the question I asked.

Thank you for your patient.

Hi Simbalee,

The fixed interface of fcn function is need only for djacobi routine. OSR are RCI solvers, it means you can use your own functions for fvec and Jacobi matrix calculation.

But you plan to use the MKL routines for Jacobi matrix calculation, you have two options:

Use the RCI solvers for Jacobi matrix calculation;

Use dajcobix routine. This routine was added into MKL 10.2 release. This routine presents an alternative interface for the djacobi function that enables the user to pass additional data into the user's objective function fcn.

res = djacobix(fcn, n, m, fjac, x, jac_eps, user_data), where user_data Pointer to void (for FORTRAN, integer user_data(*)). Input parameter. Contains additional user's data, if any. Otherwise, a dummy argument.

[ArturGuzik]

Quoting - simbalee

Thank you Artur.
(1) I have come to understand that the input "fvec" is yi - fi(x) and the meaning of f(x) is changed from an approximation function itself to its resual error within 3 pages. I trust you that your clarification is an official comment on this issue.

(2) The task I am dealing with is a calibration work of camera. That is to say, I have two kinds data,
A. 3D positions (xw, yw, zw ) of some hundreds points in 3d space
B. 2D posititons ( XI, yI ) of projection on a image of those points in 3D space.

What I am going to do is to estabished a way ( f_cam ) to project points in 3D space to 2D plane,

f_cam = ( xw, yw, zw ) -> ( XI, yI )

To use dtrnlsp, I need to define a external function to calculate fvec with a fixed prototype as you and the manual described. The prototype is defined as this way:
fcn( int *m, int *n, double *x, double *f )
here:
m should be totale number of all points times 2 for both XI and YI.
n is number of camera parameters.
x is the array of camera parameters.
f is output, which should be y - f(x). And here y is ( XI, YI ), fx is calculated projected position of ( xw, yw, zw ).

You see, there is no place to input either ( XI, YI ) or ( xw, yw, zw ) in the interface of fcn.
That is the question I asked.

Thank you for your patient.

(1) my comment is by no means official. I'm not working for Intel MKL. It's just ...my comment.

(2) I might be wrong, but fixed prototype (interface in Fortran) is required only if you want/need to use djacobi function which approximates the Jacobian matrix using central difference method. If you know your Jacobian (can get it analytically) then you don't need it and interface can be any (as you want/need).

dtrnlspbc_SOLVE requires only array (Jacobian) and vector (residuals). How you going to get them it's on your side.

If you need djacobi, then you need to obtain reference to your (XI, YI ) or (xw, yw, zw) data inside the function (as interface/prototype is indeed fixed). In Fortran it would usually be achieved by declaring module, and then USE statement inside, so data is accessed inside function.

A.

EDIT: I was writing at the same time as Nikita (didn't see his entry). His comment would be official.

[simbalee]

Quoting - Nikita Shustrov (Intel)

Hi Simbalee,

The fixed interface of fcn function is need only for djacobi routine. OSR are RCI solvers, it means you can use your own functions for fvec and Jacobi matrix calculation.

But you plan to use the MKL routines for Jacobi matrix calculation, you have two options:

Use the RCI solvers for Jacobi matrix calculation;

Use dajcobix routine. This routine was added into MKL 10.2 release. This routine presents an alternative interface for the djacobi function that enables the user to pass additional data into the user's objective function fcn.

res = djacobix(fcn, n, m, fjac, x, jac_eps, user_data), where user_data Pointer to void (for FORTRAN, integer user_data(*)). Input parameter. Contains additional user's data, if any. Otherwise, a dummy argument.

Thank you Nikita.
I am trying to pass user_data to djacobi by attaching user_data to the end fvec.
That is to say, allocating a big memory, leaving space for fvec from begining of the memory and saving user_data behind the mem of fvec. Shift of pointer to the user_data can be calculated from input value m.
It seems to be possible for me. I wil try MKL10.2 if it failed.

----- I found that MKL10.2 has not been released yet. @_@

[simbalee]

Quoting - ArturGuzik

(1) my comment is by no means official. I'm not working for Intel MKL. It's just ...my comment.

(2) I might be wrong, but fixed prototype (interface in Fortran) is required only if you want/need to use djacobi function which approximates the Jacobian matrix using central difference method. If you know your Jacobian (can get it analytically) then you don't need it and interface can be any (as you want/need).

dtrnlspbc_SOLVE requires only array (Jacobian) and vector (residuals). How you going to get them it's on your side.

If you need djacobi, then you need to obtain reference to your (XI, YI ) or (xw, yw, zw) data inside the function (as interface/prototype is indeed fixed). In Fortran it would usually be achieved by declaring module, and then USE statement inside, so data is accessed inside function.

A.

EDIT: I was writing at the same time as Nikita (didn't see his entry). His comment would be official.

Well, your answer is quite professional, so I assumed that you are Intel guy. Sorry for what I said.
Your comments are pretty trustworthy for me. Thank you very much.

[Gennady_F_Intel]

Quoting - simbalee

Thank you Nikita.
I am trying to pass user_data to djacobi by attaching user_data to the end fvec.
That is to say, allocating a big memory, leaving space for fvec from begining of the memory and saving user_data behind the mem of fvec. Shift of pointer to the user_data can be calculated from input value m.
It seems to be possible for me. I wil try MKL10.2 if it failed.

----- I found that MKL10.2 has not been released yet. @_@

Simbalee,
We are planning to release MKL 10.2 in the end of the next month.
Please stay tuned and watch the announcement on this forum page.
--Gennady

[ArturGuzik]

So the question is whether you can wait until the end of June. the idea about attaching data to the end of fvec doesn't sound good to me, as:

(a) it relies on assumption that MKL internally uses only m elements, as specified (what if some internal code uses size() instead of m?

(b) if uses size() you'll need to add zeros to Jacobian matrix (to modify your LS problem).


So you have an option to make the data global (I know that C-people hate that) or get your Jacobian analytically. The latter is by all means the best solution. Anyway, I don't quite understand why you need that in the first place. If you make projection from 3D to plane you need transformation matrix, only. On the other hand if this a curved element/surface then issue is more intricate.

A.

[Nikita_S_Intel]

No, we don't use the size() in OSR implementation. But Artur is right, it's not a good way, so if you decide to use this solution, you should be very careful. But any way, I don't recommend to use this approach. You can use the RCI version of djacobi routine. In this case you can use any implenetaion of fcn() for Jacobi matrix calculation.

I agree with Artur, in general, if you get your Jacobi matrix analytically, it will be faster than central difference algorithm.

--Nikita

[simbalee]

Quoting - ArturGuzik

So the question is whether you can wait until the end of June. the idea about attaching data to the end of fvec doesn't sound good to me, as:

(a) it relies on assumption that MKL internally uses only m elements, as specified (what if some internal code uses size() instead of m?

(b) if uses size() you'll need to add zeros to Jacobian matrix (to modify your LS problem).


So you have an option to make the data global (I know that C-people hate that) or get your Jacobian analytically. The latter is by all means the best solution. Anyway, I don't quite understand why you need that in the first place. If you make projection from 3D to plane you need transformation matrix, only. On the other hand if this a curved element/surface then issue is more intricate.

A.

The projection function I am dealing with is highly nonlinear and the function itself is non-analytical. Generally there is no analytical way to get Jacobian matrix. That is why I choose MKL to do this work.
Some C-people instinctively avoid to use global data. In fact I never thought about it untill Artur pointed out. Before trying global data, I may want to use RCI version of djacobi in case that my way failed. I will report result of this work.
I assume that MKL should not internally access over m data in fvec. If memory of m data is not enough for computation, buffer should be allocated inside the djacobi routines, instead of relying on unpredictable size of input vector. But, it would be better to check if the assumption is true.

[simbalee]

Quoting - Nikita Shustrov (Intel)

No, we don't use the size() in OSR implementation. But Artur is right, it's not a good way, so if you decide to use this solution, you should be very careful. But any way, I don't recommend to use this approach. You can use the RCI version of djacobi routine. In this case you can use any implenetaion of fcn() for Jacobi matrix calculation.

I agree with Artur, in general, if you get your Jacobi matrix analytically, it will be faster than central difference algorithm.

--Nikita

Hi Nikita.
I got a new question about RCI version of djacobi routine.
According to document of djacobi_solve, parameter f1 contains updated function values at X + eps, and f2 for X - eps. I understand that here eps means a small enough value such as 1e-10, and by updating I should recalcaulate function values with all variables pluse or minus eps. Right?
However, in the example 14-6 "djacobi_solve usage in C", both f1 and f2 are updated with exactly the same X. I am confused by this example. What should I do? There must be a reason to do so in the example. Would you please to tell me the reason?

[Nikita_S_Intel]

Hi Simbalee,

According to document of djacobi_solve, parameter f1 contains updated function values at X + eps, and f2 for X - eps. I understand that here eps means a small enough value such as 1e-10, and by updating I should recalcaulate function values with all variables pluse or minus eps. Right?

Yes, you are right. This is a general schema of central difference algorithm. DJACOBI routines use this approachfor Jacobi matrix calculation.

However, in the example 14-6 "djacobi_solve usage in C", both f1 and f2 are updated with exactly the same X. I am confused by this example. What should I do? There must be a reason to do so in the example. Would you please to tell me the reason?

The values of X change indirectlyin DJACOBI_SOLVE rouitne (Xi -> Xi+EPS or Xi-EPS). Soyou should calculate the value of fcn() on data X returned by DJACOBI_SOLVE. I mean when DJACOBI_SOLVE askes for fcn ( ... X ... ), X is changed (on Xi+EPS or Xi-EPS) already.

--Nikita

[simbalee]

Quoting - Nikita Shustrov (Intel)

Hi Simbalee,

According to document of djacobi_solve, parameter f1 contains updated function values at X + eps, and f2 for X - eps. I understand that here eps means a small enough value such as 1e-10, and by updating I should recalcaulate function values with all variables pluse or minus eps. Right?

Yes, you are right. This is a general schema of central difference algorithm. DJACOBI routines use this approachfor Jacobi matrix calculation.

However, in the example 14-6 "djacobi_solve usage in C", both f1 and f2 are updated with exactly the same X. I am confused by this example. What should I do? There must be a reason to do so in the example. Would you please to tell me the reason?

The values of X change indirectlyin DJACOBI_SOLVE rouitne (Xi -> Xi+EPS or Xi-EPS). Soyou should calculate the value of fcn() on data X returned by DJACOBI_SOLVE. I mean when DJACOBI_SOLVE askes for fcn ( ... X ... ), X is changed (on Xi+EPS or Xi-EPS) already.

--Nikita

Thanks you very much. Yourr reply is so quick and clear.

[Nikita_S_Intel]

Quoting - simbalee

Thanks you very much. Yourr reply is so quick and clear.

Thank you :)
--Nikita

[simbalee]

Quoting - simbalee

The projection function I am dealing with is highly nonlinear and the function itself is non-analytical. Generally there is no analytical way to get Jacobian matrix. That is why I choose MKL to do this work.
Some C-people instinctively avoid to use global data. In fact I never thought about it untill Artur pointed out. Before trying global data, I may want to use RCI version of djacobi in case that my way failed. I will report result of this work.
I assume that MKL should not internally access over m data in fvec. If memory of m data is not enough for computation, buffer should be allocated inside the djacobi routines, instead of relying on unpredictable size of input vector. But, it would be better to check if the assumption is true.

Well, it turns out that the data-attaching method failed.
Seems that MLK internally copies data to somewhere for further processing. And of coz the attached data is ignored.
I am going to try the RCI version of djacobi.

[simbalee]

I found a new problem of dtrnlsp routines.

dtrnlsp_get is used to retrieve results, and the values are:
iter: 0
st_cr: 0
r1: a big number
r2: 0

The results are wrong, since r2 can not be 0. But it is not the point.
The point is st_cr equals with 0, which is not listed in the document of dtrnlsp_get.

[ArturGuzik]

Quoting - simbalee

I found a new problem of dtrnlsp routines.

dtrnlsp_get is used to retrieve results, and the values are:
iter: 0
st_cr: 0
r1: a big number
r2: 0

The results are wrong, since r2 can not be 0. But it is not the point.
The point is st_cr equals with 0, which is not listed in the document of dtrnlsp_get.

Hi,

looks as initialization/interface problem (to me). You use C(++). What about arguments? Can you post relevant portion of the code, containing calling of MKL routines? Nikita may come up with a better tip/suggestion. I understand that djacobi routine works fine.

A.

[Nikita_S_Intel]

Quoting - ArturGuzik

Hi,

looks as initialization/interface problem (to me). You use C(++). What about arguments? Can you post relevant portion of the code, containing calling of MKL routines? Nikita may come up with a better tip/suggestion. I understand that djacobi routine works fine.

A.

Hi All,

Simbalee, could you pls post your RCI cycle?

--Nikita

[simbalee]

Quoting - ArturGuzik

Hi,

looks as initialization/interface problem (to me). You use C(++). What about arguments? Can you post relevant portion of the code, containing calling of MKL routines? Nikita may come up with a better tip/suggestion. I understand that djacobi routine works fine.

A.

Hi Artur.

Relevant code is posted here.
----------------------------------------------------------------------------------------------------

[cpp]RCI_Req = 10;
success = 0;
loop = 0;
maxloop = MIN_MKL_SOLV_LOOP > iternum ? MIN_MKL_SOLV_LOOP : iternum;

// solve thr problem with MKL nonlinear optimization
while ( 0 == success && loop < maxloop  ) { 
	if ( TR_SUCCESS != dtrnlsp_solve( &handle, pfvec, pfjac, &RCI_Req ) ) {
		fprintf( stderr, "Failed to solve the problem with dtrnlsp_solve!n");
		ret = ERR_MKL_SOVFAIL;
		goto DONE;
	}

	if (  -1 == RCI_Req || -2 == RCI_Req || 
		  -3 == RCI_Req || -4 == RCI_Req || -5 == RCI_Req || -6 == RCI_Req )
		 success = 1;

	if ( 1 == RCI_Req ) {	// recalculate fvec
		calc_proj_pnts_mkl( pxs, pys, px3d, py3d, pz3d, dotnum, camnum, camtmp, pfvec );
	}
	if ( 2 == RCI_Req ) {	// compute Jacobi matrix		
		if ( ERR_NONE != get_jacobi( pxs, pys, px3d, py3d, pz3d, dotnum, camnum, n, camtmp, pfjac, pf1, pf2 ) ) {
			fprintf( stderr, "Failed to get Jacobian matrix with djacobi!n");
			ret = ERR_MKL_JACFAIL;
			goto DONE;
		}
	}

	loop += 1;
}

// maximum loop number over
if( maxloop == loop && 0 == success ) {
	fprintf( stderr, "The optimization does not converge!n");
	ret = ERR_MKL_OVERLOOP;
	goto DONE;
}

// get result
iter = stop_cr = -100;  
r1 = r2 = -100.0;
if ( TR_SUCCESS != dtrnlsp_get( &handle, &iter, &stop_cr, &r1, &r2 ) ) {
	fprintf( stderr, "Failed to get results of optimization!n");
	ret = ERR_MKL_GETFAIL;
	goto DONE;
}
if ( 1 == stop_cr ) { //indicates that the algorithm has exceeded the maximal number of iterations
	fprintf( stderr, "The optimization does not converge!n");
	ret = ERR_MKL_OVERLOOP;
	goto DONE;
}

// delete handles
if ( TR_SUCCESS != dtrnlsp_delete( &handle ) ) {
	fprintf( stderr, "Failed to delete handle of optimization!n");
	ret = ERR_MKL_DELFAIL;
	goto DONE;
}
handle = NULL;

// get residual errors
*err = r2;[/cpp]

----------------------------------------------------------------------------------------------------

The weird thing is that the problem occurs when RCI_Req of dtrnlsp_solve equals with 0.
According to the document, RCI_Req = 0 indicates successful completion of the task.
However, if I use RCI_Req == 0 as a condition for quit from the loop, actually nothing happened.
The "camtmp" that is supposed to be optimized is the same with the one before optimization.
And I got
iter: 0
st_cr: 0
r1: a big number
r2: 0

But if "RCI_Req == 0" is removed from jump out conditions, as I did above, the problem disappeared.

I checked the C example of dtrnlsp again. And I found that "RCI_Req == 0" is also removed from the list of jump out coditions. They may has a reason to do so.

At this stage, optimized results can be obtained with MKL nonliear optimization routines. Although the residual error is quite big. I guess it is due to my own code of calculation.
It will be much better if an improved example of nonliear optimization routines are available.

[simbalee]

Quoting - Nikita Shustrov (Intel)

Hi All,

Simbalee, could you pls post your RCI cycle?

--Nikita

Hi Nikita. I posted some code. Sorry that all spaces at the begining of each line are lost.