Faster way to find longest chain of inequalities

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.

@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]])]

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]