Gaussian Mixture Clustering for Static SD-WAN Controller Placement

3.1. Data Preparation and Haversine Distancing

TATA, BICS and BESTEL were taken from the Internet Topology Zoo [40]. Nodes without coordinates and invalid edges were removed. Pairwise geographic separation was computed with the Haversine formula, Eq. (1) and Fig. (1) [41], which works on the sphere and so avoids the error of a planar distance over continental spans.

Dᵢⱼ = 2R arcsin √[sin²(Δφᵢⱼ/2) + cos φᵢ cos φⱼ sin²(Δλᵢⱼ/2)]     (1)

Eq. (1) yields a length in the unit of R, not a delay in μs. Tables 3–5 keep the microsecond values of the original record, but that record does not name the timing instrument, the path aggregation or the distance-to-delay conversion behind these placement-stage values. They are treated as reported latency summaries, not verified WAN propagation measurements. The emulator RTTs in Table 7 are a separate reported outcome.

In Eq. (1), D_ij is the great-circle distance between nodes i and j, R the Earth radius, φ latitude and λ longitude, with Δφ_ij and Δλ_ij the differences between the two nodes. The clustering and assignment stages take these distances as their spatial basis.

3.2. GMM Mechanism

In a GMM, an unsupervised probabilistic model, each Gaussian component is a candidate controller domain. K-Means gives every observation one hard centroid and implicitly spherical clusters; a GMM estimates μₖ, Σₖ and πₖ per component, so domains may be anisotropic and may overlap. This suits SD-WAN sites that are unevenly spread or lack sharp geographic boundaries. The mixture density is Eq. (2).

p(x) = Σₖ₌₁ᴷ πₖ N(x | μₖ, Σₖ)     (2)

In (2), πₖ is the mixing coefficient of component k and μₖ, Σₖ its mean and covariance, with πₖ ≥ 0 and Σₖ πₖ = 1. EM (Section 3.3) estimates the parameters, BIC (Eq. (3)) scores each candidate K, and the lowest-BIC model is kept.

BIC(K) = −2ℓ(θ̂K; X) + pK ln N     (3)

For fixed K, the fitted criterion is ℓ(θ; X) = Σᵢ log[Σₖ πₖ N(xᵢ | μₖ, Σₖ)], subject to πₖ ≥ 0, Σₖ πₖ = 1 and positive-definite covariance matrices. EM seeks a local solution for θ = {πₖ, μₖ, Σₖ}. BIC(K) = −2ℓ(θ̂K; X) + pK log N then selects K̂ = arg minK∈𝒦 BIC(K), where 𝒦 is the candidate set. For unrestricted full covariance in d dimensions, pK = (K − 1) + Kd + Kd(d + 1)/2. With d = 2, pK = 6K − 1. This parameter count is conditional on full covariance; the actual covariance setting and candidate range must be recovered from the implementation. BIC balances density fit and statistical model complexity. Its penalty is not an operating-cost model and does not jointly minimise latency, imbalance or deployment cost.

3.3. Expectation-Maximization (EM)

Component membership is latent, so the parameters are estimated with EM. From initial μₖ, Σₖ and πₖ, EM alternates an E-step and an M-step until the log-likelihood gain drops below a tolerance or an iteration cap is reached.

Expectation: each node xn receives a responsibility for every component, computed from the current values of:

1. Mean (μk): mean location of the component’s nodes.
2. Covariance (Σk): spread and orientation of the nodes about μk.
3. Mixing coefficient (πk): fraction of the network assigned to the component.

The responsibility γ(znk) of component k for node xn is given by Eq. (4):

γᵢₖ = πₖ N(xᵢ | μₖ, Σₖ) / [Σⱼ₌₁ᴷ πⱼ N(xᵢ | μⱼ, Σⱼ)]     (4)

N(·) is the Gaussian density, and γ is the posterior probability that x_n belongs to component k under the current estimates.

Maximization: After assigning probabilities in the E-step, we update the clusters’ parameters to better fit the data:

1. New mean (μ_k^new): Calculated as the weighted average of all nodes, where the weights are the probabilities from the E-step.
2. New covariance (Σ_k^new): Updated to reflect how spread out the nodes are around the new mean.
3. New mixing coefficient (π_k^new): adjusted to reflect the proportion of all data points that are likely to be in cluster k.

The formulas (Eq. 5-7) for these updates are:

   (5)

   (6)

   (7)

Nₖ = Σᵢ γᵢₖ is the effective size of component k. EM does not decrease the likelihood between iterations but can stop at a local maximum or a degenerate covariance, which is why reproducible implementations add covariance regularisation and several initialisations. Restart count, tolerance and regularisation value are absent from the retained record and are not asserted here. Fig. (2) and Algorithm 1 state the intended procedure, which the manuscript alone does not allow to be fully reconstructed.

A fitted mean is generally not a node. Each mean is moved to the geographically nearest topology node, and each node is then assigned to the nearest selected controller, so the final domains need not match the posterior-membership domains. Two means can also map to one node.

3.4. Clustering Metrics

ACL, WCL, ICL and NDR describe different aspects of a configuration and are not combined into an optimisation objective. Eq. (8-11) below retain the pairwise intra-domain interpretation of the original equations. In particular, intra-domain node-pair latency must not be relabelled as node-to-controller latency.

Average-case latency within domain k is the mean of the ordered node-pair delays in that domain. Let Vₖ contain nₖ nodes and let δᵢⱼ denote latency, rather than Haversine distance. For nₖ ≥ 2:

ACLₖ = [1 / (nₖ(nₖ − 1))] Σᵢ∈Vₖ Σⱼ∈Vₖ,ⱼ≠ᵢ δᵢⱼ.     (8)

Worst-case latency within domain k is the largest node-pair latency for nₖ ≥ 2:

WCLₖ = maxᵢ,ⱼ∈Vₖ,ᵢ≠ⱼ δᵢⱼ.     (9)

For distinct selected controller nodes cₖ and cₗ, the pairwise inter-controller latency is:

ICLₖₗ = δ(cₖ, cₗ), k ≠ l.     (10)

A network-level ACL requires an explicit weighting across domains, and a network-level ICL requires a specified mean, sum or maximum across controller pairs. Singleton-domain handling must also be specified. These conventions cannot be inferred from the existing table entries. For hard assignments, NDR is defined by the maximum domain size relative to the mean:

NDR = maxₖ nₖ / (N/K) = K maxₖ nₖ / N.     (11)

3.5. Computational Complexity and Scalability

Let N be the number of nodes, K the number of components, d the feature dimension, I the EM iterations and S the independent initialisations. A full-covariance EM fit costs O(I(NKd² + Kd³)), including covariance factorisation, with working memory O(NK + Nd + Kd²). Selecting among 𝒦 therefore costs O(S ΣK∈𝒦 IK(NKd² + Kd³)), rather than the cost of a single fit. With d = 2 and a bounded K, each iteration is approximately linear in N. However, materialising all pairwise geographic distances costs O(N²) time and memory. Mapping representatives and assigning nodes each require O(NK) direct distance evaluations. K-Means has an approximate per-fit cost O(INKd); covariance modelling adds work, and no measured runtime advantage is claimed for GMM.

With 33–145 nodes, the tested networks do not form a scaling benchmark. Runtime, peak memory, iteration counts, model-selection cost and coordination overhead were not recorded, so the orders above are analytical bounds only. A scaling study would need a controlled sequence of N on fixed hardware, a fixed K range and restart budget, and separate timings for preprocessing, fitting, selection and assignment. On-demand distance evaluation and a narrower 𝒦 could cut the O(N²) and per-K fitting costs; neither was benchmarked.

4. COMPARATIVE ANALYSIS

GMM is compared with K-Means, AHC, AP and DP on the same three topologies and metrics, to locate where GMM gives better placement figures and where a baseline does; no single clustering model is assumed best on every objective.

4.1. Proposed Method Implementation

All five methods were implemented in Python 3, take the same geographic topology input and are scored with the same placement metrics (Section 3.4). Fig. (3) summarises the pipeline.

Inputs are shared, but recorded controller counts are not (Table 2). Tables 3–9 therefore compare complete configurations in which algorithm, K and assignment act together, and shared inputs alone do not make this a controlled algorithm comparison. Fig. (3) shows common evaluation stages, not equal controller resources.

The benchmark has three Internet Topology Zoo networks [40]: TATA with 145 nodes and 194 edges, BICS with 84 nodes and 101 edges, and BESTEL with 33 nodes and 48 edges (Table 1, Fig. 4). Three instances make this a case study. The networks are not claimed to be the largest available, a representative sample or evidence of Internet-scale behaviour, and their layouts carry no production traffic or measured control-channel delay.