6
$\begingroup$

Let me start with an example. Let

$$\mathbf{A}=\begin{bmatrix}3&1\\2&3\\1&5\end{bmatrix},$$

and let $Q=\{\mathbf{q}\vert\mathbf{q}\in\Bbb R^3_+ \land \sum_i^n q_i=1\}$ and $\alpha=\bigl(\frac{1}{2},\frac{1}{2},\frac{1}{2}\bigr)$. Abusing notation, let $\mathbf{q}^{\frac{1}{\alpha}}=\bigl(q_1^{\frac{1}{\alpha_1}},q_2^{\frac{1}{\alpha_2}},q_3^{\frac{1}{\alpha_3}}\bigr)=\bigl(q_1^2,q_2^2,q_3^2\bigr).$ Consider the set $K=\{\mathbf{k}\in\Bbb R^2_+\vert \mathbf{k}=\mathbf{A}^{\top}\mathbf{q}^{\frac{1}{\alpha}}\land \mathbf{q}\in Q\}.$ There exist two members of $K$ such that the $j$-th component of $\underline{\mathbf{k}}^j$ is lower than the $j$-th component of any other $\mathbf{k}\in K$. These can be easily identified using calculus. In this case, $\underline{\mathbf{k}}^1=\bigl(\frac{6}{11},\frac{211}{121}\bigr),$ that corresponds to $\mathbf{q}_{\underline{\mathbf{k}}^1}=\bigl(\frac{2}{11},\frac{3}{11},\frac{6}{11}\bigr)$ is the $\mathbf{k}$ with the lowest first component and $\underline{\mathbf{k}}^2=\bigl(\frac{734}{529},\frac{15}{23}\bigr),$ corresponding to $\mathbf{q}_{\underline{\mathbf{k}}^2}=\bigl(\frac{15}{23},\frac{5}{23},\frac{3}{23}\bigr)$, the one with the lowest second component. $K$, (blue) $\underline{\mathbf{k}}^1$ and $\underline{\mathbf{k}}^2$ (red) are shown in the following figure.

figure

I need to pinpoint those $\mathbf{k}$ that are on the convex part of the frontier of $K$ between $\underline{\mathbf{k}}^1$ and $\underline{\mathbf{k}}^2$, i.e., the yellow dots. In order to identify them, I

  1. Choose $\beta\in[0,1]$ and compute $\hat{\mathbf{k}} =\beta\underline{\mathbf{k}}_1+(1-\beta)\underline{\mathbf{k}}_2.$ This yields a point in the straight segment that connects $\underline{\mathbf{k}}_1$ and $\underline{\mathbf{k}}_2,$ i.e., the green dots.
  2. Determine the proportion $r=\frac{\hat k_1}{\hat k_2}.$ By construction $\frac{\underline{k}_2^2}{\underline{k}_1^2}\leq r \leq\frac{\underline{k}_2^1}{\underline{k}_1^1}.$
  3. Minimize $k_1$ subject to $\mathbf{k}\in K$ and $k_2=r k_1.$ Let $k^*_1$ be the solution to this problem and $\mathbf{k}^*=(k^*_1,rk^*_1).$

For every $\beta \in[0,1]$ this procedure should yield a $\mathbf{k}^*$ that is on the frontier. These $\mathbf{k}^*$ have an additional property: the solution to the problem $\max\sum_{i=1}^n q_i, s.t.\mathbf{A}^{\top}\mathbf{q}^{\alpha} \leq \mathbf{k}^* \land \mathbf{q} \geq 0,$ denoted $\mathbf{q}^*$, is such that $\sum_i q_i^*=1.$

However, when I solve this maximization problem, only a handful of the $\mathbf{k}^*$ are actually such that $\sum q^*_i=1$ (purple dots). This is very unfortunate, because what I'm really interested in are the KKT multipliers in this last problem, in particular those associated with the $\mathbf{A}^{\top}\mathbf{q}^{\alpha} \leq \mathbf{k}^*$ constraints.

I'm guessing that the problem can lie in the fact that $\mathbf{k}^*$ identified in the minimization step are just approximations to the true $\mathbf{k}^*$ on the frontier. However, increasing WorkingPrecision to, e.g., 30, only makes matters worse, as fewer $\mathbf{k}^*$ result in $\sum q^*_i=1$.

Edit: Do you believe this can be the case? If so, can you propose a workaround or a way of identifying the true $\mathbf{k}^*$?

Ideally this should work for any ${n\times J}$ matrix $\mathbf{A}$ with $a_{ij}\geq0$ (and $a_{ij}>0$ for most $ij$, there cannot be any rows or columns that consist only of $0$s), $Q=\{\mathbf{q}\vert \mathbf{q}\in\Bbb R^n_+\land\sum_i^n q_i=1\}$ and an $\alpha\in\Bbb R^n_+$ such that $\alpha_i\in]0,1[$ $\forall i=1,\dots,n$.

My code follows:

A = {{3, 1}, {2, 3}, {1, 5}};
t = Dimensions[A][[1]];
f = Dimensions[A][[2]];
qVec = Array[q, t];
αVec = ConstantArray[1/2, t];
onesVec = ConstantArray[1, t];
zerosVec = ConstantArray[0, t];
needs = Transpose[A].qVec^(1/αVec);
kcrit = 
  Table[
    needs /. 
      Minimize[
        {needs[[i]], qVec.onesVec == 1, Thread[0 <= qVec <= 1]}, 
        qVec, 
        Reals][[2]], 
    {i, 1, f}];
step = 1/100;
grid = Flatten[Permutations /@ IntegerPartitions[1, {t}, Range[0, 1, step]], 1];
lg = Length[grid];
plotK = 
  ListPlot[
    Table[needs /. Thread[qVec -> grid[[i]]], {i, 1, lg}], 
    AspectRatio -> 1, PlotRange -> {{0, 3.2}, {0, 5.2}}];
plotkcrit = ListPlot[kcrit, AspectRatio -> 1, PlotStyle -> Red];
grid = Flatten[Permutations /@ IntegerPartitions[1, {f}, Range[0, 1, step]], 1];
lg = Length[grid];
line = Table[Transpose[kcrit].grid[[i]], {i, 1, lg}];
plotline = ListPlot[line, PlotStyle -> Green];

frontier = {};
frontier2 = {};
success = 0;
Do[
  r = Table[line[[g, i]]/line[[g, 1]], {i, 1, f}];
  k = 
    needs /. 
      Minimize[
        Flatten[
          {needs[[1]], 
           qVec.onesVec == 1, 
           Table[needs[[i]] == r[[i]] needs[[1]], {i, 2, f}], 
           onesVec >= qVec >= zerosVec}], 
       qVec, Reals][[2]];
  AppendTo[frontier, k];
  temp = 
    NMaximize[{qVec.onesVec, onesVec >= qVec >= zerosVec, needs <= k}, qVec, Reals];
  If[temp[[1]] == 1, 
    success++; AppendTo[frontier2, needs /. temp[[2]]];];, 
  {g, 1, lg}];

plotfrontier = ListPlot[frontier, PlotStyle -> Yellow];
plotfrontier2 = ListPlot[frontier2, PlotStyle -> Purple];

Show[
  {plotK, plotline, plotfrontier, plotkcrit, plotfrontier2, 
   Graphics[{Text["k1", {0.48, 1.75}], Text["k2", {1.4, 0.55}]}]}]

Print["Purple over Yellow: ", success/lg]

Follow-up: I've tried Chris K.'s suggestion. His code is much faster than mine, but, again, is tantamount to identifying the boundary in yellow in my original post. If I check which of his $\mathbf{k}^*$ result in $\sum q^*_i=1$, I obtain a very similar picture as above (orange is $\mathbf{k}^*$ obtained with Chris's method, purple those $\mathbf{k}^*$ that result in $\sum q^*_i=1$). It's striking that the number of purple points in both pictures is the same, 2% of the yellow/orange points, though the specific $\mathbf{k}^*$ for which this happens are different.

figure

Also, I don't see how Chris's method would generalize to $n>3$, e.g., to

$$\mathbf{A}=\begin{bmatrix}3&1\\2&3\\1&5\\5&\frac{1}{2}\end{bmatrix},$$

$\endgroup$
  • $\begingroup$ I don't see any clear-cut question this post; i.e., a question for which the form and content of an acceptable answer is apparent. $\endgroup$ – m_goldberg Mar 19 at 15:23
  • $\begingroup$ @m_goldberg, You're right, I didn't ask any specific question. I've edited my post to add it. Thanks for pointing it out. $\endgroup$ – Patricio Mar 19 at 15:37
5
$\begingroup$

If I understand your problem correctly, I think you can derive the envelope using calculus. Basically you have a function n (your needs) of two variables q1 and q2 (q3 == 1 - q1 - q2 since you write that the q's sum to one).

n[q1_, q2_] := Transpose[A].{q1^2, q2^2, (1 - q1 - q2)^2}

First, replicating your point plot:

points = ListPlot[Flatten[
  Table[n[q1, q2], {q2, 0, 1, 0.01}, {q1, 0, 1 - q2, 0.01}]
  , 1]]

Mathematica graphics

To find the envelope, set the determinant of the Jacobian matrix of the transformation n equal to zero and solve for one of the q's as a function of the other:

env = Solve[Det[D[n[q1, q2], {{q1, q2}}]] == 0, q2]

Mathematica graphics

Then you can inject q2 /. env into the definition of n and ParametricPlot (the first solution turns out to correspond to the lower envelope):

pp = ParametricPlot[n[q1, q2 /. env[[1]]], {q1, 0, 1}, PlotStyle -> Orange];
Show[points, pp]

Mathematica graphics

You already knew the turning points (how exactly did you find those??), which can be used to restrict q1 to the part you want:

pp = ParametricPlot[n[q1, q2 /. env[[1]]], {q1, 2/11, 15/23}, PlotStyle -> Orange];
Show[points, pp]

Mathematica graphics

$\endgroup$
  • $\begingroup$ The turning points are obtained looking for the global minimum of each coordinate, with kcrit = Table[ needs /. Minimize[ {needs[[i]], qVec.onesVec == 1, Thread[0 <= qVec <= 1]}, qVec, Reals][[2]], {i, 1, f}]; $\endgroup$ – Patricio Mar 19 at 18:30
  • $\begingroup$ I need to look at your answer more closely, but do you think that this would work for $\mathbf{A}_{n\times J}?.$ I think it will for $n>3,$ but I need to think a bit more for $J>2.$ $\endgroup$ – Patricio Mar 19 at 18:37
  • $\begingroup$ @Patricio Maybe, but I'd need a concrete example to think about. If you've got more than 3 q's, the envelope / frontier would be a surface, right? $\endgroup$ – Chris K Mar 19 at 19:05
  • $\begingroup$ @Patricio Check out books.google.fr/… $\endgroup$ – Chris K Mar 19 at 19:10
  • $\begingroup$ Not really. $n>3$ would still be a line, because what we're plotting is $\mathbf{A}^{\top}\mathbf{q}^{\frac{1}{\alpha}}$, which is dimension $(J\times n) \times (n \times 1)=J\times 1$. A surface would arise with $J=3$. I'll look at Gibson's text. $\endgroup$ – Patricio Mar 20 at 8:05
5
$\begingroup$

You can also use Internal`ListMin as follows:

ClearAll[f, A1, A2]
f[a_][q__] /; Length[a] == Length[{q}] + 1 := Transpose[a].{q, 1 - Plus[q]}^2

A1 = {{3, 1}, {2, 3}, {1, 5}};
pnts1 = Flatten[Table[f[A1][q1, q2], {q2, 0, 1, 0.01}, {q1, 0, 1 - q2, 0.01}], 1];

frontier1 = SortBy[First]@Internal`ListMin[pnts1];
Show[ListPlot[pnts1], 
 ConvexHullMesh[frontier1, MeshCellStyle -> {{1, All} -> Gray, 2 -> Opacity[.5, Yellow]}], 
 ListLinePlot[frontier1, PlotStyle -> Directive[Red, Thick]], 
 ImageSize -> Large]

enter image description here

A2 = {{3, 1}, {2, 3}, {1, 5}, {5, 1/2}};
pnts2 = Flatten[Table[f[A2][q1, q2, q3], {q3, 0, 1, 0.02}, {q2, 0, 1, 0.02}, 
  {q1, 0, 1 - q2 - q3, 0.02}], 2];

frontier2 = SortBy[First]@Internal`ListMin[pnts2];
Show[ListPlot[pnts2], 
 ConvexHullMesh[frontier2, MeshCellStyle -> {{1, All} -> Gray, 2 -> Opacity[.3, Yellow]}], 
 ListLinePlot[frontier2, PlotStyle -> Directive[Red, Thick]], 
 ImageSize -> Large]

enter image description here

$\endgroup$
3
$\begingroup$

Since the frontier is convex, we can just walk our way around it with whatever resolution you need by choosing weights for linear weighted sums and optimizing. In this case the weights are {surf,1-surf}.

frontierRules = Table[Minimize[{{surf, 1 - surf}.needs, qVec.onesVec == 1, 
 Thread[0 <= qVec <= 1]}, qVec, Reals][[2]], {surf, 0, 1, .25}];

They all add to one.

 qVec.onesVec /. frontierRule

{1,1,1,1,1}

The points themselves are on the frontier (purple points).

kPoints= needs /. frontier;
frontierPlot = ListPlot[kPoints, AspectRatio -> 1, PlotStyle -> {PointSize[Large], Purple}]

enter image description here

The approach is readily extensible to arbitrary dimensions. You could just generate weights either randomly or structured to cover the space in whatever resolution you'd want.

$\endgroup$
  • $\begingroup$ It works fine, but I can't see why. Can you elaborate or provide a reference? $\endgroup$ – Patricio Mar 21 at 12:41
  • $\begingroup$ I think I see what you're doing. If I'm right you're obtaining my yellow dots in a slightly different (and much more efficient) manner. I still need to check whether this results in more purple dots. $\endgroup$ – Patricio Mar 21 at 12:52
  • $\begingroup$ Since it is a convex hull, you can just randomly choose weights that are positive and sum to 1, and then optimize against that criteria. Your two end points are just optimizing with weights={1,0} and {0,1}. It extends to higher dimensions directly. $\endgroup$ – MikeY Mar 21 at 16:00
2
$\begingroup$

This is a full revamp of my previous answer (see the edit history for the old version.)

Let

mat = {{3, 1}, {2, 3}, {1, 5}};
α = 1/2;

The idea is to construct a parametrization of $\mathbf q$ that automatically satisfies the constraints $\mathbf{q}\in\Bbb R^3_+$ and $\|\mathbf q\|_1=1$. An immediate choice is to use the parametrization of a superellipsoid; in particular, the positivity constraint (i.e., keep only the first octant) allows a significant simplification of the parametric equations:

vecs[u_, v_] := {(Cos[u] Cos[v])^2, (Cos[u] Sin[v])^2, Sin[u]^2}

Visualize:

ParametricPlot3D[vecs[u, v], {u, 0, π/2}, {v, 0, π/2}, ViewPoint -> {1.3, 2.4, 2.}]

octant of a "flat" superellipsoid

Thus, we can visualize the region $K$ like so:

ParametricPlot[(vecs[u, v]^(1/α)).mat, {u, 0, π/2}, {v, 0, π/2}]

region K

Let us look at the lines of constant u and v:

mesh lines

Unfortunately, it seems the frontier desired is not part of the mesh lines, so we need to do some more work. One might notice that the frontier can be thought of as the envelope of the mesh lines, so we can try deriving the envelope equation.

Det[D[(vecs[u, v]^(1/α)).mat, {{u, v}}]] // FullSimplify
   7/8 Cos[u]^5 (-11 Sin[3 u] Sin[2 v] + 48 Sin[u]^3 Sin[4 v] +
                 Sin[u] (21 Sin[2 v] + 4 Cos[u]^2 Sin[6 v]))

It seems simpler to solve for v, so let's take the factor containing it:

fac = -11 Sin[3 u] Sin[2 v] + 48 Sin[u]^3 Sin[4 v] +
      Sin[u] (21 Sin[2 v] + 4 Cos[u]^2 Sin[6 v]);

and make a substitution:

TrigExpand[fac /. v -> ArcSin[h]] // FullSimplify
   64 h Sqrt[1 - h^2] (2 - 4 h^2 + h^4 + (-2 + 2 h^2 + h^4) Cos[2 u]) Sin[u]

We can work further with the factor that is quartic in h; let us find its roots (and impose constraints as well):

ArcSin[h] /. (Solve[2 - 4 h^2 + h^4 + (-2 + 2 h^2 + h^4) cos == 0 &&
                    -1 < cos < 1 && -1 < h < 1, h, Reals] /. cos -> Cos[2 u])
   {ConditionalExpression[
    ArcSin[Root[2 - 2 Cos[2 u] + (-4 + 2 Cos[2 u]) #1^2 + (1 + Cos[2 u]) #1^4 &, 2]],
    -1 < Cos[2 u] < 1], 
    ConditionalExpression[ArcSin[Root[
    2 - 2 Cos[2 u] + (-4 + 2 Cos[2 u]) #1^2 + (1 + Cos[2 u]) #1^4 &, 3]],
    -1 < Cos[2 u] < 1]}

As it turns out, it is the second root that yields the required frontier:

Show[ParametricPlot[(vecs[u, v]^(1/α)).mat, {u, 0, π/2}, {v, 0, π/2}], 
     ParametricPlot[With[{v = ArcSin[Root[2 - 2 Cos[2 u] + (-4 + 2 Cos[2 u]) #1^2 +
                                          (1 + Cos[2 u]) #1^4 &, 3]]},
                         (vecs[u, v]^(1/α)).mat], {u, 0, π/2}, 
                    PlotStyle -> Directive[AbsoluteThickness[4],
                                           ColorData[97, 3]]]]

frontier

One can derive a more conventional parametric equation for the frontier, if desired:

With[{v = ArcSin[Root[2 - 2 Cos[2 u] + (-4 + 2 Cos[2 u]) #1^2 +
                      (1 + Cos[2 u]) #1^4 &, 3]]},
     (vecs[u, v]^(1/α)).mat] // ToRadicals // FullSimplify
   {1/4 (37 + 15 Cos[4 u] + 7 Sqrt[2] Sqrt[7 - 8 Cos[2 u] + 3 Cos[4 u]] - 
    2 Cos[2 u] (21 + 4 Sqrt[2] Sqrt[7 - 8 Cos[2 u] + 3 Cos[4 u]])), 
    1/4 (38 + 12 Cos[4 u] + 7 Sqrt[2] Sqrt[7 - 8 Cos[2 u] + 3 Cos[4 u]] -
    Cos[2 u] (42 + 5 Sqrt[2] Sqrt[7 - 8 Cos[2 u] + 3 Cos[4 u]]))}

which can then be converted to an implicit Cartesian equation:

eq = First[GroebnerBasis[Append[MapAll[TrigExpand, Thread[{x, y} == %]], 
                                Cos[u]^2 + Sin[u]^2 == 1],
                         {x, y}, {Cos[u], Sin[u]}]]
   -2401 + 4116 x - 6321 x^2 + 2828 x^3 - 1587 x^4 + 2058 y + 2940 x y -
   1050 x^2 y + 5520 x^3 y - 2352 y^2 - 84 x y^2 - 6318 x^2 y^2 +
   574 y^3 + 2640 x y^3 - 363 y^4

You can then obtain your $\underline{\mathbf{k}}^j$ like so:

ArgMin[{x, eq == 0}, {x, y}]
   {6/11, 211/121}

ArgMin[{y, eq == 0}, {x, y}]
   {734/529, 15/23}

As I had noted, for a general $n\times j$ matrix, one would now need to use modified hyperspherical coordinates; e.g.

vec[u_, v_, w_] = CoordinateTransformData[{"Hyperspherical", 4} -> "Cartesian",
                                          "Mapping", {1, u, v, w}]^2

for a four-dimensional parametrization. This is a bit more elaborate to do, so I'll leave that for the next edit of this answer.

$\endgroup$
  • $\begingroup$ If I combine your two plots, the second does not correspond to the lower enveloppe. I'm guessing obtaining it would requiere a non-zero value for v, but the ones I've tried ($v=\pi/2$, or setting $u=0$ or $u=\pi/2$ and letting v vary) don't work either. I know nothing about hyperspherical coordinates, so I'm really lost here. I've taken a look at the link your provided, but the material there is way over my math abilities. $\endgroup$ – Patricio Mar 25 at 13:33
  • $\begingroup$ Ah, sorry about that; I had not properly checked the curve. I'm not at a computer right now, but if you are willing to wait for a day, I can edit this solution the next day. $\endgroup$ – J. M. is away Mar 25 at 14:05
  • $\begingroup$ I'll also see if I can write a quick summary of hyperspherical coordinates for your case. $\endgroup$ – J. M. is away Mar 25 at 14:21
  • $\begingroup$ I'm not only willing, but also extremely grateful $\endgroup$ – Patricio Mar 25 at 15:29
  • $\begingroup$ If you had the time to check what the appropriate v is, it'd help me a lot. Thank you. $\endgroup$ – Patricio Mar 27 at 9:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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