Low-rank Distributional Matrix Completion

# Low-rank Distributional Matrix Completion
Source: [https://arxiv.org/html/2606.04176](https://arxiv.org/html/2606.04176)
###### Abstract

We study a distributional generalization of the matrix completion problem in which each entry of the target matrix is a probability distribution rather than a scalar\. In this setting, only a subset of matrix entries is observed, and even for observed entries, the underlying distributions are not directly accessible; instead, we observe finitely many samples drawn from them\. To represent distributional entries, we employ kernel mean embeddings and introduce a notion of Tucker rank for distribution\-valued matrices to capture their low\-rank structure\. The infinite\-dimensional nature of kernel embeddings poses significant methodological challenges\. To address this, we introduce functional unfolding operators that link the proposed distributional low\-rank structure to the classical Tucker rank for finite\-dimensional tensors\. Based on this framework, we propose a novel estimator for distributional matrix completion\. We establish non\-asymptotic error bounds that characterize the statistical performance of the estimator\. Extensive experiments on synthetic data and a real\-world application demonstrate the effectiveness of the proposed method\.

## 1Introduction

Low\-rank matrix completion addresses the problem of recovering missing entries from a partially observed and potentially contaminated data matrix, where the low\-rank assumption serves as a crucial structural regularizer to facilitate identification\. Due to its flexibility and effectiveness, low\-rank matrix completion has been widely applied in a broad range of domains, including recommendation systems\(Kang et al\.,[2016](https://arxiv.org/html/2606.04176#bib.bib17); Gurini et al\.,[2018](https://arxiv.org/html/2606.04176#bib.bib14); Chen et al\.,[2021](https://arxiv.org/html/2606.04176#bib.bib8); Zhang et al\.,[2021](https://arxiv.org/html/2606.04176#bib.bib40)\), signal processing\(Weng and Wang,[2012](https://arxiv.org/html/2606.04176#bib.bib36); Zhang and Zhang,[2020](https://arxiv.org/html/2606.04176#bib.bib41); Chen et al\.,[2023](https://arxiv.org/html/2606.04176#bib.bib5); Yuchi et al\.,[2023](https://arxiv.org/html/2606.04176#bib.bib39)\), image recovery\(Changjun et al\.,[2012](https://arxiv.org/html/2606.04176#bib.bib4); Cao et al\.,[2014](https://arxiv.org/html/2606.04176#bib.bib3); Zheng et al\.,[2024](https://arxiv.org/html/2606.04176#bib.bib43)\)and data imputation\(Chen et al\.,[2017](https://arxiv.org/html/2606.04176#bib.bib6),[2020](https://arxiv.org/html/2606.04176#bib.bib7); Xu et al\.,[2023](https://arxiv.org/html/2606.04176#bib.bib38)\)\. A substantial body of work has been devoted to the theoretical and algorithmic development of matrix completion\. However, most existing approaches are tailored to scalar\-valued matrices and do not directly extend to target matrices with non\-scalar entries\. In this work, we study matrix completion problems in which the entries are probability distributions\.

As a motivating example, consider modeling the distribution of daily taxi trip counts between pairs of origin and destination locations \(see Section[7](https://arxiv.org/html/2606.04176#S7)\) in order to capture both the volume and variability of traffic flows\. Another example arises in the analysis of quarterly earnings estimates for public companies, where the goal is to study the distribution of forecasts issued by different banks and traders\. In this setting, companies index the rows and fiscal quarters across multiple years index the columns; see Section 1\.1 ofFeitelberg et al\. \([2024](https://arxiv.org/html/2606.04176#bib.bib10)\)for a detailed discussion\.

For observed entries, the underlying distributions can be inferred from the corresponding samples\. In practice, however, some entries may be entirely unobserved\. In the taxi\-trip example, missingness may result from incomplete trip records or limited GPS coverage in certain locations\. In the earnings\-estimates example, firms may release information irregularly, leading to missing entries for particular quarters or time periods\. These challenges highlight the need for efficient methods to impute missing distributional entries\.

To the best of our knowledge, there is only one existing work that explicitly addresses distributional matrix completion\. InFeitelberg et al\. \([2024](https://arxiv.org/html/2606.04176#bib.bib10)\), the authors propose to impute missing entries using a nearest\-neighbor\-based approach\. Specifically, for a missing entry located at\(i,j\)\(i,j\), similar rows are identified by using averaged empirical22\-Wasserstein distances over the overlapping columns with available observations\. The missing distribution is then estimated by taking the empirical Wasserstein barycenter of the distributions associated with these neighboring rows in columnjj\. Their method relies on a simplified observation setting in which columnjjis observed for all selected neighbors, and each neighbor has the same number of observationsNNin columnjj\. Moreover, the final estimate is represented byNNsynthetic samples generated from the estimated distribution, which requiresNNto be sufficiently large in order to accurately recover the underlying distribution\. In addition, due to the computational challenges associated with evaluating22\-Wasserstein distances and computing Wasserstein barycenters for multivariate distributions, the proposed algorithm does not readily extend to multidimensional distributional settings\.

In contrast, our approach accommodates a more general observation pattern, allowing the number of samples to vary across observed entries, and naturally extends to multidimensional distributional data\. We represent distributional entries using kernel mean embeddings \(KMEs\), thereby embedding probability distributions into a Hilbert space and circumventing explicit manifold constraints\. Unlike the nearest\-neighbor strategy inFeitelberg et al\. \([2024](https://arxiv.org/html/2606.04176#bib.bib10)\), our method is formulated within an empirical risk minimization framework and enforces a global low\-rank structure on the distributional matrix through carefully designed regularization\.

The proposed low\-rank structure is motivated by Tucker decompositionTucker \([1966](https://arxiv.org/html/2606.04176#bib.bib32)\)in tensor analysis, but operates in a more challenging setting in which one mode is infinite\-dimensional\. To tackle this, we introduce unfolding operators that map the distributional matrix into tensor products of two spaces \(see Definition[3\.1](https://arxiv.org/html/2606.04176#S3.Thmtheorem1)\), and impose nuclear\-norm regularization on the corresponding induced linear operators to promote low\-rankness\. This modeling strategy simultaneously encourages intrinsic dimension reduction in the functional component of the kernel mean embeddings and induces a low\-rank structure in the associated expectation functional matrices \(see definition in \([5](https://arxiv.org/html/2606.04176#S3.E5)\)\)\. Although the infinite\-dimensional nature of the distributional representations poses optimization challenges, we show that, with appropriate regularization and by leveraging the reproducing property of the underlying kernel, the solution admits a finite\-dimensional representation\. As a result, the proposed optimization problem can be solved efficiently using convex algorithms\. On the theoretical side, we establish error bounds for the aggregate squared maximum mean discrepancy \(MMD\) across all matrix entries\. We show that the resulting rate matches that of classical nuclear\-norm\-regularized matrix completion for \(approximately\) low\-rank matrices, despite the additional technical challenges arising from replacing scalar entries with functional representations\. Notably, unlikeFeitelberg et al\. \([2024](https://arxiv.org/html/2606.04176#bib.bib10)\), our consistency guarantees do not necessarily require the number of samples for each observed entry to diverge\.

## 2Problem Setup

Let𝒳⊆ℝd\\mathcal\{X\}\\subseteq\\mathbb\{R\}^\{d\}, and denote by𝒫​\(𝒳\)\\mathcal\{P\}\(\\mathcal\{X\}\)the collection of probability distributions of random variables taking values in𝒳\\mathcal\{X\}\. We consider a distribution\-valued matrixρ∗∈𝒫​\(𝒳\)m1×m2\\rho^\{\*\}\\in\\mathcal\{P\}\(\\mathcal\{X\}\)^\{m\_\{1\}\\times m\_\{2\}\}, whose\(i,j\)\(i,j\)\-th entryρ∗​\(i,j\)\\rho^\{\*\}\(i,j\)is a probability distribution on𝒳\\mathcal\{X\}\. The goal of this paper is to recoverρ∗\\rho^\{\*\}from partial and indirect observations\. Specifically, observations are available only for a subset of the entries ofρ∗\\rho^\{\*\}, and even for these observed entries, the underlying distribution is not directly observed\. Instead, for each observed entry\(i,j\)\(i,j\), we collect i\.i\.d\. samplesX1\(i​j\),…,XNi​j\(i​j\)∼ρ∗​\(i,j\)X^\{\(ij\)\}\_\{1\},\\dots,X^\{\(ij\)\}\_\{N\_\{ij\}\}\\sim\\rho^\{\*\}\(i,j\), whereNi​jN\_\{ij\}denotes the sample size for entry\(i,j\)\(i,j\)\. LetT​\(i,j\)∈\{0,1\}T\(i,j\)\\in\\\{0,1\\\}be an indicator of whether the distributionρ∗​\(i,j\)\\rho^\{\*\}\(i,j\)is observed:T​\(i,j\)=1T\(i,j\)=1if samples are collected fromρ∗​\(i,j\)\\rho^\{\*\}\(i,j\), andT​\(i,j\)=0T\(i,j\)=0otherwise\. The observed data thus consist of\{X1\(i​j\),…,XNi​j\(i​j\):T​\(i,j\)=1,i=1,…,m1,j=1,…,m2\}\\\{X^\{\(ij\)\}\_\{1\},\\dots,X^\{\(ij\)\}\_\{N\_\{ij\}\}:T\(i,j\)=1,i=1,\\dots,m\_\{1\},j=1,\\dots,m\_\{2\}\\\}\. Throughout the paper, we assume that the collections\{\(T​\(i,j\),X1\(i​j\),…,XNi​j\(i​j\)\):i=1,…,m1,j=1,…,m2\}\\\{\(T\(i,j\),X^\{\(ij\)\}\_\{1\},\\dots,X^\{\(ij\)\}\_\{N\_\{ij\}\}\):i=1,\\dots,m\_\{1\},j=1,\\dots,m\_\{2\}\\\}are mutually independent\.

There are various common ways to represent a distribution on𝒳\\mathcal\{X\}\. For examples, probability mass/density functions and cumulative distribution function\. However, many of these classical functional representations exhibits manifold constraints\. To avoid manifold constraints, kernel mean embedding \(KME\) has become a popular tool for Hilbertian embedding of probability distributions\. It allows distributions to be treated as elements of a reproducing kernel Hilbert space \(RKHS\)\(e\.g\. Smola et al\.,[2007](https://arxiv.org/html/2606.04176#bib.bib30); Muandet et al\.,[2017](https://arxiv.org/html/2606.04176#bib.bib27)\)\. Moreover, KME supports efficient empirical estimation from samples, and admits a natural metric—the Maximum Mean Discrepancy \(MMD\)—that facilitates distribution comparison and learning in a wide range of machine learning tasks\(Gretton et al\.,[2012](https://arxiv.org/html/2606.04176#bib.bib13)\)\. More specifically, the MMD is defined as

MMD​\(ρ1,ρ2;ℋ\)\\displaystyle\\mathrm\{MMD\}\(\\rho\_\{1\},\\rho\_\{2\};\\mathcal\{H\}\):=sup‖f‖ℋ≤1\(𝔼X∼ρ1​\[f​\(X\)\]−𝔼Y∼ρ2​\[f​\(Y\)\]\)\\displaystyle:=\\sup\_\{\\begin\{subarray\}\{c\}\\\|f\\\|\_\{\\mathcal\{H\}\}\\leq 1\\end\{subarray\}\}\\left\(\\mathbb\{E\}\_\{X\\sim\\rho\_\{1\}\}\[f\(X\)\]\-\\mathbb\{E\}\_\{Y\\sim\\rho\_\{2\}\}\[f\(Y\)\]\\right\)=‖𝔼X∼ρ1​K​\(X,⋅\)−𝔼X∼ρ2​K​\(X,⋅\)‖ℋ,\\displaystyle=\\\|\\mathbb\{E\}\_\{X\\sim\\rho\_\{1\}\}\\,K\(X,\\cdot\)\-\\mathbb\{E\}\_\{X\\sim\\rho\_\{2\}\}\\,K\(X,\\cdot\)\\\|\_\{\\mathcal\{H\}\},whereK:𝒳×𝒳→ℝK:\\mathcal\{X\}\\times\\mathcal\{X\}\\to\\mathbb\{R\}is a reproducing kernel andℋ\\mathcal\{H\}denotes the associated RKHS\. Throughout this work, we represent the entries of the target distribution\-valued matrixρ∗∈𝒫​\(𝒳\)m1×m2\\rho^\{\*\}\\in\\mathcal\{P\}\(\\mathcal\{X\}\)^\{m\_\{1\}\\times m\_\{2\}\}by their corresponding KMEs\. Specifically, for each\(i,j\)\(i,j\), we define

μ∗​\(i,j,⋅\)=𝔼X∼ρ∗​\(i,j\)​K​\(X,⋅\)∈ℋ\.\\mu^\{\*\}\(i,j,\\cdot\)\\;=\\;\\mathbb\{E\}\_\{X\\sim\\rho^\{\*\}\(i,j\)\}\\,K\(X,\\cdot\)\\;\\in\\;\\mathcal\{H\}\.\(1\)When the kernelKKis characteristic \(for example, Gaussian and Laplacian kernels\), the KME is injective, implying that the underlying distributionρ∗​\(i,j\)\\rho^\{\*\}\(i,j\)is uniquely determined byμ∗​\(i,j,⋅\)\\mu^\{\*\}\(i,j,\\cdot\)\(Sriperumbudur et al\.,[2010](https://arxiv.org/html/2606.04176#bib.bib31)\)\. Consequently, recovering the distribution\-valued matrixρ∗\\rho^\{\*\}is equivalent to recovering its KME matrixμ∗\\mu^\{\*\}\. In the following, we therefore focus on the problem of completing the KME matrixμ∗∈ℋm1×m2\.\\mu^\{\*\}\\in\\mathcal\{H\}^\{m\_\{1\}\\times m\_\{2\}\}\.

## 3Low Rank Modeling

Following the classical literature on noisy matrix completion\(e\.g\., Klopp,[2014](https://arxiv.org/html/2606.04176#bib.bib18); Candes and Plan,[2010](https://arxiv.org/html/2606.04176#bib.bib2)\), we aim to impose a suitable “low\-rank” structural assumption onμ∗\\mu^\{\*\}to facilitate recovery from partial observations\. We begin by observing that the KME matrixμ∗\\mu^\{\*\}can be naturally viewed as a three\-way tensor taking values in the tensor product spaceℝm1⊗ℝm2⊗ℋ\.\\mathbb\{R\}^\{m\_\{1\}\}\\otimes\\mathbb\{R\}^\{m\_\{2\}\}\\otimes\\mathcal\{H\}\.Our goal is to impose low\-rank structure onμ∗\\mu^\{\*\}in the Tucker sense\. Note that, unlike standard tensor learning settings, the last spaceℋ\\mathcal\{H\}is often infinite dimensional, resulting in a partly infinite\-dimensional tensor with finite marginal dimensions for the first two modes and an infinite marginal dimension for the last mode\. To deal with this type of tensors, we need more involved mathematical construct for the investigation of multilinear rank\. Motivated by the notion of functional unfolding operators\(Wang et al\.,[2022](https://arxiv.org/html/2606.04176#bib.bib35)\), we define three unfolding operators that mapμ∗\\mu^\{\*\}into tensor products of two spaces, thereby inducing linear operators whose ranks can be characterized\.

Notice that the tensor product spaceℝm1×ℝm2×ℋ\\mathbb\{R\}^\{m\_\{1\}\}\\times\\mathbb\{R\}^\{m\_\{2\}\}\\times\\mathcal\{H\}is the completion of the linear span of all elementary tensors of the forma1⊗a2⊗ha\_\{1\}\\otimes a\_\{2\}\\otimes hunder the inner product⟨a1⊗a2⊗h,a1′⊗a2′⊗h′⟩ℝm1×ℝm2×ℋ=\(a1⊺​a1′\)​\(a2⊺​a2′\)​⟨h,h′⟩ℋ\\langle a\_\{1\}\\otimes a\_\{2\}\\otimes h,a^\{\\prime\}\_\{1\}\\otimes a^\{\\prime\}\_\{2\}\\otimes h^\{\\prime\}\\rangle\_\{\\mathbb\{R\}^\{m\_\{1\}\}\\times\\mathbb\{R\}^\{m\_\{2\}\}\\times\\mathcal\{H\}\}=\(a\_\{1\}^\{\\intercal\}a\_\{1\}^\{\\prime\}\)\(a\_\{2\}^\{\\intercal\}a\_\{2\}^\{\\prime\}\)\\langle h,h^\{\\prime\}\\rangle\_\{\\mathcal\{H\}\}, wherea1,a1′∈ℝm1a\_\{1\},a\_\{1\}^\{\\prime\}\\in\\mathbb\{R\}^\{m\_\{1\}\},a2,a2′∈ℝm2a\_\{2\},a\_\{2\}^\{\\prime\}\\in\\mathbb\{R\}^\{m\_\{2\}\}, andh,h′∈ℋh,h^\{\\prime\}\\in\\mathcal\{H\}\.

###### Definition 3\.1\.

The unfolding operators𝒰j\\mathcal\{U\}\_\{j\},j=1,2,3j=1,2,3, for any elementary tensora1⊗a2⊗ha\_\{1\}\\otimes a\_\{2\}\\otimes h, are defined as follows:

𝒰1\\displaystyle\\mathcal\{U\}\_\{1\}:ℝm1⊗ℝm2⊗ℋ→ℝm1⊗\(ℝm2⊗ℋ\),\\displaystyle:\\mathbb\{R\}^\{m\_\{1\}\}\\otimes\\mathbb\{R\}^\{m\_\{2\}\}\\otimes\\mathcal\{H\}\\to\\mathbb\{R\}^\{m\_\{1\}\}\\otimes\(\\mathbb\{R\}^\{m\_\{2\}\}\\otimes\\mathcal\{H\}\),𝒰1​\(a1⊗a2⊗h\)=a1⊗\(a2⊗h\);\\displaystyle\\mathcal\{U\}\_\{1\}\(a\_\{1\}\\otimes a\_\{2\}\\otimes h\)=a\_\{1\}\\otimes\(a\_\{2\}\\otimes h\);𝒰2\\displaystyle\\mathcal\{U\}\_\{2\}:ℝm1⊗ℝm2⊗ℋ→ℝm2⊗\(ℝm1⊗ℋ\),\\displaystyle:\\mathbb\{R\}^\{m\_\{1\}\}\\otimes\\mathbb\{R\}^\{m\_\{2\}\}\\otimes\\mathcal\{H\}\\to\\mathbb\{R\}^\{m\_\{2\}\}\\otimes\(\\mathbb\{R\}^\{m\_\{1\}\}\\otimes\\mathcal\{H\}\),𝒰2​\(a1⊗a2⊗h\)=a2⊗\(a1⊗h\);\\displaystyle\\mathcal\{U\}\_\{2\}\(a\_\{1\}\\otimes a\_\{2\}\\otimes h\)=a\_\{2\}\\otimes\(a\_\{1\}\\otimes h\);𝒰3\\displaystyle\\mathcal\{U\}\_\{3\}:ℝm1⊗ℝm2⊗ℋ→ℋ⊗\(ℝm1⊗ℝm2\),\\displaystyle:\\mathbb\{R\}^\{m\_\{1\}\}\\otimes\\mathbb\{R\}^\{m\_\{2\}\}\\otimes\\mathcal\{H\}\\to\\mathcal\{H\}\\otimes\(\\mathbb\{R\}^\{m\_\{1\}\}\\otimes\\mathbb\{R\}^\{m\_\{2\}\}\),𝒰3​\(a1⊗a2⊗h\)=h⊗\(a1⊗a2\)\.\\displaystyle\\mathcal\{U\}\_\{3\}\(a\_\{1\}\\otimes a\_\{2\}\\otimes h\)=h\\otimes\(a\_\{1\}\\otimes a\_\{2\}\)\.These definitions extend to any elementμ∈ℝm1⊗ℝm2⊗ℋ\\mu\\in\\mathbb\{R\}^\{m\_\{1\}\}\\otimes\\mathbb\{R\}^\{m\_\{2\}\}\\otimes\{\\mathcal\{H\}\}by linearity\.

For a finite\-dimensional tensor, the unfolding operators defined above are equivalent to the standard matricization of the tensor defined in Definition[A\.1](https://arxiv.org/html/2606.04176#A1.Thmtheorem1)in Appendix[A](https://arxiv.org/html/2606.04176#A1)\. Any elementf∈ℱ1⊗ℱ2f\\in\\mathcal\{F\}\_\{1\}\\otimes\\mathcal\{F\}\_\{2\}can be naturally identified with a linear operator mappingℱ2\\mathcal\{F\}\_\{2\}toℱ1\\mathcal\{F\}\_\{1\}, defined byg∈ℱ2⟼⟨f​\(⋆,⋅\),g​\(⋅\)⟩ℱ2\.g\\in\\mathcal\{F\}\_\{2\}\\;\\longmapsto\\;\\left\\langle f\(\\star,\\cdot\),\\,g\(\\cdot\)\\right\\rangle\_\{\\mathcal\{F\}\_\{2\}\}\.With a slight abuse of notation, we useffto denote both the tensor element and its induced linear operator\. Under this operator interpretation, an element𝒰1​μ∈ℝm1⊗\(ℝm2⊗ℋ\)\\mathcal\{U\}\_\{1\}\\mu\\in\\mathbb\{R\}^\{m\_\{1\}\}\\otimes\\bigl\(\\mathbb\{R\}^\{m\_\{2\}\}\\otimes\\mathcal\{H\}\\bigr\)can be viewed as a linear operator mapping the spaceℝm2⊗ℋ\\mathbb\{R\}^\{m\_\{2\}\}\\otimes\\mathcal\{H\}intoℝm1\\mathbb\{R\}^\{m\_\{1\}\}\. Similar operator interpretations apply to𝒰2​μ\\mathcal\{U\}\_\{2\}\\muand𝒰3​μ\\mathcal\{U\}\_\{3\}\\mu\. And we define the multilinear rank ofμ∈ℝm1⊗ℝm2⊗ℋ\\mu\\in\\mathbb\{R\}^\{m\_\{1\}\}\\otimes\\mathbb\{R\}^\{m\_\{2\}\}\\otimes\{\\mathcal\{H\}\}as

\(rank​\(𝒰1​μ\),rank​\(𝒰2​μ\),rank​\(𝒰3​μ\)\),\\bigl\(\\mathrm\{rank\}\(\\mathcal\{U\}\_\{1\}\\mu\),\\mathrm\{rank\}\(\\mathcal\{U\}\_\{2\}\\mu\),\\mathrm\{rank\}\(\\mathcal\{U\}\_\{3\}\\mu\)\\bigr\),where the rank is understood as the rank of the corresponding induced linear operator\. When the Hilbert spaceℋ\\mathcal\{H\}is finite\-dimensional, this definition coincides with the classical Tucker rank of a finite\-dimensional third\-order tensor\(Tucker,[1966](https://arxiv.org/html/2606.04176#bib.bib32); Kolda and Bader,[2009](https://arxiv.org/html/2606.04176#bib.bib19)\)\.

Next, we illustrate how a small multilinear rank ofμ∗\\mu^\{\\ast\}leads to meaningful structural understanding\. Assume thatℋ\\mathcal\{H\}is a separable Hilbert space\. Let\{ek\}k=1q\\\{e\_\{k\}\\\}\_\{k=1\}^\{q\}denote an orthonormal basis of the Hilbert spaceℋ\\mathcal\{H\}, whereqqcan be infinite\. Then, for each\(i,j\)\(i,j\), the KME admits the expansion

μ∗​\(i,j,⋅\)=∑k=1qB​\(i,j,k\)​ek​\(⋅\),\\displaystyle\\mu^\{\*\}\(i,j,\\cdot\)\\;=\\;\\sum\_\{k=1\}^\{q\}B\(i,j,k\)\\,e\_\{k\}\(\\cdot\),\(2\)whereB​\(i,j,k\)∈ℝB\(i,j,k\)\\in\\mathbb\{R\},k=1,…,qk=1,\\dots,q, are the corresponding expansion coefficients\.

Suppose thatμ∗\\mu^\{\*\}has multilinear rank\(r1,r2,r3\)\(r\_\{1\},r\_\{2\},r\_\{3\}\), withr1<m1r\_\{1\}<m\_\{1\},r2<m2r\_\{2\}<m\_\{2\}, andr3<qr\_\{3\}<q\. Then there exist a tensorG∈ℝr1×r2×r3G\\in\\mathbb\{R\}^\{r\_\{1\}\\times r\_\{2\}\\times r\_\{3\}\}and factor matricesU1∈ℝm1×r1U\_\{1\}\\in\\mathbb\{R\}^\{m\_\{1\}\\times r\_\{1\}\},U2∈ℝm2×r2U\_\{2\}\\in\\mathbb\{R\}^\{m\_\{2\}\\times r\_\{2\}\}, andU3∈ℝq×r3U\_\{3\}\\in\\mathbb\{R\}^\{q\\times r\_\{3\}\}, such that the coefficient tensorBBadmits the Tucker decomposition\(Tucker,[1966](https://arxiv.org/html/2606.04176#bib.bib32)\)

B=G×1U1×2U2×3U3,\\displaystyle B\\;=\\;G\\times\_\{1\}U\_\{1\}\\times\_\{2\}U\_\{2\}\\times\_\{3\}U\_\{3\},\(3\)where×j\\times\_\{j\}denotes the mode\-jjtensor–matrix product, defined in Definition[A\.2](https://arxiv.org/html/2606.04176#A1.Thmtheorem2)in Appendix[A](https://arxiv.org/html/2606.04176#A1)\. Here,GGis referred to as the core tensor, and the matricesUjU\_\{j\},j=1,2,3j=1,2,3, are often assumed to have orthonormal columns\. An illustration of the Tucker decomposition is provided in Figure[3](https://arxiv.org/html/2606.04176#A1.F3)in Appendix[A](https://arxiv.org/html/2606.04176#A1)\.

With such decomposition \([3](https://arxiv.org/html/2606.04176#S3.E3)\), we have

μ∗​\(i,j,⋅\)\\displaystyle\\mu^\{\*\}\(i,j,\\cdot\)=∑l1=1r1∑l2=1r2∑l3=1r3G​\(l1,l2,l3\)​U1​\(i,l1\)​U2​\(j,l2\)​sl3​\(⋅\),\\displaystyle=\\sum\_\{l\_\{1\}=1\}^\{r\_\{1\}\}\\sum\_\{l\_\{2\}=1\}^\{r\_\{2\}\}\\sum\_\{l\_\{3\}=1\}^\{r\_\{3\}\}G\(l\_\{1\},l\_\{2\},l\_\{3\}\)U\_\{1\}\(i,l\_\{1\}\)U\_\{2\}\(j,l\_\{2\}\)s\_\{l\_\{3\}\}\(\\cdot\),\(4\)wheresl3​\(⋅\)=∑k=1qU3​\(k,l3\)​ek​\(⋅\)s\_\{l\_\{3\}\}\(\\cdot\)=\\sum\_\{k=1\}^\{q\}U\_\{3\}\(k,l\_\{3\}\)e\_\{k\}\(\\cdot\),k=1,…,r3k=1,\\dots,r\_\{3\}\. This shows that onlyr3r\_\{3\}basis functions\{sl3\}l3=1r3\\\{s\_\{l\_\{3\}\}\\\}\_\{l\_\{3\}=1\}^\{r\_\{3\}\}are sufficient to characterize the functional component ofμ∗\\mu^\{\*\}\.

By the property of KMEs, for anyf∈ℋf\\in\\mathcal\{H\},𝔼X∼ρ​\(i,j\)​\[f​\(X\)\]=⟨f,μ∗​\(i,j,⋅\)⟩ℋ\.\\mathbb\{E\}\_\{X\\sim\\rho\(i,j\)\}\[f\(X\)\]=\\langle f,\\mu^\{\*\}\(i,j,\\cdot\)\\rangle\_\{\\mathcal\{H\}\}\.Define the matrixA∈ℝm1×m2A\\in\\mathbb\{R\}^\{m\_\{1\}\\times m\_\{2\}\}by

Af​\(i,j\):=𝔼X∼ρ​\(i,j\)​\[f​\(X\)\],f∈ℋ\.\\displaystyle A\_\{f\}\(i,j\):=\\mathbb\{E\}\_\{X\\sim\\rho\(i,j\)\}\[f\(X\)\],\\qquad f\\in\{\\mathcal\{H\}\}\.\(5\)For example, iff​\(t\)=tf\(t\)=t,t∈𝒳t\\in\{\\mathcal\{X\}\}\(i\.e\., an identity function\), thenAAis the mean matrix\. Substituting \([3](https://arxiv.org/html/2606.04176#S3.Ex11)\) into the above expression yields

Af​\(i,j\)\\displaystyle A\_\{f\}\(i,j\)=∑l1=1r1∑l2=1r2∑l3=1r3G​\(l1,l2,l3\)​U1​\(i,l1\)​U2​\(j,l2\)​⟨f,sl3⟩ℋ\\displaystyle=\\sum\_\{l\_\{1\}=1\}^\{r\_\{1\}\}\\sum\_\{l\_\{2\}=1\}^\{r\_\{2\}\}\\sum\_\{l\_\{3\}=1\}^\{r\_\{3\}\}G\(l\_\{1\},l\_\{2\},l\_\{3\}\)U\_\{1\}\(i,l\_\{1\}\)U\_\{2\}\(j,l\_\{2\}\)\\langle f,s\_\{l\_\{3\}\}\\rangle\_\{\\mathcal\{H\}\}=∑l1=1r1∑l2=1r2U1​\(i,l1\)​U2​\(j,l2\)​Cf​\(l1,l2\),\\displaystyle=\\sum\_\{l\_\{1\}=1\}^\{r\_\{1\}\}\\sum\_\{l\_\{2\}=1\}^\{r\_\{2\}\}U\_\{1\}\(i,l\_\{1\}\)U\_\{2\}\(j,l\_\{2\}\)C\_\{f\}\(l\_\{1\},l\_\{2\}\),\(6\)whereCf​\(l1,l2\):=∑l3=1r3G​\(l1,l2,l3\)​⟨f,sl3⟩ℋ\.C\_\{f\}\(l\_\{1\},l\_\{2\}\):=\\sum\_\{l\_\{3\}=1\}^\{r\_\{3\}\}G\(l\_\{1\},l\_\{2\},l\_\{3\}\)\\langle f,s\_\{l\_\{3\}\}\\rangle\_\{\\mathcal\{H\}\}\.Consequently, the matrixAfA\_\{f\}admits the factorizationAf=U1​Cf​U2⊺,Cf∈ℝr1×r2,A\_\{f\}=U\_\{1\}C\_\{f\}U\_\{2\}^\{\\intercal\},C\_\{f\}\\in\\mathbb\{R\}^\{r\_\{1\}\\times r\_\{2\}\},and therefore has rank at mostmin⁡\{r1,r2\}\\min\\\{r\_\{1\},r\_\{2\}\\\}\. By imposing a small multilinear rank onμ∗\\mu^\{\*\}, we promote intrinsic dimension reduction in the functional component of the KME, thereby facilitating downstream analysis of the associated distributions, and at the same time induce a low\-rank structure for the expectation functional matrixAfA\_\{f\}\.

## 4Estimation

We adopt a regularizedMM\-estimation approach to construct the proposed estimatorμ^\\hat\{\\mu\}\. Specifically,

μ^=argminμ∈ℝm1⊗ℝm2⊗ℋ:supi,j‖μ​\(i,j,⋅\)‖ℋ≤α\\displaystyle\\hat\{\\mu\}=\\operatorname\*\{argmin\}\_\{\\mu\\in\\mathbb\{R\}^\{m\_\{1\}\}\\otimes\\mathbb\{R\}^\{m\_\{2\}\}\\otimes\{\\mathcal\{H\}\}\\;:\\;\\sup\_\{i,j\}\\\|\\mu\(i,j,\\cdot\)\\\|\_\{\\mathcal\{H\}\}\\leq\\alpha\}\[1m1​m2∑\(i,j\):T​\(i,j\)=1∥μ\(i,j,⋅\)−1Ni​j∑l=1Ni​jK\(Xl\(i​j\),⋅\)∥ℋ2\\displaystyle\\Bigg\[\\frac\{1\}\{m\_\{1\}m\_\{2\}\}\\sum\_\{\(i,j\):T\(i,j\)=1\}\\Big\\\|\\mu\(i,j,\\cdot\)\-\\frac\{1\}\{N\_\{ij\}\}\\sum\_\{l=1\}^\{N\_\{ij\}\}K\(X^\{\(ij\)\}\_\{l\},\\cdot\)\\Big\\\|^\{2\}\_\{\\mathcal\{H\}\}\(7\)\+λΨ\(μ\)\],\\displaystyle\\qquad\+\\lambda\\Psi\(\\mu\)\\Bigg\],\(8\)where the term in \([7](https://arxiv.org/html/2606.04176#S4.E7)\) involves the squared maximum mean discrepancies \(MMD\) between the estimated KME and the empirical KME at the observed entries\. The functionalΨ​\(⋅\)\\Psi\(\\cdot\)is a regularization term designed to promote low\-rank structure inμ\\muwhile simultaneously ensuring the validity of the representer theorem \(see Theorem[4\.1](https://arxiv.org/html/2606.04176#S4.Thmtheorem1)\)\. Here,λ\>0\\lambda\>0is a tuning parameter controlling the strength of regularization, andα\\alphais an upper bound on the RKHS norm of each kernel mean embeddingμ​\(i,j,⋅\)\\mu\(i,j,\\cdot\)\. In the implementation, we simply takeα\\alphato be large value \(e\.g\.,maxi,j⁡10​‖∑l=1Ni​jK​\(Xl\(i​j\),⋅\)/Ni​j‖ℋ\\max\_\{i,j\}10\\\|\\sum\_\{l=1\}^\{N\_\{ij\}\}K\(X^\{\(ij\)\}\_\{l\},\\cdot\)/\{N\_\{ij\}\}\\\|\_\{\\mathcal\{H\}\}\)\.

For the regularization termΨ​\(μ\)\\Psi\(\\mu\), we adopt a nuclear\-norm\-based penalty as a convex surrogate to promote low\-rank structure\. Specifically, we define

Ψ​\(μ\)=∑j=13βj​‖𝒰j​μ‖∗,\\displaystyle\\Psi\(\\mu\)=\\sum\_\{j=1\}^\{3\}\\beta\_\{j\}\\\|\{\\mathcal\{U\}\}\_\{j\}\\mu\\\|\_\{\*\},\(9\)where∥⋅∥∗\\\|\\cdot\\\|\_\{\*\}denotes the nuclear norm of the induced operator corresponding to each unfolding\. The weightsβj\>0\\beta\_\{j\}\>0\(such that∑j=13βj=1\\sum\_\{j=1\}^\{3\}\\beta\_\{j\}=1\) control the relative degree of low\-rank regularization imposed on the different unfoldings ofμ\\mu\. By default, we setβj=1/3\\beta\_\{j\}=1/3forj=1,2,3j=1,2,3\.

Sinceμ^∈ℝm1⊗ℝm2⊗ℋ\\hat\{\\mu\}\\in\\mathbb\{R\}^\{m\_\{1\}\}\\otimes\\mathbb\{R\}^\{m\_\{2\}\}\\otimes\\mathcal\{H\}, whereℋ\\mathcal\{H\}may be infinite\-dimensional, the optimization problem in \([8](https://arxiv.org/html/2606.04176#S4.E8)\) appears challenging due to its infinite\-dimensional nature\. However, by the following representer theorem, one can verify that the solution to \([8](https://arxiv.org/html/2606.04176#S4.E8)\) admits a finite\-dimensional representation\.

###### Theorem 4\.1\(Representer Theorem\)\.

For any\(i,j\)\(i,j\),μ^​\(i,j,…\)\\hat\{\\mu\}\(i,j,\\dots\)lies in the following finite\-dimensional space

𝒦=span\{K\(Xl\(i​j\),⋅\),l=1,…,Ni​j,Ti​j=1\}\.\\displaystyle\{\\mathcal\{K\}\}=\\mathrm\{span\}\\\{K\(X\_\{l\}^\{\(ij\)\},\\cdot\),l=1,\\dots,N\_\{ij\},T\_\{ij\}=1\\\}\.\(10\)

Specifically, the target optimization \([8](https://arxiv.org/html/2606.04176#S4.E8)\) can be recast as a finite\-dimensional optimization over parameterA∈ℝm1×m2×N~A\\in\\mathbb\{R\}^\{m\_\{1\}\\times m\_\{2\}\\times\\tilde\{N\}\}such thatμ​\(i,j,⋅\)=∑l=1N~A​\(i,j,l\)​K​\(𝒙l,⋅\)\{\\mu\}\(i,j,\\cdot\)=\\sum\_\{l=1\}^\{\\tilde\{N\}\}A\(i,j,l\)K\(\\bm\{x\}\_\{l\},\\cdot\), whereN~=∑i​jTi​j​Ni​j\\tilde\{N\}=\\sum\_\{ij\}T\_\{ij\}N\_\{ij\}and𝒙=\[x\(i,j\)\]\(i,j\):Ti,j=1⊺∈ℝN~\\bm\{x\}=\[x^\{\(i,j\)\}\]\_\{\(i,j\):T\_\{i,j\}=1\}^\{\\intercal\}\\in\\mathbb\{R\}^\{\\tilde\{N\}\}withx\(i,j\)=\[X1\(i,j\),…,XNi​j\(i,j\)\]⊺∈ℝNi​jx^\{\(i,j\)\}=\[X^\{\(i,j\)\}\_\{1\},\\dots,X^\{\(i,j\)\}\_\{N\_\{ij\}\}\]^\{\\intercal\}\\in\\mathbb\{R\}^\{N\_\{ij\}\}\. It is easy to verify that

‖μ​\(i,j,⋅\)−1Ni​j​∑l=1Ni​jK​\(Xl\(i​j\),⋅\)‖ℋ2=∑l,l′=1,…,N~A​\(i,j,l\)​A​\(i,j,l′\)​K​\(𝒙l,𝒙l′\)−2​∑lN~∑l′=1Ni​jA​\(i,j,l\)/Ni​j​K​\(𝒙l,Xl′\(i​j\)\)\+1Ni​j2​∑l,l′=1,…,Ni​jK​\(𝑿l\(i​j\),Xl′\(i​j\)\)\.\\Big\\\|\\mu\(i,j,\\cdot\)\-\\frac\{1\}\{N\_\{ij\}\}\\sum\_\{l=1\}^\{N\_\{ij\}\}K\(X^\{\(ij\)\}\_\{l\},\\cdot\)\\Big\\\|^\{2\}\_\{\\mathcal\{H\}\}\\\\ =\\sum\_\{l,l^\{\\prime\}=1,\\dots,\\tilde\{N\}\}A\(i,j,l\)A\(i,j,l^\{\\prime\}\)K\(\\bm\{x\}\_\{l\},\\bm\{x\}\_\{l^\{\\prime\}\}\)\-2\\sum\_\{l\}^\{\\tilde\{N\}\}\\sum\_\{l^\{\\prime\}=1\}^\{N\_\{ij\}\}A\(i,j,l\)/N\_\{ij\}K\(\\bm\{x\}\_\{l\},X\_\{l^\{\\prime\}\}^\{\(ij\)\}\)\\\\ \+\\frac\{1\}\{N^\{2\}\_\{ij\}\}\\sum\_\{l,l^\{\\prime\}=1,\\dots,N\_\{ij\}\}K\(\\bm\{X\}\_\{l\}^\{\(ij\)\},X\_\{l^\{\\prime\}\}^\{\(ij\)\}\)\.Re\-expressing this term in vector\-matrix notation, one can show that the optimization \([7](https://arxiv.org/html/2606.04176#S4.E7)\) involves the Gram matrixKG:=\[K​\(𝒙l,𝒙m\)\]l,m=1,…,N~∈ℝN~×N~K\_\{G\}:=\[K\(\\bm\{x\}\_\{l\},\\bm\{x\}\_\{m\}\)\]\_\{l,m=1,\\dots,\\tilde\{N\}\}\\in\\mathbb\{R\}^\{\\tilde\{N\}\\times\\tilde\{N\}\}\.

In practice, however, the total number of observations∑i,jTi​j​Ni​j\\sum\_\{i,j\}T\_\{ij\}N\_\{ij\}can be very large, rendering the optimization computationally expensive\. To address this issue, one may employ standard kernel approximation techniques such as the Nyström method\(Williams and Seeger,[2001](https://arxiv.org/html/2606.04176#bib.bib37); Drineas and Mahoney,[2005](https://arxiv.org/html/2606.04176#bib.bib9)\)or Random Fourier Features\(Rahimi and Recht,[2007](https://arxiv.org/html/2606.04176#bib.bib29)\)to obtain scalable approximations of the Gram matrixKGK\_\{G\}\. More specifically, we approximateKGK\_\{G\}byZ​Z⊺ZZ^\{\\intercal\}, whereZ∈ℝ\(∑i,jTi​j​Ni​j\)×qZ\\in\\mathbb\{R\}^\{\(\\sum\_\{i,j\}T\_\{ij\}N\_\{ij\}\)\\times q\}, withq≪∑i,jTi​j​Ni​jq\\ll\\sum\_\{i,j\}T\_\{ij\}N\_\{ij\}\.

Next, we rewrite the optimization into the form of \(finite\-dimensional\) tensor optimization problem explicitly\. To this end, we construct the tensorY∈ℝm1×m2×qY\\in\\mathbb\{R\}^\{m\_\{1\}\\times m\_\{2\}\\times q\}as follows\. For each observed entry\(i,j\)\(i,j\)withTi​j=1T\_\{ij\}=1, letY​\(i,j,⋅\)∈ℝqY\(i,j,\\cdot\)\\in\\mathbb\{R\}^\{q\}be the average of the rows ofZZcorresponding to the sample vectorx\(i,j\)x^\{\(i,j\)\}; for unobserved entries withTi​j=0T\_\{ij\}=0, setY​\(i,j,⋅\)=𝟎Y\(i,j,\\cdot\)=\\bm\{0\}\. With this construction, the optimization problem in \([8](https://arxiv.org/html/2606.04176#S4.E8)\) admits the following re\-parameterization:

B^=argminB∈ℝm1×m2×q,supi,j‖B​\(i,j,⋅\)‖2≤α\[∑i​jTi​j​‖B​\(i,j,⋅\)−Y​\(i,j,⋅\)‖22\+λ​\(∑j=13βj​‖𝒰j​B‖∗\)\]\\displaystyle\\hat\{B\}=\\operatorname\*\{argmin\}\_\{B\\in\\mathbb\{R\}^\{m\_\{1\}\\times m\_\{2\}\\times q\},\\sup\_\{i,j\}\\\|B\(i,j,\\cdot\)\\\|\_\{2\}\\leq\\alpha\}\\left\[\\sum\_\{ij\}T\_\{ij\}\\\|B\(i,j,\\cdot\)\-Y\(i,j,\\cdot\)\\\|\_\{2\}^\{2\}\+\\lambda\\left\(\\sum\_\{j=1\}^\{3\}\\beta\_\{j\}\\\|\{\\mathcal\{U\}\}\_\{j\}B\\\|\_\{\*\}\\right\)\\right\]\(11\)whereB​\(i,j,⋅\),Y​\(i,j,⋅\)∈ℝqB\(i,j,\\cdot\),Y\(i,j,\\cdot\)\\in\\mathbb\{R\}^\{q\}and∥⋅∥2\\\|\\cdot\\\|\_\{2\}denotes the Euclidean norm\. WithB^\\hat\{B\}, the estimatorμ^\\hat\{\\mu\}is obtained by

μ^​\(i,j,t\)=zt⊺​B^​\(i,j,⋅\),\\hat\{\\mu\}\(i,j,t\)\\;=\\;z\_\{t\}^\{\\intercal\}\\,\\hat\{B\}\(i,j,\\cdot\),wherezt∈ℝqz\_\{t\}\\in\\mathbb\{R\}^\{q\}is the feature vector associated with the inputt∈𝒳t\\in\\mathcal\{X\}\. More specifically,ztz\_\{t\}can be constructed aszt=\[K​\(xℓ,t\)\]xℓ∈𝒙⊺​Z\+,z\_\{t\}\\;=\\;\\bigl\[\\,K\(x\_\{\\ell\},t\)\\,\\bigr\]\_\{x\_\{\\ell\}\\in\\bm\{x\}\}^\{\\intercal\}\\,Z^\{\+\},whereZ\+Z^\{\+\}denotes the Moore–Penrose pseudoinverse ofZZ\.

We solve \([11](https://arxiv.org/html/2606.04176#S4.E11)\) using an accelerated alternating direction method of multipliers \(ADMM\) algorithm\(Kadkhodaie and Montanari,[2015](https://arxiv.org/html/2606.04176#bib.bib16)\)\. The details of the algorithm can be found in Appendix[B\.1](https://arxiv.org/html/2606.04176#A2.SS1)\.

## 5Theoretical Gaurantees

In this section, we establish theoretical guarantees for the estimatorμ^\\hat\{\\mu\}\. Our analysis quantifies estimation error through the aggregate squared maximum mean discrepancy \(MMD\) across all matrix entries, both observed and unobserved:

∑i,j‖μ^​\(i,j,⋅\)−μ∗​\(i,j,⋅\)‖ℋ2\.\\sum\_\{i,j\}\\bigl\\\|\\hat\{\\mu\}\(i,j,\\cdot\)\-\\mu^\{\*\}\(i,j,\\cdot\)\\bigr\\\|\_\{\\mathcal\{H\}\}^\{2\}\.We begin by introducing notation that will be used throughout the section\. Definem=min⁡\{m1,m2\}m=\\min\\\{m\_\{1\},m\_\{2\}\\\},M=max⁡\{m1,m2\}M=\\max\\\{m\_\{1\},m\_\{2\}\\\}andd=m1\+m2d=m\_\{1\}\+m\_\{2\}\. We write≲\\lesssim\(≳\\gtrsim\) to denote inequalities holding up to an absolute positive multiplicative constant, and writea≍ba\\asymp bwhen botha≲ba\\lesssim bandb≳ab\\gtrsim ahold\. Finally, letN=mini,j:Ti​j=1⁡Ni​jN=\\min\_\{i,j:T\_\{ij\}=1\}N\_\{ij\}denote the minimum sample size among allobservedentries\. Next, we impose following technical assumptions\.

###### Assumption 5\.1\.

There exists a constantCK\>0C\_\{K\}\>0such thatsupx∈𝒳\|K​\(x,x\)\|≤CK\\sup\_\{x\\in\{\\mathcal\{X\}\}\}\|K\(x,x\)\|\\leq C\_\{K\}\.

###### Assumption 5\.2\.

The observation indicators\{Ti​j\}i=1,…,m1,j=1,…,m2\\\{T\_\{ij\}\\\}\_\{i=1,\\dots,m\_\{1\},j=1,\\dots,m\_\{2\}\}are independent Bernoulli random variables withπi​j=Pr⁡\(Ti​j=1\)\>0\\pi\_\{ij\}=\\Pr\(T\_\{ij\}=1\)\>0\.\{Xl\(i,j\)\}l=1,…,Ni​j\\\{X^\{\(i,j\)\}\_\{l\}\\\}\_\{l=1,\\dots,N\_\{ij\}\}are independent samples generated from the distributionρ​\(i,j\)\\rho\(i,j\)and they are independent of\{Ti​j\}\\\{T\_\{ij\}\\\}\. Also,\{Xl\(i,j\)\}l=1,…,Ni​j\\\{X^\{\(i,j\)\}\_\{l\}\\\}\_\{l=1,\\dots,N\_\{ij\}\}are independent across different\(i,j\)\(i,j\)\.

DefineπL=mini,j⁡πi​j\\pi\_\{L\}=\\min\_\{i,j\}\\pi\_\{ij\}andπU=maxi,j⁡πi​j\\pi\_\{U\}=\\max\_\{i,j\}\\pi\_\{ij\}\.

###### Assumption 5\.3\.

There exists a constantξ\>0\\xi\>0such thatπU/πL≤ξ\\pi\_\{U\}/\\pi\_\{L\}\\leq\\xi\. In addition,πL−1=𝒪​\(m/log⁡d\)\\pi\_\{L\}^\{\-1\}=\\scalebox\{0\.7\}\{$\\mathcal\{O\}$\}\(m/\\log d\)\.

###### Assumption 5\.4\.

supi,j‖μ∗​\(i,j,⋅\)‖ℋ≤α\\sup\_\{i,j\}\\\|\\mu^\{\*\}\(i,j,\\cdot\)\\\|\_\{\\mathcal\{H\}\}\\leq\\alpha\.

Assumption[5\.1](https://arxiv.org/html/2606.04176#S5.Thmtheorem1)imposes a mild regularity condition on the kernel function and is introduced to control the convergence of the empirical KME\. Assumption[5\.2](https://arxiv.org/html/2606.04176#S5.Thmtheorem2)specifies the data\-generating process\. The corresponding missing structure allows non\-uniform missingness and is standard in the matrix completion literature\(e\.g\. Klopp,[2014](https://arxiv.org/html/2606.04176#bib.bib18); Negahban and Wainwright,[2012](https://arxiv.org/html/2606.04176#bib.bib28); Wang et al\.,[2021](https://arxiv.org/html/2606.04176#bib.bib34); Zhao et al\.,[2025](https://arxiv.org/html/2606.04176#bib.bib42)\)\. Assumption[5\.3](https://arxiv.org/html/2606.04176#S5.Thmtheorem3)is a commonly adopted condition for nuclear\-norm\-based matrix completion estimators\(e\.g\. Klopp,[2014](https://arxiv.org/html/2606.04176#bib.bib18); Negahban and Wainwright,[2012](https://arxiv.org/html/2606.04176#bib.bib28); Koltchinskii et al\.,[2011](https://arxiv.org/html/2606.04176#bib.bib21); Wang et al\.,[2021](https://arxiv.org/html/2606.04176#bib.bib34); Li et al\.,[2024](https://arxiv.org/html/2606.04176#bib.bib23)\)\. The first part of this assumption rules out scenarios in which the observation probabilities differ in their asymptotic orders\. The second part is a mild technical requirement to help establish the consistency of the estimator\. Finally, Assumption[5\.4](https://arxiv.org/html/2606.04176#S5.Thmtheorem4)is a realization assumption, so that the ground truth satisfies the constraint in our optimization \([8](https://arxiv.org/html/2606.04176#S4.E8)\)\.

First, we derive Theorem[5\.5](https://arxiv.org/html/2606.04176#S5.Thmtheorem5)for the case whenμ∗\\mu^\{\*\}has approximately low multilinear rank, i\.e\.,‖𝒰j​μ∗‖∗≲m1​m2\\\|\{\\mathcal\{U\}\}\_\{j\}\\mu^\{\*\}\\\|\_\{\*\}\\lesssim\\sqrt\{m\_\{1\}m\_\{2\}\}\.

###### Theorem 5\.5\.

Under Assumptions[5\.1](https://arxiv.org/html/2606.04176#S5.Thmtheorem1)\-[5\.4](https://arxiv.org/html/2606.04176#S5.Thmtheorem4), if we take

\(λ​β1\+λ​β2\)−1=\\displaystyle\(\\lambda\\beta\_\{1\}\+\\lambda\\beta\_\{2\}\)^\{\-1\}=𝒪\(\[M​πU​log⁡dm1​m2max\{ΔNlog⁡d,α\}\]−1\),\\displaystyle\\qquad\\mathop\{\}\\mathopen\{\}\\mathcal\{O\}\\mathopen\{\}\\left\(\\left\[\\frac\{\\sqrt\{M\\pi\_\{U\}\\log d\}\}\{m\_\{1\}m\_\{2\}\}\\max\\left\\\{\\Delta\_\{N\}\\sqrt\{\\log d\},\\alpha\\right\\\}\\right\]^\{\-1\}\\right\),withΔN=CK/N\+log⁡d/N\\Delta\_\{N\}=C\_\{K\}/N\+\\sqrt\{\\log d/N\}, then, with probability at least1−κ/d1\-\\kappa/dfor some constantκ\>0\\kappa\>0, we have

1m1​m2​∑i,j‖μ^​\(i,j,⋅\)−μ∗​\(i,j,⋅\)‖ℋ2≲\(λ/πL\)​\(β1​‖𝒰1​\(μ∗\)‖∗\+β2​‖𝒰2​\(μ∗\)‖∗\+β3​‖𝒰3​\(μ∗\)‖∗\)\+\(2​α\)2​πU​log⁡dπL2​m1​m2\.\\frac\{1\}\{m\_\{1\}m\_\{2\}\}\\sum\_\{i,j\}\\\|\\hat\{\\mu\}\(i,j,\\cdot\)\-\{\\mu\}^\{\*\}\(i,j,\\cdot\)\\\|\_\{\\mathcal\{H\}\}^\{2\}\\lesssim\\\\ \(\{\\lambda\}/\{\\pi\_\{L\}\}\)\\left\(\\beta\_\{1\}\\\|\{\\mathcal\{U\}\}\_\{1\}\(\\mu^\{\*\}\)\\\|\_\{\*\}\+\\beta\_\{2\}\\\|\{\\mathcal\{U\}\}\_\{2\}\(\\mu^\{\*\}\)\\\|\_\{\*\}\+\\beta\_\{3\}\\\|\{\\mathcal\{U\}\}\_\{3\}\(\\mu^\{\*\}\)\\\|\_\{\*\}\\right\)\+\(2\\alpha\)^\{2\}\\sqrt\{\\frac\{\\pi\_\{U\}\\log d\}\{\\pi\_\{L\}^\{2\}m\_\{1\}m\_\{2\}\}\}\.\\\\More specifically, if we takeλ​β1,λ​β2≍M​πU​log⁡dm1​m2​max⁡\{ΔN​log⁡d,α\}\\lambda\\beta\_\{1\},\\lambda\\beta\_\{2\}\\asymp\\frac\{\\sqrt\{M\\pi\_\{U\}\\log d\}\}\{m\_\{1\}m\_\{2\}\}\\max\\left\\\{\\Delta\_\{N\}\\sqrt\{\\log d\},\\alpha\\right\\\}, then, with probability at least1−κ/d1\-\\kappa/d, we have

1m1​m2​∑i,j‖μ^​\(i,j,⋅\)−μ∗​\(i,j,⋅\)‖ℋ2≲max⁡\{α,ΔN​log⁡d\}​πU​log⁡dπL2​m​\{‖𝒰1​\(μ∗\)‖∗m1​m2\+‖𝒰2​\(μ∗\)‖∗m1​m2\}\+\(2​α\)2​πU​log⁡dπL2​m1​m2\+λ​β3πL​‖𝒰3​\(μ∗\)‖∗\.\\frac\{1\}\{m\_\{1\}m\_\{2\}\}\\sum\_\{i,j\}\\\|\\hat\{\\mu\}\(i,j,\\cdot\)\-\{\\mu\}^\{\*\}\(i,j,\\cdot\)\\\|\_\{\\mathcal\{H\}\}^\{2\}\\lesssim\\\\ \\max\\\{\\alpha,\\Delta\_\{N\}\\sqrt\{\\log d\}\\\}\\sqrt\{\\frac\{\\pi\_\{U\}\\log d\}\{\\pi\_\{L\}^\{2\}m\}\}\\left\\\{\\frac\{\\\|\{\\mathcal\{U\}\}\_\{1\}\(\\mu^\{\*\}\)\\\|\_\{\*\}\}\{\\sqrt\{m\_\{1\}m\_\{2\}\}\}\+\\right\.\\left\.\\frac\{\\\|\{\\mathcal\{U\}\}\_\{2\}\(\\mu^\{\*\}\)\\\|\_\{\*\}\}\{\\sqrt\{m\_\{1\}m\_\{2\}\}\}\\right\\\}\+\(2\\alpha\)^\{2\}\\sqrt\{\\frac\{\\pi\_\{U\}\\log d\}\{\\pi\_\{L\}^\{2\}m\_\{1\}m\_\{2\}\}\}\+\\frac\{\\lambda\\beta\_\{3\}\}\{\\pi\_\{L\}\}\\\|\{\\mathcal\{U\}\}\_\{3\}\(\\mu^\{\*\}\)\\\|\_\{\*\}\.

By choosingλ​β3\\lambda\\beta\_\{3\}to be sufficiently small, the last term in \([28](https://arxiv.org/html/2606.04176#A3.E28)\) becomes negligible\. On the other hand, assigning a larger weight toβ3\\beta\_\{3\}promotes a reduction in the number of basis functions required to represent the KMEs, which is often beneficial in terms of downstream analysis and interpretability\. When the first term in \([28](https://arxiv.org/html/2606.04176#A3.E28)\) dominates, the resulting error bound matches the rate established for classical nuclear\-norm–regularized matrix completion of approximately low\-rank matrices\(Negahban and Wainwright,[2012](https://arxiv.org/html/2606.04176#bib.bib28)\), despite the more challenging setting considered here, where each matrix entry corresponds to a probability distribution rather than a scalar\. In our technical analysis, we develop novel lemmas to tackle these unfolding operators\. See Lemma[C\.1](https://arxiv.org/html/2606.04176#A3.Thmtheorem1)and[C\.2](https://arxiv.org/html/2606.04176#A3.Thmtheorem2)for details\.

Recall that, to our knowledge,Feitelberg et al\. \([2024](https://arxiv.org/html/2606.04176#bib.bib10)\)is the only alternative method for distributional matrix completion\. Compared to the theoretical results inFeitelberg et al\. \([2024](https://arxiv.org/html/2606.04176#bib.bib10)\), our analysis differs in two important aspects\. First, the bounds inFeitelberg et al\. \([2024](https://arxiv.org/html/2606.04176#bib.bib10)\)are pointwise, in the sense that they characterize the estimation error at a fixed location\(i,j\)\(i,j\)\. Second, their consistency guarantees requireNNto diverge\. In contrast, our results provide global guarantees over the entire matrix and do not requireNNto grow\. More specifically, the termΔN\\Delta\_\{N\}captures the convergence rate of the empirical kernel mean embedding with respect to the per\-entry sample sizeNN\. Our bound indicates that consistency of the final estimator as long asmmdiverges\. This phenomenon arises from the intrinsic low\-rank structure ofμ∗\\mu^\{\\ast\}, which enables effective information sharing across observed entries\. As a result, the proposed estimator not only imputes distributions for unobserved entries but also yields more accurate estimates for observed entries with limited sample sizes\.

Next, we consider the setting whenμ∗\\mu^\{\*\}has exactly low multilinear rank\.

###### Theorem 5\.6\.

Under Assumptions[5\.1](https://arxiv.org/html/2606.04176#S5.Thmtheorem1)\-[5\.4](https://arxiv.org/html/2606.04176#S5.Thmtheorem4), suppose thatμ∗\\mu^\{\*\}has multilinear rank\(r1,r2,r3\)\(r\_\{1\},r\_\{2\},r\_\{3\}\), wherer1,r2,r3r\_\{1\},r\_\{2\},r\_\{3\}are of constant orders\. By taking

\(λβ1\+λβ2\)−1=𝒪\(\[M​πU​log⁡dm1​m2max\{ΔNlog⁡d,α\}\]−1\),\\displaystyle\(\\lambda\\beta\_\{1\}\+\\lambda\\beta\_\{2\}\)^\{\-1\}=\\mathop\{\}\\mathopen\{\}\\mathcal\{O\}\\mathopen\{\}\\left\(\\left\[\\frac\{\\sqrt\{M\\pi\_\{U\}\\log d\}\}\{m\_\{1\}m\_\{2\}\}\\max\\left\\\{\\Delta\_\{N\}\\sqrt\{\\log d\},\\alpha\\right\\\}\\right\]^\{\-1\}\\right\),we have with probability at least1−κ/d1\-\\kappa/dfor some constantκ\>0\\kappa\>0,

1m1​m2​∑i,j‖μ^​\(i,j,⋅\)−μ​\(i,j,⋅\)‖F2≲min⁡\{λ2​β12​m1​m2​r1πL2\+λ2​β22​m1​m2​r2πL2,λ2​β32​m1​m2​r3πL2\}\+\(2​α\)2​πU​log⁡dπL2​m1​m2\.\\frac\{1\}\{m\_\{1\}m\_\{2\}\}\\sum\_\{i,j\}\\\|\\hat\{\\mu\}\(i,j,\\cdot\)\-\\mu\(i,j,\\cdot\)\\\|^\{2\}\_\{F\}\\lesssim\\\\ \\min\\left\\\{\\frac\{\\lambda^\{2\}\\beta\_\{1\}^\{2\}m\_\{1\}m\_\{2\}r\_\{1\}\}\{\\pi\_\{L\}^\{2\}\}\+\\frac\{\\lambda^\{2\}\\beta\_\{2\}^\{2\}m\_\{1\}m\_\{2\}r\_\{2\}\}\{\\pi\_\{L\}^\{2\}\},\\frac\{\\lambda^\{2\}\\beta\_\{3\}^\{2\}m\_\{1\}m\_\{2\}r\_\{3\}\}\{\\pi\_\{L\}^\{2\}\}\\right\\\}\+\(2\\alpha\)^\{2\}\\sqrt\{\\frac\{\\pi\_\{U\}\\log d\}\{\\pi\_\{L\}^\{2\}m\_\{1\}m\_\{2\}\}\}\.\(12\)Furthermore, if we takeλ​β1,λ​β2,λ​β3≍M​πU​log⁡dm1​m2​max⁡\{ΔN​log⁡d,α\}\\lambda\\beta\_\{1\},\\lambda\\beta\_\{2\},\\lambda\\beta\_\{3\}\\asymp\\frac\{\\sqrt\{M\\pi\_\{U\}\\log d\}\}\{m\_\{1\}m\_\{2\}\}\\max\\left\\\{\\Delta\_\{N\}\\sqrt\{\\log d\},\\alpha\\right\\\}, then with probablity at least1−κ/d1\-\\kappa/d,we have

1m1​m2​∑i,j‖μ^​\(i,j,⋅\)−μ​\(i,j,⋅\)‖F2≲πUπL2​\[max⁡\{α,ΔN\}\]2​min⁡\{r1,r2,r3\}m​log2⁡d\+\(2​α\)2​πU​log⁡dπL2​m1​m2\.\\frac\{1\}\{m\_\{1\}m\_\{2\}\}\\sum\_\{i,j\}\\\|\\hat\{\\mu\}\(i,j,\\cdot\)\-\\mu\(i,j,\\cdot\)\\\|^\{2\}\_\{F\}\\lesssim\\\\ \\frac\{\\pi\_\{U\}\}\{\\pi\_\{L\}^\{2\}\}\\frac\{\[\\max\\\{\\alpha,\\Delta\_\{N\}\\\}\]^\{2\}\\min\\\{r\_\{1\},r\_\{2\},r\_\{3\}\\\}\}\{m\}\\log^\{2\}d\+\(2\\alpha\)^\{2\}\\sqrt\{\\frac\{\\pi\_\{U\}\\log d\}\{\\pi\_\{L\}^\{2\}m\_\{1\}m\_\{2\}\}\}\.

Theorem[5\.6](https://arxiv.org/html/2606.04176#S5.Thmtheorem6)improves upon the convergence rate established in Theorem[5\.5](https://arxiv.org/html/2606.04176#S5.Thmtheorem5)by fully exploiting the exact low\-rank structure ofμ∗\\mu^\{\*\}\. Moreover, the resulting rate matches that obtained for classical noisy low\-rank matrix completion with nuclear\-norm regularization\(Klopp,[2014](https://arxiv.org/html/2606.04176#bib.bib18); Negahban and Wainwright,[2012](https://arxiv.org/html/2606.04176#bib.bib28)\), which is known to be minimax optimal up to logarithmic factors in the scalar\-valued setting\.

## 6Simulation Studies

In this section, we present numerical experiments to illustrate the empirical performance of the proposed method\. We generate synthetic data according to the following procedure\. We first construct a matrix of distributions of dimensionm1×m2m\_\{1\}\\times m\_\{2\}\. Fixr1=r2=3r\_\{1\}=r\_\{2\}=3\. For each simulation run, we generate parametersal,vl​∼i\.i\.d\.​Unif​\[0,5\]a\_\{l\},v\_\{l\}\\overset\{\\mathrm\{i\.i\.d\.\}\}\{\\sim\}\\mathrm\{Unif\}\[0,5\]forł=1,…,r1​r2\\l =1,\\dots,r\_\{1\}r\_\{2\}\. Letg=\[0\.5,0\.5,0\.5\]⊺g=\[0\.5,0\.5,0\.5\]^\{\\intercal\}, and draw mixing weightswk​∼i\.i\.d\.​Dirichlet​\(g\)w\_\{k\}\\overset\{\\mathrm\{i\.i\.d\.\}\}\{\\sim\}\\mathrm\{Dirichlet\}\(g\)fork=1,…,m1k=1,\\dots,m\_\{1\}andwk′′​∼i\.i\.d\.​Dirichlet​\(g\)w^\{\\prime\}\_\{k^\{\\prime\}\}\\overset\{\\mathrm\{i\.i\.d\.\}\}\{\\sim\}\\mathrm\{Dirichlet\}\(g\)fork′=1,…,m2k^\{\\prime\}=1,\\dots,m\_\{2\}\. Each matrix entry\(i,j\)\(i,j\)is then assigned a Gaussian mixture distribution

∑l1=1r1∑l2=1r2wi,l1​wj,l2′​𝒩​\(al2\+r1​\(l1−1\),vl2\+r1​\(l1−1\)\)\.\\sum\_\{l\_\{1\}=1\}^\{r\_\{1\}\}\\sum\_\{l\_\{2\}=1\}^\{r\_\{2\}\}w\_\{i,l\_\{1\}\}w^\{\\prime\}\_\{j,l\_\{2\}\}\\mathcal\{N\}\(a\_\{l\_\{2\}\+r\_\{1\}\(l\_\{1\}\-1\)\},v\_\{l\_\{2\}\+r\_\{1\}\(l\_\{1\}\-1\)\}\)\.Next, we generate an observation indicator matrixTTby uniformly sampling entries with observation probabilitypp\. For each observed entry\(i,j\)\(i,j\)withT​\(i,j\)=1T\(i,j\)=1, we drawNNi\.i\.d\. samples from the corresponding distributionρ​\(i,j\)\\rho\(i,j\)\.

We compare our proposed estimator \(LRKME\) with an existing estimatorFeitelberg et al\. \([2024](https://arxiv.org/html/2606.04176#bib.bib10)\)\(NN\) using the following relative metrics defined on the whole matrix and unobserved entries:

Re​\(MMD\)=∑i,j‖μ​\(i,j,⋅\)−μ^​\(i,j,⋅\)‖ℋ2∑i,j‖μ​\(i,j,⋅\)‖ℋ2,\\displaystyle\\mathrm\{Re\}\(\\mathrm\{MMD\}\)=\\frac\{\\sum\_\{i,j\}\\\|\\mu\(i,j,\\cdot\)\-\\hat\{\\mu\}\(i,j,\\cdot\)\\\|^\{2\}\_\{\\mathcal\{H\}\}\}\{\\sum\_\{i,j\}\\\|\\mu\(i,j,\\cdot\)\\\|^\{2\}\_\{\\mathcal\{H\}\}\},Ret​\(MMD\)=∑T​\(i,j\)=0‖μ​\(i,j,⋅\)−μ^​\(i,j,⋅\)‖ℋ2∑T​\(i,j\)=0‖μ​\(i,j,⋅\)‖ℋ2,\\displaystyle\\mathrm\{Re\}\_\{t\}\(\\mathrm\{MMD\}\)=\\frac\{\\sum\_\{T\(i,j\)=0\}\\\|\\mu\(i,j,\\cdot\)\-\\hat\{\\mu\}\(i,j,\\cdot\)\\\|^\{2\}\_\{\\mathcal\{H\}\}\}\{\\sum\_\{T\(i,j\)=0\}\\\|\\mu\(i,j,\\cdot\)\\\|^\{2\}\_\{\\mathcal\{H\}\}\},whereμ^\\hat\{\\mu\}is a general notation for the estimated KME\. Beyond distributional accuracy, various summary statistics can be derived from the estimatedμ^\\hat\{\\mu\}\(orρ^\\hat\{\\rho\}\) for evaluation purpose\. In particular, we compute the estimated mean matrixM^\\hat\{M\}and variance matrixV^\\hat\{V\}, and compare them with the corresponding population mean and variance matricesMMandVV\. We also compute mean and variance estimators using classical matrix completion with nuclear\-norm regularization \(SoftImpute\), treating the sample means and sample variances as observed scalar entries\. The relative Frobenius norm errors are defined as

Re​\(M\)=‖M−M^‖F2‖M‖F2;Ret​\(M\)=‖\(1−T\)​\(M−M^\)‖F2‖\(1−T\)​M‖F2;\\displaystyle\\mathrm\{Re\}\(M\)=\\frac\{\\\|M\-\\hat\{M\}\\\|^\{2\}\_\{F\}\}\{\\\|M\\\|^\{2\}\_\{F\}\};\\mathrm\{Re\}\_\{t\}\(M\)=\\frac\{\\\|\(1\-T\)\(M\-\\hat\{M\}\)\\\|^\{2\}\_\{F\}\}\{\\\|\(1\-T\)M\\\|^\{2\}\_\{F\}\};Re​\(V\)=‖V−V^‖F2‖V‖F2;Ret​\(V\)=‖\(1−T\)​\(V−V^\)‖F2‖\(1−T\)​V‖F2;\\displaystyle\\mathrm\{Re\}\(V\)=\\frac\{\\\|V\-\\hat\{V\}\\\|^\{2\}\_\{F\}\}\{\\\|V\\\|^\{2\}\_\{F\}\};\\mathrm\{Re\}\_\{t\}\(V\)=\\frac\{\\\|\(1\-T\)\(V\-\\hat\{V\}\)\\\|^\{2\}\_\{F\}\}\{\\\|\(1\-T\)V\\\|^\{2\}\_\{F\}\};where‖A‖F\\\|A\\\|\_\{F\}refers to the Frobenius norm of matrixAA\. In our implementation, we employ the Gaussian kernel to construct KMEs and select the kernel bandwidth using the widely adopted median heuristic\(Garreau et al\.,[2017](https://arxiv.org/html/2606.04176#bib.bib12)\)\. For simplicity in tuning parameter selection, we setβj=1/3,j=1,2,3,\\beta\_\{j\}=1/3,j=1,2,3,in \([9](https://arxiv.org/html/2606.04176#S4.E9)\) and choose the regularization parameterλ\\lambdavia cross\-validation\. Cross\-validation is also used to select tuning parameters for all competing methods considered in our experiments\.

Table[1](https://arxiv.org/html/2606.04176#S6.T1)summarizes the results for the settingm1=m2=100m\_\{1\}=m\_\{2\}=100,p=0\.6p=0\.6, andN=5N=5, based on100100simulated datasets\. As shown, the proposed methodLRKMEconsistently outperformsNNacross all evaluation metrics, with particularly notable improvements on the unobserved entries\. For both mean and variance matrix estimation,LRKMEalso achieves better performance thanSoftImpute, which may be attributed to its ability to incorporate richer distributional information and to exploit additional structural properties of the data\. WhileNNperforms reasonably well in estimating the mean, it fails to deliver competitive variance estimates\.

In addition to the one\-dimensional setting, we also consider a two\-dimensional setting\. Specifically, we construct mean vectorsaℓ​∼i\.i\.d\.​Unif​\(\[0,1\]2\)a\_\{\\ell\}\\overset\{\\text\{i\.i\.d\.\}\}\{\\sim\}\\mathrm\{Unif\}\(\[0,1\]^\{2\}\)and covariance matricesVℓ=A​Aℓ⊤V\_\{\\ell\}=AA\_\{\\ell\}^\{\\top\}, whereAℓ∈ℝ2×2A\_\{\\ell\}\\in\\mathbb\{R\}^\{2\\times 2\}has entries independently generated from a standard Gaussian distribution,l=1,…,r1​r2l=1,\\dots,r\_\{1\}r\_\{2\}\. We adopt analogous relative error metrics for evaluation\. In particular, we report relative errors on the entire matrix for the first and second components of the mean matrix, denoted byRe​\(M​1\)\\mathrm\{Re\}\(M1\)andRe​\(M​2\)\\mathrm\{Re\}\(M2\), as well as for the covariance matrix, denoted byRe​\(Cov\)\\mathrm\{Re\}\(\\mathrm\{Cov\}\)\. Relative errors restricted to the unobserved entries are defined in the same manner\. For the methodNN, however, computing the two\-dimensional Wasserstein distance is considerably more challenging and does not readily extend to this setting\. As a result, we excludeNNfrom the comparison in the two\-dimensional experiments\. As Table[2](https://arxiv.org/html/2606.04176#S6.T2)shows,LRKMEperforms better thanSoftImputeby utilizing more structural information on the estimation of summary statistics\.

Table 1:Simulation results with one\-dimensional samples form1=m2=100m\_\{1\}=m\_\{2\}=100,p=0\.6p=0\.6andN=5N=5\. Values in the parentheses are standard errors\.Table 2:Simulation results with two\-dimensional samples form1=m2=50m\_\{1\}=m\_\{2\}=50,p=0\.3p=0\.3andN=80N=80\. Values in the parentheses are standard errors\.
## 7Real Data Application

In this section, we analyze the NYC taxi trip dataset \([https://www\.nyc\.gov/site/tlc/about/tlc\-trip\-record\-data\.page](https://www.nyc.gov/site/tlc/about/tlc-trip-record-data.page)\) to illustrate the performance of the proposed method\. We adopt the cleaned and organized data directly from the public repository[https://github\.com/wooner49/sofia](https://github.com/wooner49/sofia)\. The processed dataset records source locations, destination locations, and the number of taxi trips per day between each source\-destination pair over multiple days\.

We treat source and destination locations as matrix indices, yielding a265×265265\\times 265matrix, where each entry represents the distribution of daily trip counts between a given source–destination pair\. Due to data collection and preprocessing issues—such as missing or invalid location identifiers, timestamps, or corrupted logs—many entries are unobserved\. The overall missing rate exceeds39%39\\%, and the number of samples per observed entry varies widely, ranging from11to36723672\. For evaluation, we randomly select90%90\\%of the observed entries for training and reserve the remaining10%10\\%as the test set𝒯\\mathcal\{T\}\.

We compared the two distributional matrix completion methodsLRKME\(our proposal\) andNNbased on the empirical relative error

Re^t​\(MMD\)\\displaystyle\\hat\{\\mathrm\{Re\}\}\_\{t\}\(\\mathrm\{MMD\}\)=∑\(i,j\)∈𝒯‖1Ni​j​∑l=1Ni​jK​\(Xl,⋅\)−μ^​\(i,j,⋅\)‖ℋ2∑\(i,j\)∈𝒯‖∑l=1Ni​jK​\(Xl,⋅\)/Ni​j‖ℋ2,\\displaystyle=\\frac\{\\sum\_\{\(i,j\)\\in\{\\mathcal\{T\}\}\}\\\|\\frac\{1\}\{N\_\{ij\}\}\\sum\_\{l=1\}^\{N\_\{ij\}\}K\(X\_\{l\},\\cdot\)\-\\hat\{\\mu\}\(i,j,\\cdot\)\\\|\_\{\\mathcal\{H\}\}^\{2\}\}\{\\sum\_\{\(i,j\)\\in\{\\mathcal\{T\}\}\}\\\|\\sum\_\{l=1\}^\{N\_\{ij\}\}K\(X\_\{l\},\\cdot\)/N\_\{ij\}\\\|\_\{\\mathcal\{H\}\}^\{2\}\},whereμ^\\hat\{\\mu\}is a general notation for the estimated KME\. As Table[3](https://arxiv.org/html/2606.04176#S7.T3)shows,LRKMEachieves much smallerRe^t​\(MMD\)\\hat\{\\mathrm\{Re\}\}\_\{t\}\(\\mathrm\{MMD\}\)thanNN\.

Table 3:Empirical relative error for the test set\.We take a closer look at two representative test cases in which the number of observed samples exceeds 50\. The first case exhibits a relatively heavy\-tailed distribution, as shown in Figure[1](https://arxiv.org/html/2606.04176#S7.F1)\(a\), with values ranging from 1 to above 25, whereas the second case is highly concentrated around 1, as illustrated in Figure[2](https://arxiv.org/html/2606.04176#S7.F2)\(a\)\.

For test case 1, the KME estimated byLRKMEclosely matches the shape of the empirical KME\. In particular, both display relatively large values in the tail region whenxxis large, whereas the KME estimated byNNis predominantly concentrated at smaller values\. To facilitate a more direct comparison, we further transform the KME estimated byLRKMEinto a probability density function via numerical inversion \(see Appendix[B\.2](https://arxiv.org/html/2606.04176#A2.SS2)for details\)\. As shown in Figure[1](https://arxiv.org/html/2606.04176#S7.F1)\(d\), the resulting density function fromLRKMEexhibits a pattern consistent with the histogram of the observed samples\. In contrast, the histogram derived from the empirical samples used byNN\(Figure[1](https://arxiv.org/html/2606.04176#S7.F1)\(c\)\) is largely concentrated at value 1 and fails to capture the distributional mass beyond 11\.

For test case 2, the KME estimated byLRKMEexhibits the same overall trend as the empirical KME, attaining its maximum at value 1 and then decreasing\. In contrast, the KME estimated byNNreaches its maximum at a value greater than 1\.5\. This discrepancy is also evident from the histogram of the empirical samples used byNN\(Figure[2](https://arxiv.org/html/2606.04176#S7.F2)\(c\)\), which has its mode at value 2 and extends beyond value 3\. By comparison, the density estimated fromLRKME\(Figure[2](https://arxiv.org/html/2606.04176#S7.F2)\(d\)\) assigns a dominantly high density to value 1, more accurately reflecting the underlying empirical distribution\.

![Refer to caption](https://arxiv.org/html/2606.04176v1/hist1.png)\(a\)Histogram for observed samples\.
![Refer to caption](https://arxiv.org/html/2606.04176v1/KME1.png)\(b\)Empirical KME, estimated KMEs fromNNandLRKME\.
![Refer to caption](https://arxiv.org/html/2606.04176v1/NNhist1.png)\(c\)Histogram for empirical samples estimated fromNN\.
![Refer to caption](https://arxiv.org/html/2606.04176v1/den1.png)\(d\)Estimated density from estimated KME fromLRKME\.

Figure 1:Results for test case 1\.![Refer to caption](https://arxiv.org/html/2606.04176v1/hist2.png)\(a\)Histogram for observed samples\.
![Refer to caption](https://arxiv.org/html/2606.04176v1/KME2.png)\(b\)Empirical KME, estimated KMEs fromNNandLRKME\.
![Refer to caption](https://arxiv.org/html/2606.04176v1/NNhist2.png)\(c\)Histogram for empirical samples estimated fromNN\.
![Refer to caption](https://arxiv.org/html/2606.04176v1/den2.png)\(d\)Estimated density from estimated KME fromLRKME\.

Figure 2:Results for test case 2\.
## References

- Boyd et al\. \[2010\]Stephen Boyd, Neal Parikh, Eric Chu, Borja Peleato, and Jonathan Eckstein\.Distributed optimization and statistical learning via the alternating direction method of multipliers\.*Foundations and Trends® in Machine Learning*, 3\(1\):1–122, 2010\.
- Candes and Plan \[2010\]Emmanuel J Candes and Yaniv Plan\.Matrix completion with noise\.*Proceedings of the IEEE*, 98\(6\):925–936, 2010\.
- Cao et al\. \[2014\]Feilong Cao, Miaomiao Cai, and Yuanpeng Tan\.Image interpolation via low\-rank matrix completion and recovery\.*IEEE Transactions on Circuits and Systems for Video Technology*, 25\(8\):1261–1270, 2014\.
- Changjun et al\. \[2012\]Fu Changjun, Ji Xiangyang, Zhang Yongbing, and Dai Qionghai\.A single frame super\-resolution method based on matrix completion\.In*2012 Data Compression Conference*, pages 297–306\. IEEE, 2012\.
- Chen et al\. \[2023\]Junren Chen, Cheng\-Long Wang, Michael K Ng, and Di Wang\.High dimensional statistical estimation under uniformly dithered one\-bit quantization\.*IEEE Transactions on Information Theory*, 2023\.
- Chen et al\. \[2017\]Xiaobo Chen, Zhongjie Wei, Zuoyong Li, Jun Liang, Yingfeng Cai, and Bob Zhang\.Ensemble correlation\-based low\-rank matrix completion with applications to traffic data imputation\.*Knowledge\-Based Systems*, 132:249–262, 2017\.
- Chen et al\. \[2020\]Xinyu Chen, Jinming Yang, and Lijun Sun\.A nonconvex low\-rank tensor completion model for spatiotemporal traffic data imputation\.*Transportation Research Part C: Emerging Technologies*, 117:102673, 2020\.
- Chen et al\. \[2021\]Zhaoliang Chen, Wei Zhao, and Shiping Wang\.Kernel meets recommender systems: A multi\-kernel interpolation for matrix completion\.*Expert Systems with Applications*, 168:114436, 2021\.
- Drineas and Mahoney \[2005\]Petros Drineas and Michael W\. Mahoney\.On the nyström method for approximating a gram matrix for improved kernel\-based learning\.*Journal of Machine Learning Research*, 6:2153–2175, 2005\.
- Feitelberg et al\. \[2024\]Jacob Feitelberg, Kyuseong Choi, Anish Agarwal, and Raaz Dwivedi\.Distributional matrix completion via nearest neighbors in the wasserstein space\.*arXiv preprint arXiv:2410\.13112*, 2024\.
- Gandy et al\. \[2011\]Silvia Gandy, Benjamin Recht, and Isao Yamada\.Tensor completion and low\-n\-rank tensor recovery via convex optimization\.*Inverse Problems*, 27\(2\):025010, 2011\.
- Garreau et al\. \[2017\]Damien Garreau, Wittawat Jitkrittum, and Motonobu Kanagawa\.Large sample analysis of the median heuristic\.*arXiv preprint arXiv:1707\.07269*, 2017\.
- Gretton et al\. \[2012\]Arthur Gretton, Karsten M\. Borgwardt, Malte J\. Rasch, Bernhard Schölkopf, and Alex J\. Smola\.A kernel two\-sample test\.*Journal of Machine Learning Research*, 13:723–773, 2012\.
- Gurini et al\. \[2018\]Davide Feltoni Gurini, Fabio Gasparetti, Alessandro Micarelli, and Giuseppe Sansonetti\.Temporal people\-to\-people recommendation on social networks with sentiment\-based matrix factorization\.*Future Generation Computer Systems*, 78:430–439, 2018\.
- Han et al\. \[2024\]Rungang Han, Pixu Shi, and Anru R Zhang\.Guaranteed functional tensor singular value decomposition\.*Journal of the American Statistical Association*, 119\(546\):995–1007, 2024\.
- Kadkhodaie and Montanari \[2015\]Mahdi Kadkhodaie and Andrea Montanari\.Accelerated alternating direction method of multipliers\.*arXiv preprint arXiv:1505\.01883*, 2015\.
- Kang et al\. \[2016\]Zhao Kang, Chong Peng, and Qiang Cheng\.Top\-n recommender system via matrix completion\.In*Proceedings of the AAAI conference on artificial intelligence*, volume 30, 2016\.
- Klopp \[2014\]Olga Klopp\.Noisy low\-rank matrix completion with general sampling distribution\.*Bernoulli*, 2014\.
- Kolda and Bader \[2009\]Tamara G\. Kolda and Brett W\. Bader\.Tensor decompositions and applications\.*SIAM Review*, 51\(3\):455–500, 2009\.
- Koltchinskii \[2011\]Vladimir Koltchinskii\.*Oracle inequalities in empirical risk minimization and sparse recovery problems: Ecole D’Eté de Probabilités de Saint\-Flour XXXVIII\-2008*, volume 2033\.Springer, 2011\.
- Koltchinskii et al\. \[2011\]Vladimir Koltchinskii, Karim Lounici, and Alexandre B\. Tsybakov\.Nuclear\-norm penalization and optimal rates for noisy low\-rank matrix completion\.*The Annals of Statistics*, 39\(5\):2302–2329, 2011\.
- Larsen et al\. \[2024\]Brett W Larsen, Tamara G Kolda, Anru R Zhang, and Alex H Williams\.Tensor decomposition meets rkhs: Efficient algorithms for smooth and misaligned data\.*arXiv preprint arXiv:2408\.05677*, 2024\.
- Li et al\. \[2024\]Jiangyuan Li, Jiayi Wang, Raymond KW Wong, and Kwun Chuen Gary Chan\.A pairwise pseudo\-likelihood approach for matrix completion with informative missingness\.*Advances in Neural Information Processing Systems*, 37:10735–10769, 2024\.
- Liu et al\. \[2013\]Ji Liu, Przemyslaw Musialski, Peter Wonka, and Jieping Ye\.Tensor completion for estimating missing values in visual data\.*IEEE Transactions on Pattern Analysis and Machine Intelligence*, 35\(1\):208–220, 2013\.
- Mazumder et al\. \[2010\]Rahul Mazumder, Trevor Hastie, and Robert Tibshirani\.Spectral regularization algorithms for learning large incomplete matrices\.*The Journal of Machine Learning Research*, 11:2287–2322, 2010\.
- Mu et al\. \[2014\]Cun Mu, Bo Huang, John Wright, and Donald Goldfarb\.Square deal: Lower bounds and improved relaxations for tensor recovery\.In*International Conference on Machine Learning*, pages 73–81, 2014\.
- Muandet et al\. \[2017\]Krikamol Muandet, Kenji Fukumizu, Bharath K\. Sriperumbudur, and Bernhard Schölkopf\.Kernel mean embedding of distributions: A review and beyond\.*Foundations and Trends in Machine Learning*, 10\(1\-2\):1–141, 2017\.
- Negahban and Wainwright \[2012\]Sahand N\. Negahban and Martin J\. Wainwright\.Restricted strong convexity and weighted matrix completion: Optimal bounds with noise\.*Journal of Machine Learning Research*, 13:1665–1697, 2012\.
- Rahimi and Recht \[2007\]Ali Rahimi and Benjamin Recht\.Random features for large\-scale kernel machines\.In*Advances in Neural Information Processing Systems*, 2007\.
- Smola et al\. \[2007\]Alex J\. Smola, Arthur Gretton, Le Song, and Bernhard Schölkopf\.A hilbert space embedding for distributions\.*Journal of Machine Learning Research*, 7:1449–1472, 2007\.
- Sriperumbudur et al\. \[2010\]Bharath K\. Sriperumbudur, Arthur Gretton, Kenji Fukumizu, Bernhard Schölkopf, and Gert R\. G\. Lanckriet\.Hilbert space embeddings and metrics on probability measures\.*Journal of Machine Learning Research*, 11:1517–1561, 2010\.
- Tucker \[1966\]Ledyard R\. Tucker\.Some mathematical notes on three\-mode factor analysis\.*Psychometrika*, 31\(3\):279–311, 1966\.
- Vershynin \[2018\]Roman Vershynin\.*High\-dimensional probability: An introduction with applications in data science*, volume 47\.Cambridge university press, 2018\.
- Wang et al\. \[2021\]Jiayi Wang, Raymond K\. W\. Wong, Xiaojun Mao, and Kwun Chuen Gary Chan\.Matrix completion with model\-free weighting\.In*International Conference on Machine Learning*, pages 10927–10936\. PMLR, 2021\.
- Wang et al\. \[2022\]Jiayi Wang, Raymond K\. W\. Wong, and Xiaoke Zhang\.Low\-rank covariance function estimation for multidimensional functional data\.*Journal of the American statistical Association*, 117\(538\):809–822, 2022\.
- Weng and Wang \[2012\]Zhiyuan Weng and Xin Wang\.Low\-rank matrix completion for array signal processing\.In*International Conference on Acoustics, Speech and Signal Processing \(ICASSP\)*, pages 2697–2700\. IEEE, 2012\.
- Williams and Seeger \[2001\]Christopher K\. I\. Williams and Matthias Seeger\.Using the nyström method to speed up kernel machines\.*Advances in Neural Information Processing Systems*, 2001\.
- Xu et al\. \[2023\]Xiuqin Xu, Mingwei Lin, Xin Luo, and Zeshui Xu\.HRST\-LR: a hessian regularization spatio\-temporal low rank algorithm for traffic data imputation\.*IEEE Transactions on Intelligent Transportation Systems*, 24:11001–11017, 2023\.
- Yuchi et al\. \[2023\]Henry Shaowu Yuchi, Simon Mak, and Yao Xie\.Bayesian uncertainty quantification for low\-rank matrix completion\.*Bayesian Analysis*, 18\(2\):491–518, 2023\.
- Zhang et al\. \[2021\]Qian Zhang, Jie Lu, and Yaochu Jin\.Artificial intelligence in recommender systems\.*Complex & Intelligent Systems*, 7\(1\):439–457, 2021\.
- Zhang and Zhang \[2020\]Shuimei Zhang and Yimin D Zhang\.Low\-rank hankel matrix completion for robust time\-frequency analysis\.*IEEE Transactions on Signal Processing*, 68:6171–6186, 2020\.
- Zhao et al\. \[2025\]Kun Zhao, Jiayi Wang, and Yifei Lou\.Noisy low\-rank matrix completion via transformedℓ1\\ell\_\{1\}regularization and its theoretical properties\.In*Proceedings of the 28th International Conference on Artificial Intelligence and Statistics*, Proceedings of Machine Learning Research\. PMLR, 2025\.
- Zheng et al\. \[2024\]Huiwen Zheng, Yifei Lou, Guoliang Tian, and Chao Wang\.A scale\-invariant relaxation in low\-rank tensor recovery with an application to tensor completion\.*SIAM Journal on Imaging Sciences*, 17\(1\):756–783, 2024\.

## Appendix AAdditional Information Related to Tensor

###### Definition A\.1\(Matricization\)\.

Letj∈\{1,…,d\}j\\in\\\{1,\\dots,d\\\}\. Thejj\-mode matricization of a tensorB∈ℝq1×q2×⋯×qdB\\in\\mathbb\{R\}^\{q\_\{1\}\\times q\_\{2\}\\times\\cdots\\times q\_\{d\}\}, denoted byB\(j\)B\_\{\(j\)\}, rearranges the entries ofBBinto a matrix of sizeqj×\(∏i≠jqi\)q\_\{j\}\\times\\bigl\(\\prod\_\{i\\neq j\}q\_\{i\}\\bigr\)\. Specifically, the entry ofB\(j\)B\_\{\(j\)\}indexed by\(lj,k\)\(l\_\{j\},k\)corresponds to the tensor elementBl1,…,ldB\_\{l\_\{1\},\\ldots,l\_\{d\}\}, where the column indexkkis determined byk=1\+∑i=1,i≠jd\(li−1\)​\(∏m=1,m≠ji−1qm\)\.k=1\+\\sum\_\{i=1,\\,i\\neq j\}^\{d\}\(l\_\{i\}\-1\)\\Bigl\(\\prod\_\{m=1,\\,m\\neq j\}^\{i\-1\}q\_\{m\}\\Bigr\)\.111All empty products are defined to be11\. For example,∏m=ijqm=1\\prod\_\{m=i\}^\{j\}q\_\{m\}=1wheni\>ji\>j\.

###### Definition A\.2\(jj\-mode product\)\.

Letj∈\{1,…,d\}j\\in\\\{1,\\dots,d\\\}\. For an arrayB∈ℝq1×q2×⋯×qdB\\in\\mathbb\{R\}^\{q\_\{1\}\\times q\_\{2\}\\times\\cdots\\times q\_\{d\}\}and a matrixP∈ℝpj×qjP\\in\\mathbb\{R\}^\{p\_\{j\}\\times q\_\{j\}\}, thejj\-mode product ofBBandPP, denoted byB×jPB\\times\_\{j\}P, is an array of dimension

q1×⋯×qj−1×pj×qj\+1×⋯×qd\.q\_\{1\}\\times\\cdots\\times q\_\{j\-1\}\\times p\_\{j\}\\times q\_\{j\+1\}\\times\\cdots\\times q\_\{d\}\.The\(l1,…,lj−1,k,lj\+1,…,ld\)\(l\_\{1\},\\dots,l\_\{j\-1\},k,l\_\{j\+1\},\\dots,l\_\{d\}\)\-th entry ofB×jPB\\times\_\{j\}Pis defined as

\(B×jP\)l1,…,lj−1,k,lj\+1,…,ld=∑i=1qjBl1,…,lj−1,i,lj\+1,…,ld​Pk,i\.\(B\\times\_\{j\}P\)\_\{l\_\{1\},\\dots,l\_\{j\-1\},k,l\_\{j\+1\},\\dots,l\_\{d\}\}=\\sum\_\{i=1\}^\{q\_\{j\}\}B\_\{l\_\{1\},\\dots,l\_\{j\-1\},i,l\_\{j\+1\},\\dots,l\_\{d\}\}\\,P\_\{k,i\}\.

Figure[3](https://arxiv.org/html/2606.04176#A1.F3)provides a pictorial illustration of a Tucker decomposition\.

BB\(q1×q2×q3\)\(q\_\{1\}\\times q\_\{2\}\\times q\_\{3\}\)

=

U1U\_\{1\}\(q1×r1\)\(q\_\{1\}\\times r\_\{1\}\)GG\(r1×r2×r3\)\(r\_\{1\}\\times r\_\{2\}\\times r\_\{3\}\)U3U\_\{3\}\(q3×r3\)\(q\_\{3\}\\times r\_\{3\}\)U2U\_\{2\}\(q2×r2\)\(q\_\{2\}\\times r\_\{2\}\)Figure 3:Tucker decomposition of a third\-order arrayB∈ℝq1×q2×q3B\\in\\mathbb\{R\}^\{q\_\{1\}\\times q\_\{2\}\\times q\_\{3\}\}\. The values in the parentheses are dimensions for the corresponding matrices or tensors\.
## Appendix BComputational Details

### B\.1Algorithm

In this section, we introduce the accelerated ADMM algorithm for solving \([11](https://arxiv.org/html/2606.04176#S4.E11)\)\.

We begin with an alternative form of \([11](https://arxiv.org/html/2606.04176#S4.E11)\):

minB∈ℝm1×m2×q⁡\{ℓ​\(B\)\+λ​\(∑k=13βj​‖Dj‖∗\)\}\.\\displaystyle\\min\_\{B\\in\\mathbb\{R\}^\{m\_\{1\}\\times m\_\{2\}\\times q\}\}\\left\\\{\\ell\(B\)\+\\lambda\\left\(\\sum\_\{k=1\}^\{3\}\\beta\_\{j\}\\\|D\_\{j\}\\\|\_\{\*\}\\right\)\\right\\\}\.\(13\)subject to𝒰k​Bk=Dk,k=1,…,3,supi,j‖B​\(i,j,⋅\)‖≤α,\\displaystyle\\text\{subject to\}\\ \\ \{\\mathcal\{U\}\}\_\{k\}B\_\{k\}=D\_\{k\},\\quad k=1,\\dots,3,\\quad\\sup\_\{i,j\}\\\|B\(i,j,\\cdot\)\\\|\\leq\\alpha,\(14\)whereℓ​\(B\)=∑i​jTi​j​‖B​\(i,j,⋅\)−Y​\(i,j,⋅\)‖22\\ell\(B\)=\\sum\_\{ij\}T\_\{ij\}\\\|B\(i,j,\\cdot\)\-Y\(i,j,\\cdot\)\\\|\_\{2\}^\{2\},D1∈ℝm1×\(m2​q\)D\_\{1\}\\in\\mathbb\{R\}^\{m\_\{1\}\\times\(m\_\{2\}q\)\},D2∈ℝm2×\(m1​q\)D\_\{2\}\\in\\mathbb\{R\}^\{m\_\{2\}\\times\(m\_\{1\}q\)\}andD3∈ℝq×\(m​1​m​2\)D\_\{3\}\\in\\mathbb\{R\}^\{q\\times\(m1m2\)\}\.

Then a standard ADMM algorithm solves the optimization problem \([13](https://arxiv.org/html/2606.04176#A2.E13)\) by minimizing the augmented Lagrangian with respect to different variables alternatively\. More explicitly, at the\(t\+1\)\(t\+1\)\-th iteration, the following updates are implemented:

B\(t\+1\)\\displaystyle\{B\}^\{\(t\+1\)\}=argminsupi,j‖B​\(i,j,⋅\)‖≤α\{ℓ​\(B\)\+η2​∑k=13‖𝒰k​B−Dk\(t\)\+Vk\(t\)‖F2\},\\displaystyle=\\operatorname\*\{argmin\}\_\{\\sup\_\{i,j\}\\\|B\(i,j,\\cdot\)\\\|\\leq\\alpha\}\\left\\\{\{\\ell\}\(B\)\+\\frac\{\\eta\}\{2\}\\sum\_\{k=1\}^\{3\}\\left\\\|\{\\mathcal\{U\}\}\_\{k\}B\-D\_\{k\}^\{\(t\)\}\+V\_\{k\}^\{\(t\)\}\\right\\\|\_\{F\}^\{2\}\\right\\\},\(15a\)Dk\(t\+1\)\\displaystyle D\_\{k\}^\{\(t\+1\)\}=argminDk\{λ​βk​‖Dk‖∗\+η2​‖𝒰k​B\(t\+1\)−Dk\+Vk\(t\)‖F2\},k=1,2,3,\\displaystyle=\\operatorname\*\{argmin\}\_\{D\_\{k\}\}\\left\\\{\\lambda\\beta\_\{k\}\\\|D\_\{k\}\\\|\_\{\*\}\+\\frac\{\\eta\}\{2\}\\left\\\|\{\\mathcal\{U\}\}\_\{k\}B^\{\(t\+1\)\}\-D\_\{k\}\+V\_\{k\}^\{\(t\)\}\\right\\\|\_\{F\}^\{2\}\\right\\\},\\ k=1,2,3,\(15b\)Vk\(t\+1\)\\displaystyle V\_\{k\}^\{\(t\+1\)\}=Vk\(t\)\+𝒰k​B\(t\+1\)−Dk\(t\+1\),k=1,2,3,\\displaystyle=V\_\{k\}^\{\(t\)\}\+\{\\mathcal\{U\}\}\_\{k\}B^\{\(t\+1\)\}\-D\_\{k\}^\{\(t\+1\)\},\\ k=1,2,3,\(15c\)whereVkV\_\{k\}is of the same dimension ofDkD\_\{k\}, fork=1,2,3k=1,2,3, are scaled Lagrangian multipliers andη\>0\\eta\>0is an algorithmic parameter\. An adaptive strategy to tuneη\\etais provided inBoyd et al\. \[[2010](https://arxiv.org/html/2606.04176#bib.bib1)\]\. Steps \([15a](https://arxiv.org/html/2606.04176#A2.E15.1)\), \([15b](https://arxiv.org/html/2606.04176#A2.E15.2)\) involve additional optimizations\. and now we discuss how to solve them\.

The objective function of \([15a](https://arxiv.org/html/2606.04176#A2.E15.1)\) is a quadratic function, and so we can easily solve this with a closed\-form solution, and then perform a projection to let it satisfy the sup constraint, given in Algorithm[1](https://arxiv.org/html/2606.04176#alg1)\. To solve \([15b](https://arxiv.org/html/2606.04176#A2.E15.2)\), we use proximal operatorproxv\\mathrm\{prox\}\_\{v\}defined by

proxv​\(A\)\\displaystyle\\mathrm\{prox\}\_\{v\}\(A\)=argminW\{12​‖W−A‖F2\+v​‖W‖∗\}\\displaystyle=\\operatorname\*\{argmin\}\_\{W\}\\left\\\{\\frac\{1\}\{2\}\\\|W\-A\\\|\_\{F\}^\{2\}\+v\\\|W\\\|\_\{\*\}\\right\\\}\(16a\)forv≥0v\\geq 0, whereA∈ℝn1×n2A\\in\\mathbb\{R\}^\{n\_\{1\}\\times n\_\{2\}\}denotes a general matrix\. By Lemma 1 inMazumder et al\. \[[2010](https://arxiv.org/html/2606.04176#bib.bib25)\], the solutions to \([16](https://arxiv.org/html/2606.04176#A3.EGx19)\) have closed forms\. More specifically, write the singular value decomposition ofAAasU​diag​\(\(a~1,…,a~min⁡\{n1,n2\}\)\)​V⊺U\\mathrm\{diag\}\(\(\\tilde\{a\}\_\{1\},\\dots,\\tilde\{a\}\_\{\\min\\\{n\_\{1\},n\_\{2\}\\\}\}\)\)V^\{\\intercal\}, thenproxv​\(A\)=U​diag​\(c~\)​V⊺\\mathrm\{prox\}\_\{v\}\(A\)=U\\mathrm\{diag\}\(\\tilde\{c\}\)V^\{\\intercal\}wherec~=\(\(a~1−v\)\+,\(a~2−v\)\+,…,\(a~min⁡\{n1,n2\}−v\)\+\)\\tilde\{c\}=\(\(\\tilde\{a\}\_\{1\}\-v\)\_\{\+\},\(\\tilde\{a\}\_\{2\}\-v\)\_\{\+\},\\dots,\(\\tilde\{a\}\_\{\\min\\\{n\_\{1\},n\_\{2\}\\\}\}\-v\)\_\{\+\}\)\.

The details of the computational procedure are summarized in Algorithm[1](https://arxiv.org/html/2606.04176#alg1), which presents an accelerated variant of the alternating direction method of multipliers \(ADMM\) incorporating additional steps to improve convergence speed\. We denote by𝒰k−1\\mathcal\{U\}\_\{k\}^\{\-1\}the inverse of𝒰k\\mathcal\{U\}\_\{k\}\. LetT~∈ℝm1×m2×q\\tilde\{T\}\\in\\mathbb\{R\}^\{m\_\{1\}\\times m\_\{2\}\\times q\}denote the indicator tensor, whereT~​\(i,j,k\)=1\\tilde\{T\}\(i,j,k\)=1ifTi​j=1T\_\{ij\}=1andT~​\(i,j,k\)=0\\tilde\{T\}\(i,j,k\)=0otherwise\. We further defineI~∈ℝm1×m2×q\\tilde\{I\}\\in\\mathbb\{R\}^\{m\_\{1\}\\times m\_\{2\}\\times q\}as the identity tensor with all entries equal to one\.

Algorithm 1Accelerated ADMM for solving \([11](https://arxiv.org/html/2606.04176#S4.E11)\)Input:Tensor

YY, parameters

λ\>0\\lambda\>0,

β1\\beta\_\{1\},

β2\\beta\_\{2\},

β3≥0\\beta\_\{3\}\\geq 0,

η\>0\\eta\>0,

TT,

ϵ\\epsilon\.

Initialize:

V1\(0\),D1\(0\)∈ℝm1×\(m1​q\)\{V\}\_\{1\}^\{\(0\)\},\{D\}\_\{1\}^\{\(0\)\}\\in\\mathbb\{R\}^\{m\_\{1\}\\times\(m\_\{1\}q\)\},

V2\(0\),D2\(0\)∈ℝm2×\(m2​q\)\{V\}\_\{2\}^\{\(0\)\},\{D\}\_\{2\}^\{\(0\)\}\\in\\mathbb\{R\}^\{m\_\{2\}\\times\(m\_\{2\}q\)\}, and

V3\(0\),D3\(0\)∈ℝq×\(m1​m2\)\{V\}\_\{3\}^\{\(0\)\},\{D\}\_\{3\}^\{\(0\)\}\\in\\mathbb\{R\}^\{q\\times\(m\_\{1\}m\_\{2\}\)\}\.

B\(0\)∈ℝm1×m2×qB^\{\(0\)\}\\in\\mathbb\{R\}^\{m\_\{1\}\\times m\_\{2\}\\times q\}\.

for

t=1t=1to

TTdo

B~\(t\)=\[2​Y\+∑k=13\{η​\(𝒰k−1\)​Dk\(t−1\)−η​\(𝒰k−1\)​Vk\(t−1\)\}\]/\[T~\+3∗η∗I~\],\\tilde\{B\}^\{\(t\)\}=\\left\[2Y\+\\sum\_\{k=1\}^\{3\}\\\{\\eta\(\{\\mathcal\{U\}\}\_\{k\}^\{\-1\}\)D\_\{k\}^\{\(t\-1\)\}\-\\eta\(\{\\mathcal\{U\}\}\_\{k\}^\{\-1\}\)V\_\{k\}^\{\(t\-1\)\}\\\}\\right\]/\[\\tilde\{T\}\+3\*\\eta\*\\tilde\{I\}\],where the operator

//denotes element\-wise division between tensors\. Then

B\(t\)B^\{\(t\)\}is constructed as

B\(t\)​\(i,j,⋅\)=\{α∗B~​\(i,j,⋅\)/‖B~​\(i,j,⋅\)‖,if​‖B~​\(i,j,⋅\)‖\>α,B~\(t\)​\(i,j,⋅\),otherwise\.B^\{\(t\)\}\(i,j,\\cdot\)=\\begin\{cases\}\\alpha\*\\tilde\{B\}\(i,j,\\cdot\)/\\\|\\tilde\{B\}\(i,j,\\cdot\)\\\|,&\\text\{if \}\\\|\\tilde\{B\}\(i,j,\\cdot\)\\\|\>\\alpha,\\\\ \\tilde\{B\}^\{\(t\)\}\(i,j,\\cdot\),&\\text\{otherwise\}\.\\end\{cases\}Dk\(t\)=proxλ​βk/η​\(𝒰k​B\(t\)−Vk\(t−1\)\),k=1,2,3\.D\_\{k\}^\{\(t\)\}=\\mathrm\{prox\}\_\{\\lambda\\beta\_\{k\}/\\eta\}\(\{\\mathcal\{U\}\}\_\{k\}B^\{\(t\)\}\-V\_\{k\}^\{\(t\-1\)\}\),\\qquad k=1,2,3\.Vk\(t\)=Vk\(t−1\)\+𝒰k​B\(t\)−Dk\(t\),k=1,2,3\.V\_\{k\}^\{\(t\)\}=V\_\{k\}^\{\(t\-1\)\}\+\{\\mathcal\{U\}\}\_\{k\}B^\{\(t\)\}\-D\_\{k\}^\{\(t\)\},\\qquad k=1,2,3\.
if

‖B\(t\)−B\(t−1\)‖F2/‖B\(t−1\)‖F≤ϵ\\\|B^\{\(t\)\}\-B^\{\(t\-1\)\}\\\|\_\{F\}^\{2\}/\\\|B^\{\(t\-1\)\}\\\|\_\{F\}\\leq\\epsilonthen

Break:set

B\(T\)=B\(t\)B^\{\(T\)\}=B^\{\(t\)\}
endif

endfor

### B\.2KME to Density

In this section, we provide a way to numerically transform a KME functionμ\\muto the corresponding density functionρ\\rho\. Set a grid over the range of𝒳\{\\mathcal\{X\}\}asX1,…,XKX\_\{1\},\\dots,X\_\{K\}\. Note that

μ​\(t\)=𝔼X​K​\(X,t\)=∫x∈𝒳K​\(x,t\)​ρ​\(x\)​dx≈1K​∑t=1KK​\(Xt,t\)​ρ​\(Xt\)\.\\displaystyle\\mu\(t\)=\{\\mathbb\{E\}\}\_\{X\}K\(X,t\)=\\int\_\{x\\in\{\\mathcal\{X\}\}\}K\(x,t\)\\rho\(x\)\\mathrm\{d\}x\\approx\\frac\{1\}\{K\}\\sum\_\{t=1\}^\{K\}K\(X\_\{t\},t\)\\rho\(X\_\{t\}\)\.Then we could solveρ​\(Xt\),t=1,…,K\\rho\(X\_\{t\}\),t=1,\\dots,Kvia the following optimization

minρ​\(Xt\),t=1,…,K​∑t′=1T\(μ​\(Xt′\)−1K​∑t=1KK​\(Xt′,t\)​ρ​\(Xt\)\)2,\\displaystyle\\min\_\{\\rho\(X\_\{t\}\),t=1,\\dots,K\}\\sum\_\{t^\{\\prime\}=1\}^\{T\}\\left\(\\mu\(X\_\{t\}^\{\\prime\}\)\-\\frac\{1\}\{K\}\\sum\_\{t=1\}^\{K\}K\(X\_\{t^\{\\prime\}\},t\)\\rho\(X\_\{t\}\)\\right\)^\{2\},\(17\)subject​to​ρ​\(Xt\),t=1,…,K≥0\.\\displaystyle\\mathrm\{subject\\ to\}\\rho\(X\_\{t\}\),t=1,\\dots,K\\geq 0\.\(18\)

## Appendix CProof

We introduce some notations for simplicity\. We take𝒰j​μ\{\\mathcal\{U\}\}\_\{j\}\\muasμ\(j\)\\mu\_\{\(j\)\}andλj=λ​βj\\lambda\_\{j\}=\\lambda\\beta\_\{j\}, forj=1,2,3j=1,2,3\.

### C\.1Proof of Theorem[4\.1](https://arxiv.org/html/2606.04176#S4.Thmtheorem1)

For anyμ∈ℝm1⊗ℝm2⊗ℋ\\mu\\in\\mathbb\{R\}^\{m\_\{1\}\}\\otimes\\mathbb\{R\}^\{m\_\{2\}\}\\otimes\{\\mathcal\{H\}\}, we can decompose it into two orthogonal partsμ1\\mu\_\{1\}andμ2\\mu\_\{2\}such thatμ1∈𝒦\\mu\_\{1\}\\in\{\\mathcal\{K\}\}andμ2∈\(𝒦\)⊥\\mu\_\{2\}\\in\\left\(\{\\mathcal\{K\}\}\\right\)^\{\\bot\}\. Note that for anyi,ji,j, by the reproducing property of kernelKKand the definition ofμ2\\mu\_\{2\}, the loss

Ti​j​‖μ​\(i,j,⋅\)−1Ni​j​∑l=1Ni​jK​\(Xl\(i​j\),⋅\)‖ℋ2=Ti​j​‖μ1​\(i,j,⋅\)−1Ni​j​∑l=1Ni​jK​\(Xl\(i​j\),⋅\)‖ℋ2\.T\_\{ij\}\\left\\\|\\mu\(i,j,\\cdot\)\-\\frac\{1\}\{N\_\{ij\}\}\\sum\_\{l=1\}^\{N\_\{ij\}\}K\(X^\{\(ij\)\}\_\{l\},\\cdot\)\\right\\\|\_\{\\mathcal\{H\}\}^\{2\}=T\_\{ij\}\\left\\\|\\mu\_\{1\}\(i,j,\\cdot\)\-\\frac\{1\}\{N\_\{ij\}\}\\sum\_\{l=1\}^\{N\_\{ij\}\}K\(X^\{\(ij\)\}\_\{l\},\\cdot\)\\right\\\|\_\{\\mathcal\{H\}\}^\{2\}\.Then it suffices to show that‖𝒰k​μ‖∗≥‖𝒰k​μ1‖∗\\\|\{\\mathcal\{U\}\}\_\{k\}\\mu\\\|\_\{\*\}\\geq\\\|\{\\mathcal\{U\}\}\_\{k\}\\mu\_\{1\}\\\|\_\{\*\}fork=1,2,3k=1,2,3\. LetP𝒦P\_\{\{\\mathcal\{K\}\}\}be the projection operator to space𝒦\{\\mathcal\{K\}\}and denote the adjoint operator ofAAbyA∗A^\{\*\}\. Denoteσk​\(A\)\\sigma\_\{k\}\(A\)askk\-th singular value of the operatorAA\. We have

σk​\(𝒰j​μ1\)=σk​\(𝒰j​μ​P𝒦\)≤σk​\(𝒰j​μ\),j=1,2\\displaystyle\\sigma\_\{k\}\(\{\\mathcal\{U\}\}\_\{j\}\\mu\_\{1\}\)=\\sigma\_\{k\}\(\{\\mathcal\{U\}\}\_\{j\}\\mu P\_\{\{\\mathcal\{K\}\}\}\)\\leq\\sigma\_\{k\}\(\{\\mathcal\{U\}\}\_\{j\}\\mu\),\\qquad j=1,2σk​\(𝒰3​μ1\)=σk​\(P𝒦​𝒰3​μ\)≤σk​\(𝒰3​μ\)\\displaystyle\\sigma\_\{k\}\(\{\\mathcal\{U\}\}\_\{3\}\\mu\_\{1\}\)=\\sigma\_\{k\}\(P\_\{\{\\mathcal\{K\}\}\}\{\\mathcal\{U\}\}\_\{3\}\\mu\)\\leq\\sigma\_\{k\}\(\{\\mathcal\{U\}\}\_\{3\}\\mu\)Therefore,‖𝒰k​μ‖∗≥‖𝒰k​μ1‖∗\\\|\{\\mathcal\{U\}\}\_\{k\}\\mu\\\|\_\{\*\}\\geq\\\|\{\\mathcal\{U\}\}\_\{k\}\\mu\_\{1\}\\\|\_\{\*\}fork=1,2,3k=1,2,3\.

### C\.2Proof of Theorem[5\.5](https://arxiv.org/html/2606.04176#S5.Thmtheorem5)and[5\.6](https://arxiv.org/html/2606.04176#S5.Thmtheorem6)

###### Proof of Theorem[5\.5](https://arxiv.org/html/2606.04176#S5.Thmtheorem5)\.

From the basic inequalityQ​\(μ^\)≤Q​\(μ∗\)Q\(\\hat\{\\mu\}\)\\leq Q\(\\mu^\{\*\}\), we have

1m1​m2​∑\(i,j\)T​\(i,j\)​‖μ^​\(i,j,⋅\)−μ∗​\(i,j,⋅\)‖ℋ2\+λ1​‖μ^\(1\)‖∗\+λ2​‖μ^\(2\)‖∗\+λ3​‖μ^\(3\)‖∗≤\+λ1​‖μ\(1\)∗‖∗\+λ2​‖μ\(2\)∗‖∗\+λ3​‖μ\(3\)∗‖∗\+2m1​m2​\|∑\(i,j\)T​\(i,j\)​⟨μ^​\(i,j,⋅\)−μ∗​\(i,j,⋅\),1Ni​j​∑l=1Ni​jK​\(Xl\(i​j\),⋅\)−μ∗​\(i,j,⋅\)⟩ℋ\|\.\\frac\{1\}\{m\_\{1\}m\_\{2\}\}\\sum\_\{\(i,j\)\}T\(i,j\)\\left\\\|\\hat\{\\mu\}\(i,j,\\cdot\)\-\\mu^\{\*\}\(i,j,\\cdot\)\\right\\\|^\{2\}\_\{\\mathcal\{H\}\}\+\\lambda\_\{1\}\\\|\\hat\{\\mu\}\_\{\(1\)\}\\\|\_\{\*\}\+\\lambda\_\{2\}\\\|\\hat\{\\mu\}\_\{\(2\)\}\\\|\_\{\*\}\+\\lambda\_\{3\}\\\|\\hat\{\\mu\}\_\{\(3\)\}\\\|\_\{\*\}\\leq\\\\ \+\\lambda\_\{1\}\\\|\{\\mu\}^\{\*\}\_\{\(1\)\}\\\|\_\{\*\}\+\\lambda\_\{2\}\\\|\{\\mu\}^\{\*\}\_\{\(2\)\}\\\|\_\{\*\}\+\\lambda\_\{3\}\\\|\{\\mu\}^\{\*\}\_\{\(3\)\}\\\|\_\{\*\}\\\\ \+\\frac\{2\}\{m\_\{1\}m\_\{2\}\}\\left\|\\sum\_\{\(i,j\)\}T\(i,j\)\\left\\langle\\hat\{\\mu\}\(i,j,\\cdot\)\-\\mu^\{\*\}\(i,j,\\cdot\),\\frac\{1\}\{N\_\{ij\}\}\\sum\_\{l=1\}^\{N\_\{ij\}\}K\(X^\{\(ij\)\}\_\{l\},\\cdot\)\-\\mu^\{\*\}\(i,j,\\cdot\)\\right\\rangle\_\{\\mathcal\{H\}\}\\right\|\.
DenoteΔ=μ^−μ∗\\Delta=\\hat\{\\mu\}\-\\mu^\{\*\}\. By Lemma[C\.2](https://arxiv.org/html/2606.04176#A3.Thmtheorem2), we have

πL​1m1​m2​∑i,j‖Δ​\(i,j,⋅\)‖ℋ2≲1m1​m2​∑i,jTi,j​‖Δ​\(i,j,⋅\)‖ℋ2\+\(2​α\)2​πU​log⁡dπL2​m1​m2\\displaystyle\\pi\_\{L\}\\frac\{1\}\{m\_\{1\}m\_\{2\}\}\\sum\_\{i,j\}\\\|\\Delta\(i,j,\\cdot\)\\\|\_\{\\mathcal\{H\}\}^\{2\}\\lesssim\\frac\{1\}\{m\_\{1\}m\_\{2\}\}\\sum\_\{i,j\}T\_\{i,j\}\\\|\\Delta\(i,j,\\cdot\)\\\|\_\{\\mathcal\{H\}\}^\{2\}\+\(2\\alpha\)^\{2\}\\sqrt\{\\frac\{\\pi\_\{U\}\\log d\}\{\\pi^\{2\}\_\{L\}m\_\{1\}m\_\{2\}\}\}\+\(2​α\)​πU​log⁡dm​1m1​m2​\(λ1λ1\+λ2​\[‖μ^\(1\)‖\+‖μ\(1\)∗‖\]\+λ2λ1\+λ2​\[‖μ^\(2\)‖\+‖μ\(2\)∗‖\]\)\\displaystyle\+\(2\\alpha\)\\sqrt\{\\frac\{\\pi\_\{U\}\\log d\}\{m\}\}\\frac\{1\}\{\\sqrt\{m\_\{1\}m\_\{2\}\}\}\\left\(\\frac\{\\lambda\_\{1\}\}\{\\lambda\_\{1\}\+\\lambda\_\{2\}\}\\left\[\\\|\\hat\{\\mu\}\_\{\(1\)\}\\\|\+\\\|\{\\mu\}^\{\*\}\_\{\(1\)\}\\\|\\right\]\+\\frac\{\\lambda\_\{2\}\}\{\\lambda\_\{1\}\+\\lambda\_\{2\}\}\\left\[\\\|\\hat\{\\mu\}\_\{\(2\)\}\\\|\+\\\|\{\\mu\}^\{\*\}\_\{\(2\)\}\\\|\\right\]\\right\)≲\|∑i,j⟨Δ​\(i,j,⋅\),T​\(i,j\)​ϵ​\(i,j,⋅\)⟩ℋ\|\+λ1​‖μ\(1\)∗‖∗\+λ2​‖μ\(2\)∗‖∗\+λ3​‖μ\(3\)∗‖∗−λ2​‖μ^\(2\)‖∗−λ1​‖μ^\(1\)‖∗−λ3​‖μ^\(3\)‖∗\\displaystyle\\lesssim\\left\|\\sum\_\{i,j\}\\langle\\Delta\(i,j,\\cdot\),T\(i,j\)\\epsilon\(i,j,\\cdot\)\\rangle\_\{\\mathcal\{H\}\}\\right\|\+\\lambda\_\{1\}\\\|\{\\mu\}^\{\*\}\_\{\(1\)\}\\\|\_\{\*\}\+\\lambda\_\{2\}\\\|\{\\mu\}^\{\*\}\_\{\(2\)\}\\\|\_\{\*\}\+\\lambda\_\{3\}\\\|\{\\mu\}^\{\*\}\_\{\(3\)\}\\\|\_\{\*\}\-\\lambda\_\{2\}\\\|\\hat\{\\mu\}\_\{\(2\)\}\\\|\_\{\*\}\-\\lambda\_\{1\}\\\|\\hat\{\\mu\}\_\{\(1\)\}\\\|\_\{\*\}\-\\lambda\_\{3\}\\\|\\hat\{\\mu\}\_\{\(3\)\}\\\|\_\{\*\}\+\(2​α\)2​πU​log⁡dm1​m2\+\(2​α\)​πU​log⁡dm​1m1​m2​\(λ1λ1\+λ2​\[‖μ^\(1\)‖\+‖μ\(1\)∗‖\]\+λ2λ1\+λ2​\[‖μ^\(2\)‖\+‖μ\(2\)∗‖\]\)\.\\displaystyle\+\(2\\alpha\)^\{2\}\\sqrt\{\\frac\{\\pi\_\{U\}\\log d\}\{m\_\{1\}m\_\{2\}\}\}\+\(2\\alpha\)\\sqrt\{\\frac\{\\pi\_\{U\}\\log d\}\{m\}\}\\frac\{1\}\{\\sqrt\{m\_\{1\}m\_\{2\}\}\}\\left\(\\frac\{\\lambda\_\{1\}\}\{\\lambda\_\{1\}\+\\lambda\_\{2\}\}\\left\[\\\|\\hat\{\\mu\}\_\{\(1\)\}\\\|\+\\\|\{\\mu\}^\{\*\}\_\{\(1\)\}\\\|\\right\]\+\\frac\{\\lambda\_\{2\}\}\{\\lambda\_\{1\}\+\\lambda\_\{2\}\}\\left\[\\\|\\hat\{\\mu\}\_\{\(2\)\}\\\|\+\\\|\{\\mu\}^\{\*\}\_\{\(2\)\}\\\|\\right\]\\right\)\.Apply the bound in Lemma[C\.1](https://arxiv.org/html/2606.04176#A3.Thmtheorem1), the above inequality can be further bounded by

≲ΔN​πU​log⁡m1​m2m​1m1​m2​\(λ1λ1\+λ2​\[‖μ^\(1\)‖\+‖μ\(1\)∗‖\]\+λ2λ1\+λ2​\[‖μ^\(2\)‖\+‖μ\(2\)∗‖\]\)\\displaystyle\\lesssim\\Delta\_\{N\}\\sqrt\{\\frac\{\\pi\_\{U\}\\log m\_\{1\}m\_\{2\}\}\{m\}\}\\frac\{1\}\{\\sqrt\{m\_\{1\}m\_\{2\}\}\}\\left\(\\frac\{\\lambda\_\{1\}\}\{\\lambda\_\{1\}\+\\lambda\_\{2\}\}\\left\[\\\|\\hat\{\\mu\}\_\{\(1\)\}\\\|\+\\\|\{\\mu\}^\{\*\}\_\{\(1\)\}\\\|\\right\]\+\\frac\{\\lambda\_\{2\}\}\{\\lambda\_\{1\}\+\\lambda\_\{2\}\}\\left\[\\\|\\hat\{\\mu\}\_\{\(2\)\}\\\|\+\\\|\{\\mu\}^\{\*\}\_\{\(2\)\}\\\|\\right\]\\right\)\+λ1​‖μ\(1\)∗‖∗\+λ2​‖μ\(2\)∗‖∗\+λ3​‖μ\(3\)∗‖∗−λ2​‖μ^\(2\)‖∗−λ1​‖μ^\(1\)‖∗−λ3​‖μ^\(3\)‖∗\\displaystyle\+\\lambda\_\{1\}\\\|\{\\mu\}^\{\*\}\_\{\(1\)\}\\\|\_\{\*\}\+\\lambda\_\{2\}\\\|\{\\mu\}^\{\*\}\_\{\(2\)\}\\\|\_\{\*\}\+\\lambda\_\{3\}\\\|\{\\mu\}^\{\*\}\_\{\(3\)\}\\\|\_\{\*\}\-\\lambda\_\{2\}\\\|\\hat\{\\mu\}\_\{\(2\)\}\\\|\_\{\*\}\-\\lambda\_\{1\}\\\|\\hat\{\\mu\}\_\{\(1\)\}\\\|\_\{\*\}\-\\lambda\_\{3\}\\\|\\hat\{\\mu\}\_\{\(3\)\}\\\|\_\{\*\}\+\(2​α\)2​πU​log⁡dm1​m2\+\(2​α\)​πU​log⁡dm​1m1​m2​\(λ1λ1\+λ2​\[‖μ^\(1\)‖\+‖μ\(1\)∗‖\]\+λ2λ1\+λ2​\[‖μ^\(2\)‖\+‖μ\(2\)∗‖\]\)\.\\displaystyle\+\(2\\alpha\)^\{2\}\\sqrt\{\\frac\{\\pi\_\{U\}\\log d\}\{m\_\{1\}m\_\{2\}\}\}\+\(2\\alpha\)\\sqrt\{\\frac\{\\pi\_\{U\}\\log d\}\{m\}\}\\frac\{1\}\{\\sqrt\{m\_\{1\}m\_\{2\}\}\}\\left\(\\frac\{\\lambda\_\{1\}\}\{\\lambda\_\{1\}\+\\lambda\_\{2\}\}\\left\[\\\|\\hat\{\\mu\}\_\{\(1\)\}\\\|\+\\\|\{\\mu\}^\{\*\}\_\{\(1\)\}\\\|\\right\]\+\\frac\{\\lambda\_\{2\}\}\{\\lambda\_\{1\}\+\\lambda\_\{2\}\}\\left\[\\\|\\hat\{\\mu\}\_\{\(2\)\}\\\|\+\\\|\{\\mu\}^\{\*\}\_\{\(2\)\}\\\|\\right\]\\right\)\.By taking

\(λ1\+λ2\)−1=𝒪\(\[1m1​m2max\{\(M​πUΔNlog\(m1m2\)\),\(α\)πU​M​log⁡d\}\]−1\),\\displaystyle\(\\lambda\_\{1\}\+\\lambda\_\{2\}\)^\{\-1\}=\\mathop\{\}\\mathopen\{\}\\mathcal\{O\}\\mathopen\{\}\\left\(\\left\[\\frac\{1\}\{m\_\{1\}m\_\{2\}\}\\max\\left\\\{\\left\(\\sqrt\{M\\pi\_\{U\}\}\\Delta\_\{N\}\\log\(m\_\{1\}m\_\{2\}\)\\right\),\(\\alpha\)\\sqrt\{\{\\pi\_\{U\}M\\log d\}\}\\right\\\}\\right\]^\{\-1\}\\right\),we have

1m1​m2​∑i,j‖Δ​\(i,j,⋅\)‖ℋ2≲\(2​α\)2​πU​log⁡dπL2​m1​m2\+1πL​\(λ1​‖μ\(1\)∗‖∗\+λ2​‖μ\(2\)∗‖∗\+λ3​‖μ\(3\)∗‖∗\)\.\\displaystyle\\frac\{1\}\{m\_\{1\}m\_\{2\}\}\\sum\_\{i,j\}\\\|\\Delta\(i,j,\\cdot\)\\\|\_\{\\mathcal\{H\}\}^\{2\}\\lesssim\(2\\alpha\)^\{2\}\\sqrt\{\\frac\{\\pi\_\{U\}\\log d\}\{\\pi\_\{L\}^\{2\}m\_\{1\}m\_\{2\}\}\}\+\\frac\{1\}\{\\pi\_\{L\}\}\\left\(\\lambda\_\{1\}\\\|\{\\mu\}^\{\*\}\_\{\(1\)\}\\\|\_\{\*\}\+\\lambda\_\{2\}\\\|\{\\mu\}^\{\*\}\_\{\(2\)\}\\\|\_\{\*\}\+\\lambda\_\{3\}\\\|\{\\mu\}^\{\*\}\_\{\(3\)\}\\\|\_\{\*\}\\right\)\.
∎

###### Proof of Theorem[5\.6](https://arxiv.org/html/2606.04176#S5.Thmtheorem6)\.

For a matrixB∈ℝk1×k2B\\in\\mathbb\{R\}^\{k\_\{1\}\\times k\_\{2\}\}, takeUBU\_\{B\}andVBV\_\{B\}as the left and right singular matrixces ofBBaccordingly\. DenoteSU​\(B\)S\_\{U\}\(B\)andSV​\(B\)S\_\{V\}\(B\)as the linear spaces spanned by column vectors ofUBU\_\{B\}andVBV\_\{B\}respectively\. DenoteSU⟂⟂​\(B\)S^\{\\perp\\\!\\\!\\\!\\\!\\perp\}\_\{U\}\(B\)andSV⟂⟂​\(B\)S^\{\\perp\\\!\\\!\\\!\\\!\\perp\}\_\{V\}\(B\)as the orthogonal components ofSU​\(B\)S\_\{U\}\(B\)andSV​\(B\)S\_\{V\}\(B\)respectively\. Define

PB⟂⟂​B′=𝐏𝐒𝐔⟂⟂​\(𝐁\)​B′​𝐏𝐒𝐕⟂⟂​\(𝐁\),PB​B′=B′−PB⟂⟂​B′,\\displaystyle P^\{\\perp\\\!\\\!\\\!\\\!\\perp\}\_\{B\}B^\{\\prime\}=\\mathbf\{\{P\}\_\{S^\{\\perp\\\!\\\!\\\!\\\!\\perp\}\_\{U\}\(B\)\}\}B^\{\\prime\}\\mathbf\{\{P\}\_\{S^\{\\perp\\\!\\\!\\\!\\\!\\perp\}\_\{V\}\(B\)\}\},\\qquad P\_\{B\}B^\{\\prime\}=B^\{\\prime\}\-P^\{\\perp\\\!\\\!\\\!\\\!\\perp\}\_\{B\}B^\{\\prime\},\(19\)where𝐏𝐒\\mathbf\{\{P\}\_\{S\}\}is the projection matrix on the linear spaceSS\. And we have the following inequality

‖μ^\(1\)‖∗=‖Pμ\(1\)∗​\(μ^\(1\)−μ\(1\)∗\)\+Pμ\(1\)∗⟂⟂​\(μ^\(1\)−μ\(1\)∗\)\+μ\(1\)∗‖∗\\displaystyle\\left\\\|\\hat\{\\mu\}\_\{\(1\)\}\\right\\\|\_\{\*\}=\\left\\\|P\_\{\\mu^\{\*\}\_\{\(1\)\}\}\(\\hat\{\\mu\}\_\{\(1\)\}\-\\mu^\{\*\}\_\{\(1\)\}\)\+P^\{\\perp\\\!\\\!\\\!\\\!\\perp\}\_\{\\mu^\{\*\}\_\{\(1\)\}\}\(\\hat\{\\mu\}\_\{\(1\)\}\-\\mu^\{\*\}\_\{\(1\)\}\)\+\\mu^\{\*\}\_\{\(1\)\}\\right\\\|\_\{\*\}≥‖Pμ\(1\)∗⟂⟂​\(μ^\(1\)−μ\(1\)∗\)\+μ\(1\)∗‖∗−‖Pμ\(1\)∗​\(μ^\(1\)−μ\(1\)∗\)‖∗\\displaystyle\\geq\\left\\\|P^\{\\perp\\\!\\\!\\\!\\\!\\perp\}\_\{\\mu^\{\*\}\_\{\(1\)\}\}\(\\hat\{\\mu\}\_\{\(1\)\}\-\\mu^\{\*\}\_\{\(1\)\}\)\+\\mu^\{\*\}\_\{\(1\)\}\\right\\\|\_\{\*\}\-\\left\\\|P\_\{\\mu^\{\*\}\_\{\(1\)\}\}\(\\hat\{\\mu\}\_\{\(1\)\}\-\\mu^\{\*\}\_\{\(1\)\}\)\\right\\\|\_\{\*\}=‖Pμ\(1\)∗⟂⟂​\(μ^\(1\)−μ\(1\)∗\)‖∗\+‖μ\(1\)∗‖∗−‖Pμ\(1\)∗​\(μ^\(1\)−μ\(1\)∗\)‖∗\.\\displaystyle=\\left\\\|P^\{\\perp\\\!\\\!\\\!\\\!\\perp\}\_\{\\mu^\{\*\}\_\{\(1\)\}\}\(\\hat\{\\mu\}\_\{\(1\)\}\-\\mu^\{\*\}\_\{\(1\)\}\)\\right\\\|\_\{\*\}\+\\left\\\|\\mu^\{\*\}\_\{\(1\)\}\\right\\\|\_\{\*\}\-\\left\\\|P\_\{\\mu^\{\*\}\_\{\(1\)\}\}\(\\hat\{\\mu\}\_\{\(1\)\}\-\\mu^\{\*\}\_\{\(1\)\}\)\\right\\\|\_\{\*\}\.And

‖μ\(1\)∗‖∗−‖μ^\(1\)‖∗≤‖Pμ\(1\)∗​\(μ^\(1\)−μ\(1\)∗\)‖∗−‖Pμ\(1\)∗⟂⟂​\(μ^\(1\)−μ\(1\)∗\)‖∗\.\\displaystyle\\left\\\|\\mu^\{\*\}\_\{\(1\)\}\\right\\\|\_\{\*\}\-\\left\\\|\\hat\{\\mu\}\_\{\(1\)\}\\right\\\|\_\{\*\}\\leq\\left\\\|P\_\{\\mu^\{\*\}\_\{\(1\)\}\}\(\\hat\{\\mu\}\_\{\(1\)\}\-\\mu^\{\*\}\_\{\(1\)\}\)\\right\\\|\_\{\*\}\-\\left\\\|P^\{\\perp\\\!\\\!\\\!\\\!\\perp\}\_\{\\mu^\{\*\}\_\{\(1\)\}\}\(\\hat\{\\mu\}\_\{\(1\)\}\-\\mu^\{\*\}\_\{\(1\)\}\)\\right\\\|\_\{\*\}\.\(20\)Similar inequalities can be obtained for the second and the third mode unfolding\.

Next, by Lemma[C\.2](https://arxiv.org/html/2606.04176#A3.Thmtheorem2), we have

πL​1m1​m2​∑i,j‖Δ​\(i,j,⋅\)‖ℋ2≤𝔼​1m1​m2​∑i,jTi,j​‖Δ​\(i,j,⋅\)‖ℋ2≲1m1​m2​∑i,jTi,j​‖Δ​\(i,j,⋅\)‖ℋ2\+\(2​α\)2​πU​log⁡dπL2​m1​m2\\displaystyle\\pi\_\{L\}\\frac\{1\}\{m\_\{1\}m\_\{2\}\}\\sum\_\{i,j\}\\\|\\Delta\(i,j,\\cdot\)\\\|\_\{\\mathcal\{H\}\}^\{2\}\\leq\{\\mathbb\{E\}\}\\frac\{1\}\{m\_\{1\}m\_\{2\}\}\\sum\_\{i,j\}T\_\{i,j\}\\\|\\Delta\(i,j,\\cdot\)\\\|\_\{\\mathcal\{H\}\}^\{2\}\\lesssim\\frac\{1\}\{m\_\{1\}m\_\{2\}\}\\sum\_\{i,j\}T\_\{i,j\}\\\|\\Delta\(i,j,\\cdot\)\\\|\_\{\\mathcal\{H\}\}^\{2\}\+\(2\\alpha\)^\{2\}\\sqrt\{\\frac\{\\pi\_\{U\}\\log d\}\{\\pi^\{2\}\_\{L\}m\_\{1\}m\_\{2\}\}\}\+απU​log⁡dm1m1​m2\(λ1λ1\+λ2\[∥Pμ\(1\)∗⟂⟂\(μ^\(1\)−μ\(1\)∗\)∥∗\+∥Pμ\(1\)∗\(μ^\(1\)−μ\(1\)∗\)∥∗\]\\displaystyle\+\\alpha\\sqrt\{\\frac\{\\pi\_\{U\}\\log d\}\{m\}\}\\frac\{1\}\{\\sqrt\{m\_\{1\}m\_\{2\}\}\}\\left\(\\frac\{\\lambda\_\{1\}\}\{\\lambda\_\{1\}\+\\lambda\_\{2\}\}\\left\[\\left\\\|P^\{\\perp\\\!\\\!\\\!\\\!\\perp\}\_\{\\mu^\{\*\}\_\{\(1\)\}\}\(\\hat\{\\mu\}\_\{\(1\)\}\-\\mu^\{\*\}\_\{\(1\)\}\)\\right\\\|\_\{\*\}\+\\left\\\|P\_\{\\mu^\{\*\}\_\{\(1\)\}\}\(\\hat\{\\mu\}\_\{\(1\)\}\-\\mu^\{\*\}\_\{\(1\)\}\)\\right\\\|\_\{\*\}\\right\]\\right\.\+λ2λ1\+λ2\[∥Pμ\(2\)∗⟂⟂\(μ^\(2\)−μ\(2\)∗\)∥∗\+∥Pμ\(2\)∗\(μ^\(2\)−μ\(2\)∗\)∥∗\]\)\\displaystyle\\left\.\+\\frac\{\\lambda\_\{2\}\}\{\\lambda\_\{1\}\+\\lambda\_\{2\}\}\\left\[\\left\\\|P^\{\\perp\\\!\\\!\\\!\\\!\\perp\}\_\{\\mu^\{\*\}\_\{\(2\)\}\}\(\\hat\{\\mu\}\_\{\(2\)\}\-\\mu^\{\*\}\_\{\(2\)\}\)\\right\\\|\_\{\*\}\+\\left\\\|P\_\{\\mu^\{\*\}\_\{\(2\)\}\}\(\\hat\{\\mu\}\_\{\(2\)\}\-\\mu^\{\*\}\_\{\(2\)\}\)\\right\\\|\_\{\*\}\\right\]\\right\)By the basic inquality \([C\.2](https://arxiv.org/html/2606.04176#A3.SS2)\):≤\|∑i,j⟨Δ​\(i,j,⋅\),T​\(i,j\)​ϵ​\(i,j,⋅\)⟩ℋ\|m1​m2\+λ1​‖μ\(1\)∗‖∗\+λ2​‖μ\(2\)∗‖∗\+λ3​‖μ\(3\)∗‖∗−λ2​‖μ^\(2\)‖∗−λ1​‖μ^\(1\)‖∗−λ3​‖μ^\(3\)‖∗\\displaystyle\\leq\\frac\{\\left\|\\sum\_\{i,j\}\\langle\\Delta\(i,j,\\cdot\),T\(i,j\)\\epsilon\(i,j,\\cdot\)\\rangle\_\{\\mathcal\{H\}\}\\right\|\}\{m\_\{1\}m\_\{2\}\}\+\\lambda\_\{1\}\\\|\{\\mu\}^\{\*\}\_\{\(1\)\}\\\|\_\{\*\}\+\\lambda\_\{2\}\\\|\{\\mu\}^\{\*\}\_\{\(2\)\}\\\|\_\{\*\}\+\\lambda\_\{3\}\\\|\{\\mu\}^\{\*\}\_\{\(3\)\}\\\|\_\{\*\}\-\\lambda\_\{2\}\\\|\\hat\{\\mu\}\_\{\(2\)\}\\\|\_\{\*\}\-\\lambda\_\{1\}\\\|\\hat\{\\mu\}\_\{\(1\)\}\\\|\_\{\*\}\-\\lambda\_\{3\}\\\|\\hat\{\\mu\}\_\{\(3\)\}\\\|\_\{\*\}\+\(2α\)2πU​log⁡dm1​m2\+\(2α\)πU​log⁡dm1m1​m2\(λ1λ1\+λ2\[∥Pμ\(1\)∗⟂⟂\(μ^\(1\)−μ\(1\)∗\)∥∗\+∥Pμ\(1\)∗\(μ^\(1\)−μ\(1\)∗\)∥∗\]\\displaystyle\+\(2\\alpha\)^\{2\}\\sqrt\{\\frac\{\\pi\_\{U\}\\log d\}\{m\_\{1\}m\_\{2\}\}\}\+\(2\\alpha\)\\sqrt\{\\frac\{\\pi\_\{U\}\\log d\}\{m\}\}\\frac\{1\}\{\\sqrt\{m\_\{1\}m\_\{2\}\}\}\\left\(\\frac\{\\lambda\_\{1\}\}\{\\lambda\_\{1\}\+\\lambda\_\{2\}\}\\left\[\\left\\\|P^\{\\perp\\\!\\\!\\\!\\\!\\perp\}\_\{\\mu^\{\*\}\_\{\(1\)\}\}\(\\hat\{\\mu\}\_\{\(1\)\}\-\\mu^\{\*\}\_\{\(1\)\}\)\\right\\\|\_\{\*\}\+\\left\\\|P\_\{\\mu^\{\*\}\_\{\(1\)\}\}\(\\hat\{\\mu\}\_\{\(1\)\}\-\\mu^\{\*\}\_\{\(1\)\}\)\\right\\\|\_\{\*\}\\right\]\\right\.\+λ2λ1\+λ2\[∥Pμ\(1\)∗⟂⟂\(μ^\(2\)−μ\(2\)∗\)∥∗\+∥Pμ\(2\)∗\(μ^\(2\)−μ\(2\)∗\)∥∗\]\)\\displaystyle\\left\.\+\\frac\{\\lambda\_\{2\}\}\{\\lambda\_\{1\}\+\\lambda\_\{2\}\}\\left\[\\left\\\|P^\{\\perp\\\!\\\!\\\!\\\!\\perp\}\_\{\\mu^\{\*\}\_\{\(1\)\}\}\(\\hat\{\\mu\}\_\{\(2\)\}\-\\mu^\{\*\}\_\{\(2\)\}\)\\right\\\|\_\{\*\}\+\\left\\\|P\_\{\\mu^\{\*\}\_\{\(2\)\}\}\(\\hat\{\\mu\}\_\{\(2\)\}\-\\mu^\{\*\}\_\{\(2\)\}\)\\right\\\|\_\{\*\}\\right\]\\right\)By Lemma[C\.1](https://arxiv.org/html/2606.04176#A3.Thmtheorem1):≲2m1​m2​\(M​πU​ΔN​log⁡\(m1​m2\)\)​\(λ1λ1\+λ2​‖Δ\(1\)‖∗\+λ2λ1\+λ2​‖Δ\(2\)‖∗\)\\displaystyle\\lesssim\\frac\{2\}\{m\_\{1\}m\_\{2\}\}\\left\(\\sqrt\{M\\pi\_\{U\}\}\\Delta\_\{N\}\\log\(m\_\{1\}m\_\{2\}\)\\right\)\\left\(\\frac\{\\lambda\_\{1\}\}\{\\lambda\_\{1\}\+\\lambda\_\{2\}\}\\\|\\Delta\_\{\(1\)\}\\\|\_\{\*\}\+\\frac\{\\lambda\_\{2\}\}\{\\lambda\_\{1\}\+\\lambda\_\{2\}\}\\\|\\Delta\_\{\(2\)\}\\\|\_\{\*\}\\right\)\+λ1​\(Pμ\(1\)∗​‖Δ\(1\)‖∗−Pμ\(1\)∗⟂⟂​‖Δ\(1\)‖∗\)\+λ2​\(Pμ\(2\)∗​‖Δ\(2\)‖∗−Pμ\(2\)∗⟂⟂​‖Δ\(2\)‖∗\)\+λ3​\(Pμ\(3\)∗​‖Δ3\)‖∗−Pμ\(3\)∗⟂⟂​‖Δ\(3\)‖∗\)\\displaystyle\+\\lambda\_\{1\}\\left\(P\_\{\\mu^\{\*\}\_\{\(1\)\}\}\\\|\\Delta\_\{\(1\)\}\\\|\_\{\*\}\-P^\{\\perp\\\!\\\!\\\!\\\!\\perp\}\_\{\\mu^\{\*\}\_\{\(1\)\}\}\\\|\\Delta\_\{\(1\)\}\\\|\_\{\*\}\\right\)\+\\lambda\_\{2\}\\left\(P\_\{\\mu^\{\*\}\_\{\(2\)\}\}\\\|\\Delta\_\{\(2\)\}\\\|\_\{\*\}\-P^\{\\perp\\\!\\\!\\\!\\\!\\perp\}\_\{\\mu^\{\*\}\_\{\(2\)\}\}\\\|\\Delta\_\{\(2\)\}\\\|\_\{\*\}\\right\)\+\\lambda\_\{3\}\\left\(P\_\{\\mu^\{\*\}\_\{\(3\)\}\}\\\|\\Delta\_\{3\)\}\\\|\_\{\*\}\-P^\{\\perp\\\!\\\!\\\!\\\!\\perp\}\_\{\\mu^\{\*\}\_\{\(3\)\}\}\\\|\\Delta\_\{\(3\)\}\\\|\_\{\*\}\\right\)\+\(2α\)2πU​log⁡dm1​m2\+\(2α\)πU​log⁡dm\)1m1​m2\(λ1λ1\+λ2\[∥Pμ\(1\)∗⟂⟂\(μ^\(1\)−μ\(1\)∗\)∥∗\+∥Pμ\(1\)∗\(μ^\(1\)−μ\(1\)∗\)∥∗\]\\displaystyle\+\(2\\alpha\)^\{2\}\\sqrt\{\\frac\{\\pi\_\{U\}\\log d\}\{m\_\{1\}m\_\{2\}\}\}\+\(2\\alpha\)\\sqrt\{\\frac\{\\pi\_\{U\}\\log d\}\{m\)\}\}\\frac\{1\}\{\\sqrt\{m\_\{1\}m\_\{2\}\}\}\\left\(\\frac\{\\lambda\_\{1\}\}\{\\lambda\_\{1\}\+\\lambda\_\{2\}\}\\left\[\\left\\\|P^\{\\perp\\\!\\\!\\\!\\\!\\perp\}\_\{\\mu^\{\*\}\_\{\(1\)\}\}\(\\hat\{\\mu\}\_\{\(1\)\}\-\\mu^\{\*\}\_\{\(1\)\}\)\\right\\\|\_\{\*\}\+\\left\\\|P\_\{\\mu^\{\*\}\_\{\(1\)\}\}\(\\hat\{\\mu\}\_\{\(1\)\}\-\\mu^\{\*\}\_\{\(1\)\}\)\\right\\\|\_\{\*\}\\right\]\\right\.\+λ2λ1\+λ2\[∥Pμ\(1\)∗⟂⟂\(μ^\(2\)−μ\(2\)∗\)∥∗\+∥Pμ\(2\)∗\(μ^\(2\)−μ\(2\)∗\)∥∗\]\)\\displaystyle\\left\.\+\\frac\{\\lambda\_\{2\}\}\{\\lambda\_\{1\}\+\\lambda\_\{2\}\}\\left\[\\left\\\|P^\{\\perp\\\!\\\!\\\!\\\!\\perp\}\_\{\\mu^\{\*\}\_\{\(1\)\}\}\(\\hat\{\\mu\}\_\{\(2\)\}\-\\mu^\{\*\}\_\{\(2\)\}\)\\right\\\|\_\{\*\}\+\\left\\\|P\_\{\\mu^\{\*\}\_\{\(2\)\}\}\(\\hat\{\\mu\}\_\{\(2\)\}\-\\mu^\{\*\}\_\{\(2\)\}\)\\right\\\|\_\{\*\}\\right\]\\right\)If

\(λ1\+λ2\)−1=𝒪\(\[max\{2m1​m2\(M​πUΔNlog\(m1m2\)\),\(2α\)πU​M​log⁡d​M\}\]−1\),\\displaystyle\(\\lambda\_\{1\}\+\\lambda\_\{2\}\)^\{\-1\}=\\mathop\{\}\\mathopen\{\}\\mathcal\{O\}\\mathopen\{\}\\left\(\\left\[\\max\\left\\\{\\frac\{2\}\{m\_\{1\}m\_\{2\}\}\\left\(\\sqrt\{M\\pi\_\{U\}\}\\Delta\_\{N\}\\log\(m\_\{1\}m\_\{2\}\)\\right\),\(2\\alpha\)\\sqrt\{\{\\pi\_\{U\}M\\log dM\}\}\\right\\\}\\right\]^\{\-1\}\\right\),Then we have

πL​1m1​m2​‖Δ‖F2≲min⁡\{λ1​\(Pμ\(1\)∗​‖Δ\(1\)‖∗\)\+λ2​\(Pμ\(2\)∗​‖Δ\(2\)‖∗\),λ3​\(Pμ\(3\)∗​‖Δ3\)‖∗\)\}\+\(2​α\)2​πU​log⁡dm1​m2\\displaystyle\\pi\_\{L\}\\frac\{1\}\{m\_\{1\}m\_\{2\}\}\\\|\\Delta\\\|^\{2\}\_\{F\}\\lesssim\\min\\left\\\{\\lambda\_\{1\}\\left\(P\_\{\\mu^\{\*\}\_\{\(1\)\}\}\\\|\\Delta\_\{\(1\)\}\\\|\_\{\*\}\\right\)\+\\lambda\_\{2\}\\left\(P\_\{\\mu^\{\*\}\_\{\(2\)\}\}\\\|\\Delta\_\{\(2\)\}\\\|\_\{\*\}\\right\),\\lambda\_\{3\}\\left\(P\_\{\\mu^\{\*\}\_\{\(3\)\}\}\\\|\\Delta\_\{3\)\}\\\|\_\{\*\}\\right\)\\right\\\}\+\(2\\alpha\)^\{2\}\\sqrt\{\\frac\{\\pi\_\{U\}\\log d\}\{m\_\{1\}m\_\{2\}\}\}≲min⁡\{λ1​r1​‖Δ‖F\+λ2​r2​‖Δ‖F\+λ3​r3​‖Δ‖F\}\+\(2​α\)2​πU​log⁡dm1​m2\\displaystyle\\lesssim\\min\\\{\\lambda\_\{1\}\\sqrt\{r\_\{1\}\}\\\|\\Delta\\\|\_\{F\}\+\\lambda\_\{2\}\\sqrt\{r\_\{2\}\}\\\|\\Delta\\\|\_\{F\}\+\\lambda\_\{3\}\\sqrt\{r\_\{3\}\}\\\|\\Delta\\\|\_\{F\}\\\}\+\(2\\alpha\)^\{2\}\\sqrt\{\\frac\{\\pi\_\{U\}\\log d\}\{m\_\{1\}m\_\{2\}\}\}1m1​m2​‖Δ‖F2≲min⁡\{λ12​m1​m2​r1πL2\+λ22​m1​m2​r2πL2,λ32​m1​m2​r3πL2\}\+\(2​α\)2​πU​log⁡dπL2​m1​m2\\displaystyle\\frac\{1\}\{m\_\{1\}m\_\{2\}\}\\\|\\Delta\\\|^\{2\}\_\{F\}\\lesssim\\min\\left\\\{\\frac\{\\lambda\_\{1\}^\{2\}m\_\{1\}m\_\{2\}r\_\{1\}\}\{\\pi\_\{L\}^\{2\}\}\+\\frac\{\\lambda\_\{2\}^\{2\}m\_\{1\}m\_\{2\}r\_\{2\}\}\{\\pi\_\{L\}^\{2\}\},\\frac\{\\lambda\_\{3\}^\{2\}m\_\{1\}m\_\{2\}r\_\{3\}\}\{\\pi\_\{L\}^\{2\}\}\\right\\\}\+\(2\\alpha\)^\{2\}\\sqrt\{\\frac\{\\pi\_\{U\}\\log d\}\{\\pi\_\{L\}^\{2\}m\_\{1\}m\_\{2\}\}\}∎

### C\.3Auxiliary Lemmas

###### Lemma C\.1\.

TakeN=min\(i,j\):T​\(i,j\)=1⁡Ni​jN=\\min\_\{\(i,j\):T\(i,j\)=1\}N\_\{ij\}andϵ​\(i,j,⋅\)=1Ni​j​∑l=1Ni​jK​\(Xl\(i​j\),⋅\)−μ∗​\(i,j,⋅\)\\epsilon\(i,j,\\cdot\)=\\frac\{1\}\{N\_\{ij\}\}\\sum\_\{l=1\}^\{N\_\{ij\}\}K\(X^\{\(ij\)\}\_\{l\},\\cdot\)\-\\mu^\{\*\}\(i,j,\\cdot\)\. We have, for anyA∈ℋm1×m2A\\in\{\\mathcal\{H\}\}^\{m\_\{1\}\\times m\_\{2\}\}, with probability at least1−2/d−2​exp⁡\(−m1​m2​πL\)1\-2/d\-2\\exp\(\-m\_\{1\}m\_\{2\}\\pi\_\{L\}\),

\|∑i,j⟨A​\(i,j,⋅\),T​\(i,j\)​ϵ​\(i,j,⋅\)⟩ℋ\|\\displaystyle\\left\|\\sum\_\{i,j\}\\langle A\(i,j,\\cdot\),T\(i,j\)\\epsilon\(i,j,\\cdot\)\\rangle\_\{\\mathcal\{H\}\}\\right\|≲\\displaystyle\\lesssimmin⁡\{M​πU​ΔN​log⁡\(m1\)​‖A\(1\)‖∗,M​πU​ΔN​log⁡\(m2\)​‖A\(2\)‖∗,m1​m2​πU​ΔN​‖A\(3\)‖∗\}\\displaystyle\\min\\left\\\{\\sqrt\{M\\pi\_\{U\}\}\\Delta\_\{N\}\\log\(m\_\{1\}\)\\left\\\|A\_\{\(1\)\}\\right\\\|\_\{\*\},\\sqrt\{M\\pi\_\{U\}\}\\Delta\_\{N\}\\log\(m\_\{2\}\)\\left\\\|A\_\{\(2\)\}\\right\\\|\_\{\*\},\\sqrt\{m\_\{1\}m\_\{2\}\\pi\_\{U\}\}\\Delta\_\{N\}\\left\\\|A\_\{\(3\)\}\\right\\\|\_\{\*\}\\right\\\}whereΔN=log⁡m1​m2\+log⁡dN\\Delta\_\{N\}=\\sqrt\{\\frac\{\\log m\_\{1\}m\_\{2\}\+\\log d\}\{N\}\}\.

###### Proof of Lemma[C\.1](https://arxiv.org/html/2606.04176#A3.Thmtheorem1)\.

First, conditioned on the event∑i,jT​\(i,j\)=n\\sum\_\{i,j\}T\(i,j\)=n\. From Theorem 3\.4 inMuandet et al\. \[[2017](https://arxiv.org/html/2606.04176#bib.bib27)\], we have with probability at least1−δ1\-\\delta,

‖ϵ​\(i,j,⋅\)‖ℋ≤2Ni​j​CK\+log⁡\(2/δ\)2​Ni​j,\\displaystyle\\\|\\epsilon\(i,j,\\cdot\)\\\|\_\{\\mathcal\{H\}\}\\leq\{\\frac\{2\}\{N\_\{ij\}\}\}C\_\{K\}\+\{\\sqrt\{\\frac\{\\log\(2/\\delta\)\}\{2N\_\{ij\}\}\}\},whereCK=supx∈𝒳\|K​\(x,x\)\|C\_\{K\}=\\sup\_\{x\\in\\mathcal\{X\}\}\|K\(x,x\)\|\. By applying a union bound, we have, with probability at least1−1/d1\-1/d,

sup\(i,j\):T​\(i,j\)=1‖ϵ​\(i,j,⋅\)‖ℋ≲log⁡n\+log⁡dN≲log⁡m1​m2\+log⁡dN\.\\displaystyle\\sup\_\{\(i,j\):T\(i,j\)=1\}\\\|\\epsilon\(i,j,\\cdot\)\\\|\_\{\\mathcal\{H\}\}\\lesssim\\sqrt\{\\frac\{\\log n\+\\log d\}\{N\}\}\\lesssim\\sqrt\{\\frac\{\\log m\_\{1\}m\_\{2\}\+\\log d\}\{N\}\}\.By marginalizingnn, we have the above bound hold without conditioning\. Next, define the event

ℰ:\{sup\(i,j\):T​\(i,j\)=1‖ϵ​\(i,j,⋅\)∥ℋ≲log⁡n\+log⁡dN≲log⁡m1​m2\+log⁡dN\}\.\\displaystyle\\mathcal\{E\}:\\left\\\{\\sup\_\{\(i,j\):T\(i,j\)=1\}\\\|\\epsilon\(i,j,\\cdot\)\\\|\_\{\\mathcal\{H\}\}\\lesssim\\sqrt\{\\frac\{\\log n\+\\log d\}\{N\}\}\\lesssim\\sqrt\{\\frac\{\\log m\_\{1\}m\_\{2\}\+\\log d\}\{N\}\}\\right\\\}\.\(21\)
Note that

\|∑i,j⟨A​\(i,j,⋅\),T​\(i,j\)​ϵ​\(i,j,⋅\)⟩ℋ\|=⟨A\(1\),T∘1ϵ\(1\)⟩≤‖A\(1\)‖∗​‖T∘1ϵ\(1\)‖,\\displaystyle\\left\|\\sum\_\{i,j\}\\langle A\(i,j,\\cdot\),T\(i,j\)\\epsilon\(i,j,\\cdot\)\\rangle\_\{\\mathcal\{H\}\}\\right\|=\\left\\langle A\_\{\(1\)\},T\\circ\_\{1\}\\epsilon\_\{\(1\)\}\\right\\rangle\\leq\\\|A\_\{\(1\)\}\\\|\_\{\*\}\\\|T\\circ\_\{1\}\\epsilon\_\{\(1\)\}\\\|,
It remains to bound‖T∘1ϵ\(1\)‖\\\|T\\circ\_\{1\}\\epsilon\_\{\(1\)\}\\\|, whereTTandϵ\\epsilonare independent\.

TakeΔN=supi,j‖ϵ​\(i,j,⋅\)‖ℋ\\Delta\_\{N\}=\\sup\_\{i,j\}\\\|\\epsilon\(i,j,\\cdot\)\\\|\_\{\\mathcal\{H\}\}\.

‖T∘1ϵ\(1\)‖2\\displaystyle\\\|T\\circ\_\{1\}\\epsilon\_\{\(1\)\}\\\|^\{2\}≤sup‖a‖2≤1∑j=1m2‖∑i=1m1ai​Ti​j​ϵ​\(i,j,⋅\)‖ℋ2\\displaystyle\\leq\\sup\_\{\\\|a\\\|\_\{2\}\\leq 1\}\\sum\_\{j=1\}^\{m\_\{2\}\}\\left\\\|\\sum\_\{i=1\}^\{m\_\{1\}\}a\_\{i\}T\_\{ij\}\\epsilon\(i,j,\\cdot\)\\right\\\|\_\{\{\\mathcal\{H\}\}\}^\{2\}≤sup‖a‖2≤1∑j=1m2\{∑i=1m1ai2​Ti​j2‖ϵ​\(i,j,⋅\)∥ℋ2\+∑i≠i′ai​ai′​Ti​j​Ti′​j​⟨ϵ​\(i,j,⋅\),ϵ​\(i′,j,⋅\)⟩ℋ\}\\displaystyle\\leq\\sup\_\{\\\|a\\\|\_\{2\}\\leq 1\}\\sum\_\{j=1\}^\{m\_\{2\}\}\\left\\\{\\sum\_\{i=1\}^\{m\_\{1\}\}a^\{2\}\_\{i\}T^\{2\}\_\{ij\}\\\|\\epsilon\(i,j,\\cdot\)\\\|\_\{\\mathcal\{H\}\}^\{2\}\+\\sum\_\{i\\neq i^\{\\prime\}\}a\_\{i\}a\_\{i^\{\\prime\}\}T\_\{ij\}T\_\{i^\{\\prime\}j\}\\langle\\epsilon\(i,j,\\cdot\),\\epsilon\(i^\{\\prime\},j,\\cdot\)\\rangle\_\{\\mathcal\{H\}\}\\right\\\}≤‖M‖\\displaystyle\\leq\\\|M\\\|≤‖D‖\+‖M′‖,\\displaystyle\\leq\\\|D\\\|\+\\\|M^\{\\prime\}\\\|,whereM,D,M′∈ℝm1×m1M,D,M^\{\\prime\}\\in\\mathbb\{R\}^\{m\_\{1\}\\times m\_\{1\}\}are matrices such thatM​\(i,i′\)=∑j=1m2Ti​j​Ti′​j​⟨ϵ​\(i,j,⋅\),ϵ​\(i′,j,⋅\)⟩M\(i,i^\{\\prime\}\)=\\sum\_\{j=1\}^\{m\_\{2\}\}T\_\{ij\}T\_\{i^\{\\prime\}j\}\\langle\\epsilon\(i,j,\\cdot\),\\epsilon\(i^\{\\prime\},j,\\cdot\)\\rangle,DDis diagonal andD​\(i,i\)=M​\(i,i\)D\(i,i\)=M\(i,i\),M′=M−DM^\{\\prime\}=M\-D\.

First, we notice that for the diagonal matrixDD

‖D‖≤supi\(∑j=1m2\{Ti​j2‖ϵ​\(i,j,⋅\)∥ℋ2\}\)\\displaystyle\\\|D\\\|\\leq\\sup\_\{i\}\\left\(\\sum\_\{j=1\}^\{m\_\{2\}\}\\left\\\{T^\{2\}\_\{ij\}\\\|\\epsilon\(i,j,\\cdot\)\\\|\_\{\\mathcal\{H\}\}^\{2\}\\right\\\}\\right\)Conditioned on the eventℰ\{\\mathcal\{E\}\}, one can show that for everyii,∑j=1m2Ti​j2​‖ϵ​\(i,j,⋅\)‖ℋ2−𝔼​\{∑j=1m2Ti​j2‖ϵ​\(i,j,⋅\)∥ℋ2\}\\sum\_\{j=1\}^\{m\_\{2\}\}T^\{2\}\_\{ij\}\\\|\\epsilon\(i,j,\\cdot\)\\\|\_\{\\mathcal\{H\}\}^\{2\}\-\{\\mathbb\{E\}\}\\\{\\sum\_\{j=1\}^\{m\_\{2\}\}T^\{2\}\_\{ij\}\\\|\\epsilon\(i,j,\\cdot\)\\\|\_\{\\mathcal\{H\}\}^\{2\}\\\}is a subgaussian variable with variance bounded bym2​ΔN4​πUm\_\{2\}\\Delta\_\{N\}^\{4\}\\pi\_\{U\}\. Therefore, with probability at least1−1/d1\-1/d, we have

‖D‖≲m2​ΔN4​πU​log⁡n\+m2​πU​ΔN2≲m2​πU​ΔN2​log⁡d\.\\displaystyle\\\|D\\\|\\lesssim\\sqrt\{m\_\{2\}\\Delta\_\{N\}^\{4\}\\pi\_\{U\}\}\\log n\+m\_\{2\}\\pi\_\{U\}\\Delta\_\{N\}^\{2\}\\lesssim m\_\{2\}\\pi\_\{U\}\\Delta\_\{N\}^\{2\}\\log d\.\(22\)Still, we condition on the eventℰ\{\\mathcal\{E\}\},M′M^\{\\prime\}is a random matrix where every entry is mean zero and subgaussian\. For entry\(i,i′\)\(i,i^\{\\prime\}\), the variance of∑j=1m2Ti,j​Ti′​j​⟨ϵ​\(i,j,⋅\),ϵ​\(i′,j,⋅\)⟩ℋ\\sum\_\{j=1\}^\{m\_\{2\}\}T\_\{i,j\}T\_\{i^\{\\prime\}j\}\\langle\\epsilon\(i,j,\\cdot\),\\epsilon\(i^\{\\prime\},j,\\cdot\)\\rangle\_\{\\mathcal\{H\}\}is bounded bym2​πU2​ΔN4m\_\{2\}\\pi\_\{U\}^\{2\}\\Delta^\{4\}\_\{N\}\. Therefore, by Theorem 4\.4\.5 inVershynin \[[2018](https://arxiv.org/html/2606.04176#bib.bib33)\], we have with probability at least1−exp⁡\(−m1\)1\-\\exp\(\-m\_\{1\}\),

‖M′‖≲m2​πU2​ΔN4​\(m1\)=m1​m2​πU​ΔN2\\displaystyle\\\|M^\{\\prime\}\\\|\\lesssim\\sqrt\{m\_\{2\}\\pi\_\{U\}^\{2\}\\Delta^\{4\}\_\{N\}\}\(\\sqrt\{m\_\{1\}\}\)=\\sqrt\{m\_\{1\}m\_\{2\}\}\\pi\_\{U\}\\Delta\_\{N\}^\{2\}
Therefore, we obtain with probability at least1−κ/d1\-\\kappa/dfor some universal constantκ\>0\\kappa\>0,

‖T∘1ϵ\(1\)‖≲M​πU​ΔN​log⁡d\\displaystyle\\\|T\\circ\_\{1\}\\epsilon\_\{\(1\)\}\\\|\\lesssim\\sqrt\{M\\pi\_\{U\}\}\\Delta\_\{N\}\\log d\|∑i,j⟨A​\(i,j,⋅\),T​\(i,j\)​ϵ​\(i,j,⋅\)⟩ℋ\|≲M​πU​ΔN​log⁡d​‖A\(1\)‖∗\.\\displaystyle\\left\|\\sum\_\{i,j\}\\langle A\(i,j,\\cdot\),T\(i,j\)\\epsilon\(i,j,\\cdot\)\\rangle\_\{\\mathcal\{H\}\}\\right\|\\lesssim\\sqrt\{M\\pi\_\{U\}\}\\Delta\_\{N\}\\log d\\\|A\_\{\(1\)\}\\\|\_\{\*\}\.Similarly, one could prove for the second mode unfolding that with probability at least1−κ/d1\-\\kappa/dfor some universal constantκ\>0\\kappa\>0,

\|∑i,j⟨A​\(i,j,⋅\),T​\(i,j\)​ϵ​\(i,j,⋅\)⟩ℋ\|≲M​πU​ΔN​log⁡d​‖A\(2\)‖∗\.\\displaystyle\\left\|\\sum\_\{i,j\}\\langle A\(i,j,\\cdot\),T\(i,j\)\\epsilon\(i,j,\\cdot\)\\rangle\_\{\\mathcal\{H\}\}\\right\|\\lesssim\\sqrt\{M\\pi\_\{U\}\}\\Delta\_\{N\}\\log d\\\|A\_\{\(2\)\}\\\|\_\{\*\}\.
As for the third\-mode unfolding,

\|∑i,j⟨A​\(i,j,⋅\),T​\(i,j\)​ϵ​\(i,j,⋅\)⟩ℋ\|=⟨A\(3\),T∘3ϵ\(3\)⟩≤‖A\(3\)‖∗​‖T∘3ϵ\(3\)‖\\displaystyle\\left\|\\sum\_\{i,j\}\\langle A\(i,j,\\cdot\),T\(i,j\)\\epsilon\(i,j,\\cdot\)\\rangle\_\{\\mathcal\{H\}\}\\right\|=\\left\\langle A\_\{\(3\)\},T\\circ\_\{3\}\\epsilon\_\{\(3\)\}\\right\\rangle\\leq\\\|A\_\{\(3\)\}\\\|\_\{\*\}\\\|T\\circ\_\{3\}\\epsilon\_\{\(3\)\}\\\|
‖T∘ϵ\(3\)‖2\\displaystyle\\\|T\\circ\\epsilon\_\{\(3\)\}\\\|^\{2\}≤supf∈ℋ,‖f‖ℋ≤1∑i,jTi​j2​⟨ϵ​\(i,j,⋅\),f⟩ℋ2\\displaystyle\\leq\\sup\_\{f\\in\{\\mathcal\{H\}\},\\\|f\\\|\_\{\\mathcal\{H\}\}\\leq 1\}\\sum\_\{i,j\}T\_\{ij\}^\{2\}\\langle\\epsilon\(i,j,\\cdot\),f\\rangle\_\{\\mathcal\{H\}\}^\{2\}≤∑i,jTi​j2​‖ϵ​\(i,j,⋅\)‖ℋ2\\displaystyle\\leq\\sum\_\{i,j\}T\_\{ij\}^\{2\}\\\|\\epsilon\(i,j,\\cdot\)\\\|\_\{\\mathcal\{H\}\}^\{2\}Applying similar argument as above, we have with probability at least1−1/d1\-1/d,

∑i,jTi​j2​‖ϵ​\(i,j,⋅\)‖ℋ2≲m1​m2​πU​ΔN2\\displaystyle\\sum\_\{i,j\}T\_\{ij\}^\{2\}\\\|\\epsilon\(i,j,\\cdot\)\\\|\_\{\\mathcal\{H\}\}^\{2\}\\lesssim m\_\{1\}m\_\{2\}\\pi\_\{U\}\\Delta\_\{N\}^\{2\}\(23\)and therefore

\|∑i,j⟨A​\(i,j,⋅\),T​\(i,j\)​ϵ​\(i,j,⋅\)⟩ℋ\|≲m1​m2​πU​ΔN​‖A\(3\)‖∗\.\\displaystyle\\left\|\\sum\_\{i,j\}\\langle A\(i,j,\\cdot\),T\(i,j\)\\epsilon\(i,j,\\cdot\)\\rangle\_\{\\mathcal\{H\}\}\\right\|\\lesssim\\sqrt\{m\_\{1\}m\_\{2\}\\pi\_\{U\}\}\\Delta\_\{N\}\\\|A\_\{\(3\)\}\\\|\_\{\*\}\.\(24\)∎

###### Lemma C\.2\.

For anyA∈ℋm1×m2A\\in\{\\mathcal\{H\}\}^\{m\_\{1\}\\times m\_\{2\}\}, under the condition thatπL−1=𝒪​\(m1​m2\)\\pi\_\{L\}^\{\-1\}=\\scalebox\{0\.7\}\{$\\mathcal\{O\}$\}\(m\_\{1\}m\_\{2\}\), we have with probability at least1−κ/d1\-\\kappa/dfor some constantκ\>0\\kappa\>0,

1m1​m2​∑i,jTi,j​‖A​\(i,j,⋅\)‖ℋ2≥1m1​m2​∑i,jπi,j​‖A​\(i,j,⋅\)‖ℋ2−supi,j‖A​\(i,j,⋅\)‖ℋ2​πU​log⁡dm1​m2−supi,j‖A​\(i,j,⋅\)‖ℋ​‖A\(1\)‖∗m1​m2​πU​log⁡dm,∀A\.\\frac\{1\}\{m\_\{1\}m\_\{2\}\}\\sum\_\{i,j\}T\_\{i,j\}\\\|A\(i,j,\\cdot\)\\\|\_\{\\mathcal\{H\}\}^\{2\}\\geq\\frac\{1\}\{m\_\{1\}m\_\{2\}\}\\sum\_\{i,j\}\\pi\_\{i,j\}\\\|A\(i,j,\\cdot\)\\\|\_\{\\mathcal\{H\}\}^\{2\}\-\\sup\_\{i,j\}\\\|A\(i,j,\\cdot\)\\\|\_\{\\mathcal\{H\}\}^\{2\}\\sqrt\{\\frac\{\\pi\_\{U\}\\log d\}\{m\_\{1\}m\_\{2\}\}\}\\\\ \-\\sup\_\{i,j\}\\\|A\(i,j,\\cdot\)\\\|\_\{\\mathcal\{H\}\}\\frac\{\\\|A\_\{\(1\)\}\\\|\_\{\*\}\}\{\\sqrt\{m\_\{1\}m\_\{2\}\}\}\\sqrt\{\\frac\{\\pi\_\{U\}\\log d\}\{m\}\},\\qquad\\forall A\.and

1m1​m2​∑i,jTi,j​‖A​\(i,j,⋅\)‖ℋ2≥1m1​m2​∑i,jπi,j​‖A​\(i,j,⋅\)‖ℋ2−supi,j‖A​\(i,j,⋅\)‖ℋ2​πU​log⁡dm1​m2−supi,j‖A​\(i,j,⋅\)‖ℋ​‖A\(2\)‖∗m1​m2​πU​log⁡dm∀A\.\\frac\{1\}\{m\_\{1\}m\_\{2\}\}\\sum\_\{i,j\}T\_\{i,j\}\\\|A\(i,j,\\cdot\)\\\|\_\{\\mathcal\{H\}\}^\{2\}\\geq\\frac\{1\}\{m\_\{1\}m\_\{2\}\}\\sum\_\{i,j\}\\pi\_\{i,j\}\\\|A\(i,j,\\cdot\)\\\|\_\{\\mathcal\{H\}\}^\{2\}\-\\sup\_\{i,j\}\\\|A\(i,j,\\cdot\)\\\|\_\{\\mathcal\{H\}\}^\{2\}\\sqrt\{\\frac\{\\pi\_\{U\}\\log d\}\{m\_\{1\}m\_\{2\}\}\}\\\\ \-\\sup\_\{i,j\}\\\|A\(i,j,\\cdot\)\\\|\_\{\\mathcal\{H\}\}\\frac\{\\\|A\_\{\(2\)\}\\\|\_\{\*\}\}\{\\sqrt\{m\_\{1\}m\_\{2\}\}\}\\sqrt\{\\frac\{\\pi\_\{U\}\\log d\}\{m\}\}\\qquad\\forall A\.

###### Proof\.

Consider the class𝒦​\(1,α\)=\{A:A∈ℋm1×m2,supi,j‖A​\(i,j,⋅\)‖ℋ≤1,‖A\(1\)‖∗≤α\}\{\\mathcal\{K\}\}\(1,\\alpha\)=\\\{A:A\\in\{\\mathcal\{H\}\}^\{m\_\{1\}\\times m\_\{2\}\},\\sup\_\{i,j\}\\\|A\(i,j,\\cdot\)\\\|\_\{\\mathcal\{H\}\}\\leq 1,\\\|A\_\{\(1\)\}\\\|\_\{\*\}\\leq\\alpha\\\}\. Then𝒦​\(1,α\)\{\\mathcal\{K\}\}\(1,\\alpha\)can be represented as

\{α​∑l=1m1wl​a\(l\)​f\(l\),supi,j‖∑l=1m1wl​ai\(l\)​fj\(l\)∥ℋ≤1/α,∑l=1m1wl≤1,wl≥0,‖a\(l\)‖2≤1,∑j=1m2‖fj\(l\)‖ℋ2≤1\}\\displaystyle\\left\\\{\\alpha\\sum\_\{l=1\}^\{m\_\{1\}\}w\_\{l\}a^\{\(l\)\}f^\{\(l\)\},\\sup\_\{i,j\}\\left\\\|\\sum\_\{l=1\}^\{m\_\{1\}\}w\_\{l\}a^\{\(l\)\}\_\{i\}f^\{\(l\)\}\_\{j\}\\right\\\|\_\{\\mathcal\{H\}\}\\leq 1/\\alpha,\\sum\_\{l=1\}^\{m\_\{1\}\}w\_\{l\}\\leq 1,w\_\{l\}\\geq 0,\\\|a^\{\(l\)\}\\\|\_\{2\}\\leq 1,\\sum\_\{j=1\}^\{m\_\{2\}\}\\\|f^\{\(l\)\}\_\{j\}\\\|^\{2\}\_\{\\mathcal\{H\}\}\\leq 1\\right\\\}
Takeϵi,j\\epsilon\_\{i,j\}as independent Rademacher random variables\.

supA\(1\)∈𝒦​\(1,α\)\|∑i,jTi,j​ϵi,j​‖A​\(i,j,⋅\)‖ℋ2\|\\displaystyle\\sup\_\{A\_\{\(1\)\}\\in\{\\mathcal\{K\}\}\(1,\\alpha\)\}\\left\|\\sum\_\{i,j\}T\_\{i,j\}\\epsilon\_\{i,j\}\\\|A\(i,j,\\cdot\)\\\|\_\{\\mathcal\{H\}\}^\{2\}\\right\|≤\\displaystyle\\leqsupw,a,fα​\|∑i,jTi,j​ϵi,j​⟨∑l=1m1wl​ai\(l\)​fj\(l\),∑l=1m1wl​ai\(l\)​fj\(l\)⟩ℋ\|\\displaystyle\\sup\_\{w,a,f\}\\alpha\\left\|\\sum\_\{i,j\}T\_\{i,j\}\\epsilon\_\{i,j\}\\left\\langle\\sum\_\{l=1\}^\{m\_\{1\}\}w\_\{l\}a^\{\(l\)\}\_\{i\}f^\{\(l\)\}\_\{j\},\\sum\_\{l=1\}^\{m\_\{1\}\}w\_\{l\}a^\{\(l\)\}\_\{i\}f^\{\(l\)\}\_\{j\}\\right\\rangle\_\{\\mathcal\{H\}\}\\right\|≤\\displaystyle\\leqαsupw\{\|∑l=1m1wl2sup‖a\(l\)‖2≤1,‖f\(l\)‖ℋm2≤1∑i,jTi,jϵi,jai\(l\)2∥fj\(l\)∥ℋ2\|\\displaystyle\\alpha\\sup\_\{w\}\\left\\\{\\left\|\\sum\_\{l=1\}^\{m\_\{1\}\}w\_\{l\}^\{2\}\\sup\_\{\\\|a^\{\(l\)\}\\\|\_\{2\}\\leq 1,\\\|f^\{\(l\)\}\\\|\_\{\{\\mathcal\{H\}\}^\{m\_\{2\}\}\}\\leq 1\}\\sum\_\{i,j\}T\_\{i,j\}\\epsilon\_\{i,j\}\{a^\{\(l\)\}\_\{i\}\}^\{2\}\\\|f^\{\(l\)\}\_\{j\}\\\|\_\{\\mathcal\{H\}\}^\{2\}\\right\|\\right\.\+2∑l≠lwlwl′sup‖a\(l\)‖2≤1,‖f\(l\)‖ℋm2\|∑i,jTi,jϵi,jai\(l\)ai\(l′\)⟨fj\(l\),fj\(l′\)⟩ℋ\|\}\\displaystyle\\left\.\+2\\sum\_\{l\\neq l\}w\_\{l\}w\_\{l^\{\\prime\}\}\\sup\_\{\\\|a^\{\(l\)\}\\\|\_\{2\}\\leq 1,\\\|f^\{\(l\)\}\\\|\_\{\{\\mathcal\{H\}\}^\{m\_\{2\}\}\}\}\\left\|\\sum\_\{i,j\}T\_\{i,j\}\\epsilon\_\{i,j\}a^\{\(l\)\}\_\{i\}a^\{\(l^\{\\prime\}\)\}\_\{i\}\\left\\langle f^\{\(l\)\}\_\{j\},f^\{\(l^\{\\prime\}\)\}\_\{j\}\\right\\rangle\_\{\\mathcal\{H\}\}\\right\|\\right\\\}≤\\displaystyle\\leqα​supw\{\|∑l=1m1wl2​sup‖a\(l\)‖2≤1,‖f\(l\)‖ℋm2≤1∑i,jTi,j​ϵi,j​ai\(l\)2​‖fj\(l\)‖ℋ2​\|\+2​∑l≠lwl​wl′‖​T∘ϵ‖\}\\displaystyle\\alpha\\sup\_\{w\}\\left\\\{\\left\|\\sum\_\{l=1\}^\{m\_\{1\}\}w\_\{l\}^\{2\}\\sup\_\{\\\|a^\{\(l\)\}\\\|\_\{2\}\\leq 1,\\\|f^\{\(l\)\}\\\|\_\{\{\\mathcal\{H\}\}^\{m\_\{2\}\}\}\\leq 1\}\\sum\_\{i,j\}T\_\{i,j\}\\epsilon\_\{i,j\}\{a^\{\(l\)\}\_\{i\}\}^\{2\}\\\|f^\{\(l\)\}\_\{j\}\\\|\_\{\\mathcal\{H\}\}^\{2\}\\right\|\+2\\sum\_\{l\\neq l\}w\_\{l\}w\_\{l^\{\\prime\}\}\\\|T\\circ\\epsilon\\\|\\right\\\}\(25\)The second last inequality is due to\(∑i=1m1\[ai\(l\)ai\(l′\)\]2≤\(∑i=1m1\[ai\(l\)\]2\)\(∑i=1m1\[ai\(l′\)\]2\)≤1\(\\sum\_\{i=1\}^\{m\_\{1\}\}\[a^\{\(l\)\}\_\{i\}a^\{\(l^\{\\prime\}\)\}\_\{i\}\]^\{2\}\\leq\(\\sum\_\{i=1\}^\{m\_\{1\}\}\[a^\{\(l\)\}\_\{i\}\]^\{2\}\)\(\\sum\_\{i=1\}^\{m\_\{1\}\}\[a^\{\(l^\{\\prime\}\)\}\_\{i\}\]^\{2\}\)\\leq 1and\(∑i=1m2⟨fj\(l\),fj\(l′\)⟩ℋ2≤\(∑i=1m2∥fj\(l\)∥ℋ2\)\(∑i=1m2∥fj\(l′\)∥ℋ2\)≤1\(\\sum\_\{i=1\}^\{m\_\{2\}\}\\langle f\_\{j\}^\{\(l\)\},f\_\{j\}^\{\(l^\{\\prime\}\)\}\\rangle\_\{\\mathcal\{H\}\}^\{2\}\\leq\(\\sum\_\{i=1\}^\{m\_\{2\}\}\\\|f^\{\(l\)\}\_\{j\}\\\|\_\{\\mathcal\{H\}\}^\{2\}\)\(\\sum\_\{i=1\}^\{m\_\{2\}\}\\\|f^\{\(l^\{\\prime\}\)\}\_\{j\}\\\|\_\{\\mathcal\{H\}\}^\{2\}\)\\leq 1\.

Note that for anyi,ji,j,

‖∑l=1m1wl​ai\(l\)​fj\(l\)‖ℋ2\\displaystyle\\left\\\|\\sum\_\{l=1\}^\{m\_\{1\}\}w\_\{l\}a^\{\(l\)\}\_\{i\}f^\{\(l\)\}\_\{j\}\\right\\\|^\{2\}\_\{\\mathcal\{H\}\}≤\(1/α\)2\\displaystyle\\leq\(1/\\alpha\)^\{2\}∑l=1m1wl2​\[ai\(l\)\]2​‖fj\(l\)‖ℋ2\\displaystyle\\sum\_\{l=1\}^\{m\_\{1\}\}w\_\{l\}^\{2\}\[a\_\{i\}^\{\(l\)\}\]^\{2\}\\\|f\_\{j\}^\{\(l\)\}\\\|\_\{\\mathcal\{H\}\}^\{2\}≤\(1/α\)2\+\|∑l≠lwl​wl′​ai\(l\)​ai\(l′\)​⟨fj\(l\),fj\(l′\)⟩ℋ\|\\displaystyle\\leq\(1/\\alpha\)^\{2\}\+\\left\|\\sum\_\{l\\neq l\}w\_\{l\}w\_\{l^\{\\prime\}\}a^\{\(l\)\}\_\{i\}a^\{\(l^\{\\prime\}\)\}\_\{i\}\\left\\langle f^\{\(l\)\}\_\{j\},f^\{\(l^\{\\prime\}\)\}\_\{j\}\\right\\rangle\_\{\\mathcal\{H\}\}\\right\|≤\(1/α\)2\+\|∑l≠lwl​wl′\|≤\(1/α\)2\+1\\displaystyle\\leq\(1/\\alpha\)^\{2\}\+\\left\|\\sum\_\{l\\neq l\}w\_\{l\}w\_\{l^\{\\prime\}\}\\right\|\\leq\(1/\\alpha\)^\{2\}\+1suplwl2​\[ai\(l\)\]2​‖fj\(l\)‖ℋ2\\displaystyle\\sup\_\{l\}w\_\{l\}^\{2\}\[a\_\{i\}^\{\(l\)\}\]^\{2\}\\\|f\_\{j\}^\{\(l\)\}\\\|\_\{\\mathcal\{H\}\}^\{2\}≤\(1/α\)2\+1\\displaystyle\\leq\(1/\\alpha\)^\{2\}\+1
From the contraction inequality for Rademacher complexity, we have

∑l=1m1𝔼​supa,f,w\|∑i,jTi,j​ϵi,j​wl2​\[ai\(l\)\]2​‖fj\(l\)‖ℋ2\|\\displaystyle\\sum\_\{l=1\}^\{m\_\{1\}\}\\mathbb\{E\}\\sup\_\{a,f,w\}\\left\|\\sum\_\{i,j\}T\_\{i,j\}\\epsilon\_\{i,j\}w\_\{l\}^\{2\}\[a\_\{i\}^\{\(l\)\}\]^\{2\}\\\|f\_\{j\}^\{\(l\)\}\\\|\_\{\\mathcal\{H\}\}^\{2\}\\right\|≤\\displaystyle\\leq2​\(1/α\)2\+1​∑l=1m1𝔼​supa,f,wwl​\|∑i,jTi,j​ϵi,j​\[ai\(l\)\]​‖fj\(l\)‖ℋ\|\\displaystyle 2\\sqrt\{\(1/\\alpha\)^\{2\}\+1\}\\sum\_\{l=1\}^\{m\_\{1\}\}\\mathbb\{E\}\\sup\_\{a,f,w\}w\_\{l\}\\left\|\\sum\_\{i,j\}T\_\{i,j\}\\epsilon\_\{i,j\}\[a\_\{i\}^\{\(l\)\}\]\\\|f\_\{j\}^\{\(l\)\}\\\|\_\{\\mathcal\{H\}\}\\right\|≤\\displaystyle\\leq2​\(1/α\)2\+1​supw∑lwl​𝔼​sup‖a‖2≤1,∑j=1m2‖fj‖ℋ2≤1\|∑i,jTi,j​ϵi,j​ai​‖fj‖ℋ\|\\displaystyle 2\\sqrt\{\(1/\\alpha\)^\{2\}\+1\}\\sup\_\{w\}\\sum\_\{l\}w\_\{l\}\\mathbb\{E\}\\sup\_\{\\\|a\\\|\_\{2\}\\leq 1,\\sum\_\{j=1\}^\{m\_\{2\}\}\\\|f\_\{j\}\\\|^\{2\}\_\{\\mathcal\{H\}\}\\leq 1\}\\left\|\\sum\_\{i,j\}T\_\{i,j\}\\epsilon\_\{i,j\}a\_\{i\}\\\|f\_\{j\}\\\|\_\{\\mathcal\{H\}\}\\right\|≤\\displaystyle\\leq2​\(1/α\)2\+1​supw∑lwl​𝔼​‖T∘ϵ‖≤2​\(1/α\)2\+1​supw𝔼​‖T∘ϵ‖\\displaystyle 2\\sqrt\{\(1/\\alpha\)^\{2\}\+1\}\\sup\_\{w\}\\sum\_\{l\}w\_\{l\}\\mathbb\{E\}\\\|T\\circ\\epsilon\\\|\\leq 2\\sqrt\{\(1/\\alpha\)^\{2\}\+1\}\\sup\_\{w\}\\mathbb\{E\}\\\|T\\circ\\epsilon\\\|
Therefore, due to the symmetrization inequality, we have

𝔼​supA∈𝒦​\(1,α\)\|∑i,jTi,j​‖A​\(i,j,⋅\)‖ℋ2−∑i,jπi,j​‖A​\(i,j,⋅\)‖ℋ2\|\\displaystyle\\mathbb\{E\}\\sup\_\{A\\in\{\\mathcal\{K\}\}\(1,\\alpha\)\}\\left\|\\sum\_\{i,j\}T\_\{i,j\}\\\|A\(i,j,\\cdot\)\\\|\_\{\\mathcal\{H\}\}^\{2\}\-\\sum\_\{i,j\}\\pi\_\{i,j\}\\\|A\(i,j,\\cdot\)\\\|\_\{\\mathcal\{H\}\}^\{2\}\\right\|≤2​𝔼​supA\(1\)∈𝒦​\(1,α\)\|∑i,jTi,j​ϵi,j​‖A​\(i,j,⋅\)‖ℋ2\|\\displaystyle\\leq 2\\mathbb\{E\}\\sup\_\{A\_\{\(1\)\}\\in\{\\mathcal\{K\}\}\(1,\\alpha\)\}\\left\|\\sum\_\{i,j\}T\_\{i,j\}\\epsilon\_\{i,j\}\\\|A\(i,j,\\cdot\)\\\|\_\{\\mathcal\{H\}\}^\{2\}\\right\|≲α​\[\(1/α\)2\+1\+1\]​𝔼​‖T∘ϵ‖\\displaystyle\\lesssim\\alpha\[\\sqrt\{\(1/\\alpha\)^\{2\}\+1\}\+1\]\\mathbb\{E\}\\\|T\\circ\\epsilon\\\|≲α​πU​M\\displaystyle\\lesssim\\alpha\\sqrt\{\\pi\_\{U\}M\}the third inequality is due to \([25](https://arxiv.org/html/2606.04176#A3.E25)\) and the last inequality can be obtained from Lemma S7 inWang et al\. \[[2021](https://arxiv.org/html/2606.04176#bib.bib34)\]by settingW^\\hat\{W\}to be the matrix whose entries are all ones\.

Next, we adopt the Talagrand concentration inequality and the peeling argument to establish the concentration inequality\. ConsiderA∈𝒦​\(1,α\)A\\in\{\\mathcal\{K\}\}\(1,\\alpha\), takeXi,j=\(Ti,j−πi,j\)​‖A​\(i,j,⋅\)‖ℋ2X\_\{i,j\}=\(T\_\{i,j\}\-\\pi\_\{i,j\}\)\\\|A\(i,j,\\cdot\)\\\|\_\{\\mathcal\{H\}\}^\{2\}\. It is easy to see that\|Xi,j\|≤1\|X\_\{i,j\}\|\\leq 1,

supA\(1\)∈𝒦​\(1,α\)𝔼​∑i,jXi,j2\+16​supi,j\|Xi,j\|​𝔼​supA∈𝒦​\(1,α\)\|∑i,jXi,j\|≤m1​m2​πU\+16​α​πU​M\.\\displaystyle\\sup\_\{A\_\{\(1\)\}\\in\{\\mathcal\{K\}\}\(1,\\alpha\)\}\{\\mathbb\{E\}\}\\sum\_\{i,j\}X^\{2\}\_\{i,j\}\+16\\sup\_\{i,j\}\|X\_\{i,j\}\|\{\\mathbb\{E\}\}\\sup\_\{A\\in\{\\mathcal\{K\}\}\(1,\\alpha\)\}\|\\sum\_\{i,j\}X\_\{i,j\}\|\\leq m\_\{1\}m\_\{2\}\\pi\_\{U\}\+16\\alpha\\sqrt\{\\pi\_\{U\}M\}\.\(26\)Therefore, we could takeU=1U=1andV=m1​m2​πU\+16​α​πU​MV=m\_\{1\}m\_\{2\}\\pi\_\{U\}\+16\\alpha\\sqrt\{\\pi\_\{U\}M\}in Theorem 2\.6 inKoltchinskii \[[2011](https://arxiv.org/html/2606.04176#bib.bib20)\], which leads to

Pr⁡\(\[supA\(1\)∈𝒦​\(1,α\)\|∑i,jXi,j\|−𝔼​supA\(1\)∈𝒦​\(1,α\)\|∑i,jXi,j\|\]≥t\)\\displaystyle\\Pr\\left\(\\left\[\\sup\_\{A\_\{\(1\)\}\\in\{\\mathcal\{K\}\}\(1,\\alpha\)\}\\left\|\\sum\_\{i,j\}X\_\{i,j\}\\right\|\-\{\\mathbb\{E\}\}\\sup\_\{A\_\{\(1\)\}\\in\{\\mathcal\{K\}\}\(1,\\alpha\)\}\\left\|\\sum\_\{i,j\}X\_\{i,j\}\\right\|\\right\]\\geq t\\right\)≤K​exp⁡\{−1K​tU​log⁡\(1\+t​UV\)\}≤K​exp⁡\{−1K​t2m1​m2​πU\+16​α​πU​M\+t\}\\displaystyle\\leq K\\exp\\left\\\{\-\\frac\{1\}\{K\}\\frac\{t\}\{U\}\\log\\left\(1\+\\frac\{tU\}\{V\}\\right\)\\right\\\}\\leq K\\exp\\left\\\{\-\\frac\{1\}\{K\}\\frac\{t^\{2\}\}\{m\_\{1\}m\_\{2\}\\pi\_\{U\}\+16\\alpha\\sqrt\{\\pi\_\{U\}M\}\+t\}\\right\\\}\(27\)for some universal constantK\>0K\>0\.

Whenα≤m\\alpha\\leq\\sqrt\{m\}, we have

Pr⁡\(supA\(1\)∈𝒦​\(1,α\)\|∑i,jXi,j\|≥2​πU​m1​m2​log⁡d\)\\displaystyle\\Pr\\left\(\\sup\_\{A\_\{\(1\)\}\\in\{\\mathcal\{K\}\}\(1,\\alpha\)\}\\left\|\\sum\_\{i,j\}X\_\{i,j\}\\right\|\\geq 2\\sqrt\{\\pi\_\{U\}m\_\{1\}m\_\{2\}\\log d\}\\right\)≤Pr⁡\(supA\(1\)∈𝒦​\(1,m\)\|∑i,jXi,j\|≥2​πU​m1​m2​log⁡d\)≤cd\\displaystyle\\leq\\Pr\\left\(\\sup\_\{A\_\{\(1\)\}\\in\{\\mathcal\{K\}\}\(1,\\sqrt\{m\}\)\}\\left\|\\sum\_\{i,j\}X\_\{i,j\}\\right\|\\geq 2\\sqrt\{\\pi\_\{U\}m\_\{1\}m\_\{2\}\\log d\}\\right\)\\leq\\frac\{c\}\{d\}\(28\)for some constantc\>0c\>0\. Whenα≥m\\alpha\\geq\\sqrt\{m\}, we have fort′≥1/2t^\{\\prime\}\\geq 1/2,

Pr⁡\(supA\(1\)∈𝒦​\(1,α\)\|∑i,jXi,j\|≥t′​α​M​πU​log⁡d\)≤K′​exp⁡\{−1K′​t′2​α2​M​πU​log⁡dm1​m2​πU\+t′​α​M​πU​log⁡d\},\\displaystyle\\Pr\\left\(\\sup\_\{A\_\{\(1\)\}\\in\{\\mathcal\{K\}\}\(1,\\alpha\)\}\\left\|\\sum\_\{i,j\}X\_\{i,j\}\\right\|\\geq t^\{\\prime\}\\alpha\\sqrt\{M\\pi\_\{U\}\\log d\}\\right\)\\leq K^\{\\prime\}\\exp\\left\\\{\-\\frac\{1\}\{K^\{\\prime\}\}\\frac\{\{t^\{\\prime\}\}^\{2\}\\alpha^\{2\}M\\pi\_\{U\}\\log d\}\{m\_\{1\}m\_\{2\}\\pi\_\{U\}\+t^\{\\prime\}\\alpha\\sqrt\{M\\pi\_\{U\}\\log d\}\}\\right\\\},\(29\)for some universal constantK′\>0K^\{\\prime\}\>0\. Next, consider a sequence of the sets

𝒮l:=\{A∈𝒦​\(1,α\):2l−1​m≤α≤2l​m\},l=1,2,…\\displaystyle\{\\mathcal\{S\}\}\_\{l\}:=\\\{A\\in\{\\mathcal\{K\}\}\(1,\\alpha\):2^\{l\-1\}\\sqrt\{m\}\\leq\\alpha\\leq 2^\{l\}\\sqrt\{m\}\\\},\\quad l=1,2,\\dots\(30\)Then we have

Pr⁡\(supsupi,j‖A​\(i,j,⋅\)‖ℋ≤1\|∑i,jXi,j\|‖A\(1\)‖∗≥M​πU​log⁡d\)\\displaystyle\\Pr\\left\(\\sup\_\{\\sup\_\{i,j\}\\\|A\(i,j,\\cdot\)\\\|\_\{\\mathcal\{H\}\}\\leq 1\}\\frac\{\|\\sum\_\{i,j\}X\_\{i,j\}\|\}\{\\\|A\_\{\(1\)\}\\\|\_\{\*\}\}\\geq\\sqrt\{M\\pi\_\{U\}\\log d\}\\right\)≤∑l=1∞Pr⁡\(supA∈𝒮l\|∑i,jXi,j\|‖A\(1\)‖∗≥M​πU​log⁡d\)\\displaystyle\\leq\\sum\_\{l=1\}^\{\\infty\}\\Pr\\left\(\\sup\_\{A\\in\{\\mathcal\{S\}\}\_\{l\}\}\\frac\{\|\\sum\_\{i,j\}X\_\{i,j\}\|\}\{\\\|A\_\{\(1\)\}\\\|\_\{\*\}\}\\geq\\sqrt\{M\\pi\_\{U\}\\log d\}\\right\)≤∑l=1∞Pr⁡\(supA∈𝒦​\(1,2l​m\)\|∑i,jXi,j\|≥2l−1​m​M​πU​log⁡d\)\\displaystyle\\leq\\sum\_\{l=1\}^\{\\infty\}\\Pr\\left\(\\sup\_\{A\\in\{\\mathcal\{K\}\}\(1,2^\{l\}\\sqrt\{m\}\)\}\|\\sum\_\{i,j\}X\_\{i,j\}\|\\geq 2^\{l\-1\}\\sqrt\{m\}\\sqrt\{M\\pi\_\{U\}\\log d\}\\right\)≤∑l=1∞K′​exp⁡\{−1K′​22​l−2​log⁡d1\+2l−1​log⁡dm1​m2​πU\}\\displaystyle\\leq\\sum\_\{l=1\}^\{\\infty\}K^\{\\prime\}\\exp\\left\\\{\-\\frac\{1\}\{K^\{\\prime\}\}\\frac\{2^\{2l\-2\}\\log d\}\{1\+2^\{l\-1\}\\sqrt\{\\frac\{\\log d\}\{m\_\{1\}m\_\{2\}\\pi\_\{U\}\}\}\}\\right\\\}≤∑l=1∞K′​exp⁡\{−1K′​2l−1​log⁡d1\+log⁡dm1​m2​πU\}\\displaystyle\\leq\\sum\_\{l=1\}^\{\\infty\}K^\{\\prime\}\\exp\\left\\\{\-\\frac\{1\}\{K^\{\\prime\}\}\\frac\{2^\{l\-1\}\\log d\}\{1\+\\sqrt\{\\frac\{\\log d\}\{m\_\{1\}m\_\{2\}\\pi\_\{U\}\}\}\}\\right\\\}≤∑l=1∞K′​exp⁡\{−log⁡22​K′​l​log⁡d1\+log⁡dm1​m2​πU\}\\displaystyle\\leq\\sum\_\{l=1\}^\{\\infty\}K^\{\\prime\}\\exp\\left\\\{\-\\frac\{\\log 2\}\{2K^\{\\prime\}\}\\frac\{l\\log d\}\{1\+\\sqrt\{\\frac\{\\log d\}\{m\_\{1\}m\_\{2\}\\pi\_\{U\}\}\}\}\\right\\\}≤K′​exp⁡\{−log⁡22​K′​log⁡d1\+log⁡dm1​m2​πU\}1−exp⁡\{−log⁡22​K′​log⁡d1\+log⁡dm1​m2​πU\}≤cd,\\displaystyle\\leq K^\{\\prime\}\\frac\{\\exp\\left\\\{\-\\frac\{\\log 2\}\{2K^\{\\prime\}\}\\frac\{\\log d\}\{1\+\\sqrt\{\\frac\{\\log d\}\{m\_\{1\}m\_\{2\}\\pi\_\{U\}\}\}\}\\right\\\}\}\{1\-\\exp\\left\\\{\-\\frac\{\\log 2\}\{2K^\{\\prime\}\}\\frac\{\\log d\}\{1\+\\sqrt\{\\frac\{\\log d\}\{m\_\{1\}m\_\{2\}\\pi\_\{U\}\}\}\}\\right\\\}\}\\leq\\frac\{c\}\{d\},\(31\)for some constantc\>0c\>0\.

Combine \([28](https://arxiv.org/html/2606.04176#A3.E28)\) and \([31](https://arxiv.org/html/2606.04176#A3.E31)\), we have with probability at least1−κ/d1\-\\kappa/dfor some constantκ\>0\\kappa\>0and for anyAAsuch thatsupi,j‖A​\(i,j,⋅\)‖ℋ≤1\\sup\_\{i,j\}\\\|A\(i,j,\\cdot\)\\\|\_\{\\mathcal\{H\}\}\\leq 1,

\|1m1​m2​∑i,jTi,j​‖A​\(i,j,⋅\)‖ℋ2−1m1​m2​∑i,jπi,j​‖A​\(i,j,⋅\)‖ℋ2\|≤πU​log⁡dm1​m2\+‖A\(1\)‖∗​M​πU​log⁡dm1​m2\.\\displaystyle\\left\|\\frac\{1\}\{m\_\{1\}m\_\{2\}\}\\sum\_\{i,j\}T\_\{i,j\}\\\|A\(i,j,\\cdot\)\\\|\_\{\\mathcal\{H\}\}^\{2\}\-\\frac\{1\}\{m\_\{1\}m\_\{2\}\}\\sum\_\{i,j\}\\pi\_\{i,j\}\\\|A\(i,j,\\cdot\)\\\|\_\{\\mathcal\{H\}\}^\{2\}\\right\|\\leq\\sqrt\{\\frac\{\\pi\_\{U\}\\log d\}\{m\_\{1\}m\_\{2\}\}\}\+\\\|A\_\{\(1\)\}\\\|\_\{\*\}\\frac\{\\sqrt\{M\\pi\_\{U\}\\log d\}\}\{m\_\{1\}m\_\{2\}\}\.\(32\)Therefore,

1m1​m2​∑i,jTi,j​‖A​\(i,j,⋅\)‖ℋ2≥1m1​m2​∑i,jπi,j​‖A​\(i,j,⋅\)‖ℋ2−supi,j‖A​\(i,j,⋅\)‖ℋ2​πU​log⁡dm1​m2−supi,j‖A​\(i,j,⋅\)‖ℋ​‖A\(1\)‖∗m1​m2​πU​log⁡dm\\frac\{1\}\{m\_\{1\}m\_\{2\}\}\\sum\_\{i,j\}T\_\{i,j\}\\\|A\(i,j,\\cdot\)\\\|\_\{\\mathcal\{H\}\}^\{2\}\\geq\\frac\{1\}\{m\_\{1\}m\_\{2\}\}\\sum\_\{i,j\}\\pi\_\{i,j\}\\\|A\(i,j,\\cdot\)\\\|\_\{\\mathcal\{H\}\}^\{2\}\-\\sup\_\{i,j\}\\\|A\(i,j,\\cdot\)\\\|\_\{\\mathcal\{H\}\}^\{2\}\\sqrt\{\\frac\{\\pi\_\{U\}\\log d\}\{m\_\{1\}m\_\{2\}\}\}\\\\ \-\\sup\_\{i,j\}\\\|A\(i,j,\\cdot\)\\\|\_\{\\mathcal\{H\}\}\\frac\{\\\|A\_\{\(1\)\}\\\|\_\{\*\}\}\{\\sqrt\{m\_\{1\}m\_\{2\}\}\}\\sqrt\{\\frac\{\\pi\_\{U\}\\log d\}\{m\}\}Similar results can be obtained for the second\-mode unfolding\.

∎