# Removing circular rules'

Say that I have a huge list (>300.000 elements) of polynomial equations and I want to simplify these by

1. looking for some simple equations of the form x[i] == _ or _ == x[i],
2. converting these equations to a list of rules of the form x[i] -> pol[i] where pol[i] is a polynomial that can depend on x[1],...,x[n] (but typically only depends on a few variables), and
3. substituting these rules repeatedly via ReplaceRepeated in order to reduce the number of variables in the system.

The problem that typically arises is that some of these rules might be circular. For example from the system

{
x[1] == x[2] * x[3] * x[5],
x[4] == x[5] * x[6],
x[5] == x[1],
...
}


I would get a list of rules that looks like


{
x[1] -> x[2] * x[3] * x[5],
x[4] -> x[5] * x[6],
x[5] -> x[1],
...
}


Using ReplaceRepeated will run forever.

The question I have is: what would be an efficient way to select the biggest subset of a list of rules of the form x[i] -> pol[i] such that this subset is non-circular, i.e. ReplaceRepeated would not run forever? It may be assumed that any 'self-circular' rule (in which a variable appears both on the LHS and RHS) has already been removed.

I used to have an implementation that looks like

clashingRulesQ[ a_ -> b_, c_ -> d_ ] :=
With[ { vb  = Variables[b], vd = Variables[d] },
Or[ a == c, MemberQ[a] @ vd, MemberQ[c] @ vb ]
]

removeCircularRules[ ruleList_ ] :=
DeleteDuplicates[
ruleList,
clashingRulesQ
]


but this one removes more rules than necessary since it would, e.g. reduce

{ x[1] -> x[2], x[2] -> x[3] }


to

{ x[1] -> x[2] }

• Could you just remove a rule if it replaces a "later" variable with an "earlier" one? That is, x[5] gets replaced by x[1], and $5>1$ so you remove it, but x[4] gets replaced by variables x[5]*x[6], and $4 <5,6$, and hence you keep it. You could also refine this a bit further and remove, say x[5] -> x[1] only if x[1] has a replacement rule. Otherwise, leave it. By the way, your implementation of removeCircularRules seems to leave {x[1] -> x[2] x[3] x[5], x[4] -> x[5] x[6]} for me, not { x[1] -> x[2] * x[3] * x[5] } . Commented Apr 4 at 18:46
• IGraph/M has an implementation of finding minimum feedback arc sets, if that helps. Commented May 2 at 12:10

## 1 Answer

Could set it up as an integer linear programming problem. I'm short on time at the moment so I'll just post code with minimal commentary. I'll assume you have culled equations of the form var==product where the "product" could just be a single variable. The idea is to use a 0-1 valued variable for each left-hand side, find right-hand sides where each such variable appears, and constrain so the variable corresponding to each such rhs, plus the lhs variable, sum to no more than 1.

eqns = {x[1] == x[2]*x[3]*x[5], x[4] == x[5]*x[6], x[5] == x[1]};
{lhs, rhs} = Transpose[List @@@ eqns];
rhs = Map[If[Head[#] === Times, Apply[List, #], List[#]] &, rhs];
rhsindices = rhs[[All, All, 1]];
lhsindices = lhs[[All, 1]];
lhsvars = Map[c, lhsindices]
lpoly = Total[lhsvars];
members =
Table[lhsvars*
Map[Boole[MemberQ[#, lhsindex]] &, rhsindices], {lhsindex,
lhsindices}] /. 0 -> Nothing
c1 = Map[0 <= # <= 1 &, lhsvars];
c2 = Flatten[
Table[Map[lhsvars[[j]] + # <= 1 &, members[[j]]], {j,
Length[lhsvars]}]];
constraints = Join[c1, c2]

(* Out[338]= {c[1], c[4], c[5]}

Out[340]= {{c[5]}, {}, {c[1], c[4]}}

Out[343]= {0 <= c[1] <= 1, 0 <= c[4] <= 1, 0 <= c[5] <= 1,
c[1] + c[5] <= 1, c[1] + c[5] <= 1, c[4] + c[5] <= 1} *)


Now just maximize the sum of variables subject to these constraints.

Maximize[{Total[lhsvars], constraints}, lhsvars, Integers]

(* Out[344]= {2, {c[1] -> 1, c[4] -> 1, c[5] -> 0}} *)


Now post-process (exercise for reader) that you can use {x[1],x[4]}`.