Independent set problem

Problem definition

In graph theory, an independent set is a set of vertices in a graph, no two of which are adjacent.

In the following, we are going to solve the solution space properties of the independent set problem on the Petersen graph. To start, let us define a Petersen graph instance.

using GenericTensorNetworks, Graphs

graph = Graphs.smallgraph(:petersen)

{10, 15} undirected simple Int64 graph

We can visualize this graph using the show_graph function

# set the vertex locations manually instead of using the default spring layout
rot15(a, b, i::Int) = cos(2i*π/5)*a + sin(2i*π/5)*b, cos(2i*π/5)*b - sin(2i*π/5)*a
locations = [[rot15(0.0, 60.0, i) for i=0:4]..., [rot15(0.0, 30, i) for i=0:4]...]
show_graph(graph, locations; format=:svg)

The graphical display is available in the following editors

a VSCode editor,
a Jupyter notebook,
or a Pluto notebook,

Generic tensor network representation

The independent set problem can be constructed with IndependentSet type as

iset = IndependentSet(graph)

IndependentSet{UnitWeight}(Graphs.SimpleGraphs.SimpleGraph{Int64}(15, [[2, 5, 6], [1, 3, 7], [2, 4, 8], [3, 5, 9], [1, 4, 10], [1, 8, 9], [2, 9, 10], [3, 6, 10], [4, 6, 7], [5, 7, 8]]), UnitWeight())

The tensor network representation of the independent set problem can be obtained by

problem = GenericTensorNetwork(iset; optimizer=TreeSA())

GenericTensorNetwork{IndependentSet{UnitWeight}, OMEinsum.SlicedEinsum{Int64, OMEinsum.DynamicNestedEinsum{Int64}}, Int64}(IndependentSet{UnitWeight}(Graphs.SimpleGraphs.SimpleGraph{Int64}(15, [[2, 5, 6], [1, 3, 7], [2, 4, 8], [3, 5, 9], [1, 4, 10], [1, 8, 9], [2, 9, 10], [3, 6, 10], [4, 6, 7], [5, 7, 8]]), UnitWeight()), OMEinsum.SlicedEinsum{Int64, OMEinsum.DynamicNestedEinsum{Int64}}(Int64[], 1∘5, 5∘1 -> 
├─ 5, 1∘5 -> 1∘5
│  ├─ 5
│  └─ 1∘5
└─ 1∘6, 1∘5∘6 -> 5∘1
   ├─ 6, 1∘6 -> 1∘6
   │  ├─ 6
   │  └─ 1∘6
   └─ 1∘3∘7, 5∘3∘6∘7 -> 1∘5∘6
      ├─ 1∘3∘2, 2∘7 -> 1∘3∘7
      │  ├─ 1∘2, 2∘3 -> 1∘3∘2
      │  │  ├─ 2, 2∘1 -> 1∘2
      │  │  │  ⋮
      │  │  │  
      │  │  └─ 3, 2∘3 -> 2∘3
      │  │     ⋮
      │  │     
      │  └─ 7, 2∘7 -> 2∘7
      │     ├─ 7
      │     └─ 2∘7
      └─ 5∘3∘6∘7, 3∘6∘7∘5 -> 5∘3∘6∘7
         ├─ 5∘9∘3, 6∘7∘9 -> 5∘3∘6∘7
         │  ├─ 5∘9∘4, 3∘4 -> 5∘9∘3
         │  │  ⋮
         │  │  
         │  └─ 6∘9, 7∘9 -> 6∘7∘9
         │     ⋮
         │     
         └─ 3∘6∘8, 7∘8∘5 -> 3∘6∘7∘5
            ├─ 3∘8, 6∘8 -> 3∘6∘8
            │  ⋮
            │  
            └─ 7∘10, 8∘5∘10 -> 7∘8∘5
               ⋮
               
), Dict{Int64, Int64}())

where the key word argument optimizer specifies the tensor network contraction order optimizer as a local search based optimizer TreeSA.

Here, the key word argument optimizer specifies the tensor network contraction order optimizer as a local search based optimizer TreeSA. The resulting contraction order optimized tensor network is contained in the code field of problem.

Theory (can skip)

Let $G=(V, E)$ be a graph with each vertex $v\in V$ associated with a weight $w_v$. To reduce the independent set problem on it to a tensor network contraction, we first map a vertex $v\in V$ to a label $s_v \in \{0, 1\}$ of dimension $2$, where we use $0$ ($1$) to denote a vertex absent (present) in the set. For each vertex $v$, we defined a parameterized rank-one tensor indexed by $s_v$ as

\[W(x_v^{w_v}) = \left(\begin{matrix} 1 \\ x_v^{w_v} \end{matrix}\right)\]

where $x_v$ is a variable associated with $v$. Similarly, for each edge $(u, v) \in E$, we define a matrix $B$ indexed by $s_u$ and $s_v$ as

\[B = \left(\begin{matrix} 1 & 1\\ 1 & 0 \end{matrix}\right).\]

Ideally, an optimal contraction order has a space complexity $2^{{\rm tw}(G)}$, where ${\rm tw(G)}$ is the tree-width of $G$ (or graph in the code). We can check the time, space and read-write complexities by typing

contraction_complexity(problem)

Time complexity: 2^7.977279923499916
Space complexity: 2^4.0
Read-write complexity: 2^8.675957032941747

For more information about how to improve the contraction order, please check the Performance Tips.

Solution space properties

Maximum independent set size $\alpha(G)$

We can compute solution space properties with the solve function, which takes two positional arguments, the problem instance and the wanted property.

maximum_independent_set_size = solve(problem, SizeMax())[]

4.0ₜ

Here SizeMax means finding the solution with maximum set size. The return value has Tropical type. We can get its content by typing

maximum_independent_set_size.n

4.0

Counting properties

Count all solutions and best several solutions

We can count all independent sets with the CountingAll property.

count_all_independent_sets = solve(problem, CountingAll())[]

76.0

The return value has type Float64. The counting of all independent sets is equivalent to the infinite temperature partition function

solve(problem, PartitionFunction(0.0))[]

76.0

We can count the maximum independent sets with CountingMax.

count_maximum_independent_sets = solve(problem, CountingMax())[]

(4.0, 5.0)ₜ

The return value has type CountingTropical, which contains two fields. They are n being the maximum independent set size and c being the number of the maximum independent sets.

count_maximum_independent_sets.c

5.0

Similarly, we can count independent sets of sizes $\alpha(G)$ and $\alpha(G)-1$ by feeding an integer positional argument to CountingMax.

count_max2_independent_sets = solve(problem, CountingMax(2))[]

30.0*x^3 + 5.0*x^4

The return value has type TruncatedPoly, which contains two fields. They are maxorder being the maximum independent set size and coeffs being the number of independent sets having sizes $\alpha(G)-1$ and $\alpha(G)$.

count_max2_independent_sets.coeffs

(30.0, 5.0)

Find the graph polynomial

We can count the number of independent sets at any size, which is equivalent to finding the coefficients of an independence polynomial that defined as

\[I(G, x) = \sum_{k=0}^{\alpha(G)} a_k x^k,\]

where $\alpha(G)$ is the maximum independent set size, $a_k$ is the number of independent sets of size $k$. The total number of independent sets is thus equal to $I(G, 1)$. There are 3 methods to compute a graph polynomial, :finitefield, :fft and :polynomial. These methods are introduced in the docstring of GraphPolynomial.

independence_polynomial = solve(problem, GraphPolynomial(; method=:finitefield))[]

1 + 10∙x + 30∙x² + 30∙x³ + 5∙x⁴

The return type is Polynomial.

independence_polynomial.coeffs

5-element Vector{BigInt}:
  1
 10
 30
 30
  5

Configuration properties

Find one best solution

We can use the bounded or unbounded SingleConfigMax to find one of the solutions with largest size. The unbounded (default) version uses a joint type of CountingTropical and ConfigSampler in computation, where CountingTropical finds the maximum size and ConfigSampler samples one of the best solutions. The bounded version uses the binary gradient back-propagation (see our paper) to compute the gradients. It requires caching intermediate states, but is often faster (on CPU) because it can use TropicalGEMM (see Performance Tips).

max_config = solve(problem, SingleConfigMax(; bounded=false))[]

(4.0, ConfigSampler{10, 1, 1}(0100100110))ₜ

The return value has type CountingTropical with its counting field having ConfigSampler type. The data field of ConfigSampler is a bit string that corresponds to the solution

single_solution = max_config.c.data

0100100110

This bit string should be read from left to right, with the i-th bit being 1 (0) to indicate the i-th vertex is present (absent) in the set. We can visualize this MIS with the following function.

show_graph(graph, locations; format=:svg, vertex_colors=
    [iszero(single_solution[i]) ? "white" : "red" for i=1:nv(graph)])

Enumerate all solutions and best several solutions

We can use bounded or unbounded ConfigsMax to find all solutions with largest-K set sizes. In most cases, the bounded (default) version is preferred because it can reduce the memory usage significantly.

all_max_configs = solve(problem, ConfigsMax(; bounded=true))[]

(4.0, {0101010001, 1010000011, 1001001100, 0100100110, 0010111000})ₜ

The return value has type CountingTropical, while its counting field having type ConfigEnumerator. The data field of a ConfigEnumerator instance contains a vector of bit strings.

all_max_configs.c.data

5-element Vector{StaticBitVector{10, 1}}:
 0101010001
 1010000011
 1001001100
 0100100110
 0010111000

These solutions can be visualized with the show_configs function.

show_configs(graph, locations, reshape(collect(all_max_configs.c), 1, length(all_max_configs.c)); padding_left=20)

We can use ConfigsAll to enumerate all sets satisfying the independence constraint.

all_independent_sets = solve(problem, ConfigsAll())[]

{0000000000, 0100000000, 0001000000, 0101000000, 0000000010, 0100000010, 0000000001, 0100000001, 0001000001, 0101000001, 0000000011, 0100000011, 0000000100, 0100000100, 0001000100, 0101000100, 0000000110, 0100000110, 0010000000, 0010000010, 0010000001, 0010000011, 0000001000, 0001001000, 0000001100, 0001001100, 0010001000, 0000010000, 0100010000, 0001010000, 0101010000, 0000010001, 0100010001, 0001010001, 0101010001, 0010010000, 0010010001, 0000011000, 0001011000, 0010011000, 1000000000, 1001000000, 1000000010, 1000000001, 1001000001, 1000000011, 1000000100, 1001000100, 1000000110, 1010000000, 1010000010, 1010000001, 1010000011, 1000001000, 1001001000, 1000001100, 1001001100, 1010001000, 0000100000, 0100100000, 0000100010, 0100100010, 0000100100, 0100100100, 0000100110, 0100100110, 0010100000, 0010100010, 0000101000, 0000101100, 0010101000, 0000110000, 0100110000, 0010110000, 0000111000, 0010111000}

The return value has type ConfigEnumerator.

Sample solutions

It is often difficult to store all configurations in a vector. A more clever way to store the data is using the sum product tree format.

all_independent_sets_tree = solve(problem, ConfigsAll(; tree_storage=true))[]

+ (count = 76.0)
├─ + (count = 58.0)
│  ├─ + (count = 40.0)
│  │  ├─ + (count = 27.0)
│  │  │  ├─ + (count = 26.0)
│  │  │  │  ├─ + (count = 22.0)
│  │  │  │  │  ⋮
│  │  │  │  │  
│  │  │  │  └─ * (count = 4.0)
│  │  │  │     ⋮
│  │  │  │     
│  │  │  └─ * (count = 1.0)
│  │  │     ├─ * (count = 1.0)
│  │  │     │  ⋮
│  │  │     │  
│  │  │     └─ * (count = 1.0)
│  │  │        ⋮
│  │  │        
│  │  └─ * (count = 13.0)
│  │     ├─ * (count = 1.0)
│  │     │  ├─ OnehotVec{10, 2}(6, 1)
│  │     │  └─ OnehotVec{10, 2}(6, 1)
│  │     └─ + (count = 13.0)
│  │        ├─ + (count = 12.0)
│  │        │  ⋮
│  │        │  
│  │        └─ * (count = 1.0)
│  │           ⋮
│  │           
│  └─ * (count = 18.0)
│     ├─ OnehotVec{10, 2}(1, 1)
│     └─ * (count = 18.0)
│        ├─ OnehotVec{10, 2}(1, 1)
│        └─ + (count = 18.0)
│           ├─ + (count = 17.0)
│           │  ⋮
│           │  
│           └─ * (count = 1.0)
│              ⋮
│              
└─ * (count = 18.0)
   ├─ * (count = 1.0)
   │  ├─ OnehotVec{10, 2}(5, 1)
   │  └─ OnehotVec{10, 2}(5, 1)
   └─ + (count = 18.0)
      ├─ + (count = 13.0)
      │  ├─ + (count = 12.0)
      │  │  ├─ + (count = 10.0)
      │  │  │  ⋮
      │  │  │  
      │  │  └─ * (count = 2.0)
      │  │     ⋮
      │  │     
      │  └─ * (count = 1.0)
      │     ├─ * (count = 1.0)
      │     │  ⋮
      │     │  
      │     └─ * (count = 1.0)
      │        ⋮
      │        
      └─ * (count = 5.0)
         ├─ * (count = 1.0)
         │  ├─ OnehotVec{10, 2}(6, 1)
         │  └─ OnehotVec{10, 2}(6, 1)
         └─ + (count = 5.0)
            ├─ + (count = 4.0)
            │  ⋮
            │  
            └─ * (count = 1.0)
               ⋮

The return value has the SumProductTree type. Its length corresponds to the number of configurations.

length(all_independent_sets_tree)

76.0

We can use Base.collect function to create a ConfigEnumerator or use generate_samples to generate samples from it.

collect(all_independent_sets_tree)

generate_samples(all_independent_sets_tree, 10)

10-element Vector{StaticBitVector{10, 1}}:
 0001000100
 0100000110
 0010000011
 0001001000
 0000010001
 0001011000
 1001000001
 1001000100
 1000001100
 0000100110

This page was generated using Literate.jl.