arXiv CSopen access

Rethinking Traffic Matrix Completion: Estimate the Process, Not the Entries

Rethinking Traffic Matrix Completion: Estimate the Process, Not the Entries Xiyuan Liu∗

Xidian University Xi’an, China

Zihao Wang∗

[email protected] Xidian University Xi’an, China

Wenting Wei†

Xiucheng Tian

arXiv:2605.02225v1 [cs.NI] 4 May 2026

Guanzuo Liu

[email protected] Xidian University Xi’an, China

[email protected] Xidian University Xi’an, China

[email protected] Xidian University Xi’an, China

Abstract

1

Traffic matrix measurement is fundamental for datacenter operations, but obtaining complete traffic matrices at scale remains challenging due to the prohibitive cost of global fine-grained measurement and partial observations resulting from network faults. Although existing matrix completion methods (reduce cost) achieve satisfactory performance in specific scenarios, their reliance on restrictive assumptions or black-box mappings results in a lack of interpretability and an inability to characterize uncertainty. In this paper, we propose Utimac, an uncertainty-aware traffic matrix completion for data center networks. Our analysis shows that, within a locally stationary window, log-domain traffic can be decomposed into a principal statistical component and a sparse deviation component. Based on this insight, we formulate traffic matrix completion as a parameter inference problem: multiple partially observed frames within a window are used to infer shared parameters and recover missing entries. To avoid the intractability and boundary degeneracy of the original integral-form marginal likelihood, we construct a regularized surrogate objective and solve the resulting joint optimization problem with block coordinate descent. Utimac consistently outperforms all baselines on data center networks datasets in both overall and burst scenarios, with its advantage becoming more pronounced as observations grow sparser. All code is publicly available in an anonymous repository: https://anonymous.4open.science/r/Utimac-0551/

In recent years, datacenter networks have been supporting increasingly large-scale workloads [25], and applications such as distributed training, LLM online inference, and large-scale data processing have further intensified network communication [8, 15]. As these workloads continue to expand, accurate measurement of network traffic becomes important [4], because it provides critical data support for downstream tasks such as traffic engineering, fault localization, and anomaly detection [2]. Obtaining a complete traffic matrix in a large-scale datacenter is challenging. As the network grows, exhaustively measuring all source-destination pairs requires devices to maintain large flow tables, which continuously consume the limited memory and processing resources of switches [11]. Meanwhile, the available observations are often incomplete [10]. Fine-grained monitoring is typically enabled only at selected locations and during selected periods, and the collected statistics can further become incomplete due to reporting delays, record loss, or device failures [32]. Existing datacenter traffic measurement methods can be divided into direct measurement and indirect measurement. Direct measurement obtains traffic statistics by collecting packets or flows from network devices or monitoring points; representative examples include flow-recording approaches such as NetFlow/IPFIX and systems like OpenTM [5, 6, 22], which directly measure traffic matrices using counters in switch flow tables. In addition, sketch-based designs such as UnivMon and ElasticSketch maintain compact summaries to approximately record traffic features with lower memory overhead [12, 29]. Indirect measurement infers the global traffic matrix from partial observations. Some of the works rely on explicit priors such as low-rankness or sparsity, as illustrated by Xie et al.’s studies on variable-rate measurements [27], block matrix completion [28], and deep adversarial tensor completion [26]. While some adopt data-driven

Keywords Traffic Matrix Completion, Uncertainty-aware

∗

Xiyuan Liu and Zihao Wang contributed equally to this work. Corresponding author.

†

Introduction

Liu, Wang, et al.

modeling and learn the mapping from partial observations to the full traffic matrix, with representative examples including AutoTomo and Satformer [16, 17]. Unfortunately, existing methods still have several gaps. For direct measurement, the main limitation is that the measurement overhead grows with network scale and measurement granularity[9]. For example, flow-record and per-flow counting methods require devices to maintain a large amount of flow-level tables [30], while sampling and path-counting methods introduce extra overhead for path maintenance or statistics aggregation [20]. As a result, direct measurement often leads to high resource and management costs in large-scale data center networks. In contrast, indirect measurement reduces the cost of explicit monitoring, but its performance often depends on strong assumptions or the training data. Methods based on priors such as low rankness rely on strong assumptions [19], and their performance degrades when real traffic deviates from these assumptions. Deep learning based methods can learn complex mappings, but they usually lack interpretability and are sensitive to distribution shifts, which limits their generalization. Most existing indirect measurement methods focus on point estimation and cannot provide confidence ranges for the traffic matrix. However, traffic matrix imputation is inherently an inference problem under incomplete observations, and its outputs therefore contain uncertainty. Without explicitly characterizing this uncertainty, downstream systems cannot identify high-risk predictions and may be misled in subsequent decisionmaking [13]. We propose Utimac, an traffic matrix completion method that mind the gaps of existing methods: the high cost of direct measurement, the limited interpretability of black-box inference, and the lack of uncertainty characterization in point estimation. Our data analysis shows that, within a locally stationary time window, log-domain traffic can be described by the superposition of a joint Gaussian principal component and a Laplacian deviation component. Based on this statistical structure, we formulate traffic completion as a parameter inference problem driven by partial observations: we infer the shared mean, shared covariance, and sparsity parameter from multiple partially observed frames within a window and then recover the unobserved traffic entries accordingly. Our contributions are summarized as follows:

• We reformulate traffic matrix completion as a statistical inference problem. Within each locally stationary window, we represent log-domain traffic by a joint Gaussian principal component and a Laplacian

deviation component, and tie multiple partially observed frames to a shared set of parameters. Completion is then carried out by inferring the shared mean, covariance, and sparsity parameter from partial observations and recovering the missing entries. • We derive a computable and well-posed inference objective from the integral-form marginal likelihood. To avoid high-dimensional integration and eliminate the degeneracy caused by sparsity growth and covariance collapse, we introduce a surrogate objective together with stabilizing regularizers, which yields a regularized nested optimization model. • We develop a structured solver for the regularized problem. We rewrite it as a joint optimization over deviation variables, mean, sparsity, and covariance, and solve it by block coordinate descent. This procedure alternates between frame-wise deviation updates, mean estimation, sparsity update, and covariance optimization, and finally recovers the missing traffic entries together with their uncertainty. Utimac consistently outperforms all baselines on data center networks datasets in both overall and burst scenarios, with its advantage becoming more pronounced as observations grow sparser.

2

Problem Statement

This section establishes the foundation of our formulation. We first introduce the system model and assumptions, and then formalize the traffic matrix completion problem under partial observations.

2.1

System Model and Assumptions

We model data center traffic in discrete time. At the selected network layer, let {𝑢1 , 𝑢2 , … , 𝑢𝑀 } denote the source-side network units and {𝑣1 , 𝑣2 , … , 𝑣𝑁 } denote the destination-side network units. For each source-destination pair (𝑢𝑖 , 𝑣𝑗 ), let 𝑋𝑡 (𝑖, 𝑗) ∈ ℝ+ be the traffic random variable at time 𝑡 , and let 𝑥𝑡 (𝑖, 𝑗) be one realization of 𝑋𝑡 (𝑖, 𝑗). By vectorizing all source-destination flows at time 𝑡 in a fixed order, we obtain the traffic random vector: ⊤

𝑋𝑡 = [𝑋𝑡 (𝑠1 ) 𝑋𝑡 (𝑠2 ) ⋯ 𝑋𝑡 (𝑠𝑑 )] ∈ ℝ𝑑+ , 𝑑 = 𝑀𝑁 ,

(1) where 𝑆 = {𝑠1 , 𝑠2 , … , 𝑠𝑑 } is the spatial index set after vectorization, and each index 𝑠𝑘 corresponds to one unique sourcedestination pair. The corresponding realization is denoted by 𝑥𝑡 ∈ ℝ𝑑+ . We adopt a locally stationary model. Let the time axis be partitioned into local windows. 𝑛-th window is denoted by:

𝑊𝑛 = {𝑡𝑛,1 , 𝑡𝑛,2 , … , 𝑡𝑛,𝐿𝑛 },

(2)

Rethinking Traffic Matrix Completion: Estimate the Process, Not the Entries

where 𝐿𝑛 = |𝑊𝑛 |. We assume that the distribution parameters of the traffic process remain approximately constant within the same window 𝑊𝑛 , and may change across different windows. To account for non-negativity and heavy-tailed behavior, we work in the log domain. Let 𝜀 ∈ ℝ𝑑++ be a strictly positive offset vector. We define

𝑍𝑡 = log(𝑋𝑡 + 𝜀),

(3)

where the logarithm is applied elementwise. We decompose the log-domain traffic into a structured component and a deviation component:

𝑍𝑡 = 𝑈𝑡 + 𝑂𝑡 ,

𝑡 ∈ 𝑊𝑛 .

(4)

Here, 𝑈𝑡 ∈ ℝ𝑑 represents the dominant statistical structure, and 𝑂𝑡 ∈ ℝ𝑑 represents local deviations. Later empirical analysis shows that, within a local stationary window, the dominant part of the log-domain traffic exhibits an approximately joint Gaussian pattern, while the deviation part is sparse. Based on that empirical finding, we use the following parametric model in window 𝑊𝑛 :

𝑈𝑡 ∼ 𝒩(𝜇𝑛 , Σ𝑛 ), and

𝑡 ∈ 𝑊𝑛 ,

(5)

𝑑

𝜆 𝑝𝑂 (𝑜; 𝜆𝑛 ) = ( 𝑛 ) exp(−𝜆𝑛 ‖𝑜‖1 ) , (6) 2 where 𝜇𝑛 ∈ ℝ𝑑 is the shared mean vector, Σ𝑛 ∈ ℝ𝑑×𝑑 is the shared covariance matrix with Σ𝑛 ≻ 0, and 𝜆𝑛 > 0 is the shared sparsity parameter. We use 𝜃𝑛 = (𝜇𝑛 , Σ𝑛 , 𝜆𝑛 )

(7)

to denote the shared statistical parameters in window 𝑊𝑛 . We further assume conditional independence within each window. Given 𝜃𝑛 , the pairs {(𝑈𝑡 , 𝑂𝑡 ) ∶ 𝑡 ∈ 𝑊𝑛 } are independent across time. For each 𝑡 ∈ 𝑊𝑛 , 𝑈𝑡 and 𝑂𝑡 are independent. The joint density of the latent variables at time 𝑡 is therefore

𝑝𝑈 ,𝑂 (𝑢𝑡 , 𝑜𝑡 ; 𝜃𝑛 ) = 𝑝𝑈 (𝑢𝑡 ; 𝜇𝑛 , Σ𝑛 ) 𝑝𝑂 (𝑜𝑡 ; 𝜆𝑛 ) 1 ∝ exp(− (𝑢𝑡 − 𝜇𝑛 )⊤ Σ−1 𝑛 (𝑢𝑡 − 𝜇𝑛 ) − 𝜆𝑛 ‖𝑜𝑡 ‖1 ) . 2 (8) For each time 𝑡 , let 𝑏𝑡 ∈ {0, 1}𝑑 be the binary observation mask. The observed index set is Ω𝑡 = {𝑠𝑘 ∈ 𝑆 ∶ 𝑏𝑡 (𝑠𝑘 ) = 1}.

(9)

The actual observation at time 𝑡 is the subvector of 𝑥𝑡 on |Ω | (Ω ) Ω𝑡 , denoted by 𝑥𝑡 𝑡 ∈ ℝ+ 𝑡 . We treat Ω𝑡 as exogenous and known during parameter estimation. All subsequent inference is conditioned on Ω𝑡 . Under this setting, the main task in each local window is to estimate the shared parameter 𝜃𝑛 from partial observations and then characterize the conditional distribution of the complete traffic vector.

2.2

Problem Formulation

The main goal is to estimate the shared parameter 𝜃𝑛 in each local stationary window. Once 𝜃𝑛 is estimated, the distribution of the complete traffic vector in that window is specified. Given current observations, the unobserved entries are characterized by their conditional posterior distribution. For any time 𝑡 ∈ 𝑊𝑛 , define the observed log-domain subvector (Ω𝑡 )

𝑧𝑡

(Ω𝑡 )

= log(𝑥𝑡

+ 𝜀Ω𝑡 ) ∈ ℝ|Ω𝑡 | ,

(10)

where 𝜀Ω𝑡 is the subvector of 𝜀 indexed by Ω𝑡 . By (4), we have (Ω𝑡 )

𝑧𝑡

(Ω )

(Ω )

= 𝑢 𝑡 𝑡 + 𝑜𝑡 𝑡 .

(11)

Since 𝑈𝑡 ∼ 𝒩(𝜇𝑛 , Σ𝑛 ), its marginal distribution on any observed index set remains Gaussian: (Ω𝑡 )

𝑈𝑡

(Ω )

(Ω ,Ω𝑡 )

∼ 𝒩(𝜇𝑛 𝑡 , Σ𝑛 𝑡

),

(12)

(Ω )

(Ω ,Ω )

where 𝜇𝑛 𝑡 is the subvector of 𝜇𝑛 on Ω𝑡 , and Σ𝑛 𝑡 𝑡 is the principal submatrix of Σ𝑛 indexed by Ω𝑡 . The deviation term on the same subspace has density |Ω𝑡 |

𝜆 (Ω ) 𝑝𝑂 (𝑜𝑡 𝑡 ; 𝜆𝑛 ) = ( 𝑛 ) 2

(Ω𝑡 )

exp(−𝜆𝑛 ‖𝑜𝑡

‖1 ) .

(13)

Given 𝜃𝑛 , the exact conditional likelihood of the observed log-domain vector is obtained by marginalizing out the latent deviation variable: (Ω𝑡 )

𝑝(𝑧𝑡

∣ Ω𝑡 ; 𝜃𝑛 ) = ∫

ℝ|Ω𝑡 |

(Ω𝑡 )

𝑝𝑈 (𝑧𝑡

(Ω )

(Ω ,Ω𝑡 )

− 𝑜; 𝜇𝑛 𝑡 , Σ𝑛 𝑡

)

(14)

⋅ 𝑝𝑂 (𝑜; 𝜆𝑛 ) 𝑑𝑜. By conditional independence within 𝑊𝑛 , the principled observed-data maximum likelihood estimator is

𝜃𝑛̂ML = arg

max

(Ω𝑡 )

∑ log 𝑝(𝑧𝑡

𝜇𝑛 , Σ𝑛 ≻0, 𝜆𝑛 >0 𝑡∈𝑊

∣ Ω𝑡 ; 𝜇𝑛 , Σ𝑛 , 𝜆𝑛 ) .

𝑛

(15) The objective in (15) is statistically well defined. Its direct evaluation is expensive because each frame requires the computation of a high-dimensional convolution integral. (Ω ,Ω )

For a general non-diagonal covariance matrix Σ𝑛 𝑡 𝑡 , the integral in (14) does not admit a simple closed form. To obtain a computable objective, we introduce a profiled approximation. Using the identity (Ω )

(Ω𝑡 )

𝑢𝑡 𝑡 = 𝑧 𝑡

(Ω )

− 𝑜𝑡 𝑡 ,

(16)

Liu, Wang, et al.

the reduced joint log-density for a feasible decomposition can be written as

1 (Ω ) (Ω ) (Ω ) (Ω ) ⊤ (Ω ,Ω ) −1 ℒ𝑡̃ (𝜃𝑛 , 𝑜𝑡 𝑡 ) = − (𝑧𝑡 𝑡 − 𝑜𝑡 𝑡 − 𝜇𝑛 𝑡 ) (Σ𝑛 𝑡 𝑡 ) 2

1 (Ω ) (Ω ) (Ω ) (Ω ,Ω ) ⋅ (𝑧𝑡 𝑡 − 𝑜𝑡 𝑡 − 𝜇𝑛 𝑡 ) − log det Σ𝑛 𝑡 𝑡 2

(Ω ) + |Ω𝑡 | log 𝜆𝑛 − 𝜆𝑛 ‖𝑜𝑡 𝑡 ‖1 + 𝐶𝑡 ,

(17) where 𝐶𝑡 is constant optimization variables. We approximate the frame-level marginal log-likelihood by the maximum value of the reduced joint log-density. The resulting estimator is

𝜃𝑛̂ = arg

∑

max

(Ω𝑡 )

max

𝜇𝑛 , Σ𝑛 ≻0, 𝜆𝑛 >0 𝑡∈𝑊 𝑜 (Ω𝑡 ) ∈ℝ|Ω𝑡 | 𝑛 𝑡

ℒ𝑡̃ (𝜃𝑛 , 𝑜𝑡

).

(18)

Dropping constants and changing the sign gives the equivalent minimization form:

𝜃𝑛̂ = arg

∑

min

min

𝜇𝑛 , Σ𝑛 ≻0, 𝜆𝑛 >0 𝑡∈𝑊 𝑜 (Ω𝑡 ) ∈ℝ|Ω𝑡 | 𝑛 𝑡

(Ω ) { − ℒ𝑡̃ (𝜃𝑛 , 𝑜𝑡 𝑡 ) }.

(19) The objective in (19) is a bilevel optimization problem. The outer variables are the shared parameters (𝜇𝑛 , Σ𝑛 , 𝜆𝑛 ). The inner variables are the frame-specific deviation vectors (Ω )

{𝑜𝑡 𝑡 }𝑡∈𝑊𝑛 . For fixed outer parameters, each inner problem is convex because it contains a positive-definite quadratic term and an ℓ1 term. The full problem is generally nonconvex in (𝜇𝑛 , Σ𝑛 , 𝜆𝑛 ) due to the inverse covariance term and the log det term. After obtaining 𝜃𝑛̂ = (𝜇̂𝑛 , Σ̂ 𝑛 , 𝜆̂ 𝑛 ), the conditional distribution of the unobserved entries is determined. Let Ωmis = 𝑆 ∖ Ω𝑡 𝑡

(20)

be the unobserved index set at time 𝑡 . The posterior distribution in the log domain is (Ωmis 𝑡 )

𝑝(𝑧𝑡

(Ω𝑡 )

∣ 𝑧𝑡

, Ω𝑡 ; 𝜃𝑛̂ ) .

The first observation is local stationarity over short intervals. Over such intervals, active jobs or resource contention change little, so consecutive frames are driven by similar communication patterns. As a result, the empirical mean, covariance, and overall distribution shape remain stable within a proper window, which supports using a local window for parameter sharing. The second observation concerns the dominant structure in the log domain. Raw traffic spans a wide dynamic range because flows differ substantially in message size. The log transform compresses this variation and regularizes the dominant structure. After transformation, most samples cluster around a stable center within each window, source-destination pairs remain statistically dependent due to common jobs or similar communication patterns. A detailed empirical analysis demonstrating that the dominant component of log-domain traffic approximately follows a joint Gaussian distribution is detailed in Appendix A. The third observation is the presence of sparse deviations around the dominant structure. Short-lived events such as incast bursts or transient congestion can create large deviations on a small number of entries. These events do not alter the dominant structure, but they produce localized perturbations that cannot be absorbed by the principal component. This motivates a sparse additive deviation term. Empirical evidence supporting the presence of these sparse deviations is presented in Appendix B. Based on these observations, we model log-domain traffic in each local window as the sum of a dominant structured component and a sparse deviation component. We use a joint Gaussian model for the dominant component to capture the shared center and spatial dependence of regular traffic, and a Laplacian prior for the deviation component to capture sparse local perturbations. Accordingly, for each time 𝑡 ∈ 𝑊𝑛 , we write

𝑍𝑡 = log(𝑋𝑡 + 𝜀) = 𝑈𝑡 + 𝑂𝑡 ,

(22)

where 𝑈𝑡 denotes the dominant structured component and 𝑂𝑡 denotes the sparse deviation component. We then parameterize them as

we obtain the corresponding posterior distribution in the original traffic domain: (Ωmis 𝑡 )

𝑝(𝑥𝑡

3

(Ω𝑡 )

∣ 𝑥𝑡

, Ω𝑡 ; 𝜃𝑛̂ ) .

Statistical Structure Discovery

(21)

By the inverse transform

𝑋𝑡 = exp(𝑍𝑡 ) − 𝜀,

3.1

(23)

Method

This section presents the log-domain Gaussian-Laplacian decomposition model, its regularized reformulation, a block coordinate descent solver, and uncertainty quantification.

𝑈𝑡 ∼ 𝒩(𝜇𝑛 , Σ𝑛 ),

𝑡 ∈ 𝑊𝑛 ,

(24)

(25)

and 𝑑

𝜆 𝑝𝑂 (𝑜; 𝜆𝑛 ) = ( 𝑛 ) exp(−𝜆𝑛 ‖𝑜‖1 ) , 2

(26)

where 𝜇𝑛 and Σ𝑛 are the shared mean and covariance in window 𝑊𝑛 , and 𝜆𝑛 controls the sparsity of local deviations.

Rethinking Traffic Matrix Completion: Estimate the Process, Not the Entries

3.2

Analysis and Reformulation

where

The profiled objective in (19) is computable, but it is not directly suitable for optimization. To analyze its behavior, define (Ω )

Ψ(𝜇𝑛 , Σ𝑛 , 𝜆𝑛 ) ∶= ∑ min {−ℒ𝜏̃ (𝜃𝑛 , 𝑜𝜏 𝜏 )} . (Ω ) 𝜏 ∈𝑊𝑛 𝑜𝜏 𝜏

(27)

For fixed (𝜇𝑛 , Σ𝑛 , 𝜆𝑛 ), each inner problem is a convex optimization over the deviation variable. The main difficulty lies in the outer objective, which is generally nonconvex in (𝜇𝑛 , Σ𝑛 , 𝜆𝑛 ). More importantly, (27) is not lower bounded. First, since (Ω )

𝑜𝜏 𝜏 = 0 is always feasible, we have Ψ(𝜇𝑛 , Σ𝑛 , 𝜆𝑛 ) ≤ 𝐶(𝜇𝑛 , Σ𝑛 ) − ( ∑ |Ω𝜏 |) log 𝜆𝑛 ,

(28)

(Ω )

𝐹𝑛 = ∑ 𝜙𝜏 (𝜇𝑛 , Σ𝑛 , 𝜆𝑛 ; 𝑜𝜏 𝜏 ) + 𝜌 tr(Σ−1 𝑛 ) + 𝜂𝜆𝑛 .

This form exposes four variable blocks: deviation variables, mean, sparsity parameter, and covariance. We first update the deviation variables. For fixed 𝜇𝑛 , Σ𝑛 , and 𝜆𝑛 , define (Ω )

(Ω ,Ω𝜏 ) −1

(Ω )

𝑟𝜏 = 𝑧𝜏 𝜏 − 𝜇𝑛 𝜏 ,

𝑄𝜏 = (Σ𝑛 𝜏

(Ω )

(Ω )

(Ω )

𝑜𝜏 𝜏 = 𝑧 𝜏 𝜏 − 𝜇 𝑛 𝜏 ,

(Ω ,Ω𝜏 )

log det Σ𝑛 𝜏

= |Ω𝜏 | log 𝜖 → −∞.

(30)

This shows that the covariance matrix can continue to shrink and further reduce the objective through the log-determinant term. To remove these two escape directions, we introduce the regularizers 𝜌 tr(Σ−1 𝜌 > 0, (31) 𝑛 ),

𝜂 > 0.

(32)

The resulting regularized problem is

(𝜇̂𝑛 , Σ̂ 𝑛 , 𝜆̂ 𝑛 ) = arg

min

𝜇𝑛 ,Σ𝑛 ≻0,𝜆𝑛 >0

(33) Equation (33) preserves the original nested inference structure while eliminating boundary degeneracy. It remains nonconvex and nonsmooth, but it is now better posed and numerically stable.

3.3

Block Coordinate Solution

We solve (33) using block coordinate descent. By explicitly reintroducing the deviation variables, we rewrite the problem as min

(Ω ) 𝜇𝑛 ,Σ𝑛 ≻0,𝜆𝑛 >0,{𝑜𝜏 𝜏 }𝜏 ∈𝑊𝑛

(Ω )

𝐹𝑛 (𝜇𝑛 , Σ𝑛 , 𝜆𝑛 , {𝑜𝜏 𝜏 }𝜏 ∈𝑊𝑛 ) ,

(34)

(37)

After that, we update the shared mean. Let 𝑃𝜏 be the selection matrix for Ω𝜏 , and define (Ω )

(Ω )

𝑎𝜏 = 𝑧𝜏 𝜏 − 𝑜𝜏̂ 𝜏 .

(38)

Then 𝜇𝑛 is obtained from (Ω ,Ω𝜏 )

( ∑ 𝑃𝜏⊤ (Σ𝑛 𝜏

)

−1

𝜏 ∈𝑊𝑛

(Ω ,Ω𝜏 )

𝑃𝜏 ) 𝜇𝑛 = ∑ 𝑃𝜏⊤ (Σ𝑛 𝜏

)

𝜏 ∈𝑊𝑛

−1

𝑎𝜏 . (39)

The sparsity parameter is then updated by

𝜆̂ 𝑛 =

∑𝜏 ∈𝑊𝑛 |Ω𝜏 | (Ω )

∑𝜏 ∈𝑊𝑛 ‖𝑜𝜏̂ 𝜏 ‖1 + 𝜂

.

(40)

Finally, define the residual (Ω )

Ψ(𝜇𝑛 , Σ𝑛 , 𝜆𝑛 )+𝜌 tr(Σ−1 𝑛 )+𝜂𝜆𝑛 .

(36)

(Ω ) +𝜆𝑛 ‖𝑜𝜏 𝜏 ‖1 }.

and

𝜂𝜆𝑛 ,

,

⊤ 1 (Ω ) (Ω ) (Ω ) 𝑜𝜏̂ 𝜏 = arg min { (𝑜𝜏 𝜏 − 𝑟𝜏 ) 𝑄𝜏 (𝑜𝜏 𝜏 − 𝑟𝜏 ) (Ω𝜏 ) 2 𝑜𝜏

(29)

which makes the quadratic residual term vanish. If we further set Σ𝑛 = 𝜖𝐼 with 𝜖 ↓ 0, then

)

and solve, for each 𝜏 ∈ 𝑊𝑛 ,

𝜏 ∈𝑊𝑛

where 𝐶(𝜇𝑛 , Σ𝑛 ) is independent of 𝜆𝑛 . Hence the objective decreases without bound as 𝜆𝑛 → ∞. Second, for any fixed 𝜇𝑛 and 𝜆𝑛 > 0, one can choose

(35)

𝜏 ∈𝑊𝑛

(Ω )

(Ω )

(Ω )

𝑒𝜏 𝜏 = 𝑧𝜏 𝜏 − 𝑜𝜏̂ 𝜏 − 𝜇𝑛 𝜏 ,

(41)

and update the covariance by

1 (Ω ) (Ω ,Ω ) (Ω ) Σ̂ 𝑛 = arg min { ∑ [(𝑒𝜏 𝜏 )⊤ (Σ𝑛 𝜏 𝜏 )−1 𝑒𝜏 𝜏 Σ𝑛 ≻0 2 𝜏 ∈𝑊 𝑛

1 (Ω ,Ω ) + log det Σ𝑛 𝜏 𝜏 ] + 𝜌 tr(Σ−1 𝑛 )}. 2

(42)

The complete iteration is therefore (Ω )

{𝑜𝜏 𝜏 }𝜏 ∈𝑊𝑛 ⟶ 𝜇𝑛 ⟶ 𝜆𝑛 ⟶ Σ𝑛 .

(43)

This block structure makes the regularized nested inference problem tractable in practice.

Liu, Wang, et al.

Once the shared parameter estimate 𝜃𝑛̂ = (𝜇̂𝑛 , Σ̂ 𝑛 , 𝜆̂ 𝑛 ) and

(Ω ) the frame-level optimal deviations {𝑜𝑡̂ 𝑡 }𝑡∈𝑊𝑛 have been ob-

tained from the block coordinate descent procedure in (43), the conditional distribution of each unobserved log-domain entry 𝑍𝑡 (𝑗), 𝑗 ∈ Ωmis 𝑡 , follows a one-dimensional Normal– Laplace distribution arising from the Gaussian–Laplace convolution posterior, from which a closed-form marginal CDF enables efficient construction of plug-in 95% credible intervals in the original traffic domain; the complete derivation is provided in Appendix C.

Facebook-Pod-B

Uncertainty Quantification

MAE

0.120

RMSE

0.115 0.105

4.2

Evaluation Metrics

All metrics are computed exclusively on missing positions Ωmiss (see Appendix E for formal definitions). Overall impu𝑡 tation accuracy is measured by MAE, RMSE, and wMAPE across all missing entries. Burst flow detection uses Precision, Recall, and F1 to assess identification of anomalously large flows within each time slot; full per-dataset detection tables are provided in Appendix F. Burst flow imputation accuracy uses Burst-MAE, Burst-RMSE, and Burst-wMAPE to evaluate numerical recovery on true missing burst flows, with Burst-wMAPE as the headline indicator. Metric for accuracy of the estimation intervals is PICP97.5 , the fraction of true missing values falling within the closed-form 97.5% estimation interval derived.

0.26

0.15

0.100

0.24 0.14

0.095

0.3 0.4 0.5 0.6 0.7 0.8 0.9

0.10

0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.4

0.8

0.08

1.3

0.6

0.06

1.2

0.4

0.04

1.1 1.0

0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

0.3 0.4 0.5 0.6 0.7 0.8 0.9

Observation Rate

0.3 0.4 0.5 0.6 0.7 0.8 0.9

Observation Rate

ImputeFormer

Observation Rate

Diffusion-TM

Utimac

Figure 1: Overall MAE, RMSE, and wMAPE vs. 𝑝obs on Facebook-Pod-B (top) and Facebook-ToR-A (bottom). Facebook-ToR-A

Normalized Absolute Residual

Facebook-Pod-B 0.35

0.25

0.30

0.20

0.25

0.15

0.20 0.15

0.10

0.10 0.05

0.05 0.00

0.00 0.3

0.4

0.5

0.6

0.7

Observation Rate PSW-I

0.8

0.9

ImputeFormer

0.3

0.4

0.5

0.6

0.7

Observation Rate

Diffusion-TM

0.8

0.9

Utimac

Figure 2: Normalized absolute residual distributions (5th– 95th percentile) vs. 𝑝obs on Facebook-Pod-B (left) and Facebook-ToR-A (right).

Facebook-Pod-B

Burst-MAE

Burst-RMSE

Burst-wMAPE

0.26

0.30

0.32

0.24

0.28

0.30

0.22

0.28

0.26

0.20

0.26

0.24

0.18 0.3 0.4 0.5 0.6 0.7 0.8 0.9

Facebook-ToR-A

We evaluate on real-world traffic datasets. Facebook-Pod-B and Facebook-ToR-A are DCN datasets from Facebook’s production infrastructure: Pod-B captures pod-level aggregated traffic (𝑁 =8, 𝑑=64), while ToR-A captures rack-level traffic (𝑁 =155, 𝑑=24,025). DCN traffic exhibits strong burstiness and heavy-tailed distributions [2], making these the primary evaluation datasets; additional WAN results on GÉANT [23] (𝑁 =22, 𝑑=484) are reported in Appendix F. We compare against three baselines covering complementary paradigms. PSW-I [24] minimizes an optimal-transport discrepancy between time-series patches without parametric training. ImputeFormer [14] is a low-rank-induced Transformer combining matrix-completion priors with projected attention and Fourier sparsity regularization. Diffusion-TM [31] is a DDPM-based generative model using a routing-free TMC inference branch. Utimac runs on an Intel Core i9-13980HX CPU; all neural baselines run in PyTorch on an NVIDIA GeForce RTX 5090 GPU. Entries are masked under 𝑝obs ∈ {0.3, 0.4, … , 0.9}; the same fixed-seed mask is shared by all methods. Burst flows are identified with dominance multiplier 𝛼=2.0 and threshold 𝛽=0.8.

0.28

0.16

0.110

PSW-I

4 Experiments 4.1 Experimental Setup

wMAPE

0.17

0.3 0.4 0.5 0.6 0.7 0.8 0.9

Facebook-ToR-A

3.4

0.10

0.22

0.24 0.3 0.4 0.5 0.6 0.7 0.8 0.9

1.3

0.6

0.06

1.2

0.4

0.04

1.1 1.0

0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

Observation Rate PSW-I

0.3 0.4 0.5 0.6 0.7 0.8 0.9

1.4

0.8

0.08

0.22

0.3 0.4 0.5 0.6 0.7 0.8 0.9

Observation Rate

ImputeFormer

Diffusion-TM

0.3 0.4 0.5 0.6 0.7 0.8 0.9

Observation Rate

Utimac

Figure 3: Burst-MAE, Burst-RMSE, and Burst-wMAPE vs. 𝑝obs on Facebook-Pod-B and Facebook-ToR-A.

4.3

Results

Figure 1 reports overall MAE, RMSE, and wMAPE on DCN datasets. On Facebook-Pod-B, Utimac leads all baselines across the full observation range: its wMAPE advantage

Rethinking Traffic Matrix Completion: Estimate the Process, Not the Entries

over the best competing method (ImputeFormer) is 12.3% at 𝑝obs =0.3 and narrows to 2.5% at 𝑝obs =0.9. This trend directly reflects the theoretical guarantee established in Section 3: by analytically characterizing the statistical structure of DCN traffic, Utimac constructs a statistically principled completion that is most informative precisely when observations are scarce, while learning-based methods require more supervision to close the gap. On Facebook-ToR-A, the lead is decisive and stable: Utimac’s MAE is 32.9% below PSW-I at 𝑝obs =0.3 and remains 31.3% lower at 𝑝obs =0.9. Diffusion-TM exhibits severe RMSE instability on this dataset, reaching 0.89 at 𝑝obs =0.9 (9.5× Utimac). Figure 2 confirms that Utimac’s lower mean error reflects uniformly better per-entry recovery: its 95th-percentile residual on Facebook-Pod-B is 0.288 at 𝑝obs =0.3, versus 0.368 for PSW-I and 0.343 for Diffusion-TM. Figure 3 evaluates recovery on true missing burst flows (𝛼=2.0, 𝛽=0.8). On Facebook-Pod-B, Utimac’s Burst-wMAPE advantage over PSW-I is 14.5% at 𝑝obs =0.3, narrowing to 1.2% at 𝑝obs =0.9, confirming that the inference advantage is most valuable when observations are scarce. On Facebook-ToRA, Utimac maintains the lowest Burst-wMAPE throughout; ImputeFormer shows the weakest burst recovery at medium-to-high 𝑝obs on Pod-B. Accuracy of the estimation intervals is further validated via PICP97.5 in Appendix F.2.

5

Conclusion

This paper presents Utimac, an uncertainty-aware traffic matrix completion method that models log-domain traffic as a joint Gaussian principal component plus a Laplacian sparse deviation, casting completion as parameter inference solved by block coordinate descent with closed-form predictive intervals. Experiments on real-world DCN and WAN datasets show that Utimac consistently outperforms representative baselines, with its advantage most pronounced under sparse observations.

References [1] Theodore W Anderson and Donald A Darling. 1954. A test of goodness of fit. Journal of the American statistical association 49, 268 (1954), 765–769. [2] Theophilus Benson, Aditya Akella, and David A Maltz. 2010. Network traffic characteristics of data centers in the wild. In Proceedings of the 10th ACM SIGCOMM conference on Internet measurement. 267–280. [3] Christopher M Bishop and Nasser M Nasrabadi. 2006. Pattern recognition and machine learning. Vol. 4. Springer. [4] Christopher Canel, Balasubramanian Madhavan, Srikanth Sundaresan, Neil Spring, Prashanth Kannan, Ying Zhang, Kevin Lin, and Srinivasan Seshan. 2024. Understanding incast bursts in modern datacenters. In Proceedings of the 2024 ACM on Internet Measurement Conference. 674–680. [5] Benoit Claise. 2004. Cisco systems netflow services export version 9. Technical Report.

[6] Benoit Claise, Brian Trammell, and Paul Aitken. 2013. Specification of the IP flow information export (IPFIX) protocol for the exchange of flow information. Technical Report. [7] RALPH D’agostino and Egon S Pearson. 1973. Tests for departure from normality. Empirical results for the distributions of 𝑏 2 and √𝑏 . Biometrika 60, 3 (1973), 613–622. [8] Adithya Gangidi, Rui Miao, Shengbao Zheng, Sai Jayesh Bondu, Guilherme Goes, Hany Morsy, Rohit Puri, Mohammad Riftadi, Ashmitha Jeevaraj Shetty, Jingyi Yang, et al. 2024. Rdma over ethernet for distributed training at meta scale. In Proceedings of the ACM SIGCOMM 2024 Conference. 57–70. [9] Ehab Ghabashneh, Yimeng Zhao, Cristian Lumezanu, Neil Spring, Srikanth Sundaresan, and Sanjay Rao. 2022. A microscopic view of bursts, buffer contention, and loss in data centers. In Proceedings of the 22nd ACM Internet Measurement Conference. 567–580. [10] Gonca Gürsun and Mark Crovella. 2012. On traffic matrix completion in the internet. In Proceedings of the 2012 internet measurement conference. 399–412. [11] Yuliang Li, Rui Miao, Changhoon Kim, and Minlan Yu. 2016. {FlowRadar}: A better {NetFlow} for data centers. In 13th USENIX symposium on networked systems design and implementation (NSDI 16). 311–324. [12] Zaoxing Liu, Antonis Manousis, Gregory Vorsanger, Vyas Sekar, and Vladimir Braverman. 2016. One sketch to rule them all: Rethinking network flow monitoring with univmon. In Proceedings of the 2016 ACM SIGCOMM Conference. 101–114. [13] Hao Mei, Junxian Li, Zhiming Liang, Guanjie Zheng, Bin Shi, and Hua Wei. 2023. Uncertainty-aware traffic prediction under missing data. In 2023 IEEE International Conference on Data Mining (ICDM). IEEE, 1223–1228. [14] Tong Nie, Guoyang Qin, Wei Ma, Yuewen Mei, and Jian Sun. 2024. ImputeFormer: Low rankness-induced transformers for generalizable spatiotemporal imputation. In Proceedings of the 30th ACM SIGKDD conference on knowledge discovery and data mining. 2260–2271. [15] Kun Qian, Yongqing Xi, Jiamin Cao, Jiaqi Gao, Yichi Xu, Yu Guan, Binzhang Fu, Xuemei Shi, Fangbo Zhu, Rui Miao, et al. 2024. Alibaba hpn: A data center network for large language model training. In Proceedings of the ACM SIGCOMM 2024 Conference. 691–706. [16] Yan Qiao, Kui Wu, and Xinyu Yuan. 2024. AutoTomo: Learning-based traffic estimator incorporating network tomography. IEEE/ACM Transactions on Networking 32, 6 (2024), 4644–4659. [17] Liang Qin, Xiyuan Liu, Wenting Wei, Chengbin Liang, and Huaxi Gu. 2024. Satformer: Accurate and robust traffic data estimation for satellite networks. Advances in Neural Information Processing Systems 37 (2024), 47530–47558. [18] William J Reed. 2006. The normal-Laplace distribution and its relatives. In Advances in distribution theory, order statistics, and inference. Springer, 61–74. [19] Matthew Roughan, Yin Zhang, Walter Willinger, and Lili Qiu. 2011. Spatio-temporal compressive sensing and internet traffic matrices (extended version). IEEE/ACM Transactions on Networking 20, 3 (2011), 662–676. [20] Vyas Sekar, Michael K Reiter, Walter Willinger, Hui Zhang, Ramana Rao Kompella, and David G Andersen. 2008. cSamp: A system for network-wide flow monitoring. (2008). [21] Samuel Sanford Shapiro and Martin B Wilk. 1965. An analysis of variance test for normality (complete samples). Biometrika 52, 3-4 (1965), 591–611. [22] Amin Tootoonchian, Monia Ghobadi, and Yashar Ganjali. 2010. OpenTM: traffic matrix estimator for OpenFlow networks. In International Conference on Passive and Active Network Measurement. Springer, 201–210.

Liu, Wang, et al. [23] Steve Uhlig, Bruno Quoitin, Jean Lepropre, and Simon Balon. 2006. Providing public intradomain traffic matrices to the research community. ACM SIGCOMM Computer Communication Review 36, 1 (2006), 83–86. [24] Hao Wang, Haoxuan Li, Xu Chen, Mingming Gong, Zhichao Chen, et al. 2025. Optimal transport for time series imputation. In The Thirteenth International Conference on Learning Representations. [25] Wenfeng Xia, Peng Zhao, Yonggang Wen, and Haiyong Xie. 2016. A survey on data center networking (DCN): Infrastructure and operations. IEEE communications surveys & tutorials 19, 1 (2016), 640–656. [26] Kun Xie, Yudian Ouyang, Xin Wang, Gaogang Xie, Kenli Li, Wei Liang, Jiannong Cao, and Jigang Wen. 2023. Deep adversarial tensor completion for accurate network traffic measurement. IEEE/ACM Transactions on Networking 31, 5 (2023), 2101–2116. [27] Kun Xie, Can Peng, Xin Wang, Gaogang Xie, Jigang Wen, Jiannong Cao, Dafang Zhang, and Zheng Qin. 2018. Accurate recovery of internet traffic data under variable rate measurements. IEEE/ACM transactions on networking 26, 3 (2018), 1137–1150. [28] Kun Xie, Jiazheng Tian, Gaogang Xie, Guangxing Zhang, and Dafang Zhang. 2021. Low cost sparse network monitoring based on block matrix completion. In IEEE INFOCOM 2021-IEEE Conference on Computer Communications. IEEE, 1–10. [29] Tong Yang, Jie Jiang, Peng Liu, Qun Huang, Junzhi Gong, Yang Zhou, Rui Miao, Xiaoming Li, and Steve Uhlig. 2018. Elastic sketch: Adaptive and fast network-wide measurements. In Proceedings of the 2018 Conference of the ACM Special Interest Group on Data Communication. 561–575. [30] Minlan Yu, Lavanya Jose, and Rui Miao. 2013. Software {Defined}{Traffic} Measurement with {OpenSketch}. In 10th USENIX symposium on networked systems design and implementation (NSDI 13). 29–42. [31] Xinyu Yuan, Yan Qiao, Zhenchun Wei, Zeyu Zhang, Minyue Li, Pei Zhao, Rongyao Hu, and Wenjing Li. 2025. Diffusion models meet network management: Improving traffic matrix analysis with diffusionbased approach. IEEE Transactions on Network and Service Management 22, 2 (2025), 1259–1275. [32] Qiao Zhang, Vincent Liu, Hongyi Zeng, and Arvind Krishnamurthy. 2017. High-resolution measurement of data center microbursts. In Proceedings of the 2017 Internet Measurement Conference. 78–85.

A

Empirical Validation of Joint Gaussianity for Log-Domain Traffic Vectors

This appendix validates the model assumption that the principal component approximation of the decomposition of the logarithmic domain flow vector follows a Gaussian joint distribution 𝒩(𝜇 𝑛 , Σ𝑛 ), using a locally stationary window of 𝐿 = 200 frames and 𝑑 = 56 OD pairs from the Facebook pod-b dataset.

A.1 Validation Criterion and Strategy The necessary and sufficient condition for joint Gaussianity is given by the Cramér–Wold theorem. PRoposition 1 (CRamÉR–Wold: NecessaRy and Sufficient Condition). 𝑍 ∈ ℝ𝑑 follows 𝒩(𝜇, Σ) if and only if

for every nonzero 𝑎 ∈ ℝ𝑑 ,

𝑎⊤ 𝑍 ∼ 𝒩(𝑎⊤ 𝜇, 𝑎⊤ Σ𝑎).

(44)

Directly verifying Proposition 1 requires applying a onedimensional normality test to every direction in ℝ𝑑 , which is statistically infeasible. We therefore adopt a two-stage progressive strategy: first, a global screening is performed using a necessary condition for joint Gaussianity; then, Proposition 1 is approximately verified on a finite, representative set of directions. A necessary condition for joint Gaussianity follows from the theoretical distribution of the Mahalanobis distance: if 𝑍 ∼ 𝒩(𝜇, Σ), then −1

2 (𝑧) = (𝑧 − 𝜇) 𝐷𝑀 ̂ ⊤ Σ̂ (𝑧 − 𝜇)̂ ∼ 𝜒 2 (𝑑).

(45)

The ordered squared Mahalanobis distances 2 2 𝐷𝑀,(1) ≤ ⋯ ≤ 𝐷𝑀,(𝑛)

(46)

are plotted against the theoretical quantiles

𝑘 − 0.5 (47) ) , 𝑘 = 1, … , 𝑛, 𝑛 to form the QQ plot. If the point cloud systematically deviates from the reference line 𝑦 = 𝑥 , the joint Gaussian hypothesis is directly rejected; if the data points closely track the reference line, the necessary condition in (45) is supported, though this result does not imply sufficiency. Section A.3 then approximately verifies the one-dimensional Gaussianity required by (44) along 181 representative directions. 𝑞𝑘 = 𝐹𝜒−1 2 (𝑑) (

A.2

Mahalanobis Distance QQ Plot

Figure 4 compares the Mahalanobis distance QQ plots for the raw domain and the log domain. In the raw domain, the point cloud exhibits a pronounced S-shaped deviation: the lower tail falls below the reference line 𝑦 = 𝑥 , and a small number of points in the upper tail lie significantly above the reference line (observed values near 92 at the theoretical quantile around 80), reflecting the heavy-tailed, nonGaussian nature of raw traffic. In the log domain, the point cloud closely tracks the reference line across the entire value range, with only minor deviations in the extreme lower tail, indicating that the log transformation effectively corrects the distributional shape to be compatible with joint Gaussianity.

A.3

Multi-Direction Projection Tests

To approximately verify the necessary and sufficient condition in Proposition 1, projection sequences 𝑠𝑡 = 𝑎⊤ 𝑧 𝑡 are computed on standardized log-domain data along 181 representative directions and subjected to the Shapiro–Wilk

Rethinking Traffic Matrix Completion: Estimate the Process, Not the Entries

(a) Raw domain

(b) Log domain

Figure 4: Mahalanobis distance QQ plots for the raw and log domains.

B

Empirical Evidence for Sparse Deviations around the Dominant Structure

This appendix provides empirical evidence for the observation in the main text that sparse deviations exist around a dominant traffic structure. We examine the statistical behavior of the raw traffic data from two complementary perspectives: low-dimensional projection and distributional fitting. The results show that datacenter traffic is dominated by a large number of small flows and a small number of sparse but off-center flows. We conduct a FastICA experiment on Facebook-Pod-B, projecting traffic samples from all five dataset splits (three

Figure 5: Normal QQ plots for six representative projection directions in the log domain.

,&    

,&

test [21], the D’Agostino 𝐾 2 test [7], and the Anderson– Darling test [1] (𝛼 = 0.05). The tested directions cover four categories: coordinate directions 𝑒 𝑘 (56 directions), the top 5 eigenvectors of the sample covariance matrix (PCA, 5 directions), pairwise combination directions (𝑒 𝑖 ± 𝑒 𝑗 )/√2 (20 directions), and random directions sampled uniformly from the unit sphere (100 directions). As the most representative subset for Proposition 1, the PCA directions achieve a pass rate of 4/5 (80%), the random directions achieve 84% (84/100), and the two categories combined yield 88/105 = 83.8%. The coordinate directions have a lower pass rate (17/56, 30.4%), consistent with individual OD pairs exhibiting local heavy-tail perturbations caused by short-term traffic bursts; such perturbations are averaged out in PCA and randomly mixed directions, allowing the joint Gaussian structure to emerge. Figure 5 presents QQ plots for six representative passing directions, covering the PCA, coordinate, and random categories. The point clouds closely track the OLS best-fit line in each panel, consistent with the one-dimensional Gaussianity required by (44).

    

'DWDVHW



'DWDVHW









'DWDVHW





'DWDVHW



'DWDVHW

Figure 6: FastICA projection of Facebook-Pod-B traffic onto IC1 and IC2, colour-coded by dataset split (three training sets and two validation sets). The majority of samples form a dense cluster near the origin, corresponding to the dominant low-magnitude traffic structure; a small number of scattered outliers represent sparse but off-center flows. training and two validation) jointly onto the first two independent components (IC1 and IC2), colour-coded by split identity. As shown in Figure 6, samples from all five splits overlap in a compact elliptical cluster centred at the origin, confirming that the dominant traffic structure is consistent and stationary across splits. A small fraction of points appear as sparse outliers well beyond the central mass, corresponding to sparse but off-center flows that are large relative to the bulk of OD entries within the same time frame but rare in occurrence. The cross-split consistency of both the central cluster and the peripheral outliers indicates that this dominant-plus-sparse pattern is a stable property of the dataset rather than an artifact of any particular split.

Liu, Wang, et al.

'HQVLW\

3.5×10 5 3×10 5

1.5×10 5

2.5×10 5 1.5×10 5 1×10 5

5×10 6

5×10 6

î





7UDIILFYDOXH D3RGWR3RG (PSLULFDO3')



î

/RJ1RUPDO

Conditional Gaussian Distribution of the Unobserved Principal Component (𝒪)

2×10 5

1×10 5



C.1

4×10 5

2×10 5

î



7UDIILFYDOXH E3RGWR3RG

î

Given the frame-level optimal deviation term 𝑜𝑡̂ returned by the deviation subproblem (37), the additive decomposition of the log-domain model (cf. (24)) constrains the principal component on the observed positions to satisfy

1RUPDO /DSODFH

(𝒪)

𝑢̂𝑡 Figure 7: Marginal distribution fitting on two representative OD dimensions from Facebook-Pod-B. (a) Pod5→Pod7: LogNormal and Normal–Laplace curves nearly coincide, indicating negligible sparse deviation. (b) Pod6→Pod0: the empirical distribution exhibits a sharper peak and a heavier left shift; Normal–Laplace provides a closer fit than LogNormal alone by capturing the sparse-deviation component.

𝑍𝑡 = 𝑈𝑡 + 𝑂𝑡 ,

(ℳ)

𝑚𝑡 = 𝜇̂𝑛

where 𝑈𝑡 denotes the dominant log-domain traffic structure and 𝑂𝑡 denotes the local deviation. Since 𝑂𝑡 is expected to satisfy the sparse pattern that most entries are close to zero while only a few entries have large magnitudes, we impose a Laplace prior on the deviation term:

𝑉𝑡 = Σ̂ 𝑛

(49)

This prior encourages sparsity through its sharp peak at zero, while its heavy-tailed behavior allows a small number of large deviations to be represented. Figure 7 shows marginal distribution fits for two representative OD dimensions from the Facebook-Pod-B dataset, within a locally stationary window of 𝐿 = 200 frames and 𝑑 = 56 OD pairs.

C

Imputation Details

This appendix derives, for each missing traffic entry, a 95% imputation estimated interval that quantifies the uncertainty of the imputed log-domain value under the estimated statistical model. All derivations are conducted within a fixed local window 𝑊𝑛 and a fixed time instant 𝑡 ∈ 𝑊𝑛 , conditioned on the estimated shared parameters 𝜃𝑛̂ = (𝜇̂𝑛 , Σ̂ 𝑛 , 𝜆̂ 𝑛 ) and the frame-level optimal deviation (Ω ) term 𝑜𝑡̂ 𝑡 obtained from the block coordinate descent in

Section 3.3. For brevity, within this appendix we abbreviate 𝒪 ∶= Ω𝑡 and ℳ ∶= Ωmis 𝑡 ; these abbreviations are used only locally and do not redefine any symbol of the main text.

,

(50)

(𝒪)

| 𝑈𝑡

(𝒪)

= 𝑢̂𝑡

∼ 𝒩(𝑚𝑡 , 𝑉𝑡 ),

(51)

where the conditional mean 𝑚𝑡 ∈ ℝ|ℳ| and conditional covariance 𝑉𝑡 ∈ ℝ|ℳ|×|ℳ| are

(48)

𝑝(𝑂𝑡 ) ∝ exp (−𝜆‖𝑂𝑡 ‖1 ) .

(𝒪)

− 𝑜𝑡̂

which we take as the conditioning value for the missing positions. Under the model in (25), the principal component satisfies 𝑈𝑡 ∼ 𝒩(𝜇̂𝑛 , Σ̂ 𝑛 ) within 𝑊𝑛 . Applying the standard conditional distribution formula for multivariate Gaussians [3] to (50) yields

𝑈𝑡 Based on this observation, we adopt an additive decomposition in the log domain:

(𝒪)

= 𝑧𝑡

(ℳ)

(ℳ, 𝒪)

+ Σ̂ 𝑛

(ℳ, ℳ)

(𝒪) (𝒪, 𝒪) −1 (𝒪) ) (𝑢̂𝑡 − 𝜇̂𝑛 ), (Σ̂ 𝑛

(ℳ, 𝒪)

− Σ̂ 𝑛

(𝒪, 𝒪) −1 ̂ (𝒪, ℳ) . ) Σ𝑛

(Σ̂ 𝑛

(52) (53)

The subvectors and submatrices in (52)–(53) follow the notation of (12), with 𝒪 and ℳ substituted for Ω𝑡 and Ωmis 𝑡 .

C.2

Laplace Distribution for the Unobserved Deviation Component

With the distribution of the unobserved principal component established in (51), we next characterise the distribution of the unobserved deviation component. Under the model in (26), each entry of 𝑂𝑡 follows a univariate density proportional to exp(−𝜆𝑛 |𝑜|). Comparing this with the standard parameterisation of Laplace(0, 𝑏), whose 1 density is 2𝑏 exp(−|𝑥|/𝑏), identifies the scale parameter as the reciprocal of 𝜆𝑛 . Substituting the estimate 𝜆̂ 𝑛 therefore gives 1 𝑏 ∶= > 0, (54) 𝜆̂ 𝑛 so that the deviation component on the missing positions has the marginal distribution (ℳ)

𝑂𝑡

∼ ∏ Laplace(0, 𝑏).

(55)

𝑗∈ℳ

Furthermore, since 𝑈𝑡 ⟂ 𝑂𝑡 for each 𝑡 under the model in Section 2.1, this independence is preserved under conditioning: the conditional variable in (51) remains independent of (55).

Rethinking Traffic Matrix Completion: Estimate the Process, Not the Entries

C.3

Conditional estimated Density for Unobserved Positions

admits the closed form

Combining the distributions in (51) and (55) through the additive decomposition in (24), the log-domain variable on the missing positions satisfies (ℳ)

𝑍𝑡

(ℳ)

= 𝑈𝑡

(ℳ)

+ 𝑂𝑡

.

(𝒪)

(𝒪)

, 𝑜𝑡̂

1 −|𝑧𝑗 −𝑢𝑗 |/𝑏 d𝑢, ∫ |ℳ| 𝒩(𝑢 ∣ 𝑚𝑡 , 𝑉𝑡 ) ∏ 2𝑏 𝑒

2 𝑠𝑡,𝑗 𝑧 − 𝑚𝑡,𝑗 1 + exp( 2 + ) Φ(−𝛼𝑡,𝑗 (𝑧) − 𝜅𝑡,𝑗 ), (64) 2 𝑏 2𝑏

where Φ(⋅) is the standard normal CDF.

(57)

+∞

𝐹𝑡,𝑗 (𝑧) = ∫

𝑗∈ℳ

Per-Coordinate Marginal Distribution and Closed-Form CDF

Since the 95% interval is constructed separately for each missing coordinate, we reduce (57) to its one-dimensional marginals. For each 𝑗 ∈ ℳ, the marginal at coordinate 𝑗 decomposes as

𝑍𝑡 (𝑗) = 𝑈𝑡 (𝑗) + 𝑂𝑡 (𝑗),

(58)

where the two components are independent and distributed as (𝒪) (𝒪) 𝑈𝑡 (𝑗) | 𝑈𝑡 = 𝑢̂𝑡

2 ), ∼ 𝒩(𝑚𝑡,𝑗 , 𝑠𝑡,𝑗

𝑂𝑡 (𝑗) ∼ Laplace(0, 𝑏),

(59)

2 ∶= (𝑉 ) . 𝑠𝑡,𝑗 𝑡 𝑗𝑗

(60)

Consequently, 𝑍𝑡 (𝑗) follows the convolution of a Gaussian and a Laplace distribution [18], with variance 2 + 2𝑏 2 , Var(𝑍𝑡 (𝑗)) = 𝑠𝑡,𝑗

(61)

where we used the independence of 𝑈𝑡 (𝑗) and 𝑂𝑡 (𝑗) together with the identity Var(Laplace(0, 𝑏)) = 2𝑏 2 . To present the closed-form CDF compactly, introduce the auxiliary quantities 𝑧 − 𝑚𝑡,𝑗 𝑠𝑡,𝑗 . (62) 𝛼𝑡,𝑗 (𝑧) ∶= , 𝜅𝑡,𝑗 ∶= 𝑠𝑡,𝑗 𝑏 PRoposition 2 (Closed-foRm CDF). The CDF of 𝑍𝑡 (𝑗), (𝒪) (𝒪) 𝐹𝑡,𝑗 (𝑧) ∶= 𝑃(𝑍𝑡 (𝑗) ≤ 𝑧 | 𝑧𝑡 , 𝑜𝑡̂ , 𝜃𝑛̂ ) ,

(63)

𝑧 − 𝑜 − 𝑚𝑡,𝑗 𝑠𝑡,𝑗

)

1 −|𝑜|/𝑏 𝑒 d𝑜. 2𝑏

(65)

Split the integration domain at 𝑜 = 0: 𝐹𝑡,𝑗 (𝑧) = 𝐼1 (𝑧) + 𝐼2 (𝑧), where 𝐼1 and 𝐼2 denote the integrals over [0, +∞) and (−∞, 0), respectively. Evaluation of 𝐼1 . Applying integration by parts with 𝑢1 = Φ(𝛼𝑡,𝑗 (𝑧) − 𝑜/𝑠𝑡,𝑗 ) and d𝑣1 = 𝑒 −𝑜/𝑏 d𝑜 , +∞

2𝑏 𝐼1 (𝑧) = ∫

0

Φ(𝛼𝑡,𝑗 (𝑧) − 𝑜/𝑠𝑡,𝑗 ) 𝑒 −𝑜/𝑏 d𝑜

= 𝑏 Φ(𝛼𝑡,𝑗 (𝑧)) −

𝑏 𝐽 (𝑧), 𝑠𝑡,𝑗 1

(66)

where +∞

𝐽1 (𝑧) ∶= ∫

0

with the marginal parameters defined as

𝑚𝑡,𝑗 ∶= (𝑚𝑡 )𝑗 ,

Φ(

−∞

where 𝒩(𝑢 ∣ 𝑚𝑡 , 𝑉𝑡 ) denotes the Gaussian density with mean 𝑚𝑡 and covariance 𝑉𝑡 evaluated at 𝑢 .

C.4

2 𝑠𝑡,𝑗 𝑧 − 𝑚𝑡,𝑗 1 exp( 2 − ) Φ(𝛼𝑡,𝑗 (𝑧) − 𝜅𝑡,𝑗 ) 2 𝑏 2𝑏

PRoof. By the total probability formula and the Laplace density of 𝑂𝑡 (𝑗),

, 𝜃𝑛̂ ) =

ℝ

−

(56)

By the independence established in Section C.2, the joint (ℳ) conditional estimated density of 𝑍𝑡 is the convolution of the Gaussian in (51) and the product Laplace in (55):

𝑝(𝑧 (ℳ) | 𝑧𝑡

𝐹𝑡,𝑗 (𝑧) = Φ(𝛼𝑡,𝑗 (𝑧))

𝜙(𝛼𝑡,𝑗 (𝑧) − 𝑜/𝑠𝑡,𝑗 ) 𝑒 −𝑜/𝑏 d𝑜

(67)

and 𝜙(⋅) denotes the standard normal density. Completing the square in the combined exponent of 𝐽1 (𝑧) gives 2

− 12 (𝛼𝑡,𝑗 (𝑧) − 𝑜/𝑠𝑡,𝑗 ) − 𝑏𝑜 = −

2

[𝑜 − 𝑠𝑡,𝑗 (𝛼𝑡,𝑗 (𝑧) − 𝜅𝑡,𝑗 )] 2 2𝑠𝑡,𝑗

+

2 𝑠𝑡,𝑗

2𝑏

− 2

𝑧 − 𝑚𝑡,𝑗 (68) . 𝑏

Substituting 𝑤 = [𝑜 − 𝑠𝑡,𝑗 (𝛼𝑡,𝑗 (𝑧) − 𝜅𝑡,𝑗 )]/𝑠𝑡,𝑗 so that 𝑜 = 0 corresponds to 𝑤 = 𝜅𝑡,𝑗 − 𝛼𝑡,𝑗 (𝑧), and noting that the new lower limit yields the standard normal tail probability,

𝐽1 (𝑧) = exp(

2 𝑠𝑡,𝑗

2𝑏 2

−

𝑧 − 𝑚𝑡,𝑗 𝑏

) Φ(𝛼𝑡,𝑗 (𝑧) − 𝜅𝑡,𝑗 ).

(69)

Substituting (69) into (66) and simplifying,

𝐼1 (𝑧) = 21 Φ(𝛼𝑡,𝑗 (𝑧)) − 12 exp(

2 𝑠𝑡,𝑗

𝑧 − 𝑚𝑡,𝑗

2𝑏

𝑏

− 2

) Φ(𝛼𝑡,𝑗 (𝑧) − 𝜅𝑡,𝑗 ).

(70)

Liu, Wang, et al.

𝐼2 (𝑧) = 12 Φ(𝛼𝑡,𝑗 (𝑧)) 2 𝑠𝑡,𝑗 𝑧 − 𝑚𝑡,𝑗 + 12 exp( 2 + ) Φ(−𝛼𝑡,𝑗 (𝑧) − 𝜅𝑡,𝑗 ). 𝑏 2𝑏

Burst-MAE

0.4 0.2

□

PSW-I

Burst-wMAPE 1.00 0.75

0.75

Observation Rate

(71)

Burst-RMSE 1.00

0.3 0.4 0.5 0.6 0.7 0.8 0.9

Summing (70) and (71) yields (64).

C.5

0.6

GEANT

Evaluation of 𝐼2 . Substituting 𝑜 ↦ −𝑜 maps 𝐼2 to an integral over [0, +∞) whose integrand differs from that of 𝐼1 only in the sign of 𝑜/𝑠𝑡,𝑗 . Applying the same integration-by-parts, completing-the-square, and substitution procedure as for 𝐼1 , with this sign reversed throughout, yields

0.50

0.50

0.25

0.25 0.3 0.4 0.5 0.6 0.7 0.8 0.9

Observation Rate

ImputeFormer

Diffusion-TM

0.3 0.4 0.5 0.6 0.7 0.8 0.9

Observation Rate Utimac

Figure 9: Burst-MAE, Burst-RMSE, and Burst-wMAPE vs. 𝑝obs on GÉANT.

Imputation Estimated Interval

Equipped with the closed-form marginal CDF 𝐹𝑡,𝑗 from Proposition 2, we now construct a 97.5% interval for each missing coordinate in the log domain and map it to the original traffic domain. The marginal distributions in (59) are each symmetric about their respective centres. By the independence of 𝑈𝑡 (𝑗) and 𝑂𝑡 (𝑗), the sum 𝑍𝑡 (𝑗) − 𝑚𝑡,𝑗 is therefore also symmetric about zero. Consequently, the 97.5% interval in the log domain, defined as the interval between the 1.25% and 98.75% quantiles of 𝑍𝑡 (𝑗), takes the symmetric form (𝑧)

PI0.975 (𝑡, 𝑗) = [𝑚𝑡,𝑗 − 𝛿𝑡,𝑗 , 𝑚𝑡,𝑗 + 𝛿𝑡,𝑗 ],

D

Dataset Statistics

Table 1 summarises the three datasets used in our experiments. GÉANT [23] is a wide-area network (WAN) dataset from the pan-European GÉANT backbone. Facebook-PodB, and Facebook-ToR-A are datacenter network (DCN) datasets from Facebook’s production infrastructure, representing pod-level and rack-level traffic respectively.

(72)

where 𝛿𝑡,𝑗 > 0 is the unique solution to

𝐹𝑡,𝑗 (𝑚𝑡,𝑗 + 𝛿𝑡,𝑗 ) = 0.9875.

(73) Table 1: Dataset statistics († WAN; ‡ DCN).

By the strictly increasing transformation in (22), quantile ordering is preserved, and the 97.5% imputation estimated interval in the original traffic domain is (𝑥) PI0.975 (𝑡, 𝑗) = [𝑒 𝑚𝑡,𝑗 −𝛿𝑡,𝑗 − 𝜀𝑗 ,

𝑒

𝑚𝑡,𝑗 +𝛿𝑡,𝑗

− 𝜀𝑗 ],

Dataset

(74)

GÉANT† 22 ‡ Facebook-Pod-B 8 Facebook-ToR-A‡ 155

where 𝜀𝑗 is the 𝑗 -th component of the offset vector 𝜀 in (3). Equation (74) serves as the per-entry uncertainty output of Utimac and is evaluated in Section 4.

MAE(log scale)

RMSE(log scale)

GEANT

104

Observation Rate PSW-I

wMAPE

0.50

100 0.3 0.4 0.5 0.6 0.7 0.8 0.9

462 484 56 64 23,870 24,025

0.75

102

100

72 56 7,194

𝑑

1.00

104 102

𝑁 Links OD flows

0.25 0.3 0.4 0.5 0.6 0.7 0.8 0.9

Observation Rate

ImputeFormer

Diffusion-TM

0.3 0.4 0.5 0.6 0.7 0.8 0.9

Observation Rate Utimac

Figure 8: Overall MAE, RMSE, and wMAPE vs. 𝑝obs on GÉANT.

E

Evaluation Metric Definitions

All metrics are computed exclusively on missing positions Ωmiss = {𝑘 ∣ 𝐵𝑡 (𝑘)=0}, where 𝐵𝑡 is the observation mask 𝑡 at time 𝑡 . Let 𝑥̂𝑡 (𝑘) and 𝑥𝑡 (𝑘) denote the imputed and true traffic values at position 𝑘 and time 𝑡 , and let 𝜀>0 be a small constant preventing division by zero.

Rethinking Traffic Matrix Completion: Estimate the Process, Not the Entries

Normalized Absolute Residual(log scale)

Overall Imputation Accuracy ∑ ∑ |𝑥̂𝑡 (𝑘) − 𝑥𝑡 (𝑘)|

100

,

∑ |Ωmiss | 𝑡

(75)

𝑡

√ √ √ ∑ ∑ (𝑥̂𝑡 (𝑘) − 𝑥𝑡 (𝑘))2 √ √ 𝑡 𝑘∈Ωmiss √ 𝑡 √ RMSE = √ , √ √ ∑ |Ωmiss | 𝑡 𝑡 √ 𝑡 𝑘∈Ωmiss 𝑡

.

∑ ∑ |𝑥𝑡 (𝑘)|+𝜀

10 2 10 4 10 6 0.3

(76)

PSW-I

∑ ∑ |𝑥̂𝑡 (𝑘) − 𝑥𝑡 (𝑘)| wMAPE =

GEANT

MAE =

𝑡 𝑘∈Ωmiss 𝑡

(77)

𝑡 𝑘∈Ωmiss 𝑡

0.4

0.5

ImputeFormer

∑ ∑ |𝑥̂𝑡 (𝑘) − 𝑥𝑡 (𝑘)| Burst-MAE =

𝑡 𝑘∈ℬmiss 𝑡

∑ |ℬmiss | 𝑡

(78)

(𝛼,𝛽) define the missing burst set as ℬmiss = ℬ𝑡 ∩ Ωmiss , with 𝑡 𝑡 miss predicted counterpart ℬ̂ 𝑡 obtained by applying the same definition to 𝑥̂𝑡 .

Burst Flow Detection ∑ |ℬ̂ miss ∩ ℬmiss | 𝑡 𝑡 ,

,

(83)

∑ ∑ |𝑥̂𝑡 (𝑘) − 𝑥𝑡 (𝑘)| Burst-wMAPE =

𝑡 𝑘∈ℬmiss 𝑡

∑ ∑ |𝑥𝑡 (𝑘)|+𝜀

.

(84)

𝑡 𝑘∈ℬmiss 𝑡

Burst-wMAPE is the headline indicator as it directly reflects relative recovery quality on high-volume critical entries.

F

(80)

This appendix presents two sets of supplementary results omitted from the main text for space: (i) complete imputation performance on the GÉANT WAN dataset with the same three figures as Section 4.3, and (ii) burst flow detection metrics (Precision, Recall, F1) for all three datasets.

∑ |ℬ̂ miss ∩ ℬmiss | 𝑡 𝑡 ∑ |ℬmiss |+ 𝜀 𝑡

(82)

(79)

𝑡

𝑡

,

√ √ √ ∑ ∑ (𝑥̂𝑡 (𝑘) − 𝑥𝑡 (𝑘))2 √ √ 𝑡 𝑘∈ℬmiss √ 𝑡 √ Burst-RMSE = √ , √ √ ∑ |ℬmiss | 𝑡 𝑡 √

with dominance multiplier 𝛼>1 and majority threshold (𝛼,𝛽) 𝛽∈(0, 1). Let ℬ𝑡 denote the burst flow set at time 𝑡 , and

Recall =

Utimac

𝑡

1 ∑ 1[𝑥 (𝑘) ≥ 𝛼 𝑥𝑡 (ℓ)], 𝑑−1 ℓ≠𝑘 𝑡

∑ |ℬ̂ miss |+ 𝜀 𝑡

Diffusion-TM

0.9

Evaluated only on true missing burst positions ℬmiss : 𝑡

A flow 𝑘 at time 𝑡 is defined as a burst flow if its dominance ratio 𝑟𝑡 (𝑘; 𝛼) ≥ 𝛽 , where

𝑡

0.8

Figure 10: Normalized absolute residual distributions (5th– 95th percentile) vs. 𝑝obs on GÉANT.

Burst Flow Definition

Precision =

0.7

Burst Flow Imputation Accuracy

wMAPE is preferred over per-entry MAPE because it aggregates numerator and denominator separately, making it robust to near-zero flows common in heavy-tailed traffic distributions.

𝑟𝑡 (𝑘; 𝛼) =

0.6

Observation Rate

𝑡

2 ⋅ Precision ⋅ Recall F1 = . (81) Precision + Recall + 𝜀 Recall and F1 are the primary indicators: Recall measures coverage of true burst flows, while F1 balances false alarms and missed detections.

Extended Experimental Results

F.1 GÉANT Results Figures 8–9 report the same three figures as in the main text, now for the GÉANT WAN backbone dataset (𝑁 =22, 𝑑=484, 15-min intervals). On GÉANT, the ranking of methods differs from the DCN setting in an instructive way. For overall imputation,

Liu, Wang, et al.

PSW-I achieves the lowest MAE across the full observation range (e.g., 0.030 at 𝑝obs =0.3), while Utimac ranks second with MAE = 0.037, staying well below ImputeFormer (up to 0.211 at 𝑝obs =0.9) and far below Diffusion-TM, whose MAE collapses to values exceeding 104 due to numerical instability in the per-dimension min-max normalization applied to WAN-scale data. In terms of wMAPE, Utimac achieves 0.235 at 𝑝obs =0.3 and remains stable around 0.231–0.235 across the full range, while ImputeFormer degrades from 0.994 to 1.042 and Diffusion-TM stays above 1.13 throughout. These results demonstrate that Utimac generalizes well beyond the DCN setting: even on a WAN backbone with qualitatively different traffic characteristics, it produces stable and competitive completions, trailing only PSW-I (an alignment-based method that does not require traffic-specific structural assumptions) while substantially outperforming the two neural baselines. The residual distribution in Figure 10 reinforces this picture: Utimac’s 95th-percentile residual (0.164 at 𝑝obs =0.3) is close to PSW-I (0.131) and far below Diffusion-TM, whose extreme tail residuals (p95 up to 24.97) confirm the numerical instability visible in Figure 8. For burst flow imputation on GÉANT (Figure 9), Utimac achieves a stable Burst-wMAPE around 0.209–0.212 across all observation rates, ranking second behind PSW-I (0.077–0.108) and substantially ahead of ImputeFormer (0.674–0.994) and Diffusion-TM (0.991–1.022). The gap to PSW-I on GÉANT reflects the particular advantage of alignment-based methods on low-dimensional WAN traffic matrices (𝑑=484) where distributional proximity can be computed accurately with small batches, an advantage that diminishes as dimensionality grows.

Table 2: Burst flow detection on Facebook-ToR-A. Precision ≈1.000 for all methods (see text).†

F.2 Burst Flow Detection

On GÉANT (Table 4), burst detection is considerably easier than on the DCN datasets due to the lower dimensionality. PSW-I leads with F1 ≈0.92–0.94, followed closely by Utimac at F1 ≈0.90–0.91, a gap of only 1.7–4.4%. DiffusionTM achieves the weakest detection performance with F1 ≈ 0.55, consistent with its numerical instability on GÉANT observed in Section F.1.

Tables 2 and 4 report Precision, Recall, and F1 for burst flow detection (𝛼=2.0, 𝛽=0.8) on Facebook-ToR-A and GÉANT. Burst flow detection differs from burst flow imputation (Section 4.3): detection measures whether a method correctly identifies which flows are bursty based on its imputed values, independent of the numerical recovery error. On Facebook-ToR-A (Table 2), all methods obtain Precision ≈1.000, a statistical artifact of the dataset’s extreme sparsity (𝑑=24,025): at this dimensionality, burst flows constitute a negligibly small fraction of all non-zero entries, so virtually all predicted burst positions are genuine. Recall is therefore the sole discriminative indicator. At 𝑝obs =0.9, Utimac’s Recall jumps to 1.000 while PSW-I reaches 0.474 and Diffusion-TM 0.437, revealing that Utimac’s inference framework can recover the full burst-flow structure when sufficient observations are available even in a high-dimensional sparse setting.

Observation rate 𝑝obs Method

0.3

0.4

0.5

0.6

0.7

0.8

Recall PSW-I ImputeFormer Diffusion-TM Utimac

.202 .149 .196 .151

.199 .160 .205 .157

.201 .168 .217 .167

.217 .174 .234 .186

.257 .178 .260 .220

.345 .474 .182 .185 .304 .437 .282 1.000

F1 PSW-I ImputeFormer Diffusion-TM Utimac

.337 .259 .328 .262

.332 .276 .340 .271

.334 .287 .356 .287

.357 .296 .379 .314

.409 .302 .412 .361

.513 .643 .307 .312 .466 .608 .440 1.000

†

0.9

Precision ≈1.000 for all methods due to extreme burst-flow sparsity at 𝑑=24,025; Recall and F1 are the sole discriminative indicators.

Table 3: PICP97.5 (%) of Utimac vs. observation rate 𝑝obs on Facebook-Pod-B and Facebook-ToR-A. Values close to 97.5 indicate estimated intervals. Observation rate 𝑝obs Dataset

0.3

0.4

0.5

0.6

0.7

0.8

0.9

Facebook-Pod-B .968 .951 .955 .958 .959 .966 .973 Facebook-ToR-A .970 .972 .973 .974 .975 .975 .976

F.3 Estimated Interval Coverage Beyond point estimates, Utimac derives closed-form estimated intervals for every missing entry from the analytically characterised Gaussian–Laplace posterior (Section 3). We evaluate calibration via the Prediction Interval Coverage Probability at the 97.5% nominal level (PICP97.5 ), defined as the fraction of true missing values that fall within the predicted 97.5% interval. A well-calibrated method should achieve PICP97.5 ≈ 97.5% regardless of the observation rate.

Rethinking Traffic Matrix Completion: Estimate the Process, Not the Entries

Table 4: Burst flow detection on GÉANT. Observation rate 𝑝obs Method

0.3

0.4

0.5

0.6

0.7

0.8

0.9

Precision PSW-I ImputeFormer Diffusion-TM Utimac

.924 .930 .606 .639 .583 .592 .920 .923

.934 .694 .590 .928

.938 .787 .587 .930

.942 .883 .583 .932

.944 .919 .580 .935

.948 .923 .586 .934

Recall PSW-I ImputeFormer Diffusion-TM Utimac

.921 .924 .851 .905 .544 .543 .894 .891

.928 .927 .534 .889

.933 .925 .524 .886

.935 .901 .519 .883

.936 .867 .515 .880

.941 .802 .523 .875

F1 PSW-I ImputeFormer Diffusion-TM Utimac

.922 .927 .708 .749 .563 .566 .907 .907

.931 .794 .561 .908

.935 .851 .554 .907

.939 .892 .549 .907

.940 .892 .546 .907

.945 .858 .553 .903

Table 3 reports PICP97.5 on Facebook-Pod-B and Facebook-ToR-A across all observation rates. On FacebookToR-A, coverage remains tightly around 97.5% across the full range of 𝑝obs , confirming that the posterior-derived intervals faithfully reflect the true uncertainty of the imputed values. On Facebook-Pod-B, coverage increases monotonically from 96.8% at 𝑝obs =0.3 to 97.3% at 𝑝obs =0.9: at very low observation rates the posterior uncertainty is larger and the interval width grows accordingly, so the nominal level is approached from below as more observations constrain the inference. This behaviour is consistent with the theoretical guarantee that the interval endpoints are computed from the exact marginal covariance; the convergence to the nominal level as 𝑝obs increases provides direct empirical support for the correctness of the closed-form estimated distribution.