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