What is the best way to minimize a function f[x_] subject to the constraint that x is one of several discrete values using Mathematica?

For example, let x be one of the first five primes

valueConstraint = x \[Element] Prime /@ Range[5]

The above does not actually work with Minimize, so instead, let us do the following:

fOneOf[x_, y_List] := Or @@ (x == # & /@ y); (* {x == y1 || x == y2 || ...} *)
variables = Table[Symbol["x$" <> ToString[i]], {i, 5}]
values = Prime[#] & /@ Range[5]
valueConstraint = fOneOf[#, values] & /@ variables
f[{x__Integer}] := Total[{x}]
NMinimize[{f[variables], valueConstraint}, variables \[Element] Integers] // AbsoluteTiming
(*20 seconds, correct answer*)

I think this is the most natural way to do it, but the solver is exceedingly slow with discrete equality constraints,

Let's try something else...

Rather than saying x is equal to one of several values, we can instead make x an integer index for an array of possible values. In this case, we have an integer optimization problem with x being a member of some range of valid indices:

variables = Table[Symbol["x$" <> ToString[i]], {i, 5}]
values = Prime[#] & /@ Range[5]
valueConstraint = 1 <= variables <= Length[values]
f[{x__Integer}] := Total[Part[values, {x}]]
NMinimize[{f[variables], valueConstraint}, variables \[Element] Integers] // AbsoluteTiming
(* 1 second, correct answer *)

And I can now minimize this f with 100 variables and 100 possible values in about 60s.

This is much faster, but using Part[] seems like an extremely inelegant way to do this, especially if the function is complicated and requires replacing each variable in the original expression with Part[values, x] where x is now an integer array index.

Is there some better way to tell the numerical solver that 'x' can be one of several values? Or is there a more elegant optimization technique to use for such a constrained task?

Edit: To help find a better solution, it might be worth mentioning the binary case (found here). Another issue with naively using Part[] is probably that it destroys gradient information or the potential for symbolic analysis, which is why the pattern f[{x__Integer}] is used. Technically this is not needed in the equality constraint version above, but removing it does not fix the terrible performance of NMinimize with equality constraints. It is quite possible that there is simply no way to recover gradient information with finite-set-valued variables and NMinimize performance cannot be significantly improved over the Part[] implementation. I was hoping there was something we could do with matrices to maintain gradient information and nice symbolic evaluation, but I cannot think of anything right now.

variables = Array[x, 100]
valueConstraint = Thread[0 <= variables <= 1]
f[x_] := Total[x]
NMinimize[{f[variables], valueConstraint}, variables \[Element] Integers] // AbsoluteTiming
(* 0.004 seconds *)
  • $\begingroup$ The 2nd code (so the one below "The above does not actually work ..) is missing a bracket-> [...NMinimize(* here *){]. After adding it, I received an error : "NMinimize{" ""Options expected (instead of 2) beyond position 2 in \ Total[2,2,2,2,2]. An option must be a rule or a list of rules"". The code that followed worked. $\endgroup$ Aug 15, 2022 at 2:54
  • 1
    $\begingroup$ you can also use Unique[ConstantArray["x",5]] for the variable names which is a bit quicker to type in this case but each time you evaluate that code the names change. $\endgroup$ Aug 15, 2022 at 3:19
  • 1
    $\begingroup$ In the second code I think the error was because it should be Total[{x}] rather than Total[x] $\endgroup$ Aug 15, 2022 at 3:51
  • 2
    $\begingroup$ Have you seen: library.wolfram.com/infocenter/Conferences/4317 $\endgroup$ Aug 15, 2022 at 16:02
  • 1
    $\begingroup$ I saw this post mathematica.stackexchange.com/q/277948/86543 that seems related in particular the advice in the comments to consider the resource function AntColonyOptimization resources.wolframcloud.com/FunctionRepository/resources/…. I also saw something about using neural networks for discrete optimization on another site but I do not remember much anymore. $\endgroup$ Jan 3, 2023 at 13:55

1 Answer 1


For complex conditions constraints and function f,maybe use Pick or Select and use Ordering to locate the min-value index.

Clear[reg, constraints, f, constraintsReg, indexs];
reg = Tuples[{2, 3, 5, 7, 11, 13}, 6];
constraints = {x1, x2, x3, x4, x5, x6} |-> 
   And @@ {x1 + x2 + x3 - x4 - x5 - x6 > 5, x5^2 < x1 + x2};
f[x1_, x2_, x3_, x4_, x5_, x6_] = x1^2 + x2^2 - x3 - x4 - x5 - x6;
constraintsReg = Pick[reg, constraints @@@ reg];
indexs = Ordering[f @@@ constraintsReg, 1]
{constraintsReg[[indexs]], f @@@ constraintsReg[[indexs]]}

{{{2, 3, 13, 3, 2, 7}}, {-12}}

  • $\begingroup$ This will fail if you have 100 possible values unless you have a lot of RAM $\endgroup$ Aug 17, 2022 at 3:27
  • $\begingroup$ @AlecGraves this answer uses a procedural program that tries all cases while using less memory $\endgroup$ Aug 18, 2022 at 11:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.