Re: what is the structure of fjac vector for dtrnlsp_solve?

gpwr9k5 · ‎07-13-2009

Hello. I am trying to use dtrnlsp_solve function, but I do not know the format of fjac. I calculate fjac myself. According to the manual, it is a 1D array of size m*n where m is the number of function values and n is the number of variables. But, the manual does not say how the fjack elements are arranged. There are two possible ways:

for(i=0;i/dx;

or

for(j=0;j/dx;

Which one is correct?

ArturGuzik · ‎07-13-2009

Quoting - gpwr9k5

Hello. I am trying to use dtrnlsp_solve function, but I do not know the format of fjac. I calculate fjac myself. According to the manual, it is a 1D array of size m*n where m is the number of function values and n is the number of variables. But, the manual does not say how the fjack elements are arranged. There are two possible ways:

for(i=0;i/dx;

or

for(j=0;j/dx;

Which one is correct?

Hi,

The Fortran version is an array FJac(m,n), so I believe it goes the same as with other two-dimensional matrix A of size m x n stored densely in a one-dimensional array:
A = B[i*n+j] in C (i=0, ... , m-1, j=0, ... , n-1)
A(i,j) = B(j*m+i) in Fortran (i=1, ... , m, j=1, ... , n).

so this would be your first version. See Chapter 7 of User's Manual.

A.

gpwr9k5 · ‎07-13-2009

Quoting - ArturGuzik

Hi,

The Fortran version is an array FJac(m,n), so I believe it goes the same as with other two-dimensional matrix A of size m x n stored densely in a one-dimensional array:
A = B[i*n+j] in C (i=0, ... , m-1, j=0, ... , n-1)
A(i,j) = B(j*m+i) in Fortran (i=1, ... , m, j=1, ... , n).

so this would be your first version. See Chapter 7 of User's Manual.

A.

Thanks. That is what I have been using. Something does not work. The solver keeps exiting with 0 iterations and st_cr=6. I tried different values of rs and even changed the sign of jacobian, but still the same result. If I set eps[5]=0 to avoid early exit, then initialization fails. If you can advice on how to avoid early exit with st_cr=6, that would be great.
Thanks again.

ArturGuzik · ‎07-13-2009

Quoting - gpwr9k5

Thanks. That is what I have been using. Something does not work. The solver keeps exiting with 0 iterations and st_cr=6. I tried different values of rs and even changed the sign of jacobian, but still the same result. If I set eps[5]=0 to avoid early exit, then initialization fails. If you can advice on how to avoid early exit with st_cr=6, that would be great.
Thanks again.

Well, sounds as problems with initialization/settings.

(1) does MKL example (from disk) work for you?

(2) show a bit more of _init call, so we can take a look.

str_cr = 6 means that solution converged, so maybe you have a FJac zero(s) or ..... Do you call it properly (handle, pointers etc)?

A.

gpwr9k5 · ‎07-13-2009

Quoting - ArturGuzik

Well, sounds as problems with initialization/settings.

(1) does MKL example (from disk) work for you?

(2) show a bit more of _init call, so we can take a look.

str_cr = 6 means that solution converged, so maybe you have a FJac zero(s) or ..... Do you call it properly (handle, pointers etc)?

A.

The example from disk works, saying "PASS". I am using this solver to train a neural network. I just copied the example code and changed everything that was related to my task. From what I see, all handles and variables get set correctly. Jacobian is not zero. dtrnlsp_solve computes the starting and ending residuals and they are the same. I compute the same rezidual myself (as mySSE in the code below) and it is way different.

[cpp]double NN::train(double **in,double *tgt,int nep,const double maxSSE)
{
	int nit1=1000;			// maximum number of iterations
	int nit2=100;			// maximum number of iterations of calculation of trial-step
	double rs=1.0;			// initial step bound
	double eps[6];			// precisions for stop-criteria
	int RCI_Request;		// reverse communication interface variable
	int successful;			// rci cycle variable
	int nit;				// number of iterations performed
	int st_cr;				// number of stop-criterion
	double init_sse,sse;	// initial and final sse's
	_TRNSP_HANDLE_t handle;	// TR solver handle

//	set precisions for stop-criteria
	for(int i=0;i<6;i++) eps=1.0e-30;
	cout << np << endl;

//	initialize solver (allocate memory, set initial values)
	if(dtrnlsp_init(&handle,&nw,&np,wv,eps,&nit1,&nit2,&rs) != TR_SUCCESS)
	{
		cout << "error in dtrnlsp_initn";
		cin.get();
		return 0;
	}

//	set initial rci cycle variables
	RCI_Request = 0;
	successful = 0;
	double mySSE=0.0;

//	main training cycle
	while(successful == 0)
	{
	//	call tr solver, which reurns RCI_request telling us the next step
		if (dtrnlsp_solve(&handle,err,jac,&RCI_Request) != TR_SUCCESS)
		{
			cout << "error in dtrnlsp_solven";
			cin.get();
			return 0;
		}

	//	act on RCI_request
		if (RCI_Request == -1 || RCI_Request == -2 || RCI_Request == -3 || 
			RCI_Request == -4 || RCI_Request == -5 || RCI_Request == -6)
			successful = 1; // one of the eps is below maxSSE; exit rci cycle
		if (RCI_Request == 1)
		{
		//	recalculate errors
			for(int p=0;p [/cpp]

ArturGuzik · ‎07-13-2009

Quoting - gpwr9k5

[cpp]double NN::train(double **in,double *tgt,int nep,const double maxSSE)
{

//	initialize solver (allocate memory, set initial values)
	if(dtrnlsp_init(&handle,&nw,&np,wv,eps,&nit1,&nit2,&rs) != TR_SUCCESS)

[/cpp]

Quick look -- nothing suspicious. Where do you get/set nw, np, that is sizes of the problem/arrays (used in _init).

A.

Nikita_S_Intel · ‎07-14-2009

Hi all,

I would like to recommend to change the precisions: for(int i=0;i<6;i++) eps=1.0e-30;

I suggest to use 1.0e-12 or large.

Thank you!
--Nikita

gpwr9k5 · ‎07-14-2009

Quoting - ArturGuzik

Quick look -- nothing suspicious. Where do you get/set nw, np, that is sizes of the problem/arrays (used in _init).

A.

I passthem when I call my train() function. nw=71 (number of weoghts/variables), np=1500 (number of patterns/functions). So, my fjac array is huge, 71x1500=106,500. Could that be a problem? I tried to increase esp[] to 1e-12, but now it quits with st_cr=2. My biggest puzzle is why the initial and final residuals computed by the solver are the same and so different from mine. I pass the solver thevalues err

=tgt

-out

as the function values. I compute an SSE as sum(err

^2,p=0..np-1). The solver probably uses something else. For instance, solver computes 42.9975. I compute 77431.3. Even if I divide my sse by (np) or (2*np), I still get a different residual that the solver's. How does the solver computes the residual from the passed function values?

gpwr9k5 · ‎07-14-2009

Quoting - gpwr9k5

I passthem when I call my train() function. nw=71 (number of weoghts/variables), np=1500 (number of patterns/functions). So, my fjac array is huge, 71x1500=106,500. Could that be a problem? I tried to increase esp[] to 1e-12, but now it quits with st_cr=2. My biggest puzzle is why the initial and final residuals computed by the solver are the same and so different from mine. I pass the solver thevalues err
=tgt
-out
as the function values. I compute an SSE as sum(err
^2,p=0..np-1). The solver probably uses something else. For instance, solver computes 42.9975. I compute 77431.3. Even if I divide my sse by (np) or (2*np), I still get a different residual that the solver's. How does the solver computes the residual from the passed function values?

ArturGuzik · ‎07-14-2009

Quoting - gpwr9k5

Hi,

see this thread. You pass the residual. The docs are misleading.

A.

ArturGuzik · ‎07-14-2009

Quoting - gpwr9k5

And one more, your Jac is 1500 x 71, not 71 x 1500, right?

A.

gpwr9k5 · ‎07-15-2009

Quoting - ArturGuzik

And one more, your Jac is 1500 x 71, not 71 x 1500, right?

A.

Number of function values m=1500
Number of variables n=71.
My jac is a 1D array oflength 21,300 computed as jac[i*n+j]=dFunction/dVariable, i=0..m-1, j=0..n-1
I got the residuals to agree with my calculations: r1/2=sqrt(Function^2,j=0..m-1).
I changed the m from 1500 to 300 and still got stopped with st_cr=6. All eps[]=1e-30. I checked my jac values and their magnituderanges from 1e-6 to 1.
Here is the whole code of the function, that is doing all calculations (nw=n - number of variables, np=m - number of function values, err

, p=0,..,np-1 - evaluated function, which = tgt

-Function_Output)

[cpp]double NN::train(double **in,double *tgt,int nep,const double maxSSE)
{
	int nit1=1000;			// maximum number of iterations
	int nit2=100;			// maximum number of iterations of calculation of trial-step
	double rs=1.0;			// initial step bound
	double eps[6];			// precisions for stop-criteria
	int RCI_Request;		// reverse communication interface variable
	int successful;			// rci cycle variable
	int nit;				// number of iterations performed
	int st_cr;				// number of stop-criterion
	double init_sse,sse;	// initial and final sse's
	_TRNSP_HANDLE_t handle;	// TR solver handle

//	set precisions for stop-criteria
	for(int i=0;i<6;i++) eps=1.0e-50;

//	initialize solver (allocate memory, set initial values)
	if(dtrnlsp_init(&handle,&nw,&np,wv,eps,&nit1,&nit2,&rs) != TR_SUCCESS)
	{
		cout << "error in dtrnlsp_initn";
		cin.get();
		return 0;
	}

//	set initial rci cycle variables
	RCI_Request = 0;
	successful = 0;
	double mySSE;

//	main training cycle
	while(successful == 0)
	{
	//	call tr solver, which reurns RCI_request telling us the next step
		if (dtrnlsp_solve(&handle,err,jac,&RCI_Request) != TR_SUCCESS)
		{
			cout << "error in dtrnlsp_solven";
			cin.get();
			return 0;
		}

	//	act on RCI_request
		if (RCI_Request == -1 || RCI_Request == -2 || RCI_Request == -3 || 
			RCI_Request == -4 || RCI_Request == -5 || RCI_Request == -6)
			successful = 1; // one of the eps is below maxSSE; exit rci cycle
		if (RCI_Request == 1)
		{
		//	recalculate errors
			mySSE=0.0;
			for(int p=0;p < np;p++)					// for each training pattern
			{
				ffwd(in

); err

=tgt

-out[nl-1][0]; // only one neuron in the output layer mySSE+=err

*err

; } } if (RCI_Request == 2) { // compute Jacobi matrix cout << "compute Jacobiann"; for(int p=0;p < np;p++) // for each training pattern { // update outputs of each neuron ffwd(in

); // find sensitivity for the output layer i=nl-1 s[nl-1][0]=-1.0; if(oaf==1) s[nl-1][0]*=afDeriv(out[nl-1][0]); // propagate sensitivities from output layer to hidden layers for(int i=nl-2;i > 0;i--) // for each layer except input & output for(int j=0;j < ls;j++) // for each neuron in current layer { double sum=0.0; for(int k=0;k < ls[i+1];k++) // for each neuron in later layer sum+=s[i+1]*w[i+1]; s=afDeriv(out)*sum; } // compute gradients: jac

[iw]=derr

/dw=s*out[i-1] int iw=0; for(int i=1;i < nl;i++) // for each layer except input for(int j=0;j < ls;j++) // for each neuron in current layer for(int k=0;k <= ls[i-1];k++) // for each input of curr neuron incl bias { if(k==ls[i-1]) jac[p*nw+iw]=s; // bias input else jac[p*nw+iw]=s*out[i-1]; // data input if(p*nw+iw>=21100) cout << i << "," << j << "," << k << ": " << jac[p*nw+iw] << endl; iw++; } } } } // get solution statuse if (dtrnlsp_get(&handle,&nit,&st_cr,&init_sse,&sse) != TR_SUCCESS) { cout << "error in dtrnlsp_getn"; cin.get(); return 0; } // free handle memory if (dtrnlsp_delete(&handle) != TR_SUCCESS) { cout << "error in dtrnlsp_deleten"; cin.get(); return 0; } cout << "finished after " << nit << " iterationsn"; cout << "stop criterion is " << st_cr << endl; cout << "sse: " << init_sse << " -> " << sse << endl; cout << "mySSE: " << sqrt(mySSE) << endl; return(sse); }[/cpp]

gpwr9k5 · ‎07-15-2009

I think I found my mistake