optimization routine using conjugate gradient (or other) method?
- To: mathgroup at smc.vnet.net
- Subject: [mg83642] optimization routine using conjugate gradient (or other) method?
- From: dantimatter <dantimatter at gmail.com>
- Date: Tue, 27 Nov 2007 06:08:26 -0500 (EST)
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