# Using Lagrange multipliers to minimize function of six variables and an inequality

I want to form the function

h = f − λ1 g1 − λ2 g2 − λ3 g3 − λ4 g4


where f is the function to minimize subject to the identities

g1 = n1 + n2 + n3 - 3
g2 = n1 + 2 n2 + n5 - 4
g3 = 4n3 + 2n4 + 2n5 - 8
g4 = 2n6 - 2


so that I can solve for the first partial derivatives with respect to n1, n2, n3, n4, n5, n6, λ1, λ2, λ3, λ4.

Also, the solutions for n1 to n6 must be real numbers satisfying

0 ≤ n1, n2, n3, n4, n5, n6 ≤  n1 + n2 + n3 + n4 + n5 + n6


My code so far:

f[n1_, n2_, n3_, n4_, n5_, n6_] :=
-78416.91961691371  n1 + 8920.922 n1 Log[n1/(n1 + n2 + n3 + n4 + n5 + n6)] -
293404.0587042522  n2 + 8920.922 n2 Log[n2/(n1 + n2 + n3 + n4 + n5 + n6)] -
36099.19105260982  n3 + 8920.922 n3 Log[n3/(n1 + n2 + n3 + n4 + n5 + n6)]-
53730.85220261606  n4 + 8920.922 n4 Log[n4/(n1 + n2 + n3 + n4 + n5 + n6)]-
261977.41085974447 n5 + 8920.922 n5 Log[n5/(n1 + n2 + n3 + n4 + n5 + n6)] -
117715.30870437803 n6 + 8920.922 n6 Log[n6/(n1 + n2 + n3 + n4 + n5 + n6)]

g1[n1_, n2_, n3_] := n1 + n2 + n3 - 3
g2[n1_, n2_, n5_] := n1 + 2 n2 + n5 - 4
g3[n3_, n4_, n5_] := 4 n3 + 2 n4 + 2 n5 - 8
g4[n6_] := 2 n6 - 2

h[n1_, n2_, n3_, n4_, n5_, n6_, λ1_, λ2_, λ3_, λ4_] :=
f[n1, n2, n3, n4, n5, n6] -
λ1 g1[n1, n2, n3] - λ2 g2[n1, n2, n5] - λ3 g3[n3, n4, n5] - λ4 g4[n6]

Column[
pts = {n1, n2, n3, n4, n5, n6} /.
FullSimplify @
Solve[
D[h[n1, n2, n3, n4, n5, n6, λ1, λ2, λ3, λ4], #] == 0 & /@
{n1, n2, n3, n4, n5, n6, λ1, λ2, λ3, λ4},
{n1, n2, n3, n4, n5, n6, λ1, λ2, λ3, λ4}],
Frame -> All]]


I have based my code off the example in this question.

Whenever I evaluate the code, it says'Running' for a while, but does not return any results. However, I do not know why this is occuring. Also, I do not know how to code the inequality

0 ≤ n1, n2, n3, n4, n5, n6 ≤  n1 + n2 + n3 + n4 + n5 + n6


Any help would be greatly appreciated.

Edit:

The identities in this question are now different.

g1 = n1 + n2 + n3 - 432
g2 = n1 + 2 n2 + n5 - 259
g3 = 4n3 + 2n4 + 2n5 - 628
g4 = 2n6 - 4


I am using this line of code to evaluate the problem:

vars = {n1, n2, n3, n4, n5,
n6, \[Lambda]1, \[Lambda]2, \[Lambda]3, \[Lambda]4};
eq = D[h @@ vars == 0, {vars, 1}];
Column[pts = {n1, n2, n3, n4, n5, n6} /.
FindRoot[eq, Evaluate[Transpose[{vars, 1 + 0 vars}]]],
Frame -> All]], 4]


however when i evaluate it, i get this message:

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


How can i overcome this?

Another edit: I made a mistake when defining the function f. It should be as follows.

 f[n1_, n2_, n3_, n4_, n5_, n6_] :=
-78416.91961691371  n1 + 8920.922 n1 Log[n1/(n1 + n2 + n3 + n4 + n5 + n6)] -
293404.0587042522  n2 + 8920.922 n2 Log[n2/(n1 + n2 + n3 + n4 + n5 + n6)] -
36099.19105260982  n3 + 8920.922 n3 Log[n3/(n1 + n2 + n3 + n4 + n5 + n6)]-
261977.41085974447 n4 + 8920.922 n4 Log[n4/(n1 + n2 + n3 + n4 + n5 + n6)]-
53730.85220261606 n5 + 8920.922 n5 Log[n5/(n1 + n2 + n3 + n4 + n5 + n6)] -
117715.30870437803 n6 + 8920.922 n6 Log[n6/(n1 + n2 + n3 + n4 + n5 + n6)]


Now when i evaluate using FindRoot, it gives me a list of solutions in i. Why is this? I am not using any subscripts in the definition, so why is i coming up?

• Try to use FindRoot instead of Solve. – Henrik Schumacher Feb 7 '18 at 23:11
• the function f is the sum of two summations i.e – MinimizeMyFunction Feb 10 '18 at 15:25
• I suspect that Minimize might be able to use Lagrange multipliers natively. Am I right in understanding that g1==g2==g3==g4==0? – LLlAMnYP Feb 19 '18 at 13:34
• However, it so appears that your constraints are self contradicting, since calling a Reduce on them returns False. – LLlAMnYP Feb 19 '18 at 13:39
• yes all constraints equal 0. can you elaborate on how my constraints are self contradicting? – MinimizeMyFunction Feb 19 '18 at 19:15

You may use FindRoot instead of Solve to find KKT points of equality constraints.

vars = {n1, n2, n3, n4, n5, n6, λ1, λ2, λ3, λ4};
eq = D[h @@ vars, {vars, 1}];
FindRoot[eq, Evaluate[Transpose[{vars, 1 + 0 vars}]]]


{n1 -> 0.35893, n2 -> 1.36013, n3 -> 1.28094, n4 -> 0.517309, n5 -> 0.920811, n6 -> 1., λ1 -> 100438., λ2 -> -203103., λ3 -> -37358.9, λ4 -> -66411.1}

However, getting global minima or all critical points is not straight-forward...

• Ive tried to use the code above and this was my output FindRoot[eq, {{n1, 1}, {n2, 1}, {n3, 1}, {n4, 1}, {n5, 1}, {n6, 1}, {[Lambda]1, 1}, {[Lambda]2, 1}, {[Lambda]3, 1}, {[Lambda]4, 1}}] is there something im missing here? also what would getting global minima require? – MinimizeMyFunction Feb 8 '18 at 1:31
• Well, you have to evaluate your own definitions first. – Henrik Schumacher Feb 10 '18 at 17:22

I'm missing the constraints in @Henrik Schumacher answer. Here my attempt with constraints:

constr = Map[0 <= # <= n1 + n2 + n3 + n4 + n5 + n6 &, {n1, n2, n3, n4, n5, n6}]
(* the constraints could be simplified using Reduce...*)

NMinimize[{1, Join[Map[# == 0 &, eq], constr]}, vars]
(* {1.,
{n1 -> 0.35893, n2 -> 1.36013, n3 -> 1.28094, n4 -> 0.517309,n5 -> 0.920811,n6 -> 1.,
\[Lambda]1 ->100438., \[Lambda]2 -> -203103., \[Lambda]3 ->-37358.9,\\[Lambda]4 -> -66411.1}} *)


The result coincides with @Henrik Schumacher answer.

appendix

If I understand the question in the right way one could solve the problem without lagrangemultipliers

varn = {n1, n2, n3, n4, n5, n6}

NMinimize[{Apply[f, varn],
Join[{n1+n2+n3-3==0, n1+2 n2+n5-4==0, 4 n3+2 n4+2 n5-8==0, 2 n6-2==0}, constr]}
, varn]


Again the solution is the same as @Henrik Schumacher gave...

• Well, I just did not tackle inequality constraints in my post because solving the resulting KKT conditions (that also involve inequalities) is not that straight-forward by hand (NMinimize might use penalization, interior point methods or similar methods... ). That "my" minimizers coincide with yours is pure luck. – Henrik Schumacher Feb 8 '18 at 11:06

This is not going to be an answer per se as it is not going to offer a solution to the question asked.

It will simply demonstrate reasons for thinking the question itself more thoroughly.

The function in the question is a special case of $f(\vec x)=\sum_{i=1}^{k}a_ix_i+b\sum_{i=1}^{k}x_ilog(\frac{x_i}{\sum_{j=1}^{k}x_j})$, where $a_i<0$ and $b>0$.

The first derivatives of such functions are given by $\frac{\partial f}{\partial x_i}=a_i+blog(\frac{x_i}{\sum_{j=1}^{k}{x_j}})$.

The Hessian matrix ($H$) of such functions has on its main diagonal $b(\frac{1}{x_i}-\frac{1}{\sum_{j=1}^{k}{x_j}})$ and off-diagonal elements are equal to $-\frac{b}{\sum_{j=1}^{k}{x_j}}$. Note, that $Det(H)\equiv0$.

Additionally, in the question-where $k=6$-there are $4$ equality constraints, $6$ inequality constraints and $6$ non-negativity constraints.

If we assume that the proportions of each class of constraints remain constant, as the size of the problem increases, then the following code produces random instances of problems similar to the one in the question:

BlockRandom[

Block[{vars, pattn, sum, logs, alphas, betas, obj, A1, b1, A2, lambdas, mus, ypsilons, eqcs, ineqs, nncs, l},

(* 'Nn' is the number of variables and 'Ss' is the symbol to use for the variables *)
With[{Nn = 3, Ss = "n"},

(* returns a list of variables *)
vars[n_, str_] := Table[Unique[str, Temporary], n];

(* returns a list of variables' patterns *)
pattn[vars_] := Pattern[#, Blank[]] & /@ vars;

With[{vs = vars[Nn, Ss]},
With[{vspattns = pattn[vs]},

(* returns the sum of its arguments *)
sum[vspattns] := Evaluate[Total[vs]];

(* returns a list of 'Log[argument/sum[arguments]]' *)
logs[vspattns] := Evaluate[Log[vs/sum[vs]]];

With[{sm = sum[vs], lgs = logs[vs]},

With[{mu = -140223.9568567524, s = 110397.30093092154, cst = 8920.922},

(* original constraints *)
alphas = {-78416.91961691371, -293404.0587042522, -36099.19105260982, -53730.85220261606, -261977.41085974447,-117715.30870437803};
betas = {8920.922, 8920.922, 8920.922, 8920.922, 8920.922, 8920.922};

(* override these definitions to define a different (non-random) objective function *)
alphas = RandomVariate[NormalDistribution[mu, s], Nn];
betas = ConstantArray[cst, Nn];

(* define the objective function *)
obj[Sequence @@ vspattns] := Evaluate[Total[alphas vs + betas vs lgs]]

];

(* 'neqs' equals the number of (random) equality constraints to use *)
With[{neqcs = Round[Nn 2./3]},

(* original constraints *)
A1 = {{1, 1, 1, 0, 0, 0}, {1, 2, 0, 0, 1, 0}, {0, 0, 4, 2, 2, 0}, {0, 0, 0, 0, 0, 2}};
b1 = {3, 4, 8, 2};

A2 = {{0, 1, 1, 1, 1, 1}, {1, 0, 1, 1, 1, 1}, {1, 1, 0, 1, 1, 1}, {1, 1, 1, 0, 1, 1}, {1, 1, 1, 1, 0, 1}, {1, 1, 1, 1, 1, 0}};

(* override these definitions to define different (non-random) linear constraints *)
A1 = RandomVariate[DiscreteUniformDistribution[{0, 4}], {neqcs, Nn}];
b1 = RandomVariate[DiscreteUniformDistribution[{3, 8}], neqcs];

A2 = ConstantArray[1, {Nn, Nn}] - IdentityMatrix[Nn];

(* define the multipliers for the equality constraints *)
lambdas = vars[neqcs, "λ"];

(* define the multipliers for the inequality constraints *)
mus = vars[Nn, "μ"];

(* define the multipliers for the non-negativity constraints *)
ypsilons = vars[Nn, "υ"];

Block[{lambdaspattns, muspattns, ypsilonspattns},

(* define patterns associated with the multipliers of the equality constraints *)
lambdaspattns = pattn[lambdas];

(* define patterns associated with the multipliers of the inequality constraints *)
muspattns = pattn[mus];

(* define patterns associated with the multipliers of the non-negativity constraints *)
ypsilonspattns = pattn[ypsilons];

Block[{eqcsvs, ineqsvs, nncsvs, lvs, eqcspattns, ineqspattns, nncspattns, lpattns},

(* define relevant variables for the equality constraints equations *)
eqcsvs = Flatten[{vs, lambdas}];

(* define relevant variables for the inequality constraints     equations *)
ineqsvs = Flatten[{vs, mus}];

(* define relevant variables for the non-negativity constraints equations *)
nncsvs = Flatten[{vs, ypsilons}];

(* define relevant variables for the Lagrangian *)
lvs = Flatten[{vs, lambdas, mus, ypsilons}];

(* define patterns for the definition of the equality constraints equations *)
eqcspattns = Flatten[{vspattns, lambdaspattns}];

(* define patterns for the definition of the inequality constraints equations *)
ineqspattns = Flatten[{vspattns, muspattns}];

(* define patterns for the definition of the non-negativity constraints equations *)
nncspattns = Flatten[{vspattns, ypsilonspattns}];

(* define patterns for the definition of the Lagrangian *)
lpattns = Flatten[{vspattns, lambdaspattns, muspattns, ypsilonspattns}];

(* The following definitions assist in constructing the Lagrangian *)

(* definition of the equality constraints equations *)
eqcs[Sequence @@ eqcspattns] := Evaluate[Flatten[(A1.Transpose[{vs}] - Transpose[{b1}]) lambdas]];

(* definition of the inequality constraints equations *)
ineqs[Sequence @@ ineqspattns] := Evaluate[Flatten[(A2.Transpose[{vs}]) mus]];

(* definition of the non-negativity constraints equations *)
nncs[Sequence @@ nncspattns] := Evaluate[ypsilons vs];

(* definition of the Lagrangian *)
l[Sequence @@ lpattns] := Evaluate[obj[Sequence @@ vs] - Total[eqcs[Sequence @@ eqcsvs]] - Total[ineqs[Sequence @@ ineqsvs]] - Total[nncs[Sequence @@ nncsvs]]];

(* obtain relations conducive to a minimum *)
Block[{foceqs, eqeqs, csceqs1, csceqs2, sys},

(* first order conditions *)
foceqs = Thread[D[l @@ lvs, {vs}] == 0.] // Simplify;
(* equality constraints *)
eqeqs = Thread[0. == (eqcs @@ eqcsvs)/lambdas] // Simplify;
(* complementary slackness conditions wrt to inequality constraints *)
csceqs1 = Thread[0. == (ineqs @@ ineqsvs)] // Simplify;
(* complementary slackness conditions wrt to non-
negativity constraints *)
csceqs2 = Thread[0. == (nncs @@ nncsvs)];

(* all equations bundled together *)
sys = Flatten[{foceqs, eqeqs, csceqs1, csceqs2}];

sys

]

]

]

]

]

]
]

]

], RandomSeeding -> 123]


In the code, matrices A1, b1, A2 are used to reproduce the constraints; changing their relative position, can recreate the original problem (also please remember to change the number of variables Nn to 6 in that case.)

After evaluation, the code returns a list of equations that should hold true at a minimum (it does not return the inequalities that should hold true at the minimum namely, the inequality and non-negativity constraints of the original problem and the non-negativity constraints of the respective multipliers).

However, that system of equations-either for the original problem or of the random problems I tested-cannot be solved. On my machine, I get

"Solve::ratnz: Solve was unable to solve the system with inexact coefficients. The answer was obtained by solving a corresponding exact system and numericizing the result."

and then an empty list.

I do not think that this is an issue with Mathematica, rather it has to do with the problem itself; the vanishing Hessian is probably an indication that something might go wrong. Perhaps, using the Lagrangian is not the ideal way to manage this optimization problem.

ps. the relevant inequalities can be assembled using

{
}


before or after sys, in the code.

I do not know for sure, but I would guess that Minimize probably has the method of Lagrangian multipliers built in.

Your problem seems fully analytical, so Minimize would be a viable choice.

Your original constraints are given as:

constraints = {0 == n1 + n2 + n3 - 3,
0 == n1 + 2 n2 + n5 - 4,
0 == 4 n3 + 2 n4 + 2 n5 - 8,
0 == 2 n6 - 2}~Join~Thread[{n1, n2, n3, n4, n5, n6} >= 0]


and I copy-paste your function as

function =
-78416.91961691371  n1 + 8920.922 n1 Log[n1/(n1 + n2 + n3 + n4 + n5 + n6)] -
293404.0587042522  n2 + 8920.922 n2 Log[n2/(n1 + n2 + n3 + n4 + n5 + n6)] -
36099.19105260982  n3 + 8920.922 n3 Log[n3/(n1 + n2 + n3 + n4 + n5 + n6)]-
53730.85220261606  n4 + 8920.922 n4 Log[n4/(n1 + n2 + n3 + n4 + n5 + n6)]-
261977.41085974447 n5 + 8920.922 n5 Log[n5/(n1 + n2 + n3 + n4 + n5 + n6)] -
117715.30870437803 n6 + 8920.922 n6 Log[n6/(n1 + n2 + n3 + n4 + n5 + n6)]


then trivially

Minimize[{function, constraints}, {n1, n2, n3, n4, n5, n6}]

{-942790., {n1 -> 0.35893, n2 -> 1.36013, n3 -> 1.28094,
n4 -> 0.517309, n5 -> 0.920811, n6 -> 1.}}


On the other hand if I update your constraints

constraints = {0 == n1 + n2 + n3 - 432,
0 == n1 + 2 n2 + n5 - 259,
0 == 4 n3 + 2 n4 + 2 n5 - 628,
0 == 2 n6 - 4}~Join~Thread[{n1, n2, n3, n4, n5, n6} >= 0]


then

Reduce[constraints]

False


I.e. your constraints cannot be all simultaneously satisfied, which is a hint why you might not be getting a solution.