Efficient multidimensional optimization while constraining some coordinates

Question

I have a function f[x] with a high-dimensional vector input and a scalar output. I am looking to minimize this function while constraining some of the components of x. What is an efficient way of doing this? My key concern is efficiency.

---

A trivial example

FindMinimum[{Norm[x], Indexed[x, 1] == 1.0}, {x, {1,2,3}}]

Here x is a 3D vector and I constrained its first component to be 1.0. This certainly works, but when using Indexed to constrain certain components, the performance becomes very bad. Is there a way to constrain some components of the vector without sacrificing performance?

---

A more realistic example

The following is a more involved example that better represents my actual use case, and can be used to test performance. We want to find the minimal energy state of a system of interconnected springs.

Let us generate a random system of springs:

n = 20;

SeedRandom[1234];
pts = Join[
   RandomPoint[Disk[], n],
   CirclePoints[n]
   ];

mesh = DelaunayMesh[pts]

The edges of this mesh will be our springs. We will take the rest lengths of the springs from the above-plotted configuration.

restLengths = PropertyValue[{mesh, 1}, MeshCellMeasure];

edges = MeshCells[mesh, 1][[All, 1]];

For simplicity, the energy of a spring is the squared difference between its actual length and rest length.

Clear[energy]
energy[coords_ /; MatrixQ[coords, NumericQ]] :=
 
With[{diff = Subtract @@ (coords[[#]] &) /@ Transpose[edges]},
  Total[(Norm /@ diff - restLengths)^2]
]

Update: Here I meant to write Norm[#]^2& instead of Norm. I will keep the example as-is in order not to pull out the rug from under the already posted answer.

Let us now make the system 3-dimensional and choose our initial configuration:

coords = Append[#, 0.] & /@ MeshCoordinates[mesh];

With these coordinates, the energy is zero:

energy[coords] // Chop
(* 0 *)

Let us select the point closest to the centre, and pull it out of the initial plane by setting a non-zero $z$ coordinate:

v = First@Nearest[coords -> "Index", {0, 0, 0}];

coords[[v, 3]] = 0.3;

Now let us choose which points we keep fixed. Here, these are the points on the perimeter, as well as the point that we pulled out of the plane.

fixed = Append[
   Range[n + 1, 2 n],
   v
   ];

Due to pulling that single point, the energy is no longer zero:

energy[coords]
(* 0.0509467 *)

We are now almost ready to minimize the energy function, but first we must construct the appropiate constraints for FindMinimum, in order to keep some of the points fixed during minimization.

constraints = Table[
   Indexed[x, i] == coords[[i]],
   {i, fixed}
   ];

Here's the problem: Performance is really bad with constraints applied. Compare minimization with and without constraints:

In[187]:= 
result = 
   FindMinimum[{energy[x], constraints}, {x, coords}, 
    MaxIterations -> 5]; // AbsoluteTiming

During evaluation of In[187]:= FindMinimum::cvmit: Failed to converge to the requested accuracy or precision within 5 iterations.

Out[187]= {15.676, Null}

In[189]:= 
In[190]:= 
result = 
   FindMinimum[energy[x], {x, coords}, 
    MaxIterations -> 5]; // AbsoluteTiming

During evaluation of In[190]:= FindMinimum::lstol: The line search decreased the step size to within the tolerance specified by AccuracyGoal and PrecisionGoal but was unable to find a sufficient decrease in the function. You may need more than MachinePrecision digits of working precision to meet these tolerances.

Out[190]= {0.361019, Null}

Question: How can I fix some of the coordinates without compromising performance? Is there another way than using Indexed?

Note that I need the flexibility to choose which coordinates to constrain (in this case: which points to fix). This is why I do not want to implement energy in such a way that its input would consist only of non-fixed coordinates.

In my real use case, the function I am minimizing is a compiled "black box" for better evaluation performance. I would hope to work with an x that has a dimension on the order of 100–1000.

BTW, for the above example, this is an easy way to plot the final configuration:

Graphics3D@GraphicsComplex[x /. result[[2]], Line[edges]]

Answer 1

Phew, this became more complex than I thought because the objective is not quadratic (as I had first expected). Here some preparations for an efficient evaluation of the energy and its first two derivatives:

MySparseArray[X_, r_, f_: Total, OptionsPattern[{"Background" -> 0.}]] := 
  If[(Head[X] === Rule) && (X[[1]] === {}), X[[2]], 
   With[{spopt = SystemOptions["SparseArrayOptions"]}, 
    Internal`WithLocalSettings[
     SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> f}], 
     SparseArray[X, r, OptionValue["Background"]], 
     SetSystemOptions[spopt]]]];

ClearAll[cEnergy];
cEnergy[k_Integer] := cEnergy[k] = Module[{energy, PP, P, L, e},
    PP = Table[Compile`GetElement[P, Compile`GetElement[e, i], j], {i, 1, 2}, {j, 1, 3}];
    energy = (Sqrt[Dot[PP[[2]] - PP[[1]], PP[[2]] - PP[[1]]]] - L)^2;
    Print["Compiling cEnergy[", k, "]..."];
    With[{code = D[energy, {Flatten[PP], k}]}, 
     Compile[{{P, _Real, 2}, {e, _Integer, 1}, {L, _Real}}, code, 
      CompilationTarget -> "C", RuntimeAttributes -> {Listable}, 
      Parallelization -> True, RuntimeOptions -> "Speed"]]];

ClearAll[cDerivativePattern];
cDerivativePattern[tuplesize_Integer, dim_Integer, derivative_Integer] := 
  cDerivativePattern[tuplesize, dim, derivative] = Block[{t},
    Print["Compiling cDerivativePattern[", tuplesize, ",", dim, ",", derivative, "]..."]; 
    With[{code = 
       Tuples[Flatten[ Table[{dim (Compile`GetElement[t, i] - 1) + j}, {i, 1, tuplesize}, {j, 1, dim}]], derivative]}, 
     Compile[{{t, _Integer, 1}}, code, CompilationTarget -> "C", 
      RuntimeAttributes -> {Listable}, Parallelization -> True, 
      RuntimeOptions -> "Speed"]]];

For another energy, simply edit the line

energy = (Sqrt[Dot[PP[[2]] - PP[[1]], PP[[2]] - PP[[1]]]] - L)^2;

to whatever suits you.

Now for a fixed mesh, you can execute this once as preparation:

n = 2000;
d = 3;
SeedRandom[1234];
pts = Join[RandomPoint[Disk[], n], CirclePoints[n]];
mesh = DelaunayMesh[pts];

restLengths = PropertyValue[{mesh, 1}, MeshCellMeasure];
edges = MeshCells[mesh, 1][[All, 1]];

ClearAll[pat];
pat[k_] := pat[k] = Join @@ cDerivativePattern[2, 3, k][edges];

m = MeshCellCount[mesh, 0];
(*Objective*)
ClearAll[F]
F[X_?VectorQ] := Total[cEnergy[0][Partition[X, d], edges, restLengths]];
Derivative[k_Integer][F][X_?VectorQ] := MySparseArray[
 pat[k] -> Flatten[cEnergy[k][Partition[X, d], edges, restLengths]],
 ConstantArray[d m, k],
 Total
 ];

Then, for every new set of fixed coordinates, execute this:

coords = Append[#, 0.] & /@ MeshCoordinates[mesh];
v = First@Nearest[coords -> "Index", {0, 0, 0}];

coords[[v, 3]] = 0.3;
fixed = Append[Range[n + 1, 2 n], v];
fixeddofs = Join @@ Transpose[Table[d (fixed - 1) + k, {k, 1, d}]];
free = Complement[Range[1, m], fixed];
freedofs = Join @@ Transpose[Table[d (free - 1) + k, {k, 1, d}]];

X0 = Flatten[coords];
target = X0[[fixeddofs]];


(*Equality constraint*)

constraintmatrix = KroneckerProduct[
   SparseArray[
    Transpose[{Range[Length[fixed]], fixed}] -> 1.,
    {Length[fixed], m}
    ],
   IdentityMatrix[d]
   ];

\[CapitalPhi][X_?VectorQ] := X[[fixeddofs]] - target;
\[CapitalPhi]'[X_?VectorQ] = constraintmatrix;
\[CapitalPhi]''[X_?VectorQ] = SparseArray[{}, {Length[fixeddofs], d m, d m}];

Then you can run my implementation SemiSmoothNewton of Newton's method with line search from here:

result = SemiSmoothNewton[X0, F, \[CapitalPhi], None];
result // Dataset
X = result["Solution"];

It took about 2 seconds for $n = 2000$.

And finally, you can visualize the resulting surface as follows:

MeshRegion[Partition[X, d], MeshCells[mesh, 2, "Multicells" -> True]]

Because the constraints are linear and satisfied for X0, we can simply cast this into an unconstrained problem by "forgetting" the fixed degrees of freedom:

free = Complement[Range[1, m], fixed];
freedofs = Join @@ Transpose[Table[d (free - 1) + k, {k, 1, d}]];
X0 = Flatten[coords];
X = X0;
x0 = Flatten[X][[freedofs]];

f[x_?VectorQ] := (X[[freedofs]] = x; F[X]);
f'[x_?VectorQ] := (X[[freedofs]] = x; F'[X][[freedofs]]);
f''[x_?VectorQ] := (X[[freedofs]] = x; F''[X][[freedofs, freedofs]]);

result = SemiSmoothNewton[x0, f, None, None];
result // Dataset
X[[freedofs]] = result["Solution"];

The linear system in Newton's method are now considerably smaller, so this is more efficient. For $n = 2000$, this took about 1.8 seconds.