Re: Re: optimization routine using conjugate gradient (or other) method?
- To: mathgroup at smc.vnet.net
- Subject: [mg83815] Re: [mg83769] Re: optimization routine using conjugate gradient (or other) method?
- From: DrMajorBob <drmajorbob at bigfoot.com>
- Date: Sat, 1 Dec 2007 05:45:01 -0500 (EST)
- References: <9889841.1196176137975.JavaMail.root@m35> <fim9d1$s0c$1@smc.vnet.net> <24166678.1196460323231.JavaMail.root@m35>
- Reply-to: drmajorbob at bigfoot.com
> ask one further question: how does one search the solution landscape
> to determine if there are other local minima?
That's exactly what my randomStart and Table tried to accomplish, but
you'd need trillions of samples to be sure of anything, with 160
dimensions!
You can also use randomStart (or something like it) as a start for
FindMinimum -- that's what I did first, in fact, but it was very slow
(running just ONCE) and I didn't care to wait. Again, with 160 dimensions
it's hard to imagine sampling enough to be sure.
The only way to be sure, in fact, is to use calculus (if possible) to
prove there are no other local minima.
But there probably ARE, so... good luck!
Bobby
On Fri, 30 Nov 2007 04:17:03 -0600, dantimatter <google at dantimatter.com>
wrote:
>
> 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
>
>
>
-- =
DrMajorBob at bigfoot.com