# Efficient multidimensional optimization while constraining some coordinates

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


• Did I understand this correctly? You want to solve only quadratic problems (with sparse Hessian) on affine subspaces? Then rephrasing this into a linear system would probably be the most efficient approach. I can do that for you, but I first want to be sure that this is not only a toy example. Oct 18 at 11:07
• Instead of using equality constraints, why not just set those coordinates in the objective? e.g instead of FindMinimum[{f[x,y,z], x==1},...] , try FindMinimum[f[1,y,z],...] You could achieve that with a replacement while constructing the objective. Oct 18 at 11:17
• Oh, I just realized that the problem is actually not quadratic in the vertex coordinates (because the unsquared Norm occurs in the objective)... Oct 18 at 11:19
• @flinty Originally I did not consider this type of solution because I only have a single, vector-valued x variable. If I used a separate symbol for each component of x, it would kill performance. But perhaps you are right, and I should just map a "constrained version" (i.e. values already substituted) of x into the original x, using some helper functions. This is exactly what Henrik seems to be doing. Oct 18 at 13:42
• @CElliott This is new to me. I assume that this idea might come from some specific pseudorandom number generation methods requiring this (see e.g. here, but even than it would only matter if the modulus is chosen as non-prime). I assume that if Mathematica required this for some of its random number generation methods, this would be clearly documented. But there is no mention of this, and the typical example seed used in the Mathematica documentation is 1234. Oct 21 at 15:01

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"]},
InternalWithLocalSettings[
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[CompileGetElement[P, CompileGetElement[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 (CompileGetElement[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.

• Thanks so much, this is great! And thanks for pointing to your SemiSmoothNewton function. Oct 18 at 13:45
• You're welcome, Szabolcs! Oct 19 at 7:02
• @HenrikSchumacher It is very nice solution (+1). I will try to use it with my algorithm for PDEs and integral equations. Have you any publication with SemiSmoothNewton function discussion? Nov 16 at 6:07
• Hi Alex! I learnt all what I know about the algorithm from chapter 2 in these lecture notes by Michael Hintermüller. Nov 16 at 6:19
• @HenrikSchumacher Ok! Thank you very much. Nov 16 at 16:30