Re: optimization routine using conjugate gradient (or other) method?
- To: mathgroup at smc.vnet.net
- Subject: [mg83769] Re: optimization routine using conjugate gradient (or other) method?
- From: dantimatter <google at dantimatter.com>
- Date: Fri, 30 Nov 2007 05:17:03 -0500 (EST)
- References: <9889841.1196176137975.JavaMail.root@m35> <fim9d1$s0c$1@smc.vnet.net>
thanks to all for responding. this is exactly what i needed. may i
ask one further question: how does one search the solution landscape
to determine if there are other local minima?
many thanks,
dan
On Nov 29, 5:56 am, DrMajorBob <drmajor... at bigfoot.com> wrote:
> My result seems close to Daniel Lichtblau's, although I've rewritten the
> code:
>
> Clear[a, h, ahmodT, sqSum, func]
> h[t_][x_List] := h[t] /@ x
> h[t_][x_] = (1/33 UnitStep[a, 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])] //
> Rationalize // PiecewiseExpand // FullSimplify) /. a -> x;
> ahmodT[a_List][t_] := a.h[t]@Mod[t - Range@160, 160]
> sqSum[a_List] := #.# &[m - ahmodT[a] /@ time]
> func[a_] :=
> Module[{an = Normalize[a, Total]}, 50 an.Log[an] + 1/2 sqSum[a]
> ]
>
> time = {0, 12, 24, 36, 48, 60, 84, 96, 108, 120, 132, 144, 156};
> 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};
> a = Array[f, 160];
>
> The objective function is VERY complicated:
>
> func[a] // LeafCount
>
> 107403
>
> But NMinimize finishes in under 81 seconds:
>
> constraints = Thread[a > .0001];
> Timing[{obj1, rule1} = NMinimize[{func[a], constraints}, a];]
> obj1
> a1 = a /. rule1;
> {lo, hi} = Through[{Min, Max}@a1]
>
> {80.89, Null}
>
> -237.938
>
> {0.811749, 15.0187}
>
> A random search doesn't find anything better (or even nearly as good):
>
> randomStart := RandomReal[{.001, 16}, Length@a]
> Timing[Table[func[randomStart], {5000}] // Min]
>
> {142.891, -158.832}
>
> But that's a tiny sample considering the dimensions!!
>
> I also tried FindMinimum, but it seemed to take FOREVER, so I aborted.
>
> Bobby
>
> On Tue, 27 Nov 2007 05:08:26 -0600, 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, 3=
> 6,
> > 48, 60, 84, 96, 108, 120, 132, 144, 156} and m = {0.149, 0.223, 0.72=
> 1,
> > 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
>
> -- =
>
> DrMajor... at bigfoot.com