# Constrained Optimization doesn't converge?

Given a $$n\times J$$ matrix $$T$$ with strictly positive elements, I'm trying to obtain the values ($$\lambda_i$$) that minimize $$\sum T_{ij}\lambda_i^{\alpha_i}$$ with the constraints $$\lambda_i \geq 0$$ and $$\sum \lambda_i = 1$$ for each of the $$J$$ columns. Provided $$\alpha_i \geq 1$$, the problem has a unique solution such that $$\lambda_i > 0$$, so only constraint $$\sum \lambda_i = 1$$ will be binding. I have used

Alpha=4/3
T = {{2, 3}, {3, 2}, {1, 4}}
tecniques = Dimensions[T][[1]]
factors = Dimensions[T][[2]]
AlphaVec = ConstantArray[Alpha, tecniques]
LambdaVec = Array[Lambda, tecniques]
onesVec = ConstantArray[1, tecniques]
maxmin = Table[LambdaVec /. Minimize[{(Transpose[T].LambdaVec^(AlphaVec))[[i]],
onesVec.LambdaVec == 1}, LambdaVec][[2]], {i, 1, factors}]


And this works fine when $$\alpha=2$$. However, when I set $$\alpha=\frac{4}{3}$$, Mathematica fails to return an answer, just keeps running. I'm sure there is some problem in my code, but I cannot see it (also, it looks kind of awkard to me). I've computed the solutions manually and it did not take long. Of course, this is an example and I'd like the code to work for arbitrary $$n>J,\alpha_i$$ and $$T$$.

• Please check your syntax. The code does not run at all. – Henrik Schumacher Sep 25 '18 at 11:35
• You're right, I skipped the definition of $\alpha$ – Patricio Sep 25 '18 at 12:44

Replacing Minimize by NMinimize . You will see that in the process of computation at $$\alpha =4/3$$ complex numbers are encountered. It is necessary to add the restrictions of $$\lambda\ge 0$$

T = {{2, 3}, {3, 2}, {1, 4}};
Alpha = 4/3;
tecniques = Dimensions[T][[1]];
factors = Dimensions[T][[2]];
AlphaVec = ConstantArray[Alpha, tecniques];
LambdaVec = Array[Lambda, tecniques];
onesVec = ConstantArray[1, tecniques];
h = (Transpose[T].LambdaVec^(AlphaVec));
maxmin = Table[
LambdaVec /.
NMinimize[{h[[i]], {onesVec.LambdaVec == 1,
Table[Lambda[i] >= 0, {i, 1, tecniques}]}}, LambdaVec][[
2]], {i, 1, factors}] // Quiet
{{0.201207, 0.0314277, 0.767365}, {0.239077, 0.683654, 0.077269}}


This task has an exact solution, which the author indicated and which can easily be found using the Mathematica. But it does not coincide with the numerical solution that we found using NMinimize. We consider the function h={2 Lambda[1]^(4/3) + 3 Lambda[2]^(4/3) + Lambda[3]^(4/3), 3 Lambda[1]^(4/3) + 2 Lambda[2]^(4/3) + 4 Lambda[3]^(4/3)}, find its extrema, using the standard analysis. Put h[[1]]=q, h[[2]]=q1, Lambda[1]=x,Lambda[2]=y, using constraints, we find

q = 2*x^(4/3) + 3*y^(4/3) + (1 - x - y)^(4/3); q1 =
3*x^(4/3) + 2*y^(4/3) + 4*(1 - x - y)^(4/3);


Necessary conditions for an extremum

eq = {D[q, x] == 0, D[q, y] == 0} // FullSimplify
Out[]= {9 x + y == 1, x + 28 y == 1}

Solve[eq, {x, y}]

Out[]={{x -> 27/251, y -> 8/251}}
eq1 = {D[q1, x] == 0, D[q1, y] == 0} // FullSimplify

Out[]= {91 x + 64 y == 64, 8 x + 9 y == 8}

Solve[eq1, {x, y}]

Out[]= {{x -> 64/307, y -> 216/307}}


The third value is found as Lambda[3]=1-x-y. And so we got the exact solution, which the author indicated in the comments. Now we need to get a numerical solution that would not be as rude as we indicated above using NMinimize. I will indicate a simple solution to the problem. We use a special method

maxmin2 =
Table[LambdaVec /.
NMinimize[{h[[i]], {onesVec.LambdaVec == 1,
Table[Lambda[i] >= 0, {i, 1, tecniques}]}}, LambdaVec,
WorkingPrecision -> 30, MaxIterations -> 100,
Method -> "RandomSearch"][[2]], {i, 1, factors}]
{{0.10756972111553840466091432374, 0.0318725099601588814911179025068,
0.860557768924302713847967773754}, {0.20846905537459283387622149837,
0.703583061889250814332247557003, 0.0879478827361563517915309446257}}


We compare it with the exact solution

{{27/251, 8/251, 216/251}, {64/307, 216/307, 27/307}}*1.30

{{0.107569721115537848605577689243, 0.0318725099601593625498007968127,
0.860557768924302788844621513944},{0.208469055374592833876221498371,
0.703583061889250814332247557003, 0.0879478827361563517915309446254}}


And so, we numerically reproduced the exact solution with an error of $$10^{-15}$$.

• I was just about to write the same... (+1) – Henrik Schumacher Sep 25 '18 at 12:48
• It was obvious, but the code contains a typo. You need to fix the code, and then add restrictions. But still in the process of computation, complex numbers are encountered. So I added Quiet. – Alex Trounev Sep 25 '18 at 12:54
• Here's another way to get these results: NArgMin[#.(v^(4/3)), v ∈ RegionIntersection[Hyperplane[{1, 1, 1}, 1], ConicHullRegion[{{0, 0, 0}}, IdentityMatrix[3]]]] & /@ Transpose[T] – J. M.'s discontentment Sep 25 '18 at 13:02
• @J. M. is somewhat okay It's funny, your code gives a similar warning 'NArgMin::nrnum: The function value 1.30754 -0.334771 I is not a real number at {Subscript[v, 1],Subscript[v, 2],Subscript[v, 3]} = {-0.291501,0.0796861,1.21182}. – Alex Trounev Sep 25 '18 at 13:20
• Yes, it's funny that NMinimize[] (and thus NArgMin[] as well) still evaluates outside a constraint region even when it's not supposed to. I have also tried using Surd[v, 3]^4 instead, but got different results, and was not sufficiently motivated to investigate further. – J. M.'s discontentment Sep 25 '18 at 13:23