Re: optimization routine using conjugate gradient (or other) method?

*To*: mathgroup at smc.vnet.net*Subject*: [mg83677] Re: optimization routine using conjugate gradient (or other) method?*From*: danl at wolfram.com*Date*: Wed, 28 Nov 2007 05:30:04 -0500 (EST)*References*: <figu2a$fem$1@smc.vnet.net>

On Nov 27, 5:12 am, dantimatter <dantimat... at gmail.com> wrote: > Hello all, > > I have a system that I'd like to optimize in 160 dimensions. I have > all the equations I need to set up the problem, but I have no idea how > to use Mathematica to actually do the optimization. You've all been > so helpful in the past, I'm hoping maybe some of you have suggestions > as to how I should proceed. My problem is as follows: > > I have an array of 160 numbers a = Array[f, 160]. I don't know any of > the a[[j]] but I'd like to determine them through the optimization > process. The sum of all elements in the array is atot = Sum[a[[j]], > {j, 160}]. > > I also have two fixed arrays, time and m, where time = {0, 12, 24, 36, > 48, 60, 84, 96, 108, 120, 132, 144, 156} and m = {0.149, 0.223, 0.721, > 2.366, 2.580, 5.904, 7.613, 11.936, 11.943, 8.477, 6.683, 4.641, > 4.238}. > > Lastly, I have a somewhat-complicated "time shift function" h[a_, > t_] := > 1/33 UnitStep[a] UnitStep[ > 64.5 - a] Piecewise[{{UnitStep[32.25 - a] , > t < 127.75}, {0.5 (UnitStep[32.25 - a] + > UnitStep[32.25 - (t - 127.75) - a] + > UnitStep[a - 64.5 + (t - 127.75)]), 127.75 <= t < 160}}, > 0.5 (UnitStep[32.25 - a] + UnitStep[a - 32.25])] > > Now with all this stuff set up, the function I'd like to minimize is > > func = 50 Sum[(a[[j]]/atot) Log[a[[j]]/atot], {j, 160}] + 1/2 > Sum[(m[[i]] - > Sum[a[[j]] h[Mod[time[[i]] - j, 160], time[[i]]], {j, 160}])^2, > {i, Length[time]}] > > The problem is complicated by the fact that the conjugate gradient > optimization routine may go through a point where the a[[j]] have > negative values, which makes it impossible to compute Log[a[[j]]]. > So, an alternative formula is used instead of the a[[j]]/atot) > Log[a[[j]]/atot]: > > 1/atotp (eps Log[eps/atotp] - (eps - a[[j]]) (1 + Log[eps/atotp]) + > (eps - a[[j]])^2/(2 eps)) > > in this case, 'eps' is a predefined small positive number and 'atotp' > is the sum of a[[j]] such that a[[j]] is equal to or greater than > 'eps'. > > So, does anyone have any suggestions as to how I could code this in > Mathematica? Is it something that's even possible? Please let me > know if you have any ideas. > > Many thanks, > Dan I tried this, adding constraints to keep the a[[j]] all strictly positive. constraints = Thread[a > .0001]; NMinimize[{func, constraints}, a] This gives a result of around -238, with the f[j] values in the range between (approx) 1 and 14. I do not know if this is a "sensible" result. Instead of NMinimize you might try FindMinimum. Also there are option settings you can tweak, and they might have some impact on result quality or run time (what I show above takes several minutes, with default settings). Daniel Lichtblau Wolfram Research