Influence Maximization, in depth: IC and LT diffusion, submodularity, greedy with CELF, reverse influence sampling and IMM

Diffusion models: Independent Cascade and Linear Threshold

A network is a directed graph G = (V, E) with n nodes and m edges. A diffusion model says how activation spreads once a seed set S is switched on at time zero. Two models dominate the literature, and the choice between them is a modelling decision, not a detail.

Independent Cascade (IC). Every edge (u, v) carries a probability p(u, v). When u becomes active at step t it gets exactly one chance to activate each inactive out-neighbour v at step t + 1, succeeding with probability p(u, v), independently of everything else. If it fails, it never tries v again. The process stops when a step activates nobody.

Linear Threshold (LT). Every edge carries a weight w(u, v) and the incoming weights of each node sum to at most 1. Each node draws a threshold uniformly from [0, 1] once, at the start. A node activates when the summed weight of its active in-neighbours reaches its threshold, modelling peer pressure rather than independent trials.

The influence spread sigma(S) is the expected number of active nodes when the process ends, counting the seeds themselves. The problem is: maximise sigma(S) subject to |S| = k.

Both models have a live-edge view that does all the theoretical work. For IC, flip a biased coin for every edge in advance and keep the edges that come up live; the nodes activated by S are exactly the nodes reachable from S in that random subgraph. For LT, each node independently keeps at most one incoming edge, choosing edge (u, v) with probability w(u, v) and none with the leftover probability. Kempe, Kleinberg and Tardos proved both constructions give the same distribution over final active sets as the step-by-step processes. So sigma(S) is the expected size of the reachable set from S in a random graph, and reachability is something we know how to reason about.

Hardness and the submodularity rescue

Choosing the seeds is NP-hard. Under IC, set every probability to 1 on a bipartite graph of sets and elements: the spread of k set-nodes is k plus the elements they cover, so the problem contains maximum coverage.

Evaluating sigma(S) for one fixed S is #P-hard. Chen, Wang and Wang showed this for IC in 2010, and a companion result covers LT. So the oracle greedy wants to call is itself unavailable, and every implementation estimates it by sampling.

What rescues the situation is structure. Under the live-edge view, for each fixed outcome of the coins the function S to |reachable(S)| is a coverage function: adding a seed adds the nodes it reaches that nobody else already reaches, and that extra can only shrink as S grows. A non-negative weighted sum of coverage functions is monotone and submodular, and sigma is exactly such an average. Submodularity means diminishing returns: for S contained in T and v outside T, sigma(S + v) - sigma(S) is at least sigma(T + v) - sigma(T). The Nemhauser, Wolsey and Fisher theorem then says greedy selection of the node with the largest marginal gain, k times, achieves at least (1 - 1/e), about 63 percent, of the optimal spread. With sampled estimates the guarantee becomes (1 - 1/e - epsilon) with high probability, where epsilon shrinks as the sample grows. Feige's result for max coverage means no polynomial algorithm beats 1 - 1/e in general unless P = NP.

A worked example by hand

Take the nine-node IC network in the figure. Hub A points at B, C, D and H, each with probability 0.5. Node X points at the same four nodes, each with probability 0.1. A separate chain runs E to F with probability 1.0 and F to G with probability 0.5. The budget is k = 2.

Single-node spreads come straight from linearity of expectation: each node counts itself plus the probability it reaches each other node. sigma(A) = 1 + 4 x 0.5 = 3.0. sigma(X) = 1 + 4 x 0.1 = 1.4. sigma(E) = 1 + 1.0 + 1.0 x 0.5 = 2.5. sigma(F) = 1.5, and B, C, D, H, G each have spread 1.

Greedy round one picks A with 3.0. Round two computes marginal gains given A. Adding E gains its full 2.5, because nothing A reaches overlaps the chain. Adding X gains only 1 for itself plus 4 x (0.55 - 0.5) = 0.2, a total of 1.2: with both A and X seeded, each of B, C, D and H is reached with probability 1 - 0.5 x 0.9 = 0.55, barely more than A alone gave. Adding B gains 0.5, the gap between certain and half-likely. So greedy picks E and sigma({A, E}) = 5.5.

A top-degree heuristic sees A and X tied at out-degree four and picks both, for sigma({A, X}) = 1 + 1 + 4 x 0.55 = 4.2, which is 24 percent worse. Node B shows submodularity directly: its gain is 1 on the empty set and 0.5 once A is seeded.

Greedy with Monte Carlo, and the CELF shortcut

The baseline is greedy over Monte Carlo estimates: simulate the cascade R times, often 10,000, and average. Naive greedy costs k x n x R simulations, each up to O(m). The first big saving is CELF, from Leskovec and colleagues in 2007. Submodularity means a node's marginal gain computed in an earlier round is an upper bound on its gain now. Keep candidates in a max-heap keyed on their last computed gain. Pop the top; if its gain was computed this round, it beats every other upper bound, so take it. Otherwise recompute its gain, push it back, and pop again. Most candidates are never re-evaluated after round one, and published experiments report speedups of several hundred times over naive greedy, with the exact factor depending on the graph.

import heapq
import random


def simulate_ic(graph, seeds, rng):
    """One cascade under Independent Cascade. graph: {u: [(v, p), ...]}."""
    active = set(seeds)
    frontier = list(seeds)
    while frontier:
        nxt = []
        for u in frontier:
            for v, p in graph.get(u, ()):
                if v not in active and rng.random() < p:
                    active.add(v)
                    nxt.append(v)
        frontier = nxt
    return len(active)


def spread(graph, seeds, runs=10_000, seed=7):
    # A fixed seed gives common random numbers: every candidate set is scored
    # against the same coin flips, so comparisons are far less noisy.
    rng = random.Random(seed)
    return sum(simulate_ic(graph, seeds, rng) for _ in range(runs)) / runs


def celf(graph, nodes, k, runs=10_000):
    heap = [(-spread(graph, [v], runs), v, 0) for v in nodes]
    heapq.heapify(heap)
    seeds, current = [], 0.0
    while len(seeds) < k and heap:
        neg_gain, v, computed_round = heapq.heappop(heap)
        if computed_round == len(seeds):
            # Gain is fresh for this round and beats every stale upper bound.
            seeds.append(v)
            current += -neg_gain
        else:
            gain = spread(graph, seeds + [v], runs) - current
            heapq.heappush(heap, (-gain, v, len(seeds)))
    return seeds, current

Run on the example graph, CELF evaluates all nine singletons, takes A, then re-evaluates E (gain 2.5), finds it beats every other stale bound, and stops, returning ['A', 'E'] with an estimate near 5.5. Monte Carlo noise can occasionally break the upper-bound property; common random numbers, as in the code, keep that rare.

Reverse influence sampling

CELF still needs fresh simulations for every re-evaluation. The idea that made influence maximization practical on graphs with billions of edges is reverse influence sampling, introduced by Borgs, Brautbar, Chayes and Lucier in 2014.

Pick a node v uniformly at random. Run the live-edge process backwards from v: for each in-edge (u, v) keep it with probability p(u, v), and continue from every node reached. The nodes found form a reverse-reachable (RR) set: exactly the nodes that, as seeds, would have activated v in this random world. Now the key identity: for any seed set S, the probability that S intersects a random RR set equals sigma(S) / n. So if you draw theta RR sets, the fraction that S covers, times n, is an unbiased estimate of sigma(S), and it is an estimate for every S at once. Choosing seeds becomes greedy maximum coverage over the collection, which runs in time linear in the total size of the sets.

def random_rr_set(reverse, nodes, rng):
    """reverse: {v: [(u, p), ...]} listing in-edges of v."""
    root = rng.choice(nodes)
    seen, stack = {root}, [root]
    while stack:
        v = stack.pop()
        for u, p in reverse.get(v, ()):
            if u not in seen and rng.random() < p:
                seen.add(u)
                stack.append(u)
    return seen


def ris_seeds(graph, nodes, k, theta, seed=7):
    rng = random.Random(seed)
    reverse = {}
    for u, outs in graph.items():
        for v, p in outs:
            reverse.setdefault(v, []).append((u, p))
    rr_sets = [random_rr_set(reverse, nodes, rng) for _ in range(theta)]
    covers = {}
    for i, rr in enumerate(rr_sets):
        for u in rr:
            covers.setdefault(u, set()).add(i)
    covered, seeds = set(), []
    for _ in range(k):
        best = max(nodes, key=lambda u: len(covers.get(u, set()) - covered))
        seeds.append(best)
        covered |= covers.get(best, set())
    # Optimistic: the same sets chose the seeds. Re-estimate on fresh sets.
    return seeds, len(nodes) * len(covered) / theta

The open question is theta. TIM (Tang, Xiao and Shi, 2014) and then IMM (Tang, Shi and Xiao, 2015) derive it from concentration bounds and achieve (1 - 1/e - epsilon) with high probability in near-linear expected time. Ship IMM or an adaptive successor such as OPIM-C from a maintained library rather than hand-tuning theta.

Running it at scale

On a real graph the cost of RIS is dominated by RR-set memory, not CPU. Practical guidance:

• Store RR sets as flat integer arrays with an offsets array, CSR style, not as Python sets.
• Generate RR sets in parallel with independent random streams per worker.
• Max coverage over the collection is a single pass with a bucket queue keyed by current coverage count; avoid the O(n) max scan in the toy code above.
• Estimate the final spread on a fresh, independent batch of RR sets or Monte Carlo simulations. The coverage measured on the selection sample is biased upward.
• Precompute the reverse adjacency list once.

Where the edge probabilities come from

Every number above depends on edge probabilities that nobody gives you. They are learned from action logs: if v acted soon after u, that is evidence for p(u, v). But friends adopt the same things because they are similar, not only because they persuade each other, so homophily inflates learned probabilities, and the greedy guarantee is only relative to the model. Benchmark settings such as uniform p = 0.01 or p(u, v) = 1 / indegree(v) are fine for comparing algorithms and misleading as a description of any real campaign.

Failure modes

• Ranking by centrality. Degree, PageRank or betweenness ignore overlap, exactly as the worked example shows; top-k by centrality often clusters seeds in one dense community.
• Too few simulations. With small R, greedy chases noise and picks a different seed set every run. Check stability by rerunning with another random seed.
• Biased spread report. Quoting the RIS coverage estimate from the selection sample overstates the result; report a held-out estimate.
• Static graph, moving world. Networks and probabilities drift; a seed set chosen on last quarter's graph decays. Recompute on a schedule.
• Ignoring cost and fairness. Seeds differ in price and the plain objective may leave whole communities unreached. Budgeted and fairness-aware variants exist; budgeted versions need a cost-benefit greedy plus a best-single-node fallback to keep a guarantee.
• Competition. With a rival cascade running, the single-cascade model no longer describes outcomes.

Trade-offs

What to do next

1. Implement simulate_ic and reproduce the worked example: sigma({A, E}) near 5.5 and sigma({A, X}) near 4.2.
2. Add CELF and count spread evaluations against naive greedy on a few thousand nodes.
3. Implement RIS, then compare its seeds with CELF's on the same graph and check that a held-out estimate agrees.
4. Decide how your edge probabilities will be learned, and write down the confounders before trusting them.
5. For production, adopt a maintained IMM-family implementation, store RR sets in CSR form, and schedule recomputation as the graph changes.
6. Read more on greedy algorithms and their proofs, BFS over implicit graphs, union-find for reachability in sampled graphs and reservoir sampling.