Context Navigation

source: trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.5.0/ALGLIB-2.5.0/minlm.cs @ 5138

Visit:

Last change on this file since 5138 was 3839, checked in by mkommend, 15 years ago
implemented first version of LR (ticket #1012)
File size: 47.2 KB

Rev	Line
[3839]	1	/*************************************************************************
	2	Copyright (c) 2009, Sergey Bochkanov (ALGLIB project).
	3
	4	>>> SOURCE LICENSE >>>
	5	This program is free software; you can redistribute it and/or modify
	6	it under the terms of the GNU General Public License as published by
	7	the Free Software Foundation (www.fsf.org); either version 2 of the
	8	License, or (at your option) any later version.
	9
	10	This program is distributed in the hope that it will be useful,
	11	but WITHOUT ANY WARRANTY; without even the implied warranty of
	12	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
	13	GNU General Public License for more details.
	14
	15	A copy of the GNU General Public License is available at
	16	http://www.fsf.org/licensing/licenses
	17
	18	>>> END OF LICENSE >>>
	19	*************************************************************************/
	20
	21	using System;
	22
	23	namespace alglib
	24	{
	25	public class minlm
	26	{
	27	public struct minlmstate
	28	{
	29	public bool wrongparams;
	30	public int n;
	31	public int m;
	32	public double epsg;
	33	public double epsf;
	34	public double epsx;
	35	public int maxits;
	36	public bool xrep;
	37	public double stpmax;
	38	public int flags;
	39	public int usermode;
	40	public double[] x;
	41	public double f;
	42	public double[] fi;
	43	public double[,] j;
	44	public double[,] h;
	45	public double[] g;
	46	public bool needf;
	47	public bool needfg;
	48	public bool needfgh;
	49	public bool needfij;
	50	public bool xupdated;
	51	public minlbfgs.minlbfgsstate internalstate;
	52	public minlbfgs.minlbfgsreport internalrep;
	53	public double[] xprec;
	54	public double[] xbase;
	55	public double[] xdir;
	56	public double[] gbase;
	57	public double[] xprev;
	58	public double fprev;
	59	public double[,] rawmodel;
	60	public double[,] model;
	61	public double[] work;
	62	public AP.rcommstate rstate;
	63	public int repiterationscount;
	64	public int repterminationtype;
	65	public int repnfunc;
	66	public int repnjac;
	67	public int repngrad;
	68	public int repnhess;
	69	public int repncholesky;
	70	public int solverinfo;
	71	public densesolver.densesolverreport solverrep;
	72	public int invinfo;
	73	public matinv.matinvreport invrep;
	74	};
	75
	76
	77	public struct minlmreport
	78	{
	79	public int iterationscount;
	80	public int terminationtype;
	81	public int nfunc;
	82	public int njac;
	83	public int ngrad;
	84	public int nhess;
	85	public int ncholesky;
	86	};
	87
	88
	89
	90
	91	public const int lmmodefj = 0;
	92	public const int lmmodefgj = 1;
	93	public const int lmmodefgh = 2;
	94	public const int lmflagnoprelbfgs = 1;
	95	public const int lmflagnointlbfgs = 2;
	96	public const int lmprelbfgsm = 5;
	97	public const int lmintlbfgsits = 5;
	98	public const int lbfgsnorealloc = 1;
	99
	100
	101	/*************************************************************************
	102	LEVENBERG-MARQUARDT-LIKE METHOD FOR NON-LINEAR OPTIMIZATION
	103
	104	Optimization using function gradient and Hessian. Algorithm - Levenberg-
	105	Marquardt modification with L-BFGS pre-optimization and internal
	106	pre-conditioned L-BFGS optimization after each Levenberg-Marquardt step.
	107
	108	Function F has general form (not "sum-of-squares"):
	109
	110	F = F(x[0], ..., x[n-1])
	111
	112	EXAMPLE
	113
	114	See HTML-documentation.
	115
	116	INPUT PARAMETERS:
	117	N - dimension, N>1
	118	X - initial solution, array[0..N-1]
	119
	120	OUTPUT PARAMETERS:
	121	State - structure which stores algorithm state between subsequent
	122	calls of MinLMIteration. Used for reverse communication.
	123	This structure should be passed to MinLMIteration subroutine.
	124
	125	See also MinLMIteration, MinLMResults.
	126
	127	NOTES:
	128
	129	1. you may tune stopping conditions with MinLMSetCond() function
	130	2. if target function contains exp() or other fast growing functions, and
	131	optimization algorithm makes too large steps which leads to overflow,
	132	use MinLMSetStpMax() function to bound algorithm's steps.
	133
	134	-- ALGLIB --
	135	Copyright 30.03.2009 by Bochkanov Sergey
	136	*************************************************************************/
	137	public static void minlmcreatefgh(int n,
	138	ref double[] x,
	139	ref minlmstate state)
	140	{
	141	int i_ = 0;
	142
	143
	144	//
	145	// Prepare RComm
	146	//
	147	state.rstate.ia = new int[3+1];
	148	state.rstate.ba = new bool[0+1];
	149	state.rstate.ra = new double[7+1];
	150	state.rstate.stage = -1;
	151
	152	//
	153	// prepare internal structures
	154	//
	155	lmprepare(n, 0, true, ref state);
	156
	157	//
	158	// initialize, check parameters
	159	//
	160	minlmsetcond(ref state, 0, 0, 0, 0);
	161	minlmsetxrep(ref state, false);
	162	minlmsetstpmax(ref state, 0);
	163	state.n = n;
	164	state.m = 0;
	165	state.flags = 0;
	166	state.usermode = lmmodefgh;
	167	state.wrongparams = false;
	168	if( n<1 )
	169	{
	170	state.wrongparams = true;
	171	return;
	172	}
	173	for(i_=0; i_<=n-1;i_++)
	174	{
	175	state.x[i_] = x[i_];
	176	}
	177	}
	178
	179
	180	/*************************************************************************
	181	LEVENBERG-MARQUARDT-LIKE METHOD FOR NON-LINEAR OPTIMIZATION
	182
	183	Optimization using function gradient and Jacobian. Algorithm - Levenberg-
	184	Marquardt modification with L-BFGS pre-optimization and internal
	185	pre-conditioned L-BFGS optimization after each Levenberg-Marquardt step.
	186
	187	Function F is represented as sum of squares:
	188
	189	F = f[0]^2(x[0],...,x[n-1]) + ... + f[m-1]^2(x[0],...,x[n-1])
	190
	191	EXAMPLE
	192
	193	See HTML-documentation.
	194
	195	INPUT PARAMETERS:
	196	N - dimension, N>1
	197	M - number of functions f[i]
	198	X - initial solution, array[0..N-1]
	199
	200	OUTPUT PARAMETERS:
	201	State - structure which stores algorithm state between subsequent
	202	calls of MinLMIteration. Used for reverse communication.
	203	This structure should be passed to MinLMIteration subroutine.
	204
	205	See also MinLMIteration, MinLMResults.
	206
	207	NOTES:
	208
	209	1. you may tune stopping conditions with MinLMSetCond() function
	210	2. if target function contains exp() or other fast growing functions, and
	211	optimization algorithm makes too large steps which leads to overflow,
	212	use MinLMSetStpMax() function to bound algorithm's steps.
	213
	214	-- ALGLIB --
	215	Copyright 30.03.2009 by Bochkanov Sergey
	216	*************************************************************************/
	217	public static void minlmcreatefgj(int n,
	218	int m,
	219	ref double[] x,
	220	ref minlmstate state)
	221	{
	222	int i_ = 0;
	223
	224
	225	//
	226	// Prepare RComm
	227	//
	228	state.rstate.ia = new int[3+1];
	229	state.rstate.ba = new bool[0+1];
	230	state.rstate.ra = new double[7+1];
	231	state.rstate.stage = -1;
	232
	233	//
	234	// prepare internal structures
	235	//
	236	lmprepare(n, m, true, ref state);
	237
	238	//
	239	// initialize, check parameters
	240	//
	241	minlmsetcond(ref state, 0, 0, 0, 0);
	242	minlmsetxrep(ref state, false);
	243	minlmsetstpmax(ref state, 0);
	244	state.n = n;
	245	state.m = m;
	246	state.flags = 0;
	247	state.usermode = lmmodefgj;
	248	state.wrongparams = false;
	249	if( n<1 )
	250	{
	251	state.wrongparams = true;
	252	return;
	253	}
	254	for(i_=0; i_<=n-1;i_++)
	255	{
	256	state.x[i_] = x[i_];
	257	}
	258	}
	259
	260
	261	/*************************************************************************
	262	CLASSIC LEVENBERG-MARQUARDT METHOD FOR NON-LINEAR OPTIMIZATION
	263
	264	Optimization using Jacobi matrix. Algorithm - classic Levenberg-Marquardt
	265	method.
	266
	267	Function F is represented as sum of squares:
	268
	269	F = f[0]^2(x[0],...,x[n-1]) + ... + f[m-1]^2(x[0],...,x[n-1])
	270
	271	EXAMPLE
	272
	273	See HTML-documentation.
	274
	275	INPUT PARAMETERS:
	276	N - dimension, N>1
	277	M - number of functions f[i]
	278	X - initial solution, array[0..N-1]
	279
	280	OUTPUT PARAMETERS:
	281	State - structure which stores algorithm state between subsequent
	282	calls of MinLMIteration. Used for reverse communication.
	283	This structure should be passed to MinLMIteration subroutine.
	284
	285	See also MinLMIteration, MinLMResults.
	286
	287	NOTES:
	288
	289	1. you may tune stopping conditions with MinLMSetCond() function
	290	2. if target function contains exp() or other fast growing functions, and
	291	optimization algorithm makes too large steps which leads to overflow,
	292	use MinLMSetStpMax() function to bound algorithm's steps.
	293
	294	-- ALGLIB --
	295	Copyright 30.03.2009 by Bochkanov Sergey
	296	*************************************************************************/
	297	public static void minlmcreatefj(int n,
	298	int m,
	299	ref double[] x,
	300	ref minlmstate state)
	301	{
	302	int i_ = 0;
	303
	304
	305	//
	306	// Prepare RComm
	307	//
	308	state.rstate.ia = new int[3+1];
	309	state.rstate.ba = new bool[0+1];
	310	state.rstate.ra = new double[7+1];
	311	state.rstate.stage = -1;
	312
	313	//
	314	// prepare internal structures
	315	//
	316	lmprepare(n, m, true, ref state);
	317
	318	//
	319	// initialize, check parameters
	320	//
	321	minlmsetcond(ref state, 0, 0, 0, 0);
	322	minlmsetxrep(ref state, false);
	323	minlmsetstpmax(ref state, 0);
	324	state.n = n;
	325	state.m = m;
	326	state.flags = 0;
	327	state.usermode = lmmodefj;
	328	state.wrongparams = false;
	329	if( n<1 )
	330	{
	331	state.wrongparams = true;
	332	return;
	333	}
	334	for(i_=0; i_<=n-1;i_++)
	335	{
	336	state.x[i_] = x[i_];
	337	}
	338	}
	339
	340
	341	/*************************************************************************
	342	This function sets stopping conditions for Levenberg-Marquardt optimization
	343	algorithm.
	344
	345	INPUT PARAMETERS:
	346	State - structure which stores algorithm state between calls and
	347	which is used for reverse communication. Must be initialized
	348	with MinLMCreate???()
	349	EpsG - >=0
	350	The subroutine finishes its work if the condition
	351	\|\|G\|\|<EpsG is satisfied, where \|\|.\|\| means Euclidian norm,
	352	G - gradient.
	353	EpsF - >=0
	354	The subroutine finishes its work if on k+1-th iteration
	355	the condition \|F(k+1)-F(k)\|<=EpsF*max{\|F(k)\|,\|F(k+1)\|,1}
	356	is satisfied.
	357	EpsX - >=0
	358	The subroutine finishes its work if on k+1-th iteration
	359	the condition \|X(k+1)-X(k)\| <= EpsX is fulfilled.
	360	MaxIts - maximum number of iterations. If MaxIts=0, the number of
	361	iterations is unlimited. Only Levenberg-Marquardt
	362	iterations are counted (L-BFGS/CG iterations are NOT
	363	counted because their cost is very low copared to that of
	364	LM).
	365
	366	Passing EpsG=0, EpsF=0, EpsX=0 and MaxIts=0 (simultaneously) will lead to
	367	automatic stopping criterion selection (small EpsX).
	368
	369	-- ALGLIB --
	370	Copyright 02.04.2010 by Bochkanov Sergey
	371	*************************************************************************/
	372	public static void minlmsetcond(ref minlmstate state,
	373	double epsg,
	374	double epsf,
	375	double epsx,
	376	int maxits)
	377	{
	378	System.Diagnostics.Debug.Assert((double)(epsg)>=(double)(0), "MinLMSetCond: negative EpsG!");
	379	System.Diagnostics.Debug.Assert((double)(epsf)>=(double)(0), "MinLMSetCond: negative EpsF!");
	380	System.Diagnostics.Debug.Assert((double)(epsx)>=(double)(0), "MinLMSetCond: negative EpsX!");
	381	System.Diagnostics.Debug.Assert(maxits>=0, "MinLMSetCond: negative MaxIts!");
	382	if( (double)(epsg)==(double)(0) & (double)(epsf)==(double)(0) & (double)(epsx)==(double)(0) & maxits==0 )
	383	{
	384	epsx = 1.0E-6;
	385	}
	386	state.epsg = epsg;
	387	state.epsf = epsf;
	388	state.epsx = epsx;
	389	state.maxits = maxits;
	390	}
	391
	392
	393	/*************************************************************************
	394	This function turns on/off reporting.
	395
	396	INPUT PARAMETERS:
	397	State - structure which stores algorithm state between calls and
	398	which is used for reverse communication. Must be
	399	initialized with MinLMCreate???()
	400	NeedXRep- whether iteration reports are needed or not
	401
	402	Usually algorithm returns from MinLMIteration() only when it needs
	403	function/gradient/Hessian. However, with this function we can let it stop
	404	after each iteration (one iteration may include more than one function
	405	evaluation), which is indicated by XUpdated field.
	406
	407	Both Levenberg-Marquardt and L-BFGS iterations are reported.
	408
	409
	410	-- ALGLIB --
	411	Copyright 02.04.2010 by Bochkanov Sergey
	412	*************************************************************************/
	413	public static void minlmsetxrep(ref minlmstate state,
	414	bool needxrep)
	415	{
	416	state.xrep = needxrep;
	417	}
	418
	419
	420	/*************************************************************************
	421	This function sets maximum step length
	422
	423	INPUT PARAMETERS:
	424	State - structure which stores algorithm state between calls and
	425	which is used for reverse communication. Must be
	426	initialized with MinCGCreate???()
	427	StpMax - maximum step length, >=0. Set StpMax to 0.0, if you don't
	428	want to limit step length.
	429
	430	Use this subroutine when you optimize target function which contains exp()
	431	or other fast growing functions, and optimization algorithm makes too
	432	large steps which leads to overflow. This function allows us to reject
	433	steps that are too large (and therefore expose us to the possible
	434	overflow) without actually calculating function value at the x+stp*d.
	435
	436	NOTE: non-zero StpMax leads to moderate performance degradation because
	437	intermediate step of preconditioned L-BFGS optimization is incompatible
	438	with limits on step size.
	439
	440	-- ALGLIB --
	441	Copyright 02.04.2010 by Bochkanov Sergey
	442	*************************************************************************/
	443	public static void minlmsetstpmax(ref minlmstate state,
	444	double stpmax)
	445	{
	446	System.Diagnostics.Debug.Assert((double)(stpmax)>=(double)(0), "MinLMSetStpMax: StpMax<0!");
	447	state.stpmax = stpmax;
	448	}
	449
	450
	451	/*************************************************************************
	452	One Levenberg-Marquardt iteration.
	453
	454	Called after inialization of State structure with MinLMXXX subroutine.
	455	See HTML docs for examples.
	456
	457	Input parameters:
	458	State - structure which stores algorithm state between subsequent
	459	calls and which is used for reverse communication. Must be
	460	initialized with MinLMXXX call first.
	461
	462	If subroutine returned False, iterative algorithm has converged.
	463
	464	If subroutine returned True, then:
	465	* if State.NeedF=True, - function value F at State.X[0..N-1]
	466	is required
	467	* if State.NeedFG=True - function value F and gradient G
	468	are required
	469	* if State.NeedFiJ=True - function vector f[i] and Jacobi matrix J
	470	are required
	471	* if State.NeedFGH=True - function value F, gradient G and Hesian H
	472	are required
	473	* if State.XUpdated=True - algorithm reports about new iteration,
	474	State.X contains current point,
	475	State.F contains function value.
	476
	477	One and only one of this fields can be set at time.
	478
	479	Results are stored:
	480	* function value - in MinLMState.F
	481	* gradient - in MinLMState.G[0..N-1]
	482	* Jacobi matrix - in MinLMState.J[0..M-1,0..N-1]
	483	* Hessian - in MinLMState.H[0..N-1,0..N-1]
	484
	485	-- ALGLIB --
	486	Copyright 10.03.2009 by Bochkanov Sergey
	487	*************************************************************************/
	488	public static bool minlmiteration(ref minlmstate state)
	489	{
	490	bool result = new bool();
	491	int n = 0;
	492	int m = 0;
	493	int i = 0;
	494	double stepnorm = 0;
	495	bool spd = new bool();
	496	double fbase = 0;
	497	double fnew = 0;
	498	double lambda = 0;
	499	double nu = 0;
	500	double lambdaup = 0;
	501	double lambdadown = 0;
	502	int lbfgsflags = 0;
	503	double v = 0;
	504	int i_ = 0;
	505
	506
	507	//
	508	// Reverse communication preparations
	509	// I know it looks ugly, but it works the same way
	510	// anywhere from C++ to Python.
	511	//
	512	// This code initializes locals by:
	513	// * random values determined during code
	514	// generation - on first subroutine call
	515	// * values from previous call - on subsequent calls
	516	//
	517	if( state.rstate.stage>=0 )
	518	{
	519	n = state.rstate.ia[0];
	520	m = state.rstate.ia[1];
	521	i = state.rstate.ia[2];
	522	lbfgsflags = state.rstate.ia[3];
	523	spd = state.rstate.ba[0];
	524	stepnorm = state.rstate.ra[0];
	525	fbase = state.rstate.ra[1];
	526	fnew = state.rstate.ra[2];
	527	lambda = state.rstate.ra[3];
	528	nu = state.rstate.ra[4];
	529	lambdaup = state.rstate.ra[5];
	530	lambdadown = state.rstate.ra[6];
	531	v = state.rstate.ra[7];
	532	}
	533	else
	534	{
	535	n = -983;
	536	m = -989;
	537	i = -834;
	538	lbfgsflags = 900;
	539	spd = true;
	540	stepnorm = 364;
	541	fbase = 214;
	542	fnew = -338;
	543	lambda = -686;
	544	nu = 912;
	545	lambdaup = 585;
	546	lambdadown = 497;
	547	v = -271;
	548	}
	549	if( state.rstate.stage==0 )
	550	{
	551	goto lbl_0;
	552	}
	553	if( state.rstate.stage==1 )
	554	{
	555	goto lbl_1;
	556	}
	557	if( state.rstate.stage==2 )
	558	{
	559	goto lbl_2;
	560	}
	561	if( state.rstate.stage==3 )
	562	{
	563	goto lbl_3;
	564	}
	565	if( state.rstate.stage==4 )
	566	{
	567	goto lbl_4;
	568	}
	569	if( state.rstate.stage==5 )
	570	{
	571	goto lbl_5;
	572	}
	573	if( state.rstate.stage==6 )
	574	{
	575	goto lbl_6;
	576	}
	577	if( state.rstate.stage==7 )
	578	{
	579	goto lbl_7;
	580	}
	581	if( state.rstate.stage==8 )
	582	{
	583	goto lbl_8;
	584	}
	585	if( state.rstate.stage==9 )
	586	{
	587	goto lbl_9;
	588	}
	589	if( state.rstate.stage==10 )
	590	{
	591	goto lbl_10;
	592	}
	593	if( state.rstate.stage==11 )
	594	{
	595	goto lbl_11;
	596	}
	597	if( state.rstate.stage==12 )
	598	{
	599	goto lbl_12;
	600	}
	601	if( state.rstate.stage==13 )
	602	{
	603	goto lbl_13;
	604	}
	605	if( state.rstate.stage==14 )
	606	{
	607	goto lbl_14;
	608	}
	609	if( state.rstate.stage==15 )
	610	{
	611	goto lbl_15;
	612	}
	613
	614	//
	615	// Routine body
	616	//
	617	System.Diagnostics.Debug.Assert(state.usermode==lmmodefj \| state.usermode==lmmodefgj \| state.usermode==lmmodefgh, "LM: internal error");
	618	if( state.wrongparams )
	619	{
	620	state.repterminationtype = -1;
	621	result = false;
	622	return result;
	623	}
	624
	625	//
	626	// prepare params
	627	//
	628	n = state.n;
	629	m = state.m;
	630	lambdaup = 20;
	631	lambdadown = 0.5;
	632	nu = 1;
	633	lbfgsflags = 0;
	634
	635	//
	636	// if we have F and G
	637	//
	638	if( ! ((state.usermode==lmmodefgj \| state.usermode==lmmodefgh) & state.flags/lmflagnoprelbfgs%2==0) )
	639	{
	640	goto lbl_16;
	641	}
	642
	643	//
	644	// First stage of the hybrid algorithm: LBFGS
	645	//
	646	minlbfgs.minlbfgscreate(n, Math.Min(n, lmprelbfgsm), ref state.x, ref state.internalstate);
	647	minlbfgs.minlbfgssetcond(ref state.internalstate, 0, 0, 0, Math.Max(5, n));
	648	minlbfgs.minlbfgssetxrep(ref state.internalstate, state.xrep);
	649	minlbfgs.minlbfgssetstpmax(ref state.internalstate, state.stpmax);
	650	lbl_18:
	651	if( ! minlbfgs.minlbfgsiteration(ref state.internalstate) )
	652	{
	653	goto lbl_19;
	654	}
	655	if( ! state.internalstate.needfg )
	656	{
	657	goto lbl_20;
	658	}
	659
	660	//
	661	// RComm
	662	//
	663	for(i_=0; i_<=n-1;i_++)
	664	{
	665	state.x[i_] = state.internalstate.x[i_];
	666	}
	667	lmclearrequestfields(ref state);
	668	state.needfg = true;
	669	state.rstate.stage = 0;
	670	goto lbl_rcomm;
	671	lbl_0:
	672	state.repnfunc = state.repnfunc+1;
	673	state.repngrad = state.repngrad+1;
	674
	675	//
	676	// Call LBFGS
	677	//
	678	state.internalstate.f = state.f;
	679	for(i_=0; i_<=n-1;i_++)
	680	{
	681	state.internalstate.g[i_] = state.g[i_];
	682	}
	683	lbl_20:
	684	if( ! (state.internalstate.xupdated & state.xrep) )
	685	{
	686	goto lbl_22;
	687	}
	688	lmclearrequestfields(ref state);
	689	state.f = state.internalstate.f;
	690	for(i_=0; i_<=n-1;i_++)
	691	{
	692	state.x[i_] = state.internalstate.x[i_];
	693	}
	694	state.xupdated = true;
	695	state.rstate.stage = 1;
	696	goto lbl_rcomm;
	697	lbl_1:
	698	lbl_22:
	699	goto lbl_18;
	700	lbl_19:
	701	minlbfgs.minlbfgsresults(ref state.internalstate, ref state.x, ref state.internalrep);
	702	goto lbl_17;
	703	lbl_16:
	704
	705	//
	706	// No first stage.
	707	// However, we may need to report initial point
	708	//
	709	if( ! state.xrep )
	710	{
	711	goto lbl_24;
	712	}
	713	lmclearrequestfields(ref state);
	714	state.needf = true;
	715	state.rstate.stage = 2;
	716	goto lbl_rcomm;
	717	lbl_2:
	718	lmclearrequestfields(ref state);
	719	state.xupdated = true;
	720	state.rstate.stage = 3;
	721	goto lbl_rcomm;
	722	lbl_3:
	723	lbl_24:
	724	lbl_17:
	725
	726	//
	727	// Second stage of the hybrid algorithm: LM
	728	// Initialize quadratic model.
	729	//
	730	if( state.usermode!=lmmodefgh )
	731	{
	732	goto lbl_26;
	733	}
	734
	735	//
	736	// RComm
	737	//
	738	lmclearrequestfields(ref state);
	739	state.needfgh = true;
	740	state.rstate.stage = 4;
	741	goto lbl_rcomm;
	742	lbl_4:
	743	state.repnfunc = state.repnfunc+1;
	744	state.repngrad = state.repngrad+1;
	745	state.repnhess = state.repnhess+1;
	746
	747	//
	748	// generate raw quadratic model
	749	//
	750	ablas.rmatrixcopy(n, n, ref state.h, 0, 0, ref state.rawmodel, 0, 0);
	751	for(i_=0; i_<=n-1;i_++)
	752	{
	753	state.gbase[i_] = state.g[i_];
	754	}
	755	fbase = state.f;
	756	lbl_26:
	757	if( ! (state.usermode==lmmodefgj \| state.usermode==lmmodefj) )
	758	{
	759	goto lbl_28;
	760	}
	761
	762	//
	763	// RComm
	764	//
	765	lmclearrequestfields(ref state);
	766	state.needfij = true;
	767	state.rstate.stage = 5;
	768	goto lbl_rcomm;
	769	lbl_5:
	770	state.repnfunc = state.repnfunc+1;
	771	state.repnjac = state.repnjac+1;
	772
	773	//
	774	// generate raw quadratic model
	775	//
	776	ablas.rmatrixgemm(n, n, m, 2.0, ref state.j, 0, 0, 1, ref state.j, 0, 0, 0, 0.0, ref state.rawmodel, 0, 0);
	777	ablas.rmatrixmv(n, m, ref state.j, 0, 0, 1, ref state.fi, 0, ref state.gbase, 0);
	778	for(i_=0; i_<=n-1;i_++)
	779	{
	780	state.gbase[i_] = 2*state.gbase[i_];
	781	}
	782	fbase = 0.0;
	783	for(i_=0; i_<=m-1;i_++)
	784	{
	785	fbase += state.fi[i_]*state.fi[i_];
	786	}
	787	lbl_28:
	788	lambda = 0.001;
	789	lbl_30:
	790	if( false )
	791	{
	792	goto lbl_31;
	793	}
	794
	795	//
	796	// 1. Model = RawModel+lambda*I
	797	// 2. Try to solve (RawModel+LambdaI)dx = -g.
	798	// Increase lambda if left part is not positive definite.
	799	//
	800	for(i=0; i<=n-1; i++)
	801	{
	802	for(i_=0; i_<=n-1;i_++)
	803	{
	804	state.model[i,i_] = state.rawmodel[i,i_];
	805	}
	806	state.model[i,i] = state.model[i,i]+lambda;
	807	}
	808	spd = trfac.spdmatrixcholesky(ref state.model, n, true);
	809	state.repncholesky = state.repncholesky+1;
	810	if( spd )
	811	{
	812	goto lbl_32;
	813	}
	814	if( ! increaselambda(ref lambda, ref nu, lambdaup) )
	815	{
	816	goto lbl_34;
	817	}
	818	goto lbl_30;
	819	goto lbl_35;
	820	lbl_34:
	821	state.repterminationtype = 7;
	822	lmclearrequestfields(ref state);
	823	state.needf = true;
	824	state.rstate.stage = 6;
	825	goto lbl_rcomm;
	826	lbl_6:
	827	goto lbl_31;
	828	lbl_35:
	829	lbl_32:
	830	densesolver.spdmatrixcholeskysolve(ref state.model, n, true, ref state.gbase, ref state.solverinfo, ref state.solverrep, ref state.xdir);
	831	if( state.solverinfo>=0 )
	832	{
	833	goto lbl_36;
	834	}
	835	if( ! increaselambda(ref lambda, ref nu, lambdaup) )
	836	{
	837	goto lbl_38;
	838	}
	839	goto lbl_30;
	840	goto lbl_39;
	841	lbl_38:
	842	state.repterminationtype = 7;
	843	lmclearrequestfields(ref state);
	844	state.needf = true;
	845	state.rstate.stage = 7;
	846	goto lbl_rcomm;
	847	lbl_7:
	848	goto lbl_31;
	849	lbl_39:
	850	lbl_36:
	851	for(i_=0; i_<=n-1;i_++)
	852	{
	853	state.xdir[i_] = -1*state.xdir[i_];
	854	}
	855
	856	//
	857	// Candidate lambda is found.
	858	// 1. Save old w in WBase
	859	// 1. Test some stopping criterions
	860	// 2. If error(w+wdir)>error(w), increase lambda
	861	//
	862	for(i_=0; i_<=n-1;i_++)
	863	{
	864	state.xprev[i_] = state.x[i_];
	865	}
	866	state.fprev = state.f;
	867	for(i_=0; i_<=n-1;i_++)
	868	{
	869	state.xbase[i_] = state.x[i_];
	870	}
	871	for(i_=0; i_<=n-1;i_++)
	872	{
	873	state.x[i_] = state.x[i_] + state.xdir[i_];
	874	}
	875	stepnorm = 0.0;
	876	for(i_=0; i_<=n-1;i_++)
	877	{
	878	stepnorm += state.xdir[i_]*state.xdir[i_];
	879	}
	880	stepnorm = Math.Sqrt(stepnorm);
	881	if( ! ((double)(state.stpmax)>(double)(0) & (double)(stepnorm)>(double)(state.stpmax)) )
	882	{
	883	goto lbl_40;
	884	}
	885
	886	//
	887	// Step is larger than the limit,
	888	// larger lambda is needed
	889	//
	890	for(i_=0; i_<=n-1;i_++)
	891	{
	892	state.x[i_] = state.xbase[i_];
	893	}
	894	if( ! increaselambda(ref lambda, ref nu, lambdaup) )
	895	{
	896	goto lbl_42;
	897	}
	898	goto lbl_30;
	899	goto lbl_43;
	900	lbl_42:
	901	state.repterminationtype = 7;
	902	for(i_=0; i_<=n-1;i_++)
	903	{
	904	state.x[i_] = state.xprev[i_];
	905	}
	906	lmclearrequestfields(ref state);
	907	state.needf = true;
	908	state.rstate.stage = 8;
	909	goto lbl_rcomm;
	910	lbl_8:
	911	goto lbl_31;
	912	lbl_43:
	913	lbl_40:
	914	lmclearrequestfields(ref state);
	915	state.needf = true;
	916	state.rstate.stage = 9;
	917	goto lbl_rcomm;
	918	lbl_9:
	919	state.repnfunc = state.repnfunc+1;
	920	fnew = state.f;
	921	if( (double)(fnew)<=(double)(fbase) )
	922	{
	923	goto lbl_44;
	924	}
	925
	926	//
	927	// restore state and continue search for lambda
	928	//
	929	for(i_=0; i_<=n-1;i_++)
	930	{
	931	state.x[i_] = state.xbase[i_];
	932	}
	933	if( ! increaselambda(ref lambda, ref nu, lambdaup) )
	934	{
	935	goto lbl_46;
	936	}
	937	goto lbl_30;
	938	goto lbl_47;
	939	lbl_46:
	940	state.repterminationtype = 7;
	941	for(i_=0; i_<=n-1;i_++)
	942	{
	943	state.x[i_] = state.xprev[i_];
	944	}
	945	lmclearrequestfields(ref state);
	946	state.needf = true;
	947	state.rstate.stage = 10;
	948	goto lbl_rcomm;
	949	lbl_10:
	950	goto lbl_31;
	951	lbl_47:
	952	lbl_44:
	953	if( ! ((double)(state.stpmax)==(double)(0) & (state.usermode==lmmodefgj \| state.usermode==lmmodefgh) & state.flags/lmflagnointlbfgs%2==0) )
	954	{
	955	goto lbl_48;
	956	}
	957
	958	//
	959	// Optimize using LBFGS, with inv(cholesky(H)) as preconditioner.
	960	//
	961	// It is possible only when StpMax=0, because we can't guarantee
	962	// that step remains bounded when preconditioner is used (we need
	963	// SVD decomposition to do that, which is too slow).
	964	//
	965	matinv.rmatrixtrinverse(ref state.model, n, true, false, ref state.invinfo, ref state.invrep);
	966	if( state.invinfo<=0 )
	967	{
	968	goto lbl_50;
	969	}
	970
	971	//
	972	// if matrix can be inverted, use it.
	973	// just silently move to next iteration otherwise.
	974	// (will be very, very rare, mostly for specially
	975	// designed near-degenerate tasks)
	976	//
	977	for(i_=0; i_<=n-1;i_++)
	978	{
	979	state.xbase[i_] = state.x[i_];
	980	}
	981	for(i=0; i<=n-1; i++)
	982	{
	983	state.xprec[i] = 0;
	984	}
	985	minlbfgs.minlbfgscreatex(n, Math.Min(n, lmintlbfgsits), ref state.xprec, lbfgsflags, ref state.internalstate);
	986	minlbfgs.minlbfgssetcond(ref state.internalstate, 0, 0, 0, lmintlbfgsits);
	987	lbl_52:
	988	if( ! minlbfgs.minlbfgsiteration(ref state.internalstate) )
	989	{
	990	goto lbl_53;
	991	}
	992
	993	//
	994	// convert XPrec to unpreconditioned form, then call RComm.
	995	//
	996	for(i=0; i<=n-1; i++)
	997	{
	998	v = 0.0;
	999	for(i_=i; i_<=n-1;i_++)
	1000	{
	1001	v += state.internalstate.x[i_]*state.model[i,i_];
	1002	}
	1003	state.x[i] = state.xbase[i]+v;
	1004	}
	1005	lmclearrequestfields(ref state);
	1006	state.needfg = true;
	1007	state.rstate.stage = 11;
	1008	goto lbl_rcomm;
	1009	lbl_11:
	1010	state.repnfunc = state.repnfunc+1;
	1011	state.repngrad = state.repngrad+1;
	1012
	1013	//
	1014	// 1. pass State.F to State.InternalState.F
	1015	// 2. convert gradient back to preconditioned form
	1016	//
	1017	state.internalstate.f = state.f;
	1018	for(i=0; i<=n-1; i++)
	1019	{
	1020	state.internalstate.g[i] = 0;
	1021	}
	1022	for(i=0; i<=n-1; i++)
	1023	{
	1024	v = state.g[i];
	1025	for(i_=i; i_<=n-1;i_++)
	1026	{
	1027	state.internalstate.g[i_] = state.internalstate.g[i_] + v*state.model[i,i_];
	1028	}
	1029	}
	1030
	1031	//
	1032	// next iteration
	1033	//
	1034	goto lbl_52;
	1035	lbl_53:
	1036
	1037	//
	1038	// change LBFGS flags to NoRealloc.
	1039	// L-BFGS subroutine will use memory allocated from previous run.
	1040	// it is possible since all subsequent calls will be with same N/M.
	1041	//
	1042	lbfgsflags = lbfgsnorealloc;
	1043
	1044	//
	1045	// back to unpreconditioned X
	1046	//
	1047	minlbfgs.minlbfgsresults(ref state.internalstate, ref state.xprec, ref state.internalrep);
	1048	for(i=0; i<=n-1; i++)
	1049	{
	1050	v = 0.0;
	1051	for(i_=i; i_<=n-1;i_++)
	1052	{
	1053	v += state.xprec[i_]*state.model[i,i_];
	1054	}
	1055	state.x[i] = state.xbase[i]+v;
	1056	}
	1057	lbl_50:
	1058	lbl_48:
	1059
	1060	//
	1061	// Composite iteration is almost over:
	1062	// * accept new position.
	1063	// * rebuild quadratic model
	1064	//
	1065	state.repiterationscount = state.repiterationscount+1;
	1066	if( state.usermode!=lmmodefgh )
	1067	{
	1068	goto lbl_54;
	1069	}
	1070	lmclearrequestfields(ref state);
	1071	state.needfgh = true;
	1072	state.rstate.stage = 12;
	1073	goto lbl_rcomm;
	1074	lbl_12:
	1075	state.repnfunc = state.repnfunc+1;
	1076	state.repngrad = state.repngrad+1;
	1077	state.repnhess = state.repnhess+1;
	1078	ablas.rmatrixcopy(n, n, ref state.h, 0, 0, ref state.rawmodel, 0, 0);
	1079	for(i_=0; i_<=n-1;i_++)
	1080	{
	1081	state.gbase[i_] = state.g[i_];
	1082	}
	1083	fnew = state.f;
	1084	lbl_54:
	1085	if( ! (state.usermode==lmmodefgj \| state.usermode==lmmodefj) )
	1086	{
	1087	goto lbl_56;
	1088	}
	1089	lmclearrequestfields(ref state);
	1090	state.needfij = true;
	1091	state.rstate.stage = 13;
	1092	goto lbl_rcomm;
	1093	lbl_13:
	1094	state.repnfunc = state.repnfunc+1;
	1095	state.repnjac = state.repnjac+1;
	1096	ablas.rmatrixgemm(n, n, m, 2.0, ref state.j, 0, 0, 1, ref state.j, 0, 0, 0, 0.0, ref state.rawmodel, 0, 0);
	1097	ablas.rmatrixmv(n, m, ref state.j, 0, 0, 1, ref state.fi, 0, ref state.gbase, 0);
	1098	for(i_=0; i_<=n-1;i_++)
	1099	{
	1100	state.gbase[i_] = 2*state.gbase[i_];
	1101	}
	1102	fnew = 0.0;
	1103	for(i_=0; i_<=m-1;i_++)
	1104	{
	1105	fnew += state.fi[i_]*state.fi[i_];
	1106	}
	1107	lbl_56:
	1108
	1109	//
	1110	// Stopping conditions
	1111	//
	1112	for(i_=0; i_<=n-1;i_++)
	1113	{
	1114	state.work[i_] = state.xprev[i_];
	1115	}
	1116	for(i_=0; i_<=n-1;i_++)
	1117	{
	1118	state.work[i_] = state.work[i_] - state.x[i_];
	1119	}
	1120	stepnorm = 0.0;
	1121	for(i_=0; i_<=n-1;i_++)
	1122	{
	1123	stepnorm += state.work[i_]*state.work[i_];
	1124	}
	1125	stepnorm = Math.Sqrt(stepnorm);
	1126	if( (double)(stepnorm)<=(double)(state.epsx) )
	1127	{
	1128	state.repterminationtype = 2;
	1129	goto lbl_31;
	1130	}
	1131	if( state.repiterationscount>=state.maxits & state.maxits>0 )
	1132	{
	1133	state.repterminationtype = 5;
	1134	goto lbl_31;
	1135	}
	1136	v = 0.0;
	1137	for(i_=0; i_<=n-1;i_++)
	1138	{
	1139	v += state.gbase[i_]*state.gbase[i_];
	1140	}
	1141	v = Math.Sqrt(v);
	1142	if( (double)(v)<=(double)(state.epsg) )
	1143	{
	1144	state.repterminationtype = 4;
	1145	goto lbl_31;
	1146	}
	1147	if( (double)(Math.Abs(fnew-fbase))<=(double)(state.epsf*Math.Max(1, Math.Max(Math.Abs(fnew), Math.Abs(fbase)))) )
	1148	{
	1149	state.repterminationtype = 1;
	1150	goto lbl_31;
	1151	}
	1152
	1153	//
	1154	// Now, iteration is finally over:
	1155	// * update FBase
	1156	// * decrease lambda
	1157	// * report new iteration
	1158	//
	1159	if( ! state.xrep )
	1160	{
	1161	goto lbl_58;
	1162	}
	1163	lmclearrequestfields(ref state);
	1164	state.xupdated = true;
	1165	state.f = fnew;
	1166	state.rstate.stage = 14;
	1167	goto lbl_rcomm;
	1168	lbl_14:
	1169	lbl_58:
	1170	fbase = fnew;
	1171	decreaselambda(ref lambda, ref nu, lambdadown);
	1172	goto lbl_30;
	1173	lbl_31:
	1174
	1175	//
	1176	// final point is reported
	1177	//
	1178	if( ! state.xrep )
	1179	{
	1180	goto lbl_60;
	1181	}
	1182	lmclearrequestfields(ref state);
	1183	state.xupdated = true;
	1184	state.f = fnew;
	1185	state.rstate.stage = 15;
	1186	goto lbl_rcomm;
	1187	lbl_15:
	1188	lbl_60:
	1189	result = false;
	1190	return result;
	1191
	1192	//
	1193	// Saving state
	1194	//
	1195	lbl_rcomm:
	1196	result = true;
	1197	state.rstate.ia[0] = n;
	1198	state.rstate.ia[1] = m;
	1199	state.rstate.ia[2] = i;
	1200	state.rstate.ia[3] = lbfgsflags;
	1201	state.rstate.ba[0] = spd;
	1202	state.rstate.ra[0] = stepnorm;
	1203	state.rstate.ra[1] = fbase;
	1204	state.rstate.ra[2] = fnew;
	1205	state.rstate.ra[3] = lambda;
	1206	state.rstate.ra[4] = nu;
	1207	state.rstate.ra[5] = lambdaup;
	1208	state.rstate.ra[6] = lambdadown;
	1209	state.rstate.ra[7] = v;
	1210	return result;
	1211	}
	1212
	1213
	1214	/*************************************************************************
	1215	Levenberg-Marquardt algorithm results
	1216
	1217	Called after MinLMIteration returned False.
	1218
	1219	Input parameters:
	1220	State - algorithm state (used by MinLMIteration).
	1221
	1222	Output parameters:
	1223	X - array[0..N-1], solution
	1224	Rep - optimization report:
	1225	* Rep.TerminationType completetion code:
	1226	* -1 incorrect parameters were specified
	1227	* 1 relative function improvement is no more than
	1228	EpsF.
	1229	* 2 relative step is no more than EpsX.
	1230	* 4 gradient is no more than EpsG.
	1231	* 5 MaxIts steps was taken
	1232	* 7 stopping conditions are too stringent,
	1233	further improvement is impossible
	1234	* Rep.IterationsCount contains iterations count
	1235	* Rep.NFunc - number of function calculations
	1236	* Rep.NJac - number of Jacobi matrix calculations
	1237	* Rep.NGrad - number of gradient calculations
	1238	* Rep.NHess - number of Hessian calculations
	1239	* Rep.NCholesky - number of Cholesky decomposition calculations
	1240
	1241	-- ALGLIB --
	1242	Copyright 10.03.2009 by Bochkanov Sergey
	1243	*************************************************************************/
	1244	public static void minlmresults(ref minlmstate state,
	1245	ref double[] x,
	1246	ref minlmreport rep)
	1247	{
	1248	int i_ = 0;
	1249
	1250	x = new double[state.n-1+1];
	1251	for(i_=0; i_<=state.n-1;i_++)
	1252	{
	1253	x[i_] = state.x[i_];
	1254	}
	1255	rep.iterationscount = state.repiterationscount;
	1256	rep.terminationtype = state.repterminationtype;
	1257	rep.nfunc = state.repnfunc;
	1258	rep.njac = state.repnjac;
	1259	rep.ngrad = state.repngrad;
	1260	rep.nhess = state.repnhess;
	1261	rep.ncholesky = state.repncholesky;
	1262	}
	1263
	1264
	1265	/*************************************************************************
	1266	Prepare internal structures (except for RComm).
	1267
	1268	Note: M must be zero for FGH mode, non-zero for FJ/FGJ mode.
	1269	*************************************************************************/
	1270	private static void lmprepare(int n,
	1271	int m,
	1272	bool havegrad,
	1273	ref minlmstate state)
	1274	{
	1275	state.repiterationscount = 0;
	1276	state.repterminationtype = 0;
	1277	state.repnfunc = 0;
	1278	state.repnjac = 0;
	1279	state.repngrad = 0;
	1280	state.repnhess = 0;
	1281	state.repncholesky = 0;
	1282	if( n<=0 \| m<0 )
	1283	{
	1284	return;
	1285	}
	1286	if( havegrad )
	1287	{
	1288	state.g = new double[n-1+1];
	1289	}
	1290	if( m!=0 )
	1291	{
	1292	state.j = new double[m-1+1, n-1+1];
	1293	state.fi = new double[m-1+1];
	1294	state.h = new double[0+1, 0+1];
	1295	}
	1296	else
	1297	{
	1298	state.j = new double[0+1, 0+1];
	1299	state.fi = new double[0+1];
	1300	state.h = new double[n-1+1, n-1+1];
	1301	}
	1302	state.x = new double[n-1+1];
	1303	state.rawmodel = new double[n-1+1, n-1+1];
	1304	state.model = new double[n-1+1, n-1+1];
	1305	state.xbase = new double[n-1+1];
	1306	state.xprec = new double[n-1+1];
	1307	state.gbase = new double[n-1+1];
	1308	state.xdir = new double[n-1+1];
	1309	state.xprev = new double[n-1+1];
	1310	state.work = new double[Math.Max(n, m)+1];
	1311	}
	1312
	1313
	1314	/*************************************************************************
	1315	Clears request fileds (to be sure that we don't forgot to clear something)
	1316	*************************************************************************/
	1317	private static void lmclearrequestfields(ref minlmstate state)
	1318	{
	1319	state.needf = false;
	1320	state.needfg = false;
	1321	state.needfgh = false;
	1322	state.needfij = false;
	1323	state.xupdated = false;
	1324	}
	1325
	1326
	1327	/*************************************************************************
	1328	Increases lambda, returns False when there is a danger of overflow
	1329	*************************************************************************/
	1330	private static bool increaselambda(ref double lambda,
	1331	ref double nu,
	1332	double lambdaup)
	1333	{
	1334	bool result = new bool();
	1335	double lnlambda = 0;
	1336	double lnnu = 0;
	1337	double lnlambdaup = 0;
	1338	double lnmax = 0;
	1339
	1340	result = false;
	1341	lnlambda = Math.Log(lambda);
	1342	lnlambdaup = Math.Log(lambdaup);
	1343	lnnu = Math.Log(nu);
	1344	lnmax = Math.Log(AP.Math.MaxRealNumber);
	1345	if( (double)(lnlambda+lnlambdaup+lnnu)>(double)(lnmax) )
	1346	{
	1347	return result;
	1348	}
	1349	if( (double)(lnnu+Math.Log(2))>(double)(lnmax) )
	1350	{
	1351	return result;
	1352	}
	1353	lambda = lambdalambdaupnu;
	1354	nu = nu*2;
	1355	result = true;
	1356	return result;
	1357	}
	1358
	1359
	1360	/*************************************************************************
	1361	Decreases lambda, but leaves it unchanged when there is danger of underflow.
	1362	*************************************************************************/
	1363	private static void decreaselambda(ref double lambda,
	1364	ref double nu,
	1365	double lambdadown)
	1366	{
	1367	nu = 1;
	1368	if( (double)(Math.Log(lambda)+Math.Log(lambdadown))<(double)(Math.Log(AP.Math.MinRealNumber)) )
	1369	{
	1370	lambda = AP.Math.MinRealNumber;
	1371	}
	1372	else
	1373	{
	1374	lambda = lambda*lambdadown;
	1375	}
	1376	}
	1377	}
	1378	}

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Update cookies preferences