Faster way to find longest chain of inequalities

Question

I have a list of algebraic expression, all variables are positive, some of them have a specific size relationship.
Such as: $\frac{(a b+a c+b c)^2}{a+b+c}\leq \frac{1}{3} (a+b+c) (a b+a c+b c)\leq \frac{(a b+a c+b c) \left(a^2+b^2+c^2\right)}{a+b+c}\leq \frac{1}{3} (a+b+c) \left(a^2+b^2+c^2\right)\leq \frac{\left(a^2+b^2+c^2\right)^2}{a+b+c}$
I need to find chain of inequalities as long as possible. Since some algebraic relations are indeterminate, it seems not feasible to use Sort.
The algorithm I used is generate permutation first than using multiple random numbers to determine less or greater. It works slowly, Is there a better way?

sets={a^3+b^3+c^3,a b^2+a^2 c+b c^2,(a^2+b^2+c^2)^2/(a+b+c),(a b+a c+b c)^2/(a+b+c),(a b c (a+b+c)^2)/(a b+a c+b c),(a^4+b^4+c^4)^2/(a^5+b^5+c^5),(a^5+b^5+c^5)^2/(a^7+b^7+c^7),1/3 (a+b+c) (a b+a c+b c),(3 (a b^3+a^3 c+b c^3))/(a+b+c),(3 (a^3 b+b^3 c+a c^3))/(a+b+c),(3 a b c (a^2+b^2+c^2))/(a b+a c+b c),1/3 (a+b+c) (a^2+b^2+c^2),(3 (a^2 b^2+a^2 c^2+b^2 c^2))/(a+b+c),a^2 b+a b^2+a^2 c-3 a b c+b^2 c+a c^2+b c^2,1/2 (a^2 b+3 a b c+b^2 c+a c^2),1/2 (a b^2+a^2 c+3 a b c+b c^2),(a b c (a^2+b^2+c^2)^2)/(a^2 b^2+a^2 c^2+b^2 c^2),1/2 (a^3+a^2 b+b^3+b^2 c+a c^2+c^3),1/2 (a^3+a b^2+b^3+a^2 c+b c^2+c^3),((a b+a c+b c) (a^2+b^2+c^2))/(a+b+c)};
(*Sort[sets,With[{c=1},Resolve[ForAll[{a,b},And[a>0,b>0],#<#2]]]&]*)

n=4;
subsets=Permutations[sets,{n}];

cf=Compile[{a,b,c},Boole[Less@@@subsets]//Evaluate,RuntimeAttributes->{Listable},RuntimeOptions->"Speed"];

res=MinimalBy[LessEqual@@@Pick[subsets,Min/@Transpose[cf@@RandomReal[1,{3,1000}]],1],LeafCount,3]

AllTrue[res,Resolve[ForAll[{a,b,c},And[a>0,b>0,c>0],#]]&]

Answer 1

This is a surprisingly complicated problem, definitely.

Starting from the definition of sets, my first approach was to create a table of all meaningful arrangements of the expressions in sets:

tablele = 
  Table[sets[[i]] <= sets[[j]], {i, 1, Length[sets]}, {j, 1, Length[sets]}];

Much of the difficulty, however, is in determining which of these inequalities are or aren't feasible. We can generate a set of random conditions to try, and eliminate all positions in the table which fail them:

genlist[l_] := Intersection @@ Map[Position[tablele /. #, True] &, l, 1];
cer = {a -> #[[1]], b -> #[[2]], c -> #[[3]]} & /@ RandomReal[{0, 10}, {200, 3}];
pos = genlist[cer];

This will eliminate a decent number of cases from consideration with relatively little computational effort (a few seconds or so). The result stored in pos is the remaining set of candidate positions which may be globally true for the conditions given (a, b, c all positive).

Now, for each remaining position in pos, we want to attempt to find a counterexample specific to the expression and disprove it:

ie[i_, j_] := sets[[i]] <= sets[[j]]; (* Helper function to generate <= expressions *)
genspecificce = 
  Table[TimeConstrained[
    FindInstance[{Not[ie @@ position], a > 0, b > 0, c > 0}, {a, b, c}],
    2], {position, pos}];

Note the use of TimeConstrained here. FindInstance will typically either return a counter example or impossible result very quickly, or take absolutely forever. With a 2 second time constraint, we can recover the simpler solutions very quickly.

We can collect all of the impossible results ({} returned, signifying that no such instance exists) into a list of definitely true inequalities:

definite = DeleteCases[Pick[pos, # == {} & /@ genspecificce], {}];

We can also collect the specific counterexamples generated and apply them in addition to our random counterexamples:

ces = DeleteCases[Flatten[genspecificce, 1], $Aborted];
pos2 = genlist[Join[cer, ces]];

Using this cleaned up pos2 list and the definite list, we can construct a list of positions which may or may not be true:

maybe = Complement[pos2, definite]

These positions are the only expression pairs which have been excluded from consideration thus far. Feel free to expend whatever computing efforts you desire on trying to prove or disprove these expressions, but these are all of the difficult ones remaining, so it may or may not work out.

Once you are satisfied with the definite list, you can remove all the loops and convert it to a graph:

g = Graph[DirectedEdge @@ # & /@ DeleteCases[definite, {x_, x_}], VertexLabels -> "Name"]

This is fairly messy for the set of expressions listed above, but we can construct a "find longest path" method using DepthFirstScan:

findLongestPaths[graph_, start_] := 
  Module[{array = <|start -> {0, None}|>, dfs, maximal, group, paths},
   dfs = Reap[
     DepthFirstScan[graph, 
      start, {"FrontierEdge" -> (Sow[#, 0] &), 
       "ForwardEdge" -> (Sow[#, 1] &)}], {0, 1}];
   Do[array[e[[2]]] = {array[e[[1]]][[1]] + 1, e[[1]]}, {e, 
     dfs[[2, 1, 1]]}];
   If[Length[dfs[2]] > 1, 
    Do[If[array[e[[1]]][[1]] >= array[e[[2]]][[1]], 
      array[e[[2]]] = {array[e[[1]]][[1]] + 1, e[[1]]}], {e, 
      dfs[[2, 2, 1]]}]];
   group = GroupBy[array, #[[1]] &];
   (*group[Max[Keys[group]]]*)
   paths = 
    NestWhileList[array[#][[2]] &, #, NumericQ[#] &] & /@ 
     group[Max[Keys[group]]][[All, 2]];
   Table[Reverse[Join[{k}, paths[k]][[1 ;; -1]]], {k, 
      Keys[paths]}] /. None -> start
   ];

This method returns all of the longest paths found from a given starting vertex. Now we just need to select an appropriate starting vertex. I chose every vertex with an in-degree of 0 and combined the results:

candidatePaths = 
 Flatten[findLongestPaths[g, #] & /@ 
   Pick[VertexList[g], # == 0 & /@ VertexInDegree[g]], 1]

We can clean this up a bit by sorting and substituting the original expressions back in:

LessEqual @@ # & /@ 
 Map[Part[sets, #] &, Sort[candidatePaths, Length[#1] > Length[#2] &]]

The first result in this final sorted group is:

1/2 (a b^2 + a^2 c + 3 a b c + b c^2) <= 1/2 (a b^2 + a^2 c + 3 a b c + b c^2) <= a b^2 + a^2 c + b c^2 <= 1/3 (a + b + c) (a^2 + b^2 + c^2) <= 1/2 (a^3 + a b^2 + b^3 + a^2 c + b c^2 + c^3)

Strongly suggesting that the longest chain of inequalities is going to be of length 5. However, I did neglect every possible comparison in maybe, so that is the limitation of this answer.

Also note that FindLongestPaths only traces back one path from the most distant vertex. If you actually need every longest path, that becomes more complicated, but it can be modified to do that too.

Answer 2

@eyorble's solution can already yield great results, especially when dealing with inequalities. My main idea is similar to @eyorble's, but with more generalization and simplification.

Step 0: Create a fast FindInstance

The problem with FindInstance is that it will try to find exact solutions, which significantly slow down the computation and might leave out some solutions. So here we first test the inequality by substituting variables with random numbers. In this case, because all functions are of the same order, we can simply use var = RandomReal[{0, 1}, Length@var];. But in other use cases, you might want to tune the random function for better performance.

myFindInstance[eqn_, var_, dom_] :=
 Catch[Block[var,
   Do[var = RandomReal[{0, 1}, Length@var];
     If[eqn, Throw[1]], {10000}];];
  TimeConstrained[
   Throw@Length@
     FindInstance[eqn && (And @@ Thread[0. < var]), var, dom, 
      WorkingPrecision -> 15, RandomSeeding -> Automatic]
   , 2, Throw@0]
  ]

myFindInstance returns 0 if no instance is found, and returns 1 otherwise.

Step 1: Test inequalities and determine equations' relationship

First we shall define function edge as follows:

edge[{0, 1, 0}, i_, j_] := {Labeled[DirectedEdge[i, j], Equal], 
   Labeled[DirectedEdge[j, i], Equal]};
edge[{1, 0, 0}, i_, j_] := Labeled[DirectedEdge[j, i], Greater];
edge[{0, 0, 1}, i_, j_] := Labeled[DirectedEdge[i, j], Greater];
edge[{1, 1, 0}, i_, j_] := Labeled[DirectedEdge[j, i], GreaterEqual];
edge[{0, 1, 1}, i_, j_] := Labeled[DirectedEdge[i, j], GreaterEqual];
edge[___] := Nothing;

where the first parameters are the result of myFindInstance with (in)equality {eqn1<eqn2, eqn1==eqn2, eqn1>eqn2} and the second are the id of these two equations.

Then, we try to construct a relation graph between these equations:

gsss = Block[{g = Graph[Range@Length@sets, {}], symb, e},
  Do[
   If[Length[FindShortestPath[g, i, j]] == 
     Length[FindShortestPath[g, j, i]] == 0,
    e = edge[
      myFindInstance[#[sets[[i]], sets[[j]]], 
         DeleteDuplicates@Cases[sets[[{i, j}]], _Symbol, Infinity], 
         Reals] & /@ {Less, Equal, Greater}, i, j];
    If[e =!= Nothing, g = EdgeAdd[g, e[[1]]]; 
     PropertyValue[{g, e[[1]]}, EdgeLabels] = e[[2]]; 
     PropertyValue[{g, e[[1]]}, EdgeWeight] = -1]
    ], {i, Length[sets] - 1}, {j, i + 1, Length@sets}]; g]

Here we can reduce the computation by utilizing the fact that if a>=b and b>=c, then a>=c is automatically guaranteed (the FindShortestPath part).

Step 2: Find the longest possible chain

We can use a trick in this step: if we set EdgeWeight of each vertex to -1 then the shortest path will actually be the longest chain! so the code for finding the longest chain is simple:

FindShortestPath[gsss, ##] & @@@ With[{dm = GraphDistanceMatrix[gsss]}, Position[dm, Min@dm]]

Step 3: Visualization

No explanation.

Column[Inequality @@ (Riffle[sets[[#]], 
     MovingMap[
      PropertyValue[{gsss, DirectedEdge[#[[1]], #[[2]]]}, 
        EdgeLabels] &, #, 1]]) & /@ (FindShortestPath[gsss, ##] & @@@ 
   With[{dm = GraphDistanceMatrix[gsss]}, Position[dm, Min@dm]])]

The result will be something like:

$\begin{array}{l} a^3+b^3+c^3\geq \frac{\left(a^2+b^2+c^2\right)^2}{a+b+c}\geq \frac{1}{2} \left(a^2 b+a^3+a c^2+b^2 c+b^3+c^3\right)\geq \frac{1}{3} (a+b+c) \left(a^2+b^2+c^2\right)\geq \frac{(a b+a c+b c) \left(a^2+b^2+c^2\right)}{a+b+c}\geq \frac{1}{3} (a+b+c) (a b+a c+b c)\geq \frac{(a b+a c+b c)^2}{a+b+c} \\ a^3+b^3+c^3\geq \frac{\left(a^2+b^2+c^2\right)^2}{a+b+c}\geq \frac{1}{2} \left(a^2 b+a^3+a c^2+b^2 c+b^3+c^3\right)\geq \frac{1}{3} (a+b+c) \left(a^2+b^2+c^2\right)\geq \frac{(a b+a c+b c) \left(a^2+b^2+c^2\right)}{a+b+c}\geq \frac{1}{3} (a+b+c) (a b+a c+b c)\geq \frac{a b c (a+b+c)^2}{a b+a c+b c} \\ \end{array}$

There are two possible longest chains with a length of 7. The chains are longer than @eyorble's solution, and I am not quite sure whether they are correct, but hey, at least I'm unable to find any counter-example using Mathematica.

The complete code is as follows:

sets = {a^3 + b^3 + c^3, 
   a b^2 + a^2 c + 
    b c^2, (a^2 + b^2 + c^2)^2/(a + b + c), (a b + a c + b c)^2/(a + 
      b + c), (a b c (a + b + c)^2)/(a b + a c + 
      b c), (a^4 + b^4 + c^4)^2/(a^5 + b^5 + 
      c^5), (a^5 + b^5 + c^5)^2/(a^7 + b^7 + c^7), 
   1/3 (a + b + c) (a b + a c + b c), (3 (a b^3 + a^3 c + b c^3))/(a +
       b + c), (3 (a^3 b + b^3 c + a c^3))/(a + b + 
      c), (3 a b c (a^2 + b^2 + c^2))/(a b + a c + b c), 
   1/3 (a + b + c) (a^2 + b^2 + 
      c^2), (3 (a^2 b^2 + a^2 c^2 + b^2 c^2))/(a + b + c), 
   a^2 b + a b^2 + a^2 c - 3 a b c + b^2 c + a c^2 + b c^2, 
   1/2 (a^2 b + 3 a b c + b^2 c + a c^2), 
   1/2 (a b^2 + a^2 c + 3 a b c + 
      b c^2), (a b c (a^2 + b^2 + c^2)^2)/(a^2 b^2 + a^2 c^2 + 
      b^2 c^2), 1/2 (a^3 + a^2 b + b^3 + b^2 c + a c^2 + c^3), 
   1/2 (a^3 + a b^2 + b^3 + a^2 c + b c^2 + 
      c^3), ((a b + a c + b c) (a^2 + b^2 + c^2))/(a + b + c)};

myFindInstance[eqn_, var_, dom_] :=
 Catch[Block[var,
   Do[var = RandomReal[{0, 1}, Length@var];
     If[eqn, Throw[1]], {10000}];];
  TimeConstrained[
   Throw@Length@
     FindInstance[eqn && (And @@ Thread[0. < var]), var, dom, 
      WorkingPrecision -> 15, RandomSeeding -> Automatic]
   , 2, Throw@0]
  ]

edge[{0, 1, 0}, i_, j_] := {Labeled[DirectedEdge[i, j], Equal], 
   Labeled[DirectedEdge[j, i], Equal]};
edge[{1, 0, 0}, i_, j_] := Labeled[DirectedEdge[j, i], Greater];
edge[{0, 0, 1}, i_, j_] := Labeled[DirectedEdge[i, j], Greater];
edge[{1, 1, 0}, i_, j_] := Labeled[DirectedEdge[j, i], GreaterEqual];
edge[{0, 1, 1}, i_, j_] := Labeled[DirectedEdge[i, j], GreaterEqual];
edge[___] := Nothing;

gsss = Block[{g = Graph[Range@Length@sets, {}], symb, e},
  Do[
   If[Length[FindShortestPath[g, i, j]] == 
     Length[FindShortestPath[g, j, i]] == 0,
    e = edge[
      myFindInstance[#[sets[[i]], sets[[j]]], 
         DeleteDuplicates@Cases[sets[[{i, j}]], _Symbol, Infinity], 
         Reals] & /@ {Less, Equal, Greater}, i, j];
    If[e =!= Nothing, g = EdgeAdd[g, e[[1]]]; 
     PropertyValue[{g, e[[1]]}, EdgeLabels] = e[[2]]; 
     PropertyValue[{g, e[[1]]}, EdgeWeight] = -1]
    ], {i, Length[sets] - 1}, {j, i + 1, Length@sets}]; g]

Column[Inequality @@ (Riffle[sets[[#]], 
      MovingMap[
       PropertyValue[{gsss, DirectedEdge[#[[1]], #[[2]]]}, 
         EdgeLabels] &, #, 1]]) & /@ (FindShortestPath[gsss, ##] & @@@
     With[{dm = GraphDistanceMatrix[gsss]}, Position[dm, Min@dm]])]

Answer 3

Without random numbers evaluation.

Given sets and ordering into a graph.

sets = {a^3 + b^3 + c^3, a b^2 + a^2 c + b c^2, (a^2 + b^2 + c^2)^2/(a + b + c), (a b + a c + b c)^2/(a + b + c), (a b c (a + b + c)^2)/(a b + a c + b c), (a^4 + b^4 + c^4)^2/(a^5 + b^5 + c^5), (a^5 + b^5 + c^5)^2/(a^7 + b^7 + c^7), 1/3 (a + b + c) (a b + a c + b c), (3 (a b^3 + a^3 c + b c^3))/(a + b + c), (3 (a^3 b + b^3 c + a c^3))/(a + b + c), (3 a b c (a^2 + b^2 + c^2))/(a b + a c + b c), 1/3 (a + b + c) (a^2 + b^2 + c^2), (3 (a^2 b^2 + a^2 c^2 + b^2 c^2))/(a + b + c), a^2 b + a b^2 + a^2 c - 3 a b c + b^2 c + a c^2 + b c^2, 1/2 (a^2 b + 3 a b c + b^2 c + a c^2), 1/2 (a b^2 + a^2 c + 3 a b c + b c^2), (a b c (a^2 + b^2 + c^2)^2)/(a^2 b^2 + a^2 c^2 + b^2 c^2), 1/2 (a^3 + a^2 b + b^3 + b^2 c + a c^2 + c^3), 1/2 (a^3 + a b^2 + b^3 + a^2 c + b c^2 + c^3), ((a b + a c + b c) (a^2 + b^2 + c^2))/(a + b + c)};

Clear[compare]
compare[set_, sets_] := Module[{error = 10^(-14), chain = {set}, seti, val, i, offset = 100000},
  For[i = 1, i <= Length[sets], i++,
   seti = sets[[i]];
   If[Complement[{set}, {seti}] != {},
    val = Quiet@NMinimize[{set - seti, a >= 0, b >= 0, c >= 0, 
     a^2 + b^2 + c^2 <= offset}, {a, b, c}][[1]];
    If[Abs[val] < error, AppendTo[chain, seti]]
    ]
   ];
  Return[chain]
]

Net = {};
For[k = 1, k <= Length[sets], k++,
 AppendTo[Net, compare[sets[[k]], sets]]]

GR = {};

For[i = 1, i <= Length[Net], i++, n1 = Net[[i]]; 
 If[Length[n1] > 1, AppendTo[GR, Table[n1[[1]] -> n1[[i]], {i, 2, Length[n1]}]]]]

Graph[Flatten[GR], VertexLabels -> "Name", VertexStyle -> Red,  VertexSize ->{0.2, 0.025}, ImageSize -> 1500, AspectRatio -> 1]

Answer 4

For all combinations of sets[[i]]-sets[[j]], when we Together them, the denominator is positive, so we only have to worry about the numerator. This seems to help with the solving.

Also, because all of the numerators of sets[[i]]-sets[[j]] are homogeneous functions, we can pick a variable and set it arbitrarily. We only care if sets[[i]]-sets[[j]] is positive or negative, not its absolute value. So set c = 1

To find all of the edges indicating dominance (so sets[[i]] > sets[[j]] is greater for all a, b) use the below. Note the rule for c->1.

 res = Boole@Outer[
           FindInstance[{(Numerator@Together[sets[[#1]] - sets[[#2]]] /. c -> 1) > 0, a > 0, b > 0}, {a, b}] != {} && 
           FindInstance[{(Numerator@Together[sets[[#1]] - sets[[#2]]] /. c -> 1) < 0, a > 0, b > 0}, {a, b}] == {} &, 
           Range[20], Range[20]]; // Timing
 (* 9.73 sec *)

ag = AdjacencyGraph[resBoole];

Using the brute force method from this answer with a minor mod

 allPaths = FindPath[ag, #2, #1, Infinity, All] & @@@ Tuples[Range[20], 2]//Apply[Join];
MaximalBy[allPaths, Length@Union@# &]
(* {{1, 6, 19, 12, 20, 8, 4},
    {1, 6, 18, 12, 20, 8, 4}, 
    {1, 3, 19, 12, 20, 8, 4}, 
    {1, 3, 18, 12, 20, 8, 4}, 
    {1, 6, 19, 12, 20, 11, 5}, 
    {1, 6, 19, 12, 20, 8, 5}, 
    {1, 6, 18, 12, 20, 11, 5}, 
    {1, 6, 18, 12, 20, 8, 5}, 
    {1, 3, 19, 12, 20, 11, 5}, 
    {1, 3, 19, 12, 20, 8, 5}, 
    {1, 3, 18, 12, 20, 11, 5}, 
    {1, 3, 18, 12, 20, 8, 5}}
  *)

Again, these are ordered so sets[[1]] > sets[[6]] > etc.