Finding (bottom) SCCs in many systems at once using symbolic algorithms

Samuel Pastva

A gentle reminder that we should not forget this is still happening...

Motivation

Logical Modelling in Systems Biology

Boolean Network

f_1(x) = x_1 \land x_2
f_2(x) = \neg x_1 \lor x_3
f_3(x) = (x_1 \land x_2) \lor (\neg x_1 \land \neg x_3)

n Boolean variables—genes, proteins, other biochemical entities

Update functions—one per variable

Boolean Network (with inputs)

f_1(x, {\color{red} p}) = x_1 \land ({\color{red} p_1} \Rightarrow x_2)
f_2(x, {\color{red} p}) = \neg x_1 \lor (x_3 \land \neg {\color{red} p_1})
f_3(x, {\color{red} p}) = (x_1 \land x_2) \lor (\neg x_1 \land \neg x_3 \land {\color{red} p_2})

m Boolean inputs/constants/parameters—arbitrary but fixed values

Synchronous dynamics

x_{t+1} = f(x_{t}, p)

Asynchronous dynamics

f_1(x_t, p) \mid\mid f_2(x_t, p) \mid\mid f_3(x_t, p)

and others...

What happens to the system eventually?

Bottom/Terminal strongly connected components

(attractors)

Single states, cycles, more complex sets...

Long-lived transients—other non-trivial SCCs

Different input valuations lead to different outcomes

G = (V, C, E)~\text{(coloured graph)}
V = \mathbb{B}^n~\text{(states)}
C = \mathbb{B}^m~\text{(input val.)}
E \subseteq V \times C \times V~\text{(col. edges)}

Typically at least 2^50, often 2^100 or more in total

Symbolic Algorithms

with Binary Decision Diagrams*

*Bryant, Randal E. "Graph-based algorithms for Boolean function manipulation." Computers, IEEE Transactions on 100.8 (1986): 677-691.

x_1 x_2 x_3
1 0 1
1 0 0
0 1 1
(x_1 \land \overline{x_2}) \lor (\overline{x_1} \land x_2 \land x_3)

Formula is satisfied exactly for valuations from the set

Set of states:

Boolean formula:

x_1
x_2
x_2
x_3
x_3
x_3
x_3
0
0
0
1
1
1
0
0

Binary decision tree:

x_1
x_2
x_2
x_3
x_3
x_3
x_3
0
0
1
1
0

Binary Decision Diagram

x_1
x_2
x_2
x_3
0
1
1
0

Binary Decision Diagram

x_1
x_2
x_2
x_3
0
1

Binary Decision Diagram

Canonical for each Boolean function (assuming fixed variable order)

Often exponentially smaller than the original tree or a DNF formula (not always)

Bit-vector set ~ Boolean formula ~ BDD

(or relation)

Set op. ~ Logical op. ~ Graph op.

O(n^2) in the number of BDD nodes (mostly)

Not the size of the set! (or formula/tree/etc.)

Symbolic complexity model: one unit of work = one BDD operation

Symbolic complexity model is often surprising

\mathcal{O}(n^2)~\text{symbolic operations}
\mathcal{O}(n)~\text{symbolic operations}
\mathcal{O}(n)
\mathcal{O}(n^3)

The first case is "practically" much better!

Nobody can usually give you these bounds... 

...but we know

it happens "in the wild":

Reachability with saturation

\texttt{Post} (i, X \subseteq V) = \{ y \in V \mid \exists x \in X.~x \xrightarrow{i} y \}
\text{Successors using a single}~f_i:
\text{Successors using all}~f_i:
\texttt{PostAll} (X \subseteq V) = \bigcup_{i} \texttt{Post}(i, X)
\texttt{Post} (i, X \subseteq V {\color{red} \times C}) = \{ (y, c) \mid \exists (x, c) \in X.~x \xrightarrow{i}_{\color{red} c} y \}
\text{Successors using a single}~f_i:
\text{Successors using all}~f_i:
\texttt{PostAll} (X \subseteq V {\color{red} \times C}) = \bigcup_{i} \texttt{Post}(i, X)

"Normal" reachability

def FWD(initial):
  result = initial
  next = PostAll(initial)
  while next != result:
    result = next
    next = result.union(PostAll(result))
  return result

"Saturation" reachability*

def FWD(initial):
  result = initial
  done = False
  while not done:
    done = True
    for i in [1,n]:
      next = Post(i, result)
      if not next.is_subset(result):
        result = result.union(next)
        done = False
  return result

*Ciardo, Gianfranco, Gerald Lüttgen, and Radu Siminiceanu. "Saturation: an efficient iteration strategy for symbolic state—space generation." International Conference on Tools and Algorithms for the Construction and Analysis of Systems. Springer, Berlin, Heidelberg, 2001.

A simplified case for Boolean variables... more complex domains need extra loop here

"Normal" reachability:

  • Asymptotically optimal (linear in the graph diameter)
  • Typically out of memory for n>60

"Saturation" reachability:

  • Not optimal (O(|V|log(|V|)), but depends)
  • Typically makes good progress well into n>100, often n>1000 (depends on relationships between individual update functions)

Symbolic (bottom) SCC computation

(finally)

SCC closed = If x is in the set, the whole SCC of x is in the set as well.

F = FWD(set) is SCC closed.

B = BWD(set) is SCC closed.

Intersection of F and B is SCC closed.

set

FWD(set)

BWD(set)

If set = {x}, then this intersection is exactly the SCC of x.

If FWD(set) minus set is empty, then the SCC is a BSCC

Forward-Backward SCC detection*

def SCC(universe):
  if universe.is_empty():
    return
  pivot = Pick(universe)
  F = FWD(pivot)
  B = BWD(pivot)
  W = F.intersect(B)
  report(W)

  SCC(F.minus(W))
  SCC(B.minus(W))
  SCC(universe.minus(F).minus(B))

*Xie, Aiguo, and Peter A. Beerel. "Implicit enumeration of strongly connected components." 1999 IEEE/ACM International Conference on Computer-Aided Design. Digest of Technical Papers (Cat. No. 99CH37051). IEEE, 1999.

Bottom SCC detection

def BSCC(universe):
  if universe.is_empty():
    return
  pivot = Pick(universe)
  B = BWD(pivot)
  
  F = pivot
  while F.is_subset(B):
    update = PostAll(F)
    if update.is_subset(F):
      break
    F = F.union(update) 
  if F.is_subset(B):
    report(F)
  BSCC(universe.minus(B))

B is computed fully, but for F, we only check if it is contained in B.

Ideally, this is using saturation instead

Bottom SCC detection

Coloured

Pick one state per colour

Partition colours to two sets: one escapes B, the other trapped in B

def BSCC(universe):
  if universe.is_empty():
    return
  pivot = Pick(universe)
  B = BWD(pivot)
  
  F = pivot
  while F.is_subset(B):
    update = PostAll(F)
    if update.is_subset(F):
      break
    F = F.union(update) 
  if F.is_subset(B):
    report(F)
  BSCC(universe.minus(B))

Overall, this works better, but still surprisingly poorly...

Why?

1. Poor pivot selection

1. Poor pivot selection

1. Poor pivot selection

1. Poor pivot selection

2. Small BSCC, large graph

2. Small BSCC, large graph

Even if BSCC is small, BWD is still limited by graph diameter.

Pivot selection problem: If the chance of hitting a bottom SCC is so small, what if we just selected more pivots at once?

def Reduce(universe):
  pivots = ???
  
  B = BWD(pivots)
  F = FWD(pivots)
  universe = universe.minus(B.minus(F))
  
  S = F.minus(B)
  B_S = BWD(S)
  universe = universe.minus(B_S.minus(S))
  
  return universe

How do we pick the pivots?

All states that can update variable x:

  • (usually) easy to obtain the symbolic representation
  • (usually) evenly distributed across the state space
  • BSCCs (usually) do not update all variables

(Transition Guided Reduction*)

*Beneš, Nikola, et al. "Computing bottom SCCs symbolically using transition guided reduction." International Conference on Computer Aided Verification. Springer, Cham, 2021.

Now we can continue reducing based on some other update function...

This also helps with the diameter problem:

This also helps with the diameter problem:

This also helps with the diameter problem:

BWD — one step

This also helps with the diameter problem:

FWD — two steps

This also helps with the diameter problem:

FWD — two steps

(Usually) not limited by graph diameter ⇒ distance between updates to the same variable instead.

Each reduction makes the next reductions easier ⇒ no need to reconsider eliminated states.

Any sate that can reach a single-state BSCC will be eventually eliminated.

Any state that can update a variable that is constant in some reachable BSCC will be eventually eliminated.

Transition Guided Reduction

...but not all functions are created equal. Some reductions are still (symbolically) harder than others.

Can we somehow guess which reductions will be faster/better?

No? At least not really...

But we can just, well, run them all at once.

Interleaved

Transition Guided Reduction

1. Run all reductions for each variable, one step at a time.

2. Always advance the reduction with the smallest intermediate BDD.

4. Profit

3. If one reduction eliminates some states, remove them from all reductions.

Are the colours actually useful?

Extrapolating from 1024 single-coloured samples over a smaller subset of models...

Anywhere from 100x to 5000x speed-up compared to a parameter scan.

Conclusions

  • Symbolic algorithms can be surprising.
  • Internal "parallelism" of the symbolic representation can often compute "more" faster than "less".
  • Transition guided reduction: Filter states based on transitions that do not appear in (some) BSCCs.
  • Interleaving smooths-out the edge cases.