Stochastic Shape Models

This tutorial introduces Stochastic Shape Models (SSMs), an extension of Stochastic Block Models that reduces the number of parameters while improving the flexibility to model complex network structures.

SSMs address a fundamental challenge in network modeling: standard SBMs with k blocks require k² parameters (one for each block pair), which grows quickly. SSMs reduce this by using a small number of "shapes" (connectivity patterns) that are shared across multiple block pairs.

What is a Stochastic Shape Model?

A Stochastic Shape Model is defined by three components:

1. Shape parameters θ: A vector of S edge probabilities (or distributions), where S << k²
2. Shape assignment matrix Ψ: A k×k matrix that assigns each block pair (i,j) to one of the S shapes
3. Block sizes: The proportion of nodes in each block (same as in SBMs)

The key insight: many block pairs can share the same connectivity pattern!

Mathematical Formulation

Definition 2.1 (Stochastic Shape Model $(S S M)$ ). Assume we have defined $s \in \mathbb{N}^{+}$nonintersecting closed regions in $\mathcal{T}=[0,1]^2 \cap\{x \leq y\}$, let us call them $S_c$ for $c \in[s]$. Define $S_s=\mathcal{T} \backslash\left\{\cup_{c<s} S_c\right\}$. We can then define the function $f$ for the $s$ constants $0<\theta_c<1$, for $c \in[s]$ to be

\[f(x, y)=\left\{\begin{array}{lll} \theta_c & \text { if } & (x, y) \in S_c \\ \theta_c & \text { if } & (y, x) \in S_c \end{array} .\right.\]

For nodes i and j with latent positions $ξ_i$ and $ξ_j$

\[\mathbb{P}(A_{i,j} = 1) = f(ξ_i, ξ_j) = \theta_{c} \quad \text{if } (ξ_i, ξ_j) \in S_c\]

Instead of k² unique probabilities, we only need S parameters where S << k².

Motivation: Parameter Efficiency

Consider a k=10 block network:

• Standard SBM: Requires 100 parameters (10×10 matrix)
• SSM with 5 shapes: Requires only 5 parameters + assignment matrix

This makes SSMs particularly useful for:

• Large networks with many blocks but limited data
• Hierarchical structures where similar patterns repeat
• Model selection with limited samples to avoid overfitting

Setup

Load the packages we'll need:

using Graphons
using Random
using CairoMakie
using Distributions

Random.seed!(42)

TaskLocalRNG()

Example 1: Simple Bipartite Structure

Let's start with a simple example: a network with 3 blocks where we want:

• High connectivity within blocks (assortative)
• Low connectivity between blocks (sparse)

Using only 2 shapes instead of 9 SBM parameters:

Define the shapes:

θ = [
    0.8,  # Shape 1: high probability (within-block)
    0.1,
]  # Shape 2: low probability (between-block)

2-element Vector{Float64}:
 0.8
 0.1

Assign shapes to block pairs:

block_pair_to_shape = [
    1 2 2;   # Block 1: high within, low to others
    2 1 2;   # Block 2: high within, low to others
    2 2 1
]   # Block 3: high within, low to others

3×3 Matrix{Int64}:
 1  2  2
 2  1  2
 2  2  1

Block sizes:

sizes = [0.3, 0.4, 0.3]

3-element Vector{Float64}:
 0.3
 0.4
 0.3

Create the SSM:

ssm_simple = SSM(θ, block_pair_to_shape, sizes, cumsum(sizes))

SSM{Bool, BitMatrix, Vector{Float64}, Matrix{Int64}, Vector{Float64}, Vector{Float64}}([0.8, 0.1], [1 2 2; 2 1 2; 2 2 1], [0.3, 0.4, 0.3], [0.3, 0.7, 1.0])

Let's visualize both the shape assignment and a sampled graph:

fig = Figure(size = (900, 400))

ax1 = Axis(
    fig[1, 1],
    title = "Shape Assignment Matrix\n(which shape for each block pair)",
    xlabel = "Block",
    ylabel = "Block",
    aspect = 1,
)

ax2 = Axis(
    fig[1, 2],
    title = "Equivalent Block Probability Matrix\n(θ values)",
    xlabel = "Block",
    ylabel = "Block",
    aspect = 1,
)

ax3 = Axis(fig[1, 3], title = "Sampled Graph (n=200)", aspect = 1)

Axis with 0 plots:

Visualize shape assignments

heatmap!(ax1, block_pair_to_shape)

Makie.Heatmap{Tuple{Makie.EndPoints{Float32}, Makie.EndPoints{Float32}, Matrix{Float32}}}

Visualize the equivalent probability matrix using direct plotting

hm = heatmap!(ax2, ssm_simple, colormap = :binary, colorrange = (0, 1))

Makie.Heatmap{Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Matrix{Float64}}}

Sample and visualize graph

hidedecorations!(ax3)
A_simple = sample_graph(ssm_simple, 200)
heatmap!(ax3, A_simple, colormap = :binary)

Colorbar(fig[1, 4], hm, label = "Edge probability")

fig

Key observation: We achieved the same result as a 3-block SBM but using only 2 parameters instead of 9! The shape assignment matrix shows that diagonal blocks use shape 1 (high connectivity) while off-diagonal use shape 2 (low connectivity).

Example 2: Hierarchical Core-Periphery

SSMs excel at modeling hierarchical structures. Let's create a 6-block network with:

• A core of 2 blocks (densely connected)
• A periphery of 4 blocks (sparsely connected)
• Medium connectivity between core and periphery

We'll use only 3 shapes:

θ_hier = [
    0.9,  # Shape 1: dense (core-core)
    0.5,  # Shape 2: medium (core-periphery)
    0.1,
]  # Shape 3: sparse (periphery-periphery)

3-element Vector{Float64}:
 0.9
 0.5
 0.1

Create the shape assignment matrix for hierarchical structure:

block_pair_to_shape_hier = [
    1 1 2 2 2 2;   # Core block 1
    1 1 2 2 2 2;   # Core block 2
    2 2 3 3 3 3;   # Periphery block 1
    2 2 3 3 3 3;   # Periphery block 2
    2 2 3 3 3 3;   # Periphery block 3
    2 2 3 3 3 3
]

sizes_hier = [0.15, 0.15, 0.175, 0.175, 0.175, 0.175]  # Core smaller than periphery

ssm_hier = SSM(θ_hier, block_pair_to_shape_hier, sizes_hier, cumsum(sizes_hier))

SSM{Bool, BitMatrix, Vector{Float64}, Matrix{Int64}, Vector{Float64}, Vector{Float64}}([0.9, 0.5, 0.1], [1 1 … 2 2; 1 1 … 2 2; … ; 2 2 … 3 3; 2 2 … 3 3], [0.15, 0.15, 0.175, 0.175, 0.175, 0.175], [0.15, 0.3, 0.475, 0.6499999999999999, 0.825, 0.9999999999999999])

Visualize the hierarchical structure:

fig = Figure(size = (1100, 400))

ax1 = Axis(
    fig[1, 1],
    title = "Shape Assignment Matrix\n(6 blocks, 3 shapes)",
    xlabel = "Block",
    ylabel = "Block",
    aspect = 1,
)

ax2 = Axis(
    fig[1, 2],
    title = "Equivalent 6×6 Probability Matrix",
    xlabel = "Block",
    ylabel = "Block",
    aspect = 1,
)

ax3 = Axis(fig[1, 3], title = "Sampled Graph (n=300)", aspect = 1)

heatmap!(ax1, block_pair_to_shape_hier, colormap = :tab10, colorrange = (1, 3))

hm = heatmap!(ax2, ssm_hier, colormap = :binary, colorrange = (0, 1))

hidedecorations!(ax3)
A_hier = sample_graph(ssm_hier, 300)
heatmap!(ax3, A_hier, colormap = :binary)

Colorbar(fig[1, 4], hm, label = "Edge probability")

fig

Key advantage: This 6-block hierarchical model uses only 3 parameters instead of 36! The clear two-level hierarchy (dense core in top-left, sparse periphery in bottom-right) is efficiently captured by shape sharing.

Example 3: Modular Networks with Weak Ties

Real networks often have multiple dense communities connected by sparse "weak ties". Let's model a 4-community network where:

• Communities 1-2 form one meta-module
• Communities 3-4 form another meta-module
• Within-module connections are stronger than between-module

Using 4 shapes:

θ_modular = [
    0.9,  # Shape 1: within-community
    0.6,  # Shape 2: within-module, between-community
    0.2,  # Shape 3: between-module
    0.05,
] # Shape 4: no connection (rare edges)

block_pair_to_shape_modular = [
    1 2 3 4;   # Community 1
    2 1 4 3;   # Community 2
    3 4 1 2;   # Community 3
    4 3 2 1
]

sizes_modular = [0.25, 0.25, 0.25, 0.25]

ssm_modular =
    SSM(θ_modular, block_pair_to_shape_modular, sizes_modular, cumsum(sizes_modular))

SSM{Bool, BitMatrix, Vector{Float64}, Matrix{Int64}, Vector{Float64}, Vector{Float64}}([0.9, 0.6, 0.2, 0.05], [1 2 3 4; 2 1 4 3; 3 4 1 2; 4 3 2 1], [0.25, 0.25, 0.25, 0.25], [0.25, 0.5, 0.75, 1.0])

Visualize:

fig = Figure(size = (1100, 400))

ax1 = Axis(
    fig[1, 1],
    title = "Shape Assignment\n(4 communities, 4 shapes)",
    xlabel = "Block",
    ylabel = "Block",
    aspect = 1,
)

ax2 = Axis(
    fig[1, 2],
    title = "Probability Matrix",
    xlabel = "Block",
    ylabel = "Block",
    aspect = 1,
)

ax3 = Axis(fig[1, 3], title = "Sampled Graph (n=200)", aspect = 1)

heatmap!(ax1, block_pair_to_shape_modular, colormap = :tab10, colorrange = (1, 4))

hm = heatmap!(ax2, ssm_modular, colormap = :binary, colorrange = (0, 1))

hidedecorations!(ax3)
A_modular = sample_graph(ssm_modular, 200)
heatmap!(ax3, A_modular, colormap = :binary)

Colorbar(fig[1, 4], hm, label = "Edge probability")

fig

Notice the symmetric checkerboard pattern! This reflects the two meta-modules. The SSM uses only 4 parameters to capture this complex structure, compared to 16 parameters for an equivalent SBM.

Converting SSMs to SBMs

Every SSM can be converted to an equivalent SBM by expanding the shape assignments:

sbm_from_ssm = Graphons.convert_to_sbm(ssm_modular)

println("SSM type: ", typeof(ssm_modular))
println("SSM parameters: ", length(ssm_modular.θ), " shapes")
println("\nConverted SBM type: ", typeof(sbm_from_ssm))
println("SBM matrix size: ", size(sbm_from_ssm.θ))
println("SBM parameters: ", length(sbm_from_ssm.θ), " entries")

SSM type: SSM{Bool, BitMatrix, Vector{Float64}, Matrix{Int64}, Vector{Float64}, Vector{Float64}}
SSM parameters: 4 shapes

Converted SBM type: SBM{Matrix{Float64}, Vector{Float64}, Vector{Float64}}
SBM matrix size: (4, 4)
SBM parameters: 16 entries

The SSM and SBM are mathematically equivalent, but the SSM representation is more compact and interpretable when many block pairs share patterns.

Decorated SSMs: Adding Edge Attributes

Just like SBMs can be decorated with distributions, SSMs can too! This combines parameter efficiency with rich edge attributes.

Let's create an SSM with weighted edges:

Define shape distributions (instead of probabilities):

θ_weighted = [
    Normal(5.0, 0.5),   # Shape 1: strong positive connections
    Normal(0.0, 0.3),   # Shape 2: weak/zero connections
    Normal(-2.0, 0.5),   # Shape 3: negative/inhibitory connections
]

3-element Vector{Distributions.Normal{Float64}}:
 Distributions.Normal{Float64}(μ=5.0, σ=0.5)
 Distributions.Normal{Float64}(μ=0.0, σ=0.3)
 Distributions.Normal{Float64}(μ=-2.0, σ=0.5)

3-block network with signed edges:

block_pair_to_shape_weighted = [
    1 2 3;   # Block 1: positive within, weak to 2, negative to 3
    2 1 2;   # Block 2: positive within, weak elsewhere
    3 2 1
]

sizes_weighted = [0.35, 0.3, 0.35]

ssm_weighted =
    SSM(θ_weighted, block_pair_to_shape_weighted, sizes_weighted, cumsum(sizes_weighted))

SSM{Float64, Matrix{Float64}, Vector{Distributions.Normal{Float64}}, Matrix{Int64}, Vector{Float64}, Vector{Float64}}(Distributions.Normal{Float64}[Distributions.Normal{Float64}(μ=5.0, σ=0.5), Distributions.Normal{Float64}(μ=0.0, σ=0.3), Distributions.Normal{Float64}(μ=-2.0, σ=0.5)], [1 2 3; 2 1 2; 3 2 1], [0.35, 0.3, 0.35], [0.35, 0.6499999999999999, 0.9999999999999999])

Sample a weighted graph:

A_weighted = sample_graph(ssm_weighted, 150)

println("Weighted matrix type: ", typeof(A_weighted))
println("Matrix size: ", size(A_weighted))
println("Value range: [", minimum(A_weighted), ", ", maximum(A_weighted), "]")

Weighted matrix type: Matrix{Float64}
Matrix size: (150, 150)
Value range: [-3.6264762062670908, 6.883325873524579]

Visualize the weighted network:

fig = Figure(size = (1100, 400))

ax1 = Axis(
    fig[1, 1],
    title = "Shape Assignment",
    xlabel = "Block",
    ylabel = "Block",
    aspect = 1,
)

ax2 = Axis(
    fig[1, 2],
    title = "Shape Distributions\n(means)",
    xlabel = "Block",
    ylabel = "Block",
    aspect = 1,
)

ax3 = Axis(fig[1, 3], title = "Sampled Weighted Graph (n=150)", aspect = 1)

heatmap!(ax1, block_pair_to_shape_weighted, colormap = :tab10, colorrange = (1, 3))

Makie.Heatmap{Tuple{Makie.EndPoints{Float32}, Makie.EndPoints{Float32}, Matrix{Float32}}}

Show mean of each distribution

θ_means = [mean(d) for d in θ_weighted]
θ_matrix_weighted = [θ_means[block_pair_to_shape_weighted[i, j]] for i = 1:3, j = 1:3]
hm_means = heatmap!(ax2, θ_matrix_weighted, colormap = :RdBu, colorrange = (-3, 6))

hidedecorations!(ax3)
hm_weights = heatmap!(ax3, A_weighted, colormap = :RdBu, colorrange = (-4, 7))

Colorbar(fig[1, 4], hm_weights, label = "Edge weight")

fig

Interpretation: The SSM creates three distinct types of connections:

• Diagonal (shape 1): Strong positive weights (red)
• Off-diagonal corners (shape 3): Negative weights (blue)
• Other pairs (shape 2): Weak weights (white)

Example 4: Degree-Corrected SSM

We can create more realistic networks by varying block sizes to simulate degree heterogeneity while keeping the shape structure:

Same shapes as before, but with very unequal block sizes:

θ_dc = [0.9, 0.3, 0.05]
block_pair_to_shape_dc = [
    1 2 3;
    2 1 2;
    3 2 1
]
sizes_dc = [0.6, 0.3, 0.1]  # Highly unequal blocks

ssm_dc = SSM(θ_dc, block_pair_to_shape_dc, sizes_dc, cumsum(sizes_dc))

SSM{Bool, BitMatrix, Vector{Float64}, Matrix{Int64}, Vector{Float64}, Vector{Float64}}([0.9, 0.3, 0.05], [1 2 3; 2 1 2; 3 2 1], [0.6, 0.3, 0.1], [0.6, 0.8999999999999999, 1.0])

Compare equal vs unequal sizes:

fig = Figure(size = (900, 400))

ax1 = Axis(fig[1, 1], title = "Equal Block Sizes\n[0.33, 0.33, 0.33]", aspect = 1)

ax2 = Axis(fig[1, 2], title = "Unequal Block Sizes\n[0.6, 0.3, 0.1]", aspect = 1)

hidedecorations!.([ax1, ax2])

2-element BitVector:
 0
 0

Sample with equal sizes

ssm_equal =
    SSM(θ_dc, block_pair_to_shape_dc, [0.33, 0.33, 0.34], cumsum([0.33, 0.33, 0.34]))
A_equal = sample_graph(ssm_equal, 300)
heatmap!(ax1, A_equal, colormap = :binary)

Makie.Heatmap{Tuple{Makie.EndPoints{Float32}, Makie.EndPoints{Float32}, Matrix{Float32}}}

Sample with unequal sizes

A_unequal = sample_graph(ssm_dc, 300)
heatmap!(ax2, A_unequal, colormap = :binary)

fig

Analyzing Shape Assignments

Let's analyze which block pairs are assigned to which shapes:

function analyze_shape_assignments(ssm)
    K = length(ssm.size)
    S = length(ssm.θ)

    println("Number of blocks: ", K)
    println("Number of shapes: ", S)
    println("\nShape assignments:")

    for s = 1:S
        pairs = [(i, j) for i = 1:K for j = 1:K if ssm.block_pair_to_shape[i, j] == s]
        n_pairs = length(pairs)
        println("\nShape $s (θ = $(ssm.θ[s])):")
        println("  Used by $n_pairs block pairs ($(round(100*n_pairs/K^2, digits=1))%)")
        if n_pairs ≤ 6
            println("  Pairs: ", pairs)
        end
    end
end

analyze_shape_assignments(ssm_modular)

Number of blocks: 4
Number of shapes: 4

Shape assignments:

Shape 1 (θ = 0.9):
  Used by 4 block pairs (25.0%)
  Pairs: [(1, 1), (2, 2), (3, 3), (4, 4)]

Shape 2 (θ = 0.6):
  Used by 4 block pairs (25.0%)
  Pairs: [(1, 2), (2, 1), (3, 4), (4, 3)]

Shape 3 (θ = 0.2):
  Used by 4 block pairs (25.0%)
  Pairs: [(1, 3), (2, 4), (3, 1), (4, 2)]

Shape 4 (θ = 0.05):
  Used by 4 block pairs (25.0%)
  Pairs: [(1, 4), (2, 3), (3, 2), (4, 1)]

This helps understand how parameter sharing works in the model.

References

Stochastic Shape Models were introduced in (Verdeyme and Olhede, 2024)

The paper provides:

• Theoretical properties (consistency, rates of convergence)
• Model selection procedures for choosing S
• Applications to real network data
• Comparisons with standard SBMs and other graphon models

Real-World Example: Multiplex Network from Research

Let's recreate an example from recent research on decorated graphon estimation (Dufour & Olhede, 2024). This example demonstrates how SSMs can provide smoother approximations than standard SBMs for complex multiplex networks.

The graphon w3 models a 4-category multiplex network where edge probabilities are determined by a softmax transformation of four distinct spatial patterns:

using LogExpFunctions: softmax!
function w3(x, y)
    tabulation = zeros(4)
    tabulation[1] = 3 * x * y                                      # Linear interaction
    tabulation[2] = 3 * sin(2 * π * x) * sin(2 * π * y)           # Periodic pattern
    tabulation[3] = exp(-3 * ((x - 0.5)^2 + (y - 0.5)^2))        # Gaussian bump
    tabulation[4] = 2 - 3 * (x + y)                               # Decreasing pattern
    softmax!(tabulation)
    return DiscreteNonParametric(0:3, tabulation)
end

w3 (generic function with 1 method)

Create the decorated graphon:

graphon_w3 = DecoratedGraphon(w3)

DecoratedGraphon{Int64, Matrix{Int64}, typeof(Main.w3), Distributions.DiscreteNonParametric{Int64, Float64, UnitRange{Int64}, Vector{Float64}}}(Main.w3)

Visualize the four category probabilities:

fig = Figure(size = (1000, 250))

for i = 1:4
    ax = Axis(fig[1, i], title = "Category $i", aspect = 1)
    hidedecorations!(ax)
    hm = heatmap!(ax, graphon_w3, k = i, colormap = :binary, colorrange = (0, 1))
end

Colorbar(
    fig[1, 5],
    colormap = :binary,
    colorrange = (0, 1),
    label = "Probability",
    height = Relative(0.8),
)

fig

Comparing SBM vs SSM Approximations

k = 10

10

Create a standard SBM approximation with k=10 blocks:

sbm_w3 = discretized_graphon(graphon_w3, k)

DecoratedSBM{Distributions.DiscreteNonParametric{Int64, Float64, UnitRange{Int64}, Vector{Float64}}, Matrix{Int64}, Matrix{Distributions.DiscreteNonParametric{Int64, Float64, UnitRange{Int64}, Vector{Float64}}}, Vector{Float64}, Vector{Float64}}(Distributions.DiscreteNonParametric{Int64, Float64, UnitRange{Int64}, Vector{Float64}}[Distributions.DiscreteNonParametric{Int64, Float64, UnitRange{Int64}, Vector{Float64}}(support=0:3, p=[0.09399344867091115, 0.09399344867091115, 0.11749023749685668, 0.694522865161321]) Distributions.DiscreteNonParametric{Int64, Float64, UnitRange{Int64}, Vector{Float64}}(support=0:3, p=[0.11568105451872908, 0.11568105451872908, 0.15616569878524947, 0.6124721921772922]) … Distributions.DiscreteNonParametric{Int64, Float64, UnitRange{Int64}, Vector{Float64}}(support=0:3, p=[0.2588403746197946, 0.2588403746197946, 0.349426171333799, 0.1328930794266118]) Distributions.DiscreteNonParametric{Int64, Float64, UnitRange{Int64}, Vector{Float64}}(support=0:3, p=[0.2764062879879321, 0.2764062879879321, 0.3455032332628625, 0.1016841907612732]); Distributions.DiscreteNonParametric{Int64, Float64, UnitRange{Int64}, Vector{Float64}}(support=0:3, p=[0.11568105451872908, 0.11568105451872908, 0.15616569878524947, 0.6124721921772922]) Distributions.DiscreteNonParametric{Int64, Float64, UnitRange{Int64}, Vector{Float64}}(support=0:3, p=[0.10607996027902702, 0.35307618262717355, 0.15304397130799247, 0.38779988578580693]) … Distributions.DiscreteNonParametric{Int64, Float64, UnitRange{Int64}, Vector{Float64}}(support=0:3, p=[0.38431100461229367, 0.08273378991798222, 0.4278296043585822, 0.10512560111114182]) Distributions.DiscreteNonParametric{Int64, Float64, UnitRange{Int64}, Vector{Float64}}(support=0:3, p=[0.34810443597479224, 0.24942772772556748, 0.3367194011044718, 0.06574843519516833]); … ; Distributions.DiscreteNonParametric{Int64, Float64, UnitRange{Int64}, Vector{Float64}}(support=0:3, p=[0.2588403746197946, 0.2588403746197946, 0.349426171333799, 0.1328930794266118]) Distributions.DiscreteNonParametric{Int64, Float64, UnitRange{Int64}, Vector{Float64}}(support=0:3, p=[0.38431100461229367, 0.08273378991798223, 0.4278296043585822, 0.10512560111114182]) … Distributions.DiscreteNonParametric{Int64, Float64, UnitRange{Int64}, Vector{Float64}}(support=0:3, p=[0.682128970276014, 0.22016474545155387, 0.0954323419246286, 0.0022739423478034767]) Distributions.DiscreteNonParametric{Int64, Float64, UnitRange{Int64}, Vector{Float64}}(support=0:3, p=[0.8583249103622071, 0.05963937704246795, 0.08051123867863832, 0.0015244739166864606]); Distributions.DiscreteNonParametric{Int64, Float64, UnitRange{Int64}, Vector{Float64}}(support=0:3, p=[0.2764062879879321, 0.2764062879879321, 0.3455032332628625, 0.1016841907612732]) Distributions.DiscreteNonParametric{Int64, Float64, UnitRange{Int64}, Vector{Float64}}(support=0:3, p=[0.34810443597479224, 0.24942772772556748, 0.3367194011044718, 0.06574843519516833]) … Distributions.DiscreteNonParametric{Int64, Float64, UnitRange{Int64}, Vector{Float64}}(support=0:3, p=[0.8583249103622071, 0.05963937704246795, 0.08051123867863832, 0.0015244739166864606]) Distributions.DiscreteNonParametric{Int64, Float64, UnitRange{Int64}, Vector{Float64}}(support=0:3, p=[0.8985275315618875, 0.04473505164427972, 0.05591806574224748, 0.0008193510515855016])], [0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.10000000000000012], [0.1, 0.2, 0.30000000000000004, 0.4, 0.5, 0.6, 0.7, 0.7999999999999999, 0.8999999999999999, 1.0])

For the SSM, we'll use fewer shapes (s=30) than the SBM parameters (k(k+1)/2=120): This demonstrates the parameter efficiency of SSMs. In the Graphons package, to automatically create an SSM from an SBM, we can load the Clustering package and use the SSM constructor:

using Clustering

s = 10
ssm_w3 = SSM(sbm_w3, s)

SSM{Int64, Matrix{Int64}, Vector{Distributions.DiscreteNonParametric{Int64, Float64, UnitRange{Int64}, Vector{Float64}}}, Matrix{Int64}, Vector{Float64}, Vector{Float64}}(Distributions.DiscreteNonParametric{Int64, Float64, UnitRange{Int64}, Vector{Float64}}[Distributions.DiscreteNonParametric{Int64, Float64, UnitRange{Int64}, Vector{Float64}}(support=0:3, p=[0.37701308300949044, 0.057300190569311134, 0.4561196455570192, 0.10956708086417917]), Distributions.DiscreteNonParametric{Int64, Float64, UnitRange{Int64}, Vector{Float64}}(support=0:3, p=[0.4524869259329671, 0.37634821303851523, 0.16321442440576264, 0.007950436622755045]), Distributions.DiscreteNonParametric{Int64, Float64, UnitRange{Int64}, Vector{Float64}}(support=0:3, p=[0.7969230246086652, 0.10056338589251304, 0.10003532323851312, 0.002478266260308711]), Distributions.DiscreteNonParametric{Int64, Float64, UnitRange{Int64}, Vector{Float64}}(support=0:3, p=[0.07289706037797436, 0.7287557657458588, 0.1206950205933855, 0.07765215328278134]), Distributions.DiscreteNonParametric{Int64, Float64, UnitRange{Int64}, Vector{Float64}}(support=0:3, p=[0.34512956263361255, 0.24503678429958642, 0.35147280795456637, 0.05836084511223467]), Distributions.DiscreteNonParametric{Int64, Float64, UnitRange{Int64}, Vector{Float64}}(support=0:3, p=[0.16160447134760791, 0.22335352693464455, 0.251871544036227, 0.3631704576815205]), Distributions.DiscreteNonParametric{Int64, Float64, UnitRange{Int64}, Vector{Float64}}(support=0:3, p=[0.2170192049020255, 0.2895571048610721, 0.32996016789075067, 0.16346352234615177]), Distributions.DiscreteNonParametric{Int64, Float64, UnitRange{Int64}, Vector{Float64}}(support=0:3, p=[0.18215213268469374, 0.6042038835195751, 0.12730502413200678, 0.08633895966372432]), Distributions.DiscreteNonParametric{Int64, Float64, UnitRange{Int64}, Vector{Float64}}(support=0:3, p=[0.5545682296850535, 0.1276056596438109, 0.2951546201342859, 0.022671490536849744]), Distributions.DiscreteNonParametric{Int64, Float64, UnitRange{Int64}, Vector{Float64}}(support=0:3, p=[0.1158792730696676, 0.1158792730696676, 0.15811410053889863, 0.610127353321766])], [10 10 … 7 5; 10 6 … 1 5; … ; 7 1 … 3 3; 5 5 … 3 3], [0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.10000000000000012], [0.1, 0.2, 0.30000000000000004, 0.4, 0.5, 0.6, 0.7, 0.7999999999999999, 0.8999999999999999, 1.0])

We can now visualize the difference between the SBM and SSM approximations:

fig = Figure(size = (1000, 500))

for i = 1:4
    ax = Axis(fig[1:2, i], title = "Category $i", aspect = 1)
    hidedecorations!(ax)
    hm = heatmap!(ax, sbm_w3, k = i, colormap = :binary, colorrange = (0, 1))
    ax2 = Axis(fig[3:4, i], aspect = 1)
    hidedecorations!(ax2)
    hm = heatmap!(ax2, ssm_w3, k = i, colormap = :binary, colorrange = (0, 1))
end

Colorbar(
    fig[2:3, 5],
    colormap = :binary,
    colorrange = (0, 1),
    label = "Probability",
    height = Relative(0.8),
)

fig

We can also automatically estimate the optimal number of shapes for the SSM. First let's build a larger SBM approximation to have a better starting point:

k_big = 22
sbm_big = discretized_graphon(graphon_w3, k_big)

DecoratedSBM{Distributions.DiscreteNonParametric{Int64, Float64, UnitRange{Int64}, Vector{Float64}}, Matrix{Int64}, Matrix{Distributions.DiscreteNonParametric{Int64, Float64, UnitRange{Int64}, Vector{Float64}}}, Vector{Float64}, Vector{Float64}}(Distributions.DiscreteNonParametric{Int64, Float64, UnitRange{Int64}, Vector{Float64}}[Distributions.DiscreteNonParametric{Int64, Float64, UnitRange{Int64}, Vector{Float64}}(support=0:3, p=[0.09399344867091115, 0.09399344867091115, 0.11749023749685668, 0.694522865161321]) Distributions.DiscreteNonParametric{Int64, Float64, UnitRange{Int64}, Vector{Float64}}(support=0:3, p=[0.10312775922129815, 0.10312775922129815, 0.13316895992199376, 0.6605755216354099]) … Distributions.DiscreteNonParametric{Int64, Float64, UnitRange{Int64}, Vector{Float64}}(support=0:3, p=[0.2691302009549929, 0.2691302009549929, 0.3475280488531343, 0.11421154923687978]) Distributions.DiscreteNonParametric{Int64, Float64, UnitRange{Int64}, Vector{Float64}}(support=0:3, p=[0.2764062879879321, 0.2764062879879321, 0.3455032332628625, 0.1016841907612732]); Distributions.DiscreteNonParametric{Int64, Float64, UnitRange{Int64}, Vector{Float64}}(support=0:3, p=[0.10312775922129815, 0.10312775922129815, 0.13316895992199376, 0.6605755216354099]) Distributions.DiscreteNonParametric{Int64, Float64, UnitRange{Int64}, Vector{Float64}}(support=0:3, p=[0.10946593871355442, 0.14109769973160852, 0.14572486320310038, 0.6037114983517368]) … Distributions.DiscreteNonParametric{Int64, Float64, UnitRange{Int64}, Vector{Float64}}(support=0:3, p=[0.3161107879176657, 0.2125966028102219, 0.36979468730413045, 0.10149792196798188]) Distributions.DiscreteNonParametric{Int64, Float64, UnitRange{Int64}, Vector{Float64}}(support=0:3, p=[0.3064917006623227, 0.2656908817610156, 0.3430868531618862, 0.08473056441477554]); … ; Distributions.DiscreteNonParametric{Int64, Float64, UnitRange{Int64}, Vector{Float64}}(support=0:3, p=[0.2691302009549929, 0.2691302009549929, 0.3475280488531343, 0.11421154923687978]) Distributions.DiscreteNonParametric{Int64, Float64, UnitRange{Int64}, Vector{Float64}}(support=0:3, p=[0.3161107879176657, 0.2125966028102219, 0.36979468730413045, 0.10149792196798188]) … Distributions.DiscreteNonParametric{Int64, Float64, UnitRange{Int64}, Vector{Float64}}(support=0:3, p=[0.850920495902962, 0.07266589123892432, 0.07504889931208018, 0.0013647135460335927]) Distributions.DiscreteNonParametric{Int64, Float64, UnitRange{Int64}, Vector{Float64}}(support=0:3, p=[0.8827614572097703, 0.05069930267605584, 0.06546805105745343, 0.0010711890567205327]); Distributions.DiscreteNonParametric{Int64, Float64, UnitRange{Int64}, Vector{Float64}}(support=0:3, p=[0.2764062879879321, 0.2764062879879321, 0.3455032332628625, 0.1016841907612732]) Distributions.DiscreteNonParametric{Int64, Float64, UnitRange{Int64}, Vector{Float64}}(support=0:3, p=[0.3064917006623227, 0.2656908817610156, 0.3430868531618862, 0.08473056441477554]) … Distributions.DiscreteNonParametric{Int64, Float64, UnitRange{Int64}, Vector{Float64}}(support=0:3, p=[0.8827614572097703, 0.05069930267605584, 0.06546805105745343, 0.0010711890567205327]) Distributions.DiscreteNonParametric{Int64, Float64, UnitRange{Int64}, Vector{Float64}}(support=0:3, p=[0.8985275315618875, 0.04473505164427972, 0.05591806574224748, 0.0008193510515855016])], [0.045454545454545456, 0.045454545454545456, 0.045454545454545456, 0.045454545454545456, 0.045454545454545456, 0.045454545454545456, 0.045454545454545456, 0.045454545454545456, 0.045454545454545456, 0.045454545454545456  …  0.045454545454545456, 0.045454545454545456, 0.045454545454545456, 0.045454545454545456, 0.045454545454545456, 0.045454545454545456, 0.045454545454545456, 0.045454545454545456, 0.045454545454545456, 0.04545454545454557], [0.045454545454545456, 0.09090909090909091, 0.13636363636363635, 0.18181818181818182, 0.2272727272727273, 0.27272727272727276, 0.31818181818181823, 0.3636363636363637, 0.40909090909090917, 0.45454545454545464  …  0.5909090909090909, 0.6363636363636364, 0.6818181818181818, 0.7272727272727272, 0.7727272727272726, 0.818181818181818, 0.8636363636363634, 0.9090909090909088, 0.9545454545454543, 0.9999999999999998])

Now let's sample a graph to be able to compute a criterion for model selection:

latents = collect(rand(Uniform(0, 1), 1000))
A = sample_graph(graphon_w3, latents);

We are now ready to estimate the SSM with model selection over a range of shapes:

shape_range = 1:(k_big*(k_big+1)÷2-1)
ssm_estimated, criterion_values = Graphons.estimate_ssm(sbm_big, A, latents, shape_range)
index_argmin = argmin(criterion_values)
k_opt = shape_range[index_argmin]
println("Optimal number of shapes selected by argmin: $k_opt")

Optimal number of shapes selected by argmin: 58

We can plot the BIC values to find the number of shapes that minimizes the criterion:

fig = Figure(size = (600, 200))
ax = Axis(fig[1, 1], xlabel = "Number of Shapes", ylabel = "BIC Value")
lines!(ax, shape_range, criterion_values)
scatter!(
    ax,
    [k_opt],
    [criterion_values[index_argmin]],
    marker = :rect,
    color = :red,
    label = "potential elbow",
)
fig

let's compare the knee-estimated SSM with the argmin SSM:

using Kneedle
kr = kneedle(shape_range, criterion_values, "convex_dec", 1, scan_type = :smoothing)

Kneedle.KneedleResult{Int64}([1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0  …  243.0, 244.0, 245.0, 246.0, 247.0, 248.0, 249.0, 250.0, 251.0, 252.0], [1.187336634518873e6, 1.1822161487846365e6, 1.1771899476051119e6, 1.1722835293292657e6, 1.1675223923060652e6, 1.1629320348844775e6, 1.1585379554134696e6, 1.1543656522420086e6, 1.150388012553249e6, 1.1465602696348627e6  …  1.1101587339249968e6, 1.1102076998085128e6, 1.1102567423644713e6, 1.1103058613516411e6, 1.1103550565287906e6, 1.110404327654689e6, 1.110453674488105e6, 1.1105030967878073e6, 1.1105525943125647e6, 1.1106021668213191e6], [33])

Let's extract the optimal number of shapes using the Kneedle algorithm:

k_knee = knees(kr)[1]
println("Optimal number of shapes selected: $k_knee")

Optimal number of shapes selected: 33

let's compare the knee-estimated SSM with the argmin SSM:

ssm_knee = SSM(sbm_big, k_knee)

fig = Figure(size = (1000, 500))

for i = 1:4
    ax = Axis(fig[1:2, i], title = "Category $i", aspect = 1)
    hidedecorations!(ax)
    hm = heatmap!(ax, ssm_estimated, k = i, colormap = :binary, colorrange = (0, 1))
    ax2 = Axis(fig[3:4, i], aspect = 1)
    hidedecorations!(ax2)
    hm = heatmap!(ax2, ssm_knee, k = i, colormap = :binary, colorrange = (0, 1))
end

Colorbar(
    fig[2:3, 5],
    colormap = :binary,
    colorrange = (0, 1),
    label = "Probability",
    height = Relative(0.8),
)

fig

We can also measure the approximation error using mean squared error (MSE):

sum_squared_errors(x, y) = sum((x .- y) .^ 2)
function mse(graphon1, graphon2, xs = 0:0.01:1)
    mean(
        sum_squared_errors(params(graphon1(x, y))[2], params(graphon2(x, y))[2]) for
        x in xs, y in xs
    )
end
mse_sbm = mse(graphon_w3, sbm_big)
mse_ssm_estimated = mse(graphon_w3, ssm_estimated)
mse_ssm_knee = mse(graphon_w3, ssm_knee)
println("MSE of SBM approximation: ", round(mse_sbm, digits = 5))
println("MSE of argmin SSM approximation: ", round(mse_ssm_estimated, digits = 5))
println("MSE of knee SSM approximation: ", round(mse_ssm_knee, digits = 5))

MSE of SBM approximation: 0.00351
MSE of argmin SSM approximation: 0.00463
MSE of knee SSM approximation: 0.00563

loss comparison

mses = zeros(length(shape_range) + 1)
mses[end] = mse(graphon_w3, sbm_big)
for (index, s) in enumerate(shape_range)
    ssm_temp = SSM(sbm_big, s)
    mses[index] = mse(graphon_w3, ssm_temp)
end
s_sbm = k_big * (k_big + 1) ÷ 2

fig = Figure(size = (600, 200))
ax = Axis(
    fig[1, 1],
    xlabel = "Number of Shapes",
    ylabel = "Mean Squared Error",
    yscale = log10,
    xticks = [0, 33, 58, 100, 200, s_sbm],
)
lines!(ax, 1:s_sbm, mses)
scatter!(ax, k_knee, mses[k_knee], color = :red, marker = :rect, label = "Elbow")
scatter!(ax, k_opt, mses[k_opt], color = :black, marker = :rect, label = "Argmin")
scatter!(ax, s_sbm, mses[end], marker = :rect, label = "SBM")

axislegend(ax)
fig

This page was generated using Literate.jl.