what is the structure of fjac vector for dtrnlsp_solve?

what is the structure of fjac vector for dtrnlsp_solve?

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
or

for(j=0;j
Which one is correct?

13 posts / 0 new
Last post
For more complete information about compiler optimizations, see our Optimization Notice.

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
or

for(j=0;j
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[i][j] = 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.

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[i][j] = 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.

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.

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.

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[i]=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 

Quoting - gpwr9k5

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)

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

A.

Hi all,

I would like to recommend to change the precisions: for(int i=0;i<6;i++) eps[i]=1.0e-30;
I suggest to use 1.0e-12 or large.

Thank you!
--Nikita

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[p]=tgt[p]-out[p] as the function values. I compute an SSE as sum(err[p]^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?

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[p]=tgt[p]-out[p] as the function values. I compute an SSE as sum(err[p]^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?

Quoting - gpwr9k5

Hi,

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

A.

Quoting - gpwr9k5

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

A.

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[j]/dVariable[j], i=0..m-1, j=0..n-1
I got the residuals to agree with my calculations: r1/2=sqrt(Function[j]^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], p=0,..,np-1 - evaluated function, which = tgt[p]-Function_Output)

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[i]=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[p]);
				err[p]=tgt[p]-out[nl-1][0];	// only one neuron in the output layer
				mySSE+=err[p]*err[p];
			}
		}
		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[p]);

			//	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[i];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][k]*w[i+1][k][j];
						s[i][j]=afDeriv(out[i][j])*sum;
					}

			//	compute gradients: jac[p][iw]=derr[p]/dw[i][j][k]=s[i][j]*out[i-1][k]
				int iw=0;
				for(int i=1;i < nl;i++)					// for each layer except input
					for(int j=0;j < ls[i];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[i][j];		// bias input
							else jac[p*nw+iw]=s[i][j]*out[i-1][k];	// 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);
}

I think I found my mistake

Leave a Comment

Please sign in to add a comment. Not a member? Join today