MathGroup Archive 2007

[Date Index] [Thread Index] [Author Index]

Search the Archive

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

  • To: mathgroup at smc.vnet.net
  • Subject: [mg83754] Re: [mg83642] optimization routine using conjugate gradient (or other) method?
  • From: DrMajorBob <drmajorbob at bigfoot.com>
  • Date: Thu, 29 Nov 2007 06:40:32 -0500 (EST)
  • References: <9889841.1196176137975.JavaMail.root@m35>
  • Reply-to: drmajorbob at bigfoot.com

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 <dantimatter 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
>
>



-- =

DrMajorBob at bigfoot.com


  • Prev by Date: Re: Re: Possible bug in HamiltonianCycle
  • Next by Date: Re: FindInstance puzzler
  • Previous by thread: Re: optimization routine using conjugate gradient (or other) method?
  • Next by thread: Re: optimization routine using conjugate gradient (or other) method?