# Minimizing a function of many coordinates

I'm trying to obtain a series of points on the unit sphere with a somewhat homogeneous distribution, by minimizing a function depending on distances (I took $\exp(-d)$). My points are represented by spherical angles $\theta$ and $\phi$, starting by choosing equidistributed random vectors:

pts = Apply[{2 π #1, ArcCos[2 #2 - 1]} &, RandomReal[1, {100, 2}], 1];


The energy function is defined first:

energy[p_] := Module[{cart},
cart = Apply[{Sin[#1]*Cos[#2], Sin[#1]*Sin[#2], Cos[#1]} &, p, 1];
Total[Outer[Exp[-Norm[#1 - #2]] &, cart, cart, 1], 2]
]


But now, I can’t manage to get the right routine for minimization. I tried FindMinimum, which does local minimization from a given starting point, which is what I want. But it should operate on an expression of literal variables, so I'm kind of screwed:

FindMinimum[energy[p], {p, pts}]


Outer::normal: Nonatomic expression expected at position 2 in Outer[Exp[-Norm[#1-Slot[<<1>>]]]&,p,p,1]. >>
FindMinimum::nrnum: The function value […] is not a real number at {p} = […] >>


The above obviously doesn't work, but I don't think it's wise to introduce a series of 200 literal variables. There has to be another way, hasn't it? Or is there an efficient way of introducing a lot of variables?

-
I upvoted this because it seems to be a good example for a case where the usual advice to define the function only for numeric variables (_?NumericQ) is not so good. By allowing symbolic arguments in the function to be minimized, you also allow FindMinimum to avoid taking numerical derivatives. At least that's how I approached the answer. –  Jens Apr 8 '12 at 2:15

To get arbitrarily many formal variables, you can use Array. But with those variables, your function definition won't work because of the Apply statement. So I modified your definition as follows (I reduced the point number for testing purposes):

pts = Apply[{ArcCos[2 #2 - 1], 2 \[Pi] #1} &, RandomReal[1, {10, 2}], 1];
Clear[energy];
Clear[a];
vars = Array[a, {Length[pts], 2}];
energy[p_] :=
Module[{cart},
cart = Map[{Sin[#[[1]]]*Cos[#[[2]]], Sin[#[[1]]]*Sin[#[[2]]],
Cos[#[[1]]]} &, p];
Total[Outer[Exp[-Norm[#1 - #2]] &, cart, cart, 1], 2]];
FindMinimum[energy[vars], Transpose[{Flatten@vars, Flatten@pts}]]


{32.2548, {a[1, 1] -> 1.93787, a[1, 2] -> 1.72361, a[2, 1] -> 1.11355, a[2, 2] -> 0.893035, a[3, 1] -> 6.21077, a[3, 2] -> 2.1405, a[4, 1] -> 3.06917, a[4, 2] -> 2.14062, a[5, 1] -> 1.06997, a[5, 2] -> 2.50937, a[6, 1] -> 4.21367, a[6, 2] -> 1.69561, a[7, 1] -> 5.07748, a[7, 2] -> 2.48594, a[8, 1] -> 4.31041, a[8, 2] -> 0.111206, a[9, 1] -> 4.25016, a[9, 2] -> 3.31368, a[10, 1] -> 5.11923, a[10, 2] -> 0.955784}}

The form of the array passed to energy matches the $N\times2$ dimension list that is expected by the line creating the cart variable. In the FindMinimum statement the dummy variables and initial conditions are specified as a single list of pairs by using Flatten on both.

There is the usual wrinkle that the minimization may need to be tweaked for precision, but that's a different issue.

Finally, to get the minimizing point list, you have to do

vars/.Last[%]


Edit

Depending on the function to be optimized, it's sometimes faster to avoid the use of derivatives by specifying the initial conditions for FindMinimum in the form of three numbers:

FindMinimum[energy[vars],
Transpose[{Flatten@vars, Flatten@pts, Flatten@pts - .1, Flatten@pts + .1}]]


Edit 2

I did get a significant speed-up with this for your example, but the performance depends on the random starting points (and on the choice of bracket width) so I can't say anything definitive. That seems like a topic for a different question.

Edit 3

Though I didn't look at the speed issue in detail, forcing FindMinimum to work with numerical derivatives may be the worst option here. That will happen if you define your function energy only for numerical arguments, as in

energy[p : {{_?NumericQ, _?NumericQ} ..}] :=


followed by either your own or my initial definition above. So although that's a common advice people give in these applications, it is not going to be the fastest approach here.

Edit 4

I just had another idea on how to improve the speed of my solution: the use of Norm might make it harder to estimate the Hessian for this function. And indeed, when I got rid of Norm there was a significant speed gain (note that the initial solution above is already faster than the _?NumericQ approach even when the latter is compiled while mine is not). I think this is worth adding here because Norm seems like a natural thing to use in pair potentials, even if the energy expression becomes more complicated than the one in this question.

So here is the new version, with Norm replaced by Sqrt[(#1 - #2).(#1 - #2)]. Observe that I have now put back the original particle number of 100 because on my laptop this takes less than 8 seconds to evaluate!

pts = Apply[{ArcCos[2 #2 - 1], 2 \[Pi] #1} &, RandomReal[1, {100, 2}], 1];
Clear[energy];
Clear[a];
vars = Array[a, {Length[pts], 2}];
energy[p_] :=
Module[{cart},
cart = Map[{Sin[#[[1]]]*Cos[#[[2]]], Sin[#[[1]]]*Sin[#[[2]]],
Cos[#[[1]]]} &, p];
Total[Outer[Exp[-Sqrt[(#1 - #2).(#1 - #2)]] &, cart, cart, 1], 2]];
FindMinimum[energy[vars], Transpose[{Flatten@vars, Flatten@pts}]]


Oh, and one more thing I changed (unrelated to the question), is to switch your definitions of polar and azimuthal angles in pts.

-
Thanks Jens, it works indeed. It's horribly slow, though (takes a couple of minutes for 100 points on my laptop). I wonder what alternatives exist to FindMinimum for large number of variables… –  F'x Apr 7 '12 at 15:55
Yes, that's because of the Outer in your function; so you asked for it... –  Jens Apr 7 '12 at 15:59
Well, the energy calculation is $N^2$, yes, but I'm sure the minimization itself is not efficient. I have a C program at work that does that in much less time… and because it's operating on decent-sized data sets, I expected Mathematica to be good at it. –  F'x Apr 7 '12 at 16:11
Yes, this is a general-purpose routine, and one can probably figure out some way to speed things up. There's also a commercial package called KNITRO: wolfram.com/products/applications/knitro –  Jens Apr 7 '12 at 16:37
Might try playing with the Method option to FindMinimum. Also if you modify those energies to explicit squares then "LevenbergMarquardt" could be used. Instead of Exp[-Sqrt[#.#]] each term would be something like 1/#.#. or maximize the sum of squares (instead of minimizing sum of reciprocals thereof). Then either LM or "QuadraticProgramming" can be used. Or use the L_1 distances and LinearProgramming. –  Daniel Lichtblau Apr 24 '12 at 23:34

I played around with a few plausible variants on the energy. In all cases I compared the result using the original energy function. Some things I learned:

(1) Some variants will tend to give results that do quite well when gauged via the original energy.

(2) Others (not shown below) will do poorly because they weigh the far values too heavily. This, alas, means we cannot easily use the GaussNewton (LevenbergMarquardt) method, since it is the operations of squaring that hurts us. Well, maybe there are ways around this.

(3) Summing over only distinct pairs rather than all pairs cuts the time in half. I will speculate that the bulk of time is spent in evaluating derivatives and not in the function evaluations themselves, as I am fairly certain Total[Outer[...]] will beat Sum even when the latter only need account for half or so as many pairs.

(4) For one energy variant I got a modest speed improvement using the ConjugateGradient method.

(5) Scaling appears to be quadratic in the number of points (no huge surprise, I guess).

(6) We can handle 200 points in 24 seconds on my desktop machine.

In[254]:=
pts = Apply[{ArcCos[2 #2 - 1], 2 \[Pi] #1} &, RandomReal[1, {100, 2}],
1];
Clear[a];
vars = Array[a, {Length[pts], 2}];


Here is the basic case.

In[292]:=
energy[p_] :=
Module[{cart},
cart = Map[{Sin[#[[1]]]*Cos[#[[2]]], Sin[#[[1]]]*Sin[#[[2]]],
Cos[#[[1]]]} &, p];
Total[Outer[Exp[-Sqrt[(#1 - #2).(#1 - #2)]] &, cart, cart, 1], 2]];

In[293]:=
t = Timing[{min, vals} =
Quiet[FindMinimum[energy[vars],
Transpose[{Flatten@vars, Flatten@pts}],
MaxIterations -> 1000]];];
{t, min}

Out[294]= {{14.1, Null}, 2978.01}


We use Sum over distinct pairs from here onward.

In[295]:=
energy2[p_] :=
Module[{cart},
cart = Map[{Sin[#[[1]]]*Cos[#[[2]]], Sin[#[[1]]]*Sin[#[[2]]],
Cos[#[[1]]]} &, p];
Sum[Exp[-Sqrt[(cart[[j]] - cart[[k]]).(cart[[j]] -
cart[[k]])]], {j, Length[p] - 1}, {k, j + 1, Length[p]}]];

In[296]:=
t2 = Timing[{min2, vals2} =
Quiet[FindMinimum[energy2[vars],
Transpose[{Flatten@vars, Flatten@pts}],
MaxIterations -> 1000]];];
{t2, min2, energy[vars /. vals2]}

Out[297]= {{6.58, Null}, 1439., 2978.01}


Minimize the sum of reciprocals of the pairwise distances squared.

In[298]:= energy3[p_] := Module[{cart},
cart = Map[{Sin[#[[1]]]*Cos[#[[2]]], Sin[#[[1]]]*Sin[#[[2]]],
Cos[#[[1]]]} &, p];
Sum[1/((cart[[j]] - cart[[k]]).(cart[[j]] - cart[[k]])), {j,
Length[p] - 1}, {k, j + 1, Length[p]}]
]

In[299]:=
t3 = Timing[{min3, vals3} =
Quiet[FindMinimum[energy2[vars],
Transpose[{Flatten@vars, Flatten@pts}],
MaxIterations -> 1000]];];
{t3, min3, energy[vars /. vals3]}

Out[300]= {{6.72, Null}, 1439., 2978.01}


This variant on energy happened to get a bit faster using a nondefault method setting.

In[301]:=
t3b = Timing[{min3b, vals3b} =
Quiet[FindMinimum[energy3[vars],
Transpose[{Flatten@vars, Flatten@pts}], MaxIterations -> 1000,
{t3b, min3b, energy[vars /. vals3b]}

Out[302]= {{5.23, Null}, 5340.65, 2978.01}


Maximize sum of distances. I will mention that using the sum of squares, which i would prefer to do, fails to give a useful result. It puts half the points one place and the other half at the polar opposite, I believe. That comes from the further distances getting relatively more weight in the objective function.

In[304]:= energy4[p_] := Module[{cart},
cart = Map[{Sin[#[[1]]]*Cos[#[[2]]], Sin[#[[1]]]*Sin[#[[2]]],
Cos[#[[1]]]} &, p];
Sum[Sqrt[((cart[[j]] - cart[[k]]).(cart[[j]] - cart[[k]]))], {j,
Length[p] - 1}, {k, j + 1, Length[p]}]
]

In[305]:=
t4 = Timing[{min4, vals4} =
Quiet[FindMaximum[energy4[vars],
Transpose[{Flatten@vars, Flatten@pts}],
MaxIterations -> 1000]];];
{t4, min4, energy[vars /. vals4]}

Out[306]= {{8.44, Null}, 6662.64, 2978.}


Similar to a couple of tries above, but with distances instead of squared distances.

In[308]:= energy5[p_] := Module[{cart},
cart = Map[{Sin[#[[1]]]*Cos[#[[2]]], Sin[#[[1]]]*Sin[#[[2]]],
Cos[#[[1]]]} &, p];
Sum[1/((cart[[j]] - cart[[k]]).(cart[[j]] - cart[[k]]))^(1/2), {j,
Length[p] - 1}, {k, j + 1, Length[p]}]
]

In[309]:=
t5 = Timing[{min5, vals5} =
Quiet[FindMinimum[energy5[vars],
Transpose[{Flatten@vars, Flatten@pts}],
MaxIterations -> 1000]];];
{t5, min5, energy[vars /. vals5]}

Out[310]= {{6.04, Null}, 4448.45, 2978.01}


Notice that all of these agreed fairly closely to six places in terms of the original energy function.

Now we'll go to 200 points and use the fastest variant from above.

In[319]:=
pts200 = Apply[{ArcCos[2 #2 - 1], 2 \[Pi] #1} &,
RandomReal[1, {200, 2}], 1];
vars200 = Array[a, {Length[pts200], 2}];
t200 = Timing[{min200, vals200} =
Quiet[FindMinimum[energy3[vars200],
Transpose[{Flatten@vars200, Flatten@pts200}],
MaxIterations -> 1000,
{t200, min200, energy[vars200 /. vals200]}

Out[322]= {{23.59, Null}, 24816.3, 11891.3}


Here I crib Mark McClure's code to show both the original points and the result of the optimization. The pictures will have to speak for themselves because i'm not going to speak for them.

In[323]:=
pts3D = {Sin[#[[1]]]*Cos[#[[2]]], Sin[#[[1]]]*Sin[#[[2]]],
Cos[#[[1]]]} & /@ pts200;
Graphics3D[Point[pts3D]]


In[325] := pts3Db = {Sin[#[[1]]]*Cos[#[[2]]], Sin[#[[1]]]*Sin[#[[2]]],
Cos[#[[1]]]} & /@ (vars200 /. vals200);
Graphics3D[Point[pts3Db]]


-

First, I believe that you have switched the roles of $\phi$ and $\theta$ in your first definition. Thus, a slight edit of your code yields the following

pts = Apply[{ArcCos[2 #2 - 1], 2 #1*Pi} &,
RandomReal[1, {10000, 2}], 1];
pts3D = {Sin[#[[1]]]*Cos[#[[2]]], Sin[#[[1]]]*Sin[#[[2]]],
Cos[#[[1]]]} & /@ pts;
Graphics3D[Point[pts3D]]


Now, pts3D is already a nice uniform distribution on the sphere. Is this all you want??

If you do need to go through the energy minimization, then you can use FindMinimum since it does allow the variables to be tensors, and your original FindMinimum command is fine. The problem is that you haven't restricted your definition of energy to work only with numeric values. Thus, you can minimize the energy (with compilation for speed) like so:

energy1[p:{{_?NumericQ,_?NumericQ}..}] := Module[{cart},
cart = Apply[{Sin[#1]*Cos[#2], Sin[#1]*Sin[#2], Cos[#1]} &, p, 1];
Total[Outer[Exp[-Norm[#1 - #2]] &, cart, cart, 1], 2]
];
cEnergy2 = Compile[{{p,_Real,2}},Module[{cart},
cart = Map[{Sin[#[[1]]]*Cos[#[[2]]],
Sin[#[[1]]]*Sin[#[[2]]],Cos[#[[1]]]} &, p, 1];
Sum[Exp[-Sqrt[(u-v).(u-v)]],{u,cart},{v,cart}]
], CompilationTarget -> "C", RuntimeOptions -> "Speed"];
energy2[p:{{_?NumericQ,_?NumericQ}..}] := cEnergy2[p];

SeedRandom[1];
pts=Apply[{ArcCos[2 #2-1], 2#1*Pi}&,RandomReal[1,{20,2}],1];
FindMinimum[energy2[p],{p,pts}]//AbsoluteTiming
FindMinimum[energy1[p],{p,pts}]//AbsoluteTiming


I think this produces the result that you want.

Note that I've used compilation in to speed up the code by a factor of nearly twenty on my machine. But the time complexity is such that even 100 points is out of reach.

-
Yes, his original function works fine with the pattern energy[p : {{_?NumericQ, _?NumericQ} ..}]. However, my initial solution is already 10 times faster than using his original function with the NumericQ trick for a 10 point array. –  Jens Apr 7 '12 at 18:05
Your updated solution in the compiled form is still slower than my uncompiled solution. –  Jens Apr 8 '12 at 2:09

I was planning to get to this later but I am not sure I'll have time today. To force it to work I can do

pts = Apply[{2 \[Pi] #1, ArcCos[2 #2 - 1]} &, RandomReal[1, {10, 2}],
1];
energy[p_] :=
Module[{cart},
cart = Apply[{Sin[#1]*Cos[#2], Sin[#1]*Sin[#2], Cos[#1]} &, p, 1];
Total[Outer[Exp[-Norm[#1 - #2]] &, cart, cart, 1], 2]]

fvars = Table[{x[i], y[i]}, {i, 1, Length@pts}];

NMinimize[energy[fvars], Flatten@fvars]

(*{32.2548, {x[1] \[Rule] -2.0039, y[1] \[Rule] -0.778507,
x[2] \[Rule] -1.13768, y[2] \[Rule] -1.5639, x[3] \[Rule] 1.13769,
y[3] \[Rule] 0.00689281, x[4] \[Rule] 2.00391,
y[4] \[Rule] -0.778497, x[5] \[Rule] -2.00391,
y[5] \[Rule] 0.792293, x[6] \[Rule] -9.32847\[Times]10^-6,
y[6] \[Rule] 1.20086, x[7] \[Rule] 1.1377, y[7] \[Rule] -1.5639,
x[8] \[Rule] 3.14158, y[8] \[Rule] 2.10272, x[9] \[Rule] 2.0039,
y[9] \[Rule] 0.792295, x[10] \[Rule] -1.13769,
y[10] \[Rule] 0.00689875}}
*)


But I wouldn't do it this way.

-
When I get time I'll try to make this work faster. –  acl Apr 7 '12 at 16:01