TIDE-MC: Two-Sided Interpolative Decomposition for Billion-Scale GPU Matrix Completion

TIDE-MC: Two-Sided Interpolative Decomposition for Billion-Scale GPU Matrix Completion Chengying Huan1 , Yubo Wang2 , Pinhuan Wang1 , Lizheng Chen1 , Jie Zhang3 Fangxin Liu4 , Qing Wang1 , Ruixuan Liu5 , Shaonan Ma6 , Mingxing Zhang7 Zhibin Wang1 , Rong Gu1 , Guihai Chen1 , Chen Tian1 1 Nanjing University

4 Shanghai Jiao Tong University

2 Northeastern University (China)

5 Dalian Maritime University

arXiv:2608.00977v1 [cs.DC] 2 Aug 2026

Abstract

7 Tsinghua University

with low-dimensional matrices 𝑃 and 𝑄 while optimizing only over the observed entries. In real-world scenarios, matrix completion often operates on enormous user-item matrices or high-dimensional data tables, of which only a tiny fraction of entries are observed. Such matrices are typically sparse and approximately lowrank [3, 53]. For example, users in recommendation systems rate only a small subset of items [43, 64]; in computer vision, sensor measurements may miss pixels or views; and in genomic or network analysis, only a subset of possible interactions is measured [51]. These shared characteristics make matrix completion both practically important and increasingly difficult to scale as data dimensions continue to grow. As these applications grow to industrial and scientific scales, matrix completion becomes a large-scale systems and HPC challenge in addition to an optimization problem. Modern deployments must handle matrices with billions of potential entries under strict resource constraints, which further amplifies the need for scalable techniques in largescale data processing systems [35, 40, 49, 61, 67]. For example, Netflix manages user–item rating matrices with over 100 million observed entries [5], while Yahoo! Music processes ratings from 200,000 users across 136,000 songs, amounting to over 27 billion potential entries [64]. Under such scale, practical matrix completion must operate within limited memory while sustaining high throughput and low latency for near real-time imputation. Traditional matrix completion methods largely fall into two families: nuclear norm minimization (NN)-centric approaches [24, 31] and singular value decomposition (SVD)centric approaches [2, 14, 20, 23]. NN-centric methods recover missing entries by minimizing the nuclear norm of the matrix, either directly [52] or after reducing the matrix via column selection, as in CSNN [31]. SVD-centric methods operate on low-rank factorizations (e.g., 𝐴 = 𝑈 Σ𝑉 𝑇 ) and use iterative solvers such as CSPGD and CGM [68], kernel-based approximations like NTK [48], or block-wise partitioning as in MVGMC [29] to improve scalability. Table 1 summarizes the time and space complexity of the above methods and shows that, despite these optimizations, their costs still scale unfavorably with the matrix size, limiting their applicability to large-scale matrix completion.

Matrix completion supports large-scale recommendation and scientific computing, yet existing GPU solvers commonly assume that the observed matrix or its dense factors fit in device memory. On real workloads, this assumption leads to out-of-memory failures or severe PCIe overhead under naive paging. We present TIDE-MC, a bounded-memory GPU framework built on Two-Sided Interpolative Decomposition (TSID). TSID uses a sampled template submatrix as an anchor for reconstructing the full low-rank matrix, allowing computation and storage to scale with the template and active data chunks rather than the complete matrix. TIDE-MC realizes this formulation through two execution stages. First, a conflict-free synchronization engine recovers the template using parallel factorization and hierarchical gradient aggregation. Second, a chunked reconstruction pipeline extends the recovered template to the remaining matrix while overlapping PCIe transfers with GPU computation. An asymmetric gradientclipping scheme stabilizes mixed-precision Tensor Core execution. Across 15 benchmarks, TIDE-MC completes workloads that cause existing GPU solvers to run out of memory. Compared with the evaluated state-of-the-art baselines, it achieves up to 11,647× speedup, reduces peak memory usage by up to 8.5×, and lowers reconstruction error by up to 99.7%. These results show that template-anchored decomposition and stage-specific GPU execution can scale matrix completion beyond device-memory capacity.

1

3 Peking University

6 Qiyuan Lab

Introduction

Matrix completion has become a cornerstone of modern data-driven applications, such as collaborative filtering [4, 6, 16], genomic data analysis [8, 11, 25, 37, 41], network reconstruction [50], graph neural networks [26], information processing systems [56], image reconstruction [10, 24, 31, 32, 48, 52, 58], industrial recommendation systems [5, 30], and multi-label learning [33]. Formally, matrix completion involves recovering a complete matrix 𝑅b from a sparse set of observations Ω of an underlying matrix 𝑅 ∈ R𝑚×𝑛 . In practice, this is typically achieved by exploiting low-rank structure,e.g., approximating 𝑅b or a candidate completion 1

Table 1. Comparison of time and resident space complexity for prior systems and TIDE-MC. The matrix is in R𝑛×𝑛 ; 𝛼 is the template sampling rate (0.1–0.3); 𝛿 is the observed-entry fraction(8%–20%); 𝑘 is the latent rank, approximately 𝛼𝑛; 𝑡 1, 𝑡 2 are the iteration counts of Stage I and Stage II(ranging from 100 to 10,000); and Δ𝑐 is the chunk width. For TIDEMC, space denotes the peak resident GPU working set rather than the total streamed input size. Type NNCentric

Name NN [52] CSNN [31]

Time O (𝑘𝑛 3 ) O (𝛼 2𝑛 3 )

Space O (𝑛 2 ) O (𝛼𝑛 2 )

SVDCentric

CSPGD [31] CGM [31] NTK [48] MVGMC [29]

O (𝛼𝑛 3 ) O (𝑘𝑛 2 ) O (𝑘𝑛 2 ) O (𝑘𝑛 3 (1 + 𝛿))

O (𝛼 2𝑛 3 ) O (𝑛 2 ) O (𝛿𝑛 2 ) O (𝑛 2 )

O (𝑘𝛿𝑛 2 (𝑡 2 + 𝛼𝑡 1 ))

O (𝑘 (1 + 𝛼)𝑛 + |ΩΔ𝑐 | + 𝑘Δ𝑐)

Our System TIDE-MC

trix completion. Guided by TSID, we introduce a two-stage design: Stage I recovers a template submatrix, and Stage II extends it to the full matrix under a fixed memory budget. • Bounded-Memory Two-Stage Execution Framework. We design a TSID-guided execution framework. Stage I recovers the template submatrix using a conflict-free parallel SGD-based RMF kernel with gradient aggregation, reducing atomic contention and mapping irregular updates to dense Tensor Core tiles. Stage II reconstructs the full matrix via roofline-guided chunked alternating least squares (ALS), using double buffering and autotuned chunk sizes to overlap PCIe transfers with computation and keep GPUs near the roofline (Section 5.1 & Section 5.2). • Numerical Guardrail for Mixed-Precision Execution. To accelerate calculations, we use FP16/FP32 mixed precision in both stages, and for safe exploitation, we utilize asymmetric gradient clipping as a numerical guardrail. It caps gradient explosion while allowing shrinkage, preventing FP16 overflow and NaNs on power-law workloads and enabling stable Tensor Core execution(Section 5.3). • Comprehensive Evaluation at Billion Scale. We benchmark TIDE-MC against state-of-the-art matrix completion methods on 8 synthetic datasets and 7 real datasets. TIDE-MC delivers up to 11,647× speedup, 8.5× lower peak memory usage, and 99.7% lower reconstruction error on billion-scale matrices, whereas competing frameworks either time out or run out of memory (Section 6).

Limitations as Computing Systems. Although NN-centric and SVD-centric methods differ mathematically, modern implementations of both remain monolithic whole-matrix engines: they assume that the full sparse matrix, or large dense factors or Gram matrices, reside entirely in GPU highbandwidth memory (HBM) and organize computation around whole-matrix operators. This design is the main source of scalability problems, because it tightly couples the logical problem size to the physical memory capacity of the hardware. Once the matrix exceeds HBM capacity, such systems either hit an out-of-memory error or become stall-bound on PCIe under naive paging or offloading. Key Insight. Rather than proposing another optimizer, we seek a resource-reduction primitive for matrix completion. We introduce Two-Sided Interpolative Decomposition (TSID), which provides an exact-rank skeleton identity showing that a well-conditioned sampled template can serve as an algebraic anchor for reconstructing the ideal rank-𝑘 component, while practical error is governed by the rank-𝑘 residual, template recovery error, and extension error. TSID therefore reduces a large completion problem to a smaller resident subproblem and yields a template–extension structure. Our Design. Building on this abstraction, we use TSID as the foundation of a submatrix-centric scheme for matrix completion: instead of requiring the entire matrix to reside in HBM, we keep only a template submatrix resident on the GPU and stream the remaining columns (or rows) as chunks. This decouples the logical matrix size from the resident working set: the working set scales linearly with the matrix dimension, while the logical problem remains quadratic. On top of this scheme, we present TIDE-MC, a framework for matrix completion that exposes bounded-memory completion as a reusable service. TIDE-MC replaces the monolithic execution model of prior solvers with a two-stage bounded-memory pipeline [34, 35, 39, 40, 49, 61, 67] that balances compute density and memory bandwidth. Our specific contributions are: • TSID-Guided Abstraction. We elevate TSID from a mathematical result to a resource-reduction primitive for ma-

2

Background

Matrix completion aims to recover the missing entries of a partially observed matrix. Formally, we have the definition: Definition 1 (Matrix Completion). Given an incomplete matrix 𝑅 ∈ R𝑚×𝑛 with observed entries indexed by Ω, the matrix completion problem seeks to recover a full matrix 𝑅b ∈ R𝑚×𝑛 such that 𝑅b𝑖 𝑗 = 𝑅𝑖 𝑗 for all (𝑖, 𝑗) ∈ Ω, and 𝑅b accurately estimates the unobserved entries. Traditional matrix completion methods largely fall into two categories: nuclear norm minimization (NN)-centric methods and singular value decomposition (SVD)-centric methods in practice. 2.1

NN-Centric Methods

NN-centric approaches, such as nuclear norm minimization (NN) and Column-Selected Nuclear Norm minimization (CSNN), recover low-rank matrices by convex rank surrogates. Figure 1(a) illustrates a NN-centric pipeline: given 𝑅, the solver iteratively computes gradients over the observed entries and applies SVT-based shrinkage, while assuming that the full matrix resides in device memory. NN [52] relaxes rank minimization by replacing the rank with the sum of singular values. Given Ω ⊆ [𝑛 1 ] × [𝑛 2 ], NN 2

Gradient: 𝑋 𝑋← ←𝑋 𝑋− − 𝜂𝐿 𝜂𝐿 Shrink:

SVD: � = ���

Input Matrix R

Truncate:

��

𝑈Σ𝑉𝑇𝑇→→ 𝑈Σ𝑉 𝑈soft Σ,Σ,𝜏𝜏 𝑉𝑉𝑇𝑇 𝑈soft

Repeat until convergence Full matrix on device 𝑂(𝑛3 ) time, 𝑂 𝑛2 memory

flow, in which a low-rank decomposition 𝑅 ≈ 𝑈 Σ𝑉 𝑇 of the entire matrix is used to reconstruct the completed matrix or enable further refinement, again assuming full-matrix residency in device memory. CSPGD [31] accelerates CSNN by solving a proximal gradient formulation [65] on the sampled submatrix: 1 2 min R Ω𝐼 (𝑋 − 𝐶) 𝐹 + 𝜆 ∥𝑋 ∥ ∗, 𝑋 ∈R𝑛1 ×𝑑 2 with iterative updates  h1 2 i , R Ω𝐼 (𝑋 (𝑡 ) − 𝐶) 𝑋 (𝑡 +1) = SVT𝜂𝜆 𝑋 (𝑡 ) − 𝜂 ∇ 𝐹 2 where SVT𝜏 is singular value thresholding [7]. CSPGD avoids a full SDP [47] and benefits from GPU-accelerated tensor operations, but still requires repeated SVDs of 𝑛 1 × 𝑑 and sublinear convergence. CGM. The conditional gradient method [31], a Frank–Wolfe variant [15, 57], optimizes 𝑓 (𝑋 ) = 12 ∥R Ω𝐼 (𝑋 − 𝐶) ∥ 2𝐹 over a nuclear-norm ball. At each iteration, it computes the top singular pair (𝑢 (𝑡 ) , 𝑣 (𝑡 ) ) of ∇𝑓 (𝑋 (𝑡 ) ) and updates 𝑋 (𝑡 +1) =
𝑇 1−𝛾𝑡 𝑋 (𝑡 ) +𝛾𝑡 𝜏 𝑢 (𝑡 ) 𝑣 (𝑡 ) . GPU-accelerated SVD makes each step efficient, but the O (1/𝑡) rate and growing intermediate ranks can still be costly for high-precision recovery. NTKMC [48] derives a problem-specific Neural Tangent Kernel and casts completion as kernel ridge regression with
−1 closed-form predictor 𝑓 ∗ (𝑥) = 𝐾 (𝑥, 𝑋 ) 𝐾 (𝑋, 𝑋 ) + 𝜆𝐼 𝑌 , where 𝑋 and 𝑌 are coordinates and observed values in Ω. Training inverts a |Ω| × |Ω| Gram matrix (O (|Ω| 3 ) time, O (|Ω| 2 ) memory), and each prediction costs O (|Ω|). Dense kernel operations limit scalability and GPU efficiency. MVGMC [29] augments matrix completion with multi-view information, e.g., rating and item-similarity matrices, to handle sparse data and cold-start scenarios [38, 44]. The optimization couples rating reconstruction and similarity consistency:

SVD-based MC

Nuclear-norm Solver

Reconstructed Matrix �

(a) NN-centric.

Σk

= �k Σk ��� Reconstruct: � ss.

���

Full matrix on device

�(�2 �) time, �(�2 ) memory

(b) SVD-centric.

Figure 1. Workflow of NN-centric and SVD-centric matrix completion methods, both requiring the full matrix to reside in device memory during execution. solves 𝑅b = arg

min

𝑋 ∈R𝑛1 ×𝑛2

∥𝑋 ∥ ∗

s.t.

R Ω (𝑋 − 𝑅) = 0,

where ∥𝑋 ∥ ∗ denotes the nuclear norm and R Ω retains observed entries. This problem can be formulated as a semidefinite program [42, 45] and is often solved by singular value thresholding (SVT) [7], which iteratively applies gradients and singular-value shrinkage. While NN offers strong recovery guarantees, repeated full-matrix SVDs [7, 36] lead to high memory usage, limited GPU utilization, and poor scalability. CSNN [31, 32] reduces this cost by applying NN to a uniformly sampled subset of 𝑑 columns. Let 𝐼 ⊆ [𝑛 2 ] and 𝐶 = 𝑅:,𝐼 ∈ R𝑛1 ×𝑑 . CSNN solves 𝐶b = arg min ∥𝑋 ∥ ∗ 𝑋 ∈R𝑛1 ×𝑑

s.t.

R Ω𝐼 (𝑋 − 𝐶) = 0,

where Ω𝐼 = {(𝑖, 𝑗) ∈ Ω | 𝑗 ∈ 𝐼 }, and then reconstructs the full matrix via least squares [18]: 2 1 𝑍b = arg min R Ω (𝑅) − R Ω (𝐶b𝑍 ) , 𝑅b = 𝐶b 𝑍b. 𝑑 ×𝑛 𝐹 2 2 𝑍 ∈R

Restricting SDP computation [42, 45] to a tall-and-skinny block (𝑑 ≪ 𝑛 2 ) lowers runtime and memory, and enables GPU-accelerated conic solvers [68]. Limitations. From a systems viewpoint, NN and CSNN remain expensive and difficult to scale. NN incurs quadratic space and cubic time complexity, while CSNN still depends on dense subproblems with substantial memory overhead. Moreover, neither method provides a streamable or chunked execution model. Consequently, on large datasets such as Netflix and Yahoo! Music (Table 1), state-of-the-art NNcentric implementations either time out or exhaust GPU memory. 2.2

min ∥𝑅 − 𝑃𝑄𝑇 ∥ 𝐹 + 𝜆∥𝑃 ∥ 𝐹 + 𝛽 ∥𝑆 − 𝑄𝑄𝑇 ∥ 𝐹 , Ω

(1)

where ∥𝑆 − 𝑄𝑄𝑇 ∥ 𝐹 encourages item factors in 𝑄 to match the similarity matrix 𝑆, and 𝜆, 𝛽 are the regularization parameters. The similarity matrix is built from item attributes: Í 𝑆𝑖𝑖 ′ = 1𝐽 𝐽𝑗=1 I(𝑎𝑖 𝑗 = 𝑎𝑖 ′ 𝑗 ), where 𝐽 is the number of attributes, and I is the indicator function. Limitations. Despite differing mathematically, these systems behave as whole-matrix engines: they require the full matrix (or dense factors) to fit in device memory or rely on paging that is oblivious to low-rank structure. They lack a task-aware mechanism for memory-efficient execution, motivating our search for a resource-reduction primitive that enables a submatrix-centric design with low space cost.

SVD-Centric Methods

SVD-centric methods trade some accuracy for higher computational efficiency. Representative examples include columnselected proximal gradient descent (CSPGD), conditional gradient methods (CGM), NTK-based matrix completion (NTKMC), and multi-view graph-based matrix completion (MVGMC). Figure 1(b) illustrates a typical SVD-centric work-

3

Design Principle

As discussed in Section 2, traditional matrix completion systems face a fundamental scalability bottleneck. In real3

Table 2. Illustrative comparison between the full matrix and an arbitrary subset on a synthetic low-rank matrix completion workload. (𝛼 = 0.01) Metric Shape Top-𝑘 energy retained Subspace similarity to 𝑅 Resident working set

Full matrix 𝑅

Arbitrary subset 𝑅𝑠

𝑚 ×𝑛 1.00 1.00 high

𝛼 (𝑚 × 𝑛) 0.90 0.80 low

We write the ideal rank-𝑘 component as 𝑅𝑘 = 𝑈 1 Σ1𝑉1𝑇 . Select index sets 𝐼 ⊆ {1, . . . , 𝑚} and 𝐽 ⊆ {1, . . . , 𝑛}, each of size 𝑘, such that the ideal core 𝑅𝑘∗ = (𝑅𝑘 )𝐼,𝐽 = 𝑈 1 (𝐼, :) Σ1 𝑉1 (𝐽, :)𝑇 is nonsingular. Define the left and right interpolation matrices
−1
−1 𝑊𝐿 = 𝑈 1 𝑈 1 (𝐼, :) ∈ R𝑚×𝑘 , 𝑊𝑅 = 𝑉1 𝑉1 (𝐽, :) ∈ R𝑛×𝑘 .

world scenarios, the target matrix is typically very sparse and approximately low-rank [3, 53]. As matrix dimensions grow, storage and intermediate structures often scale as O (𝑛 2 ), while the associated computation ranges from O (𝑛 2 ) to O (𝑛 3 ). At scales of hundreds of thousands of rows or columns, even storing dense factors or Gram matrices can exceed GPU memory. The core systems challenge is thus to carry out large scale matrix completion without assuming that the full matrix fit in accelerator memory. A tempting approach is to treat matrix completion as a generic sparse computation and rely on paging or swapping: keep most data in CPU memory or on disk and fault in segments on demand. However, without semantic knowledge of the completion task, such mechanisms move arbitrary subsets of the matrix at the page granularity. They neither exploit low-rank structure nor respect the algebraic dependencies between columns; thus, they spend most of their time shuffling data over PCIe rather than computing. Therefore, we algebraically reduce the completion problem such that only a small, well-chosen part of the matrix needs to reside on the GPU. Observation. Classical interpolative decomposition (ID) results show that low-rank matrices can be approximated by subsets of their rows or columns. Table 2 provides empirical evidence that a suitably chosen submatrix can preserve the dominant low-rank structure of 𝑅, suggesting that recovery of the full matrix may be approximated by first recovering a smaller representative submatrix and then extending it. This motivates us to develop a Two-Sided Interpolative Decomposition (TSID) for matrix completion. Under mild assumptions, TSID states that if sufficiently informative subsets of rows and columns are selected from 𝑅, then the resulting submatrix 𝑅 ∗ both preserves the low-rank structure of 𝑅 and serves as an effective template; furthermore, recovering 𝑅 ∗ from the observed entries is, up to standard lowrank approximation error, equivalent to recovering 𝑅. 3.1

The superscript ∗ denotes the sampled template. In the ideal theorem below it refers to the fully specified rank-𝑘 core 𝑅𝑘∗ ; in the practical matrix-completion algorithm, the corresponding observed template is only partially observed and must be estimated from data. Theorem and Proof of TSID. The theorem below is an exact-rank skeleton/CUR identity. Its role is to justify why a well-conditioned template can serve as an algebraic anchor for bounded-memory extension. Theorem 1 (Ideal TSID skeleton identity). If 𝑅𝑘∗ = (𝑅𝑘 )𝐼,𝐽 is nonsingular, then the interpolation induced by the selected rows and columns exactly reconstructs the ideal rank-𝑘 component: 𝑅e𝑘 = 𝑊𝐿 𝑅 ∗ 𝑊 𝑇 = 𝑅𝑘 . 𝑘

Proof. By the definitions above, 𝑅e𝑘 = 𝑊𝐿 𝑅𝑘∗ 𝑊𝑅𝑇 = 𝑈 1 𝑈 1 (𝐼, :)

3.2

Σ2

𝑉2𝑇

From TSID to a Two-Stage Framework

Given the TSID reduction in Theorem 1, matrix completion is reorganized into two bounded-memory stages: (i) recovering a compact template submatrix and (ii) extending this template to reconstruct the full target matrix. In Stage I (template recovery), a sampled template submatrix 𝑅 ∗ is factorized using regularized matrix factorization (RMF). Let Ω ∗ ⊆ Ω denote the observed entries that fall in-

matrix, and let

#" # 𝑉1𝑇

−𝑇 𝑈 1 (𝐼, :)Σ1𝑉1 (𝐽, :)𝑇 𝑉1 (𝐽, :) 𝑉1𝑇

This proves the exact identity for the ideal rank-𝑘 component. If 𝑅 is not exactly rank 𝑘, then 𝑅 = 𝑅𝑘 + (𝑅 − 𝑅𝑘 ), and the identity above does not reconstruct the residual 𝑅 − 𝑅𝑘 . Therefore, TSID should be interpreted as a skeleton-based abstraction for reducing the resident problem size, not as an end-to-end proof of optimal sparse matrix completion. □

Notation and Setup. Let 𝑅 ∈ R𝑚×𝑛 denote the full target 0

−1

= 𝑈 1 Σ1𝑉1𝑇 = 𝑅𝑘 .

TSID Theoretical Basis

" Σ1   𝑅 = 𝑈 Σ 𝑉 𝑇 = 𝑈1 𝑈2 0

𝑅

For a matrix 𝑅 with nonzero residual 𝑅 − 𝑅𝑘 , this identity applies only to the ideal component 𝑅𝑘 . The resulting error depends on row/column sampling, core conditioning, template recovery, and extension error.

be its singular value decomposition, where Σ1 ∈ R𝑘 ×𝑘 contains the top 𝑘 singular values, 𝑈 1 ∈ R𝑚×𝑘 , and 𝑉1 ∈ R𝑛×𝑘 . 4

Figure 2. Overall workflow of TIDE-MC: submatrix sampling and recovery (Stage I) and chunk-based ALS reconstruction (Stage II). side the sampled template. The submatrix is parameterized as 𝑅b∗ ≈ 𝑃𝑄 ⊤ , with 𝑃 ∈ R𝑚×𝑘 and 𝑄 ∈ R𝛼𝑛×𝑘 , by solving ∑︁
2
min 𝑅 ∗𝑗𝑖 − 𝑝 ⊤𝑗 𝑞𝑖 + 𝜆 ∥𝑃 ∥ 2𝐹 + ∥𝑄 ∥ 2𝐹 , (2) 𝑃,𝑄

In practical sparse workloads, TSID should be viewed as a systems abstraction for reducing resident memory, not as a standalone guarantee of optimal matrix completion. 3.3

( 𝑗,𝑖 ) ∈Ω∗

where 𝑝 𝑗 ,𝑞𝑖 are rows of 𝑃,𝑄, and 𝜆 controls ℓ2 regularization. This objective is optimized by stochastic gradient descent. In Stage II (extension), the recovered template 𝑅b∗ ≈ 𝑃𝑄 ⊤ is kept fixed, and the remaining columns of the original matrix are reconstructed by solving regularized least-squares problems conditioned on this template. The objective decomposes over columns, which makes it amenable to chunked, pipelined execution under a fixed memory budget; the concrete execution engine is described in Sections 4.1 and 5. Practical Approximation Error. The identity above is lossless only for the ideal component 𝑅𝑘 . To see what changes in the practical setting, write 𝑅 = 𝑅𝑘 +𝐸, where 𝐸 is the residual outside the rank-𝑘 component. The sampled template of the full matrix is then 𝑅 ∗ = 𝑅𝐼,𝐽 = 𝑅𝑘∗ + 𝐸𝐼,𝐽 . If Stage I recovers an estimated template 𝑅b∗ and Stage II forms 𝑅b = 𝑊𝐿 𝑅b∗𝑊𝑅𝑇 , then

Complexity

Let 𝑅 be the partially observed input matrix and let Ω be its observed-entry set. For an 𝑛 × 𝑛 matrix with observed-entry density 𝛿 and template sampling rate 𝛼, we have |Ω| = 𝛿𝑛 2 and |Ω ∗ | ≈ 𝛼𝛿𝑛 2 . Let 𝑘 be the latent rank, 𝑡 1, 𝑡 2 the iteration counts of Stage I and Stage II, and Δ𝑐 the chunk width. The two stages require 𝑇I = 𝑂 (𝑡 1𝑘 |Ω∗ |) = 𝑂 (𝑡 1𝑘𝛼𝛿𝑛 2 ), 𝑇II = 𝑂 (𝑡 2𝑘 |Ω|) = 𝑂 (𝑡 2𝑘𝛿𝑛 2 ). Thus the total arithmetic work is 𝑂 (𝑘𝛿𝑛 2 (𝑡 2 + 𝛼𝑡 1 )). This is still linear in the observed entries; TIDE-MC does not claim sublinear input processing. Its reduction is in peak resident GPU memory: at any time, the system keeps only the template factors and one active chunk, 𝑆 resident = 𝑂 (𝑘 (1 + 𝛼)𝑛 + |ΩΔ𝑐 | + 𝑘Δ𝑐) , where ΩΔ𝑐 is the observed-entry set in the active chunk. Therefore, the resident memory is controlled by the template size and chunk width rather than by the full matrix.

𝑅b − 𝑅 = 𝑊𝐿 (𝑅b∗ − 𝑅𝑘∗ )𝑊𝑅𝑇 − 𝐸,   ∥ 𝑅b − 𝑅∥ 𝐹 ≤ ∥𝐸 ∥ 𝐹 + ∥𝑊𝐿 ∥ 2 ∥𝑊𝑅 ∥ 2 ∥ 𝑅b∗ − 𝑅 ∗ ∥ 𝐹 + ∥𝐸𝐼,𝐽 ∥ 𝐹 . Thus TSID does not imply error-free submatrix reconstruction in general. It shows lossless recovery only under the exact-rank, fully observed, well-conditioned-template assumption. In practical sparse matrix completion, the reconstruction quality depends on the rank-𝑘 residual, the conditioning of the selected rows and columns, the sampling pattern, and the template recovery error. This formulation also explains why the implementation does not simply set 𝑃 = 𝑅 ∗ and 𝑄 = 𝐼 . That choice is only a degenerate factorization when the 𝑘 × 𝑘 template is fully observed, noise-free, and already known. In sparse matrix completion, the sampled template 𝑅 ∗ is only partially observed and may be noisy. Stage I therefore uses regularized matrix factorization to estimate a low-rank template representation 𝑅b∗ ≈ 𝑃𝑄𝑇 from the observed entries. Stage II then uses the recovered template representation in a template-conditioned least-squares extension. Thus RMF/ALS-style updates are optimization mechanisms for the practical sparse setting, not mathematical requirements of the exact skeleton identity.

4

TIDE-MC Framework Overview

4.1

Overview

Section 3 introduced a TSID-guided Template-Extension decomposition of matrix completion at the algorithmic level. Figure 2 shows how TIDE-MC realizes this decomposition as a bounded-memory execution workflow. The key idea is to break the dependence of the resident working set on the logical matrix size: at any time, the GPU holds only a small submatrix or a streaming chunk, rather than the full matrix. The workflow proceeds in two phases. TIDE-MC first samples a subset of columns from the input sparse matrix 𝑅 to form a compact submatrix 𝑅 ∗ , and recovers 𝑅 ∗ as a low-rank template. This template captures the dominant structure of 𝑅 while remaining small enough to stay resident in GPU memory. TIDE-MC then treats the recovered template as a fixed anchor and reconstructs the remaining columns chunk by chunk, streaming them through the GPU under a fixed memory budget. In this way, TSID provides the algebraic basis for reducing full-matrix completion to template recovery 5

Table 3. Submatrix sampling overhead.

followed by template-guided extension, while the system realizes this reduction as a submatrix-centric dataflow. This execution model has three structural properties. First, it bounds the resident working set to the template and one chunk at a time, enabling completion beyond GPU memory capacity. Second, it exposes both phases through a common template-guided workflow that can be mapped to GPU linearalgebra kernels. Third, it separates the Template-Extension abstraction from the stage-specific execution mechanisms, allowing each phase to be optimized independently. 4.2

(a) Synthetic 10K×10K. Time(s) ΔAcc. 0.0073 0 33.7 0.00017% 9.86 0.0008%

Strategy Random Rank-Inc. Similarity

5

(b) Netflix Prize. Strategy Random Rank-Inc. Similarity

Trial(s) Total(s) 0.002 0.007 50.34 151 39.93 119.8

TIDE-MC System

The Template–Extension execution model in Section 4.1 is instantiated in TIDE-MC as a two-module execution engine. The Template module is a conflict-aware, hierarchical update engine that recovers the TSID template under massive GPU parallelism, and the Extension module is a rooflinedriven, pipelined dataflow that reconstructs the full matrix under a bounded GPU memory budget. Both modules share a common precision-management layer, i.e., an asymmetric gradient guardrail that makes aggressive mixed-precision (FP16/FP32) execution numerically stable.

Main Challenges

Although the workflow in Section 4.1 establishes a boundedmemory execution path, a straightforward implementation is still insufficient for performance. A direct TSID-guided realization would sample a submatrix 𝑅 ∗ , apply an off-theshelf SGD-based RMF solver to recover a template, and then run a standard ALS solver to reconstruct the full matrix conditioned on that template. In practice, this naive realization still encounters three major bottlenecks. Challenge I: Atomic Synchronization Overhead (Stage I). During template recovery, RMF must frequently and efficiently update shared latent factors. A naive GPU implementation assigns threads to individual ratings ( 𝑗, 𝑖) and updates user and item factors with atomic operations. On power-law datasets, popular items become hotspots, and thousands of updates contend for the same memory locations. In our profiling, this naive RMF kernel spends over 40% of its cycles stalled on memory locks, making Stage I synchronizationbound and unable to fully utilize the GPU. Challenge II: Stall-Bound Streaming Pipeline (Stage II). Global reconstruction requires streaming large portions of the matrix that do not fit in GPU memory. A naive implementation performs a sequence of blocking steps: copy a chunk over PCIe, compute ALS updates on it, copy it back, and then proceed to the next chunk. On our Stage II matrix completion microbenchmark, This exposes PCIe latency and limited bandwidth directly on the critical path: GPU SMs spend up to 65% of runtime idle, waiting for data. Without a dedicated latency-hiding dataflow, Stage II quickly becomes strictly bandwidth-bound. Challenge III: Numerical Instability Under Mixed Precision. To fully exploit Tensor Cores, both RMF and ALS should run in FP16 or mixed precision. However, naive SGDbased RMF is numerically brittle on sparse, skewed data. The limited dynamic range of FP16, combined with large gradients in early epochs, leads to overflow and divergence. Falling back to FP32 erodes the benefits of specialized hardware, leaving a gap between efficiency and stability. These challenges motivate the system mechanisms of TIDE-MC. The next section presents how we address them with a conflict-aware synchronization engine, a pipelined chunked reconstruction dataflow, and a numerical guardrail for stable mixed-precision execution.

5.1

Conflict-Free Synchronization Engine

Stage I is the compute-intensive core of TIDE-MC: it repeatedly applies regularized matrix factorization (RMF) to the sampled submatrix 𝑅 ∗ until convergence. This module implements the Template phase of the Template–Extension execution model in Section 4.1. From a systems perspective, the main difficulty is designing a synchronization engine that enforces these updates on GPUs without global write conflicts or excessive atomic operations while keeping Tensor Cores highly utilized. SGD Updates for RMF. We adopt the standard RMF objective in Equation 2 and optimize it with stochastic gradient descent. For an observed entry ( 𝑗, 𝑖) ∈ Ω, the per-entry gradients are:
for 𝑝 𝑗 : ∇𝑓𝑝 𝑗 (𝑝 𝑗 , 𝑞𝑖 ) = −2 𝑅 𝑗𝑖 − 𝑝 𝑗 𝑞𝑖⊤ 𝑞𝑖 + 2𝜆𝑝 𝑗 , (3)
for 𝑞𝑖 : ∇𝑓𝑞𝑖 (𝑝 𝑗 , 𝑞𝑖 ) = −2 𝑅 𝑗𝑖 − 𝑝 𝑗 𝑞𝑖⊤ 𝑝 𝑗 + 2𝜆𝑞𝑖 . (4) In our multi-view extension, the item-factor update incorporates the similarity regularizer from Equation 1, yielding 𝑛 ∑︁
∇𝑓𝑞𝑖 (𝑝 𝑗 , 𝑞𝑖 ) = −2 (𝑅 𝑗𝑖 − 𝑝 ⊤𝑗 𝑞𝑖 ) 𝑝 𝑗 + 2𝛽 𝑆𝑖𝑖 ′ 𝑞𝑖 − 𝑆𝑖𝑖 ′ 𝑞𝑖 ′ , 𝑖 ′ =1

where 𝑖 ′ 𝑆 (𝑞𝑖 −𝑞 ) enforces similarity consistency among items. These formulas are standard from an optimization viewpoint; rest of this subsection focuses on how TIDE-MC turns them into a high-throughput, conflict-aware GPU kernel. Sampling the TSID Template Submatrix. Stage I does not operate on the full matrix. Instead, it only sees a submatrix 𝑅 ∗ obtained by sampling columns (or rows) from the original 𝑅: 𝑅 ∗ = {𝑅𝑖 𝑗 | 𝑗 ∈ C, (𝑖, 𝑗) ∈ Ω}, where C is the sampled column set with |C| = 𝑑 = ⌈𝛼𝑛⌉, and 𝑘 is the latent rank with 𝑘 ≪ 𝑑 ≤ 𝑛. Under the TSID view, a well-sampled template provides the low-rank anchor used by Stage II, so the sampling strategy must be both cheap and robust. We therefore use simple random column sampling. Table 3 comÍ

6

𝑖𝑖 ′

𝑖′

pares this choice with two heavier alternatives, rank increment and similarity-based selection [62]. Random sampling has negligible overhead on both workloads, while the heavier strategies are orders of magnitude slower and improve accuracy by less than 10−3 on the synthetic workload. Thus, random sampling keeps the Stage I startup cost negligible without measurable accuracy loss in our setting. Conflict-Free Parallel SGD. Naively parallelizing Equations 3–4 on GPUs creates severe write conflicts on shared item factors 𝑞𝑖 : many threads process different users 𝑗 that rated the same item 𝑖, and all attempt to update 𝑞𝑖 concurrently. Using global-memory atomics (e.g., atomicAdd) for every gradient component serializes these updates and saturates the memory controller; ignoring conflicts, in turn, harms convergence on power-law data with “hot” items. TIDE-MC resolves this by turning the per-entry updates into a hierarchical reduction over the GPU memory hierarchy. Each thread first computes a local gradient contribution for its own copy of 𝑞𝑖 in registers. Within a warp, the local contributions that target the same 𝑖 are combined using shuffle instructions, producing one partial sum per warp. Warp partial sums are then accumulated in shared memory to obtain a block-level delta for each 𝑞𝑖 . Only this aggregated delta is finally applied to global memory, and we use at most one atomic operation per block to commit it. Conflicts are absorbed at the register and shared-memory levels, and global memory only sees a small number of coarse-grained, commutative updates. Formally, letting ∇𝑞𝑖 𝑓 (𝑝 𝑗 , 𝑞𝑖 ) denote the per-thread gra
dient, the aggregated update is 𝑞𝑖′ =

Í𝑘

𝑗 =1

𝑞𝑖 +𝜂 ∇𝑞𝑖 𝑓 (𝑝 𝑗 ,𝑞𝑖 ) 𝑛

Stage I Execution Sampled Submatrix R* Conflict-Free Parallel SGD

𝑝 ′𝑗

= 𝑝 𝑗 + 𝜂 ∇𝑝 𝑗 𝑓

𝑝 𝑗 , 𝑞𝑖 , 𝑞𝑖′

CUDA Core

CUDA Core Tensor Core Large dense RMF updates

GEMM GPU

Figure 3. Stage I as a conflict-free engine: sparse RMF updates are aggregated and mapped onto Tensor Core tiles. the dominant computations as matrix multiplications Δ𝐽 = MatMul(𝑃, 𝑄). Tensor Cores perform FP16 multiplications with FP32 accumulation, providing high throughput without sacrificing accuracy. We implement this kernel in PyTorch using Automatic Mixed Precision (AMP), which manages loss scaling and precision transitions, while the sparse access pattern is confined to lightweight gather/scatter around the GEMMs. Figure 3 illustrates the Stage I workflow as a GPU-resident aggregation-and-compute engine: sparse RMF accesses are handled by efficient sparse-matrix primitives, while the dominant factor computations are mapped to dense Tensor Core kernels for Δ𝐽 = MatMul(𝑃, 𝑄). Stage I Execution Example. Given a sparse input matrix 𝑅 ∈ R𝑛1 ×𝑛2 , Stage I first samples subcolumns according to 𝛼 to form 𝑅 ∗ ∈ R𝑛1 ×𝑑 with 𝑑 = 𝛼𝑛 2 . We then initialize 𝑃 ∈ R𝑛1 ×𝑘 and 𝑄 ∈ R𝑑 ×𝑘 and run the GPU-parallel SGD kernel based on (5) until 𝑅 ∗ ≈ 𝑃𝑄 ⊤ converges. This kernel is fully GPU-resident and is the main consumer of GPU FLOPs in Stage I; hierarchical aggregation keeps global memory traffic and atomic operations low, allowing Tensor Cores to run in a compute-bound regime.

=

𝑗 =1 ∇𝑞𝑖 𝑓 (𝑝 𝑗 ,𝑞𝑖 ) . The corresponding parallel update rules 𝑛

PCIe

CUDA Core

Í𝑘

𝑞𝑖 +𝜂 are:

CPU

∇𝑞𝑖 𝑓 𝑝 𝑗 , 𝑞𝑖 , (5) = 𝑞𝑖 + 𝜂 𝑛

where 1/𝑛 factor normalizes the step size for the number of contributing threads. We refer to this scheme as conflict-free hierarchical synchronization: fine-grained, high-contention updates never directly contend in global memory, and the remaining global atomics operate only on aggregated deltas. Regularizing Sparse RMF Updates for Tensor Cores. The SGD updates in RMF are irregular in the sense that each rating touches a different pair of user and item factors, so threads access non-contiguous rows of 𝑃 and 𝑄 and issue fine-grained scatter/gather operations instead of reading or writing contiguous tiles. Even with conflict-aware aggregation, such sparse, index-based updates would keep the kernel memory-bound unless we can reshape them into dense matrix operations that fully utilize Tensor Cores. TIDE-MC therefore batches RMF updates into small matrix–matrix multiplications over the low-rank factors rather than over the sparse rating matrix itself. In each SGD iteration, we group a mini-batch of users and items into dense tiles and express

5.2

Roofline-Driven Pipelined Dataflow

Once Stage I has recovered the TSID template 𝑅b∗ and its factors, Stage II must reconstruct the full matrix under tight GPU-memory constraints. From a systems perspective, Stage II is a pipelined dataflow engine whose goal is to stream chunks through the GPU at a rate determined by a roofline-style balance between computation and PCIe bandwidth. Chunk-Based ALS Kernel. We use a chunked ALS-style reconstruction. The full matrix 𝑅b ∈ R𝑚×𝑛 is recovered from a restored submatrix 𝑅b∗ ∈ R𝑚×𝛼𝑛 (the template) using the objective in Equation (6): 2 ∑︁  𝑅𝑖 𝑗 − 𝑅b∗𝑖 𝐶 𝑗 + 𝜆∥𝐶 ∥ 2𝐹 , (6) L (𝐶) = (𝑖,𝑗 ) ∈Ω

where 𝐶 ∈ R𝛼𝑛×𝑛 is the extension coefficient matrix whose 7

𝑗-th column 𝐶 𝑗 represents column 𝑗 in the fixed template basis 𝑅b∗ , and Ω indexes observed entries. We update 𝐶 using gradient  descent with the Adam optimizer [28]: 𝐶 ← 𝐶 − 𝜂 · Adam ∇𝐶 L (𝐶) , where 𝜂 is the learning rate. To avoid materializing all columns of 𝐶 on the GPU, we partition 𝑅 (and thus 𝐶) into column chunks of width Δ𝑐, and process each chunk independently given the fixed template. Double-Buffered Pipelining and Iteration Cost. We implement Stage II as a double-buffered pipeline that overlaps GPU computation with data transfers (see Figure 4). While the GPU processes the ALS gradient step on chunk 𝑖, chunk 𝑖 + 1 is asynchronously prefetched via PCIe, and the results of chunk 𝑖 − 1 are written back to CPU memory. Semantic Memory Management. Chunking provides a tunable memory footprint, which allows TIDE-MC to operate on matrices that are far larger than GPU memory. The full matrix and 𝐶 reside in CPU memory, and only one or two chunks are resident on the GPU at a time. The runtime uses multiple CUDA streams to overlap three activities: (1) hostto-device prefetch of the next chunk into a pre-allocated device buffer, (2) ALS computation on the current chunk, and (3) device-to-host flush of the previous chunk into a ring buffer allocated in pinned host memory. In contrast, the template 𝑅b∗ and its factors (𝑃, 𝑄) are read-mostly and reused across all chunks, so they are pinned in GPU HBM in Stage II. This semantic split, e.g., template vs. streaming chunks, treats different data structures according to their access patterns rather than managing all memory uniformly. Tensor Core Acceleration in Stage II. The ALS update repeatedly multiplies 𝑅b∗ with column blocks of 𝐶. We express these operations as batched matrix multiplications that run on Tensor Cores in mixed precision. PyTorch AMP handles FP16/FP32 transitions and dynamic loss scaling. Compared to a pure CUDA-core implementation, this reduces the ALS compute time by up to an order of magnitude on our largest datasets without loss in reconstruction accuracy. Stage II Execution Example. As shown in Figure 4, after Stage I the submatrix 𝑅 ∗ is fixed and serves as the template for Stage II. To reconstruct the full matrix, we partition the original matrix into blocks that fit in GPU memory, then load each block onto the GPU in turn. For each block, we run ALS with 𝑅 ∗ fixed while updating the factor matrix 𝐶 ∈ R𝛼𝑛×𝑛 as in Equation 6, repeating until target accuracy or a maximum iteration count is reached. Once a block is completed, it is offloaded to CPU memory and the next block starts, as illustrated in Figure 4. 5.3

R

1 i

i i

i + 1 + i + 1

Figure 4. Stage II as a roofline-driven pipelined dataflow engine: ALS compute on chunk 𝑖 overlaps with prefetch of 𝑖 + 1 and flush of 𝑖 − 1. dient clipping as a numerical stabilization guardrail that enables aggressive mixed-precision computation without sacrificing the convergence of the whole algorithm. Asymmetric Gradient Clipping. To ensure stable convergence and enhance accuracy, we refine the SGD updates in Equation 5 by rewriting them as:  
𝑝 ′𝑗 = 𝑝 𝑗 + 𝜂 min max 𝜉, −∇𝑝 𝑗 𝑓 (𝑝 𝑗 , 𝑞𝑖 ) , ∇𝑝 𝑗 𝑓 (𝑝 𝑗 , 𝑞𝑖 ) ,  
𝜂 min max 𝜉, −∇𝑞𝑖 𝑓 (𝑝 𝑗 , 𝑞𝑖 ) , ∇𝑞𝑖 𝑓 (𝑝 𝑗 , 𝑞𝑖 ) , 𝑞𝑖′ = 𝑞𝑖 + 𝑛 where 𝜉 > 0 is the clipping threshold. For each coordinate, positive gradients larger than 𝜉 are truncated to 𝜉, while negative gradients are left unchanged. This asymmetric clipping lets parameters shrink freely but caps explosive growth. Formally, if ∇𝑓 > 𝜉, the update is forced to use 𝜉; multiplying by the step size 𝜂 then guarantees ∥Δ𝑝 ∥ ∞ ≤ 𝜂 𝜉 per coordinate, preventing exploding updates and numerical overflow. In FP16, where the representable range is limited, this bound is critical: it ensures that even in early epochs with a large reconstruction error, the RMF kernel stays within a safe dynamic range. Combined with AMP’s loss scaling, this guardrail allows TIDE-MC to double effective memory bandwidth and compute throughput compared to a pure FP32 implementation while maintaining stable convergence. Effect on Convergence. Figure 5 compares the loss over iterations with and without clipping. Without clipping (blue), the loss oscillates and fails to converge, as large gradients repeatedly overshoot the optimum. With clipping (red), the loss descends smoothly and stabilizes near the minimum. In our end-to-end experiments, clipping consistently improves accuracy and reduces the number of iterations needed to reach a target reconstruction error, with negligible runtime overhead. From a system viewpoint, this means the mixedprecision engine can be run in a “fail-safe” mode: even aggressive TSID sampling and high learning rates do not cause catastrophic divergence.

Asymmetric Gradient Guardrail

Although SGD-based RMF delivers high throughput, its Stage I updates can become numerically unstable or even diverge when gradients grow large and unbounded, especially under FP16 execution on Tensor Cores. Instead of treating this purely as an “optimization trick”, TIDE-MC integrates gra-

6

Evaluation

Implementation Details. TIDE-MC is implemented in about 8

Without Clipping

Loss

1 0.8 0.6 0.4

0

20

40

60

80 100

Iteration

Table 4. Datasets analysis.

With Clipping

1 0.8 0.6 0.4 0.2 0

0

20

40

60

Real-World Dataset Image Construction [31] Noise Clear [31] Last.fm 2K [1] MovieLens 25M [22] MovieLens Latest 33M [22] Netflix Prize [43] Yahoo! Music [64]

80 100

Iteration

Figure 5. Comparison of loss over iterations with and without gradient clipping (threshold 𝜉). 7,000 lines of Python. Using PyTorch’s device management, inputs are automatically moved to the target device, so all tensor operations in parallel SGD and ALS routines run on GPU. The code is organized as a modular, object-oriented framework that encapsulates submatrix sampling, iterative factorization, and gradient clipping, with logging and visualization utilities for debugging and performance analysis. For Tensor Core acceleration, we use PyTorch [46] with CUDA 12.4 [60], and run critical batched matrix multiplications in mixed precision to fully exploit NVIDIA Tensor Cores. 6.1

Synthetic Dataset

Users

Items

Ratings

Size

1080 771 1.9K 162K 331K 480,189 200,000 Row 1K 1K 3K 3K 10K 50K 100K 100K

1920 621 17K 59K 83K 17,770 136,000 Column 1K 3K 3K 10K 100K 100K 100K 150K

– 25M 33M 100,480,507 262,810,175 Sparsity

3.2MB 44KB 123MB 35.6GB 102.4GB 31GB 101GB Size 3.8MB 11MB 34MB 114MB 3.72GB 18.6GB 37.2GB 55.87GB

80%

resentation is space-efficient for sparse matrices. Baselines and Variants. We compare TIDE-MC against five external state-of-the-art (SOTA) matrix-completion solvers and several internal variants. For NN-centric methods, we use CSNN as the nuclear-norm SOTA baseline. For SVDcentric, kernel-based, and graph-based methods, we use CSPGD, CGM, NTK, and MVGMC as SOTA baselines on their respective tasks. Implementations of CSNN, CSPGD, and CGM follow the open-source code from [31]; NTK uses the implementation from [48]; MVGMC follows [29]. These five methods appear in Table 5 as the external SOTA baselines. Accuracy Evaluation Metrics. To measure accuracy, we adopt the normalized mean absolute error (NMAE) [12], which is widely used in matrix completion studies [31, 55, 59, 63]. Given a held-out test set of observed entries Ωtest , Í 1 b NMAE is defined as |Ωtest | (𝑚max (𝑖,𝑗 ) ∈Ωtest |𝑅𝑖 𝑗 − 𝑅𝑖 𝑗 |, −𝑚 min ) where 𝑚 max and 𝑚 min are the maximum and minimum values in 𝑅. Lower NMAE indicates higher reconstruction accuracy. Parameter Selection. We typically set the subcolumn selection ratio (i.e., the percentage of total columns selected for each submatrix) to 0.5 for small datasets to ensure high accuracy and to 0.1 for large-scale datasets to balance memory usage and computational efficiency. The number of chunks is generally chosen to be between 10 and 30 to achieve optimal performance (detailed in Section 6.8).

Experimental Setup

Test Platform. Our experiments are conducted on a server equipped with two Intel(R) Xeon(R) Gold 6426Y processors with 512GB of CPU RAM and 32 threads, running Ubuntu 20.04.6 LTS with CUDA 12.4. All methods that require a graphics processing unit are executed on four NVIDIA A100 GPUs, each equipped with 40GB of HBM. Datasets. We evaluate on both synthetic and real-world datasets (Table 4). Synthetic datasets emulate the high sparsity of real recommenders, with 80% of user–item entries missing. We generate each synthetic matrix as the product of two randomly sampled low-rank factor matrices, 𝑃𝑄 ⊤ , and then randomly mask entries to the target sparsity, ensuring that the underlying complete matrix is low-rank. They range from 1K × 1K to 100K × 150K and are used to study behavior under varying sparsity and matrix sizes. For real-world workloads, we include both image and recommendation datasets. Image Construction and Noise Clear support an image reconstruction and denoising case study, with qualitative results shown in Figure 6; both are derived from grayscale bridge images from the public PxHere repository, with pixels treated as matrix entries [31]. Last.fm 2K is a music recommendation dataset containing user–artist interactions and auxiliary tag information [1]. MovieLens 25M and MovieLens Latest 33M provide large-scale user–movie rating matrices from the GroupLens benchmark suite [22]. Netflix Prize contains 100,480,507 ratings from 480,189 users on 17,770 movies [43]. Yahoo! Music, used in the 2011 KDD Cup, has 262,810,175 ratings from 200,000 users on 136,000 music items and is a standard benchmark for large-scale recommendation [64]. Data Model. Input matrices are represented in the coordinate (COO) format [13], which stores only the row indices, column indices, and values of non-zero elements. This rep-

6.2

Overall Performance

Table 5 compares TIDE-MC with five SOTA baselines (CSNN, CSPGD, CGM, MVGMC, and NTK) on diverse matrix completion tasks using NMAE and end-to-end runtime. TIDE-MC consistently achieves the best overall efficiency while maintaining strong accuracy. On mid-sized matrices, it delivers up to 11,647× speedup over NN-centric methods with comparable accuracy. Against SVD-centric methods, it achieves up to 1,377× speedup and up to 99.7% lower NMAE. For accuracy, on a 1K×3K matrix, TIDE-MC achieves an NMAE of 0.00075, improving over CSPGD and CGM by 46.4% and 73.1%, respectively. On a 3K×3K matrix, it attains an NMAE of 0.0027, 9

Table 5. End-to-end performance comparison in terms of accuracy (NMAE) and runtime (seconds). The symbol ✗ denotes unavailable results due to excessive runtime or missing implementations. NTK requires a strict training setup, making it unavailable on most datasets. We set 𝛼 = 0.1 for all experiments below. Dataset

NMAE

CSNN

CSPGD Time (s)

CGM Time (s)

NMAE

0.0028

532.00 (202x)

0.3662

7.98 (3.03x)

✗

✗

0.0008

2.63

0.8419

487.00 (40.5x)

0.3953

34.55(2.87x)

✗

✗

0.0027

12.02

145.00 (5.8x)

0.8438

579.00 (23.1x)

0.3450

141.93 (5.67x)

✗

✗

0.0019

25.00

7.30K (1377x)

0.0672

90.00 (17x)

0.2447

19.23 (3.62x)

0.1100

21.00 (3.96x)

0.0454

5.31

0.1899 OOM OOM OOM OOM

1.24 (17x) OOM OOM OOM OOM

0.1900 OOM OOM OOM OOM

3.52 (3.62x) OOM OOM OOM OOM

0.4530 ✗ ✗ ✗ ✗

0.85 (3.96x) ✗ ✗ ✗ ✗

0.1496 0.1352 0.1385 0.1382 0.1473

0.77 2.88K 5.04K 53.91 1.01K

Time (s)

NMAE

1K×3K

1.08 × 10-6

21.40K (8136x)

0.0014

51.00 (19.3x)

3K×3K

0.0013

140.00K (11647x)

0.0029

116.00 (9.65x)

3K×10K

0.0020

260.00K (10400x)

0.0020

Image Construction

0.0668

4.62K (872x)

0.0474

Last.fm 2K MovieLens 25M MovieLens Latest 33M Netflix Prize Yahoo! Music

0.4134 OOM OOM OOM OOM

495.40 (643x) OOM OOM OOM OOM

0.2281 OOM OOM OOM OOM

43.90 OOM OOM OOM OOM

NMAE

nearly matching CSNN and outperforming CGM by 99.7%. On image reconstruction, it achieves an NMAE of 0.04543, improving over CSNN and CGM by 31.9% and 32.4%. These gains mainly come from gradient clipping, which stabilizes training and reduces reconstruction error. Real-World Large Dataset Testing. On MovieLens (25M and Latest 33M), Netflix and Yahoo! Music datasets, CSNN, CSPGD, CGM, and MVGMC fail with out-of-memory errors due to full-matrix execution. NTK is not reported because its required training setup is unavailable for these datasets. In contrast, TIDE-MC completes Netflix with NMAE 0.1382 in 53.91s and Yahoo! Music with NMAE 0.1473 in 1.01Ks. Although Yahoo! Music and MovieLens Latest 33M have comparable storage sizes, their runtimes differ substantially, indicating that end-to-end time is affected not only by data size but also by matrix shape, sparsity pattern, and chunk-level reconstruction cost. On Last.fm 2K, TIDE-MC achieves the best NMAE of 0.1496 and finishes in 0.77s, outperforming the slowest highlighted baseline CSNN by 643×. On MovieLens 25M and MovieLens Latest 33M, all baselines run out of memory or are unavailable, while TIDE-MC completes them with NMAE 0.1352 and 0.1385, respectively. Entries marked with ✗ in Table 5 denote such failures. Although NTK attains low NMAE on some small image tasks, its scalability is limited. Overall, TIDE-MC combines strong accuracy with substantially better efficiency and scalability. Noise Clear. To complement the quantitative results, we evaluate TIDE-MC on a real-world image denoising task. We compare against CSPGD [31], the strongest baseline in Table 5. Figure 6 shows the original image and the reconstructions produced by CSPGD and TIDE-MC. TIDE-MC removes noise more effectively and produces a clearer reconstruction, providing evidence of improved reconstruction quality. 6.3

MVGMC Time (s)

(a) Origin image.

NMAE

(b) CSPGD.

NTK Time (s)

TIDE-MC (Ours) NMAE Time (s)

(c) TIDE-MC (Ours).

Figure 6. Qualitative comparison for image denoising. Table 6. Profiled GPU memory usage in GB. TIDE-MC uses the chunk-based pipeline; > 100 indicates that the method exhausts the largest available GPU-memory budget (more than 100GB) and fails with OOM. Dataset 1K × 3K 3K × 3K 3K × 10K Last.fm 2K MovieLens 25M MovieLens 33M Netflix Prize Yahoo! Music

CSNN 0.84 2.15 6.74 4.2 > 100 > 100 > 100 > 100

CSPGD 0.80 0.84 0.90 0.95 > 100 > 100 > 100 > 100

CGM 0.83 1.44 4.95 3.4 > 100 > 100 > 100 > 100

MVGMC 0.76 0.78 3.22 1.6 > 100 > 100 > 100 > 100

TIDE-MC 0.30 0.47 0.79 0.85 27 38 22 23

uses 0.85GB, lower than all baselines. On larger real-world workloads, CSNN, CSPGD, CGM, and MVGMC exceed the largest available GPU-memory budget in our testbed and fail with OOM, whereas TIDE-MC completes within 27GB on MovieLens 25M, 38GB on MovieLens Latest 33M, 22GB on Netflix, and 23GB on Yahoo! Music. Overall, the chunkbased pipeline keeps the resident working set within a single 40GB A100 GPU and enables larger-scale completion than full-matrix baselines. 6.4

Efficiency of TIDE-MC’s Optimizations

We conduct an ablation study (Figure 7) to evaluate each optimization in TIDE-MC on matrix sizes from 1K×1K to 3K×10K, plus an additional wide-matrix microbenchmark of 30×10K. Across all sizes, we compare the base RMF+ALS version and different optimization combinations in execution time (left) and NMAE (right). Efficiency of Opt1. Opt1, the conflict-free synchronization optimization in Stage I (Section 5.1), consistently reduces

Memory Usage Comparison

Table 6 reports the profiled GPU memory usage on synthetic and real-world recommendation workloads. On synthetic matrices, TIDE-MC keeps the resident footprint below 0.8GB, while baselines require up to 6.74GB. On the 3K × 10K, TIDE-MC reduces memory usage from 6.74GB for CSNN and 4.95GB for CGM to only 0.79GB. On Last.fm 2K, TIDE-MC 10

Base + Opt1 Base + Opt1 + Opt2(TIDE-MC)

NMAE

Time (s)

10

0

Time (s)

1.5

20

10−2

Time (s)

CUDA Core Tensor Core

200

Time (s)

4 GPU

100 0

0 16K ×

3 6 8 16K 2K × 32K 4K × 64K 0K × 80K

Matrix Size

10K ×

1 5 1 100K 0K × 100K 00K × 100K00K × 150K

Matrix Size

(a) Tensor Core vs. CUDA Core. (b) Evaluation on multiple GPUs.

Figure 8. Performance comparison under different hardware configurations.

6.6

execution time across all matrix sizes, with larger benefits on bigger matrices. For example, it cuts runtime on a 1K × 3K matrix from 8.3s to 4.71s, with up to a 43% reduction overall. Efficiency of Opt2. Opt2, the roofline-driven pipelined dataflow in Stage II (Section 5.2), consistently outperforms the base version, reducing runtime on a 30 × 10K matrix from 6.5s to 2.2s (65%). Its chunk-based parallelism is especially effective when rows are much fewer than columns, making it well suited to wide real-world matrices. Efficiency of Opt3. Opt3, the asymmetric gradient guardrail in Section 5.3, improves convergence and accuracy. Unless otherwise specified, we set the clipping threshold to 𝜉 = 0.1 in all Opt3 experiments. As shown in Figure 7, clipping reduces NMAE from 2.0 to 1.8 at the first iteration and from 0.423 to 0.023 at the 120th iteration (94% relative improvement). It also reaches NMAE=0.1 in about 60 iterations versus 100 without clipping, giving a 40% speedup. Efficiency of Tensor Cores. Figure 8 compares the runtime with and without Tensor Core acceleration. We disable Tensor Cores with TF_DISABLE_CUBLAS_ TENSOR_OP_MATH=1, forcing all kernels onto CUDA cores. Tensor Cores consistently run faster, with the gap widening for larger matrices: for a (16K × 16K) matrix, the runtime drops from 30.7s (CUDA cores) to 22.18s (Tensor Cores), a 27.7% reduction that highlights the advantage for matrix computations. 6.5

100

Accuracy

4

400

3

300

2

200

1 10K×100K

50K×100K 100K×100K 100K×150K Matrix Size

0

runtime and accuracy under the same submatrix-recovery setting. As shown in Figure 9, our RMF-based method reduces NMAE by up to 99.6% and is 4.09× faster than SVD. It also converges more smoothly, quickly, and stably due to efficient parallel SGD-based RMF optimization. Moreover, as matrix size increases, our accuracy improves while SVD remains nearly unchanged, reflecting the suitability of RMF for submatrix recovery and the benefit of gradient clipping.

200

100

40 60 80 Matrix Size (K)

Figure 10. Evaluation on large scale matrices.

2 GPU

300

20

100

1 GPU

400

100

500

0

500

40 60 80 Matrix Size (K)

Time

50 100 Iterator

Figure 7. Ablation study of TIDE-MC configurations: Base (no optimizations), Opt1 (conflict-free synchronization RMF), Opt2 (chunk-based pipelining), and Opt3 (gradient clipping). 300

10−1

Figure 9. Time and accuracy comparison: RMF vs. SVD. 0

Matrix Size

4 0 20

1

1K×3K 3K×3K 3K×10K 30×10K

6 2

0.5

1K×1K

100

8

2

SVD

RMF (Ours)

10

NMAE

30

With Opt3 Without Opt3

NMAE

Base Base + Opt2

Multi-GPU Scalability.

TIDE-MC scales across multiple GPUs. In Stage I, the template submatrix 𝑅 ′ is recovered by SGD-based RMF on a single GPU. This phase is lightweight and memory-bound, so multi-GPU parallelism brings limited benefit. In Stage II, TIDE-MC exploits multi-GPU parallelism by distributing column chunks across GPUs in a chunk-based pipeline, where each GPU processes a disjoint subset and results are merged to reconstruct the matrix. Asynchronous prefetching and write-back overlap I/O with computation to maximize utilization. We evaluate TIDE-MC with 1, 2, and 4 GPUs (Figure 8) on multiple matrix sizes. Runtime decreases as GPUs increase. For the smallest problem, runtime drops from 36.18s on 1 GPU to 22.12s on 2 GPUs (38.9%), and to 18.90s on 4 GPUs (47.8%). For a 100K × 150K matrix, runtime falls from 507s to 340.80s (32.8%), and further to 220.84s (56.4%). These results show that TIDE-MC distributes work effectively across GPUs with modest overhead, and that its chunk-based design scales well to larger problems. 6.7

Evaluation on Large-Scale Matrix

Figure 10 shows the execution time and NMAE of TIDEMC on large matrices, where all baselines run out of memory. Runtime increases from 36.18s on 10K×100K to 507s on 100K×150K, while NMAE does not degrade monotonically: it stays around 2.4 for 10K×100K and 50K×100K, then rises to 3.2× on 100K×150K. This suggests that, under current hyperparameters, TIDE-MC performs best near the 100K×100K

Comparison of RMF with SVD In Stage I

To quantify the benefit of RMF in Stage I, we replace it with a well-optimized SVD-based implementation and compare 11

1.5

1.5

1 1 0

0.5 0 0.2 0.4 0.6 0.8 1 Subcolumn Ratio

0

0.6

1

0.4

0.5

0.2 10

20

30

Number of Chunks

0

·10−2

·10−2

0.2 0.4 0.6 0.8 1 Subcolumn Ratio

0

0

10

20

30

Figure 11. Sampling ratio and chunk number sensitivity.

0.5 1 2 3 4 5 6 7 8 910

regularization (𝜆 × 10−2 )

·10−3

2 1 0

·10−3 3 2 1 0

1 2 3 4 5 6 7 8 910

learning rate (𝜂 × 10−2 )

1 2 3 4 5 6 7 8

clipping grad(𝜆 × 10−1 )

Figure 12. NMAE sensitivity.

scale. 6.8

1

0

Number of Chunks

3

NMAE

2

NMAE

1 0.8

NMAE

2

2 NMAE

Time (s)

3

encountering severe missingness. To validate robustness, we performed 1,000 random sampling trials on a 1K × 1K dataset. The performance variance across runs was less than 0.02%, confirming the stability of our approach.

Parameter Sensitivity

Effect of Subcolumn Sampling Ratio. The subcolumn sampling ratio determines the proportion of columns selected during Stage I. As shown in Figure 11, increasing the ratio from 0.01 to 0.1 leads to a sharp drop in runtime, i.e., from 2.49s to 0.847s, because the submatrix captures enough structure to enable faster convergence. At very low ratios (e.g., 0.01), the submatrix lacks sufficient information, resulting in more iterations and higher computational costs. From 0.1 to 0.5, runtime continues to decrease slowly, reaching a minimum of 0.45s at 0.5. This range balances representativeness and overhead. However, beyond 0.5, runtime increases again, up to 1.1s at a ratio of 1.0 due to the growing cost of processing a large submatrix. While a larger submatrix reduces Stage II workload, the overhead of handling more data offsets these gains, ultimately reducing overall efficiency. Thus, for small datasets, a ratio of 0.5 typically yields the best performance. For large-scale datasets, however, higher ratios incur excessive memory usage; therefore, we adopt a default ratio of 0.1 to balance efficiency and memory consumption. Parameter in Stage I. The number of chunks determines how the matrix is partitioned in Stage II for full-matrix reconstruction using chunk-based pipelined ALS. While it does not affect accuracy (per TSID theory), it strongly impacts runtime by influencing memory usage and parallelism. As shown in Figure 12, increasing the chunk count from 5 to 10 reduces runtime from 2.02s to 0.66s. With only 5 chunks, each chunk is large, limiting parallelism and increasing memory pressure, which slows computation. At 10 chunks, the matrix is split into smaller parts, enabling better parallel processing, a reduced memory footprint per chunk, and effective I/O-computation overlap. However, further increasing the chunk count leads to an increase in runtime due to overheads such as frequent CPU-GPU data transfers and more synchronization across tasks. This diminishes the benefits of parallelism and leads to performance degradation. For the 1K × 10K dataset, 10 chunks offer the best trade-off. While the optimal count may vary with dataset size, 10 remains a robust default for balancing efficiency and scalability. Effect of Chunk Count. In the worst case, submatrix sampling may degrade performance under highly nonuniform missing patterns. However, in real-world large-scale datasets, the likelihood of this occurring is very low. Our random sampling selects 10% of columns, which reduces the chance of

6.9

Limitation and Impact

Our design assumes that the sampled template submatrix 𝑅 ∗ fits in GPU memory. For extremely large problems, even the template may exceed device capacity, and we do not yet provide an out-of-core implementation for Stage I. Extending the template construction to out-of-core or multi-GPU settings is an important direction for future work. On the positive side, by reducing matrix completion to a template-resident, streamable workflow, TIDE-MC enables billion-scale completion on a single node and improves hardware efficiency on modern accelerators. This design can benefit a range of data-driven applications, such as recommendation systems and scientific data recovery, where large sparse matrices must be processed under tight memory budgets.

7

Related Work

Matrix Completion Algorithms. Early theoretical work established near-exact recovery guarantees using convex nuclear-norm relaxations [7, 9, 31, 52] and nonconvex lowrank factorizations of the form 𝑅 ≈ 𝑃𝑄 ⊤ [27, 29, 48]. Subsequent Bayesian formulations introduced sparsity-promoting priors over factors to improve robustness and uncertainty estimation [3]. More recent work has continued to strengthen the theory of matrix completion, including nearly linear-time robust alternating minimization [19] and distribution-free uncertainty quantification via conformalized matrix completion [21]. Other Matrix Completion Systems. To handle web-scale recommendation workloads, several systems parallelize matrix factorization across clusters or heterogeneous devices. A block-partitioned distributed SGD algorithm for matrices with billions of entries [17], and distributed variants such as DALS, ASGD, and DSGD++ that exploit both data and model parallelism [54]. More recent systems optimizie execution on heterogeneous CPU–GPU platforms [66]. These systems substantially accelerate matrix factorization, but they still treat matrix completion as a whole-matrix workload whose data or factors are partitioned across workers or accelerators. Positioning of TIDE-MC. TIDE-MC shares the same lowrank factorization objective as above work but targets a different point in the design space. Rather than scaling out 12

by sharding factors across a cluster, TIDE-MC scales up by using TSID as a resource-reduction primitive: a small template submatrix is kept in GPU memory while the remaining columns are streamed and reconstructed under a fixed memory budget. This submatrix-centric design restructures matrix completion into a bounded-memory execution workflow, enabling single-node, billion-scale completion that can be reused across applications.

8

tions on Circuits and Systems for Video Technology 25, 8 (2015), 1261– 1270. https://doi.org/10.1109/TCSVT.2014.2372351 [11] Eric C Chi, Hua Zhou, Guo-Kang Chen, Diego O Del Vecchyo, and Kenneth Lange. 2013. Genotype imputation via matrix completion. Genome Research 23, 3 (2013), 509–518. https://doi.org/10.1101/gr. 145821.112 [12] Davide Chicco, Matthijs J Warrens, and Giuseppe Jurman. 2021. The coefficient of determination R-squared is more informative than SMAPE, MAE, MAPE, MSE and RMSE in regression analysis evaluation. Peerj computer science 7 (2021), e623. https://doi.org/10.7717/peerj-cs.623 [13] Stephen Chou, Fredrik Kjolstad, and Saman Amarasinghe. 2018. Format abstraction for sparse tensor algebra compilers. Proc. ACM Program. Lang. 2, OOPSLA, Article 123 (Oct. 2018), 30 pages. https: //doi.org/10.1145/3276493 [14] Evgenia Chunikhina, Raviv Raich, and Thinh Nguyen. 2014. Performance analysis for matrix completion via iterative hard-thresholded SVD. In 2014 IEEE Workshop on Statistical Signal Processing (SSP). 392– 395. https://doi.org/10.1109/SSP.2014.6884658 [15] Kenneth L. Clarkson. 2010. Coresets, sparse greedy approximation, and the Frank-Wolfe algorithm. ACM Trans. Algorithms 6, 4, Article 63 (Sept. 2010), 30 pages. https://doi.org/10.1145/1824777.1824783 [16] Chenjiao Feng, Jiye Liang, Peng Song, and Zhiqiang Wang. 2020. A fusion collaborative filtering method for sparse data in recommender systems. Inf. Sci. 521, C (June 2020), 365–379. https://doi.org/10.1016/ j.ins.2020.02.052 [17] Rainer Gemulla, Erik Nijkamp, Peter J. Haas, and Yannis Sismanis. 2011. Large-scale matrix factorization with distributed stochastic gradient descent. In Proceedings of the 17th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. ACM, 69–77. https://doi.org/10.1145/2020408.2020426 [18] Gene H. Golub and Charles F. van Loan. 1980. An Analysis of the Total Least Squares Problem. SIAM J. Numer. Anal. 17, 6 (1980), 883–893. https://doi.org/10.1137/0717073 arXiv:https://doi.org/10.1137/0717073 [19] Yuzhou Gu, Zhao Song, Junze Yin, and Licheng Zhang. 2023. Low Rank Matrix Completion via Robust Alternating Minimization in Nearly Linear Time. ArXiv abs/2302.11068 (2023). https://api.semanticscholar. org/CorpusID:257079026 [20] Xin Guan, Chang-Tsun Li, and Yu Guan. 2017. Matrix Factorization With Rating Completion: An Enhanced SVD Model for Collaborative Filtering Recommender Systems. IEEE Access 5 (2017), 27668–27678. https://doi.org/10.1109/ACCESS.2017.2772226 [21] Yu Gui, Rina Barber, and Cong Ma. 2023. Conformalized matrix completion. In Thirty-seventh Conference on Neural Information Processing Systems. https://openreview.net/forum?id=6f320HfMeS [22] F. Maxwell Harper and Joseph A. Konstan. 2015. The MovieLens Datasets: History and Context. ACM Trans. Interact. Intell. Syst. 5, 4, Article 19 (Dec. 2015), 19 pages. https://doi.org/10.1145/2827872 [23] Trevor Hastie, Rahul Mazumder, Jason D Lee, and Reza Zadeh. 2015. Matrix Completion and Low-Rank SVD via Fast Alternating Least Squares. Journal of machine learning research : JMLR 16 (2015), 3367—3402. https://europepmc.org/articles/PMC6530939 [24] Yao Hu, Debing Zhang, Jieping Ye, Xuelong Li, and Xiaofei He. 2013. Fast and Accurate Matrix Completion via Truncated Nuclear Norm Regularization. IEEE Transactions on Pattern Analysis and Machine Intelligence 35, 9 (2013), 2117–2130. https://doi.org/10.1109/TPAMI. 2012.271 [25] Arnav Kapur, Kshitij Marwah, and Gil Alterovitz. 2016. Gene expression prediction using low-rank matrix completion. BMC Bioinformatics 17, 1 (2016), 243. https://doi.org/10.1186/s12859-016-1106-6 [26] Petr Kasalicky, Antoine Ledent, and Rodrigo Alves. 2023. Uncertaintyadjusted Inductive Matrix Completion with Graph Neural Networks. In Proceedings of the 17th ACM Conference on Recommender Systems (Singapore, Singapore) (RecSys ’23). Association for Computing Machinery, New York, NY, USA, 1169–1174. https://doi.org/10.1145/

Conclusion

We presented TIDE-MC, a submatrix-centric matrix completion system that uses TSID as a resource-reduction primitive to decouple logical matrix size from the GPU working set. This enables bounded-memory completion at billion-scale and yields up to 11,647× speedup with lower peak memory and higher accuracy. We hope that viewing TSID as a reusable resource-reduction primitive will inspire future systems for large-scale numerical workloads.

References [1] 2011. HetRec ’11: Proceedings of the 2nd International Workshop on Information Heterogeneity and Fusion in Recommender Systems (Chicago, Illinois). Association for Computing Machinery, New York, NY, USA. [2] Hervé Abdi. 2006. Singular Value Decomposition ( SVD ) and Generalized Singular Value Decomposition ( GSVD ). https://api. semanticscholar.org/CorpusID:8786634 [3] S. Derin Babacan, Martin Luessi, Rafael Molina, and Aggelos K. Katsaggelos. 2012. Sparse Bayesian Methods for Low-Rank Matrix Estimation. IEEE Transactions on Signal Processing 60, 8 (2012), 3964–3977. https://doi.org/10.1109/TSP.2012.2197748 [4] Dheeraj Baby and Soumyabrata Pal. 2024. Online Matrix Completion: A Collaborative Approach with Hott Items. In Proceedings of the 41st International Conference on Machine Learning (Proceedings of Machine Learning Research, Vol. 235), Ruslan Salakhutdinov, Zico Kolter, Katherine Heller, Adrian Weller, Nuria Oliver, Jonathan Scarlett, and Felix Berkenkamp (Eds.). PMLR, 2248–2276. https://proceedings.mlr. press/v235/baby24a.html [5] James Bennett and Stan Lanning. 2007. The netflix prize. (2007). https://www.cs.uic.edu/~liub/KDD-cup-2007/proceedings/ The-Netflix-Prize-Bennett.pdf [6] Christian Borgs, Jennifer Chayes, Christina E. Lee, and Devavrat Shah. 2017. Thy friend is my friend: iterative collaborative filtering for sparse matrix estimation. In Proceedings of the 31st International Conference on Neural Information Processing Systems (Long Beach, California, USA) (NIPS’17). Curran Associates Inc., Red Hook, NY, USA, 4718–4729. https://api.semanticscholar.org/CorpusID:19440176 [7] Jian-Feng Cai, Emmanuel J. Candès, and Zuowei Shen. 2010. A Singular Value Thresholding Algorithm for Matrix Completion. SIAM Journal on Optimization 20, 4 (2010), 1956–1982. https://doi.org/10.1137/ 080738970 arXiv:https://doi.org/10.1137/080738970 [8] Tianxi Cai, T. Tony Cai, and Anru Zhang. 2016. Structured Matrix Completion with Applications to Genomic Data Integration. J. Amer. Statist. Assoc. 111, 514 (2016), 621–633. https://doi.org/10.1080/01621459.2015. 1021005 arXiv:https://doi.org/10.1080/01621459.2015.1021005 PMID: 28042188. [9] Emmanuel J. Candès and Benjamin Recht. 2009. Exact matrix completion via convex optimization. Foundations of Computational Mathematics 9, 6 (2009), 717–772. https://doi.org/10.1007/s10208-009-9045-5 [10] Feilong Cao, Miaomiao Cai, and Yuanpeng Tan. 2015. Image Interpolation via Low-Rank Matrix Completion and Recovery. IEEE Transac13

age and Analysis (St. Louis, Missouri) (SC ’21). Association for Computing Machinery, New York, NY, USA, Article 58, 15 pages. https: //doi.org/10.1145/3458817.3476209 [41] Nagarajan Natarajan and Inderjit S. Dhillon. 2014. Inductive matrix completion for predicting gene–disease associations. Bioinformatics 30, 12 (06 2014), i60–i68. https://doi.org/10.1093/bioinformatics/ btu269 arXiv:https://academic.oup.com/bioinformatics/articlepdf/30/12/i60/48927398/bioinformatics_30_12_i60.pdf [42] Yuri Nesterov, Henry Wolkowicz, and Yinyu Ye. 2000. Semidefinite Programming Relaxations of Nonconvex Quadratic Optimization. Springer US, Boston, MA, 361–419. https://doi.org/10.1007/978-1-4615-43817_13 [43] Netflix. 2006. Netflix Prize Dataset. https://www.netflixprize.com/. https://www.netflixprize.com/ Accessed: 2025-04-14. [44] Uroš Ocepek, Jože Rugelj, and Zoran Bosnić. 2015. Improving matrix factorization recommendations for examples in cold start. Expert Systems with Applications 42, 19 (2015), 6784–6794. https://doi.org/10. 1016/j.eswa.2015.04.071 [45] Pablo A. Parrilo. 2003. Semidefinite programming relaxations for semialgebraic problems. Mathematical Programming 96, 2 (May 2003), 293–320. https://doi.org/10.1007/s10107-003-0387-5 [46] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala. 2019. PyTorch: An Imperative Style, High-Performance Deep Learning Library. In Advances in Neural Information Processing Systems, H. Wallach, H. Larochelle, A. Beygelzimer, F. d'AlchéBuc, E. Fox, and R. Garnett (Eds.), Vol. 32. Curran Associates, Inc. https://proceedings.neurips.cc/paper_files/paper/2019/file/ bdbca288fee7f92f2bfa9f7012727740-Paper.pdf [47] Gábor Pataki. 2000. The Geometry of Semidefinite Programming. Springer US, Boston, MA, 29–65. https://doi.org/10.1007/978-1-46154381-7_3 [48] Adityanarayanan Radhakrishnan, George Stefanakis, Mikhail Belkin, and Caroline Uhler. 2022. Simple, fast, and flexible framework for matrix completion with infinite width neural networks. Proceedings of the National Academy of Sciences 119, 16 (2022), e2115064119. https: //doi.org/10.1073/pnas.2115064119 [49] Samyam Rajbhandari, Jeff Rasley, Olatunji Ruwase, and Yuxiong He. 2020. ZeRO: Memory optimizations Toward Training Trillion Parameter Models. In SC20: International Conference for High Performance Computing, Networking, Storage and Analysis. 1–16. https://doi.org/ 10.1109/SC41405.2020.00024 [50] Soheil Sadeghi Eshkevari and Shamim N. Pakzad. 2020. Signal Reconstruction from Mobile Sensors Network Using Matrix Completion Approach. In Topics in Modal Analysis & Testing, Volume 8, Michael L. Mains and Brandon J. Dilworth (Eds.). Springer International Publishing, Cham, 61–75. https://doi.org/10.1007/978-3-030-12684-1_7 [51] Angela Serra, Pietro Coretto, Michele Fratello, and Roberto Tagliaferri. 2017. Robust and sparse correlation matrix estimation for the analysis of high-dimensional genomics data. Bioinformatics 34, 4 (10 2017), 625–634. https://doi.org/10.1093/bioinformatics/ btx642 arXiv:https://academic.oup.com/bioinformatics/articlepdf/34/4/625/48913355/bioinformatics_34_4_625.pdf [52] Peter J. Shin, Peder E. Z. Larson, Michael A. Ohliger, Michael Elad, John M. Pauly, Daniel B. Vigneron, and Michael Lustig. 2014. Calibrationless parallel imaging reconstruction based on structured low-rank matrix completion. Magnetic Resonance in Medicine 72, 4 (2014), 959– 970. https://doi.org/10.1002/mrm.24997 [53] Min Tao and Xiaoming Yuan. 2011. Recovering Low-Rank and Sparse Components of Matrices from Incomplete and Noisy Observations. SIAM Journal on Optimization 21, 1 (2011), 57–81. https://doi.org/10.

3604915.3610654 [27] Raghunandan H. Keshavan, Andrea Montanari, and Sewoong Oh. 2010. Matrix completion from a few entries. IEEE Transactions on Information Theory 56, 6 (2010), 2980–2998. https://doi.org/10.1109/ TIT.2010.2046205 [28] Diederik P. Kingma and Jimmy Ba. 2017. Adam: A Method for Stochastic Optimization. arXiv:1412.6980 [cs.LG] https://arxiv.org/abs/ 1412.6980 [29] Arezou Koohi and Houman Homayoun. 2019. Parallel Multi-View Graph Matrix Completion for Large Input Matrix. In 2019 IEEE 9th Annual Computing and Communication Workshop and Conference (CCWC). IEEE, Las Vegas, NV, USA, 337–341. https://dblp.org/rec/ conf/ccwc/KoohiH19 [30] Yehuda Koren, Robert Bell, and Chris Volinsky. 2009. Matrix factorization techniques for recommender systems. Computer 42, 8 (2009), 30–37. https://doi.org/10.1109/MC.2009.263 [31] Antonina Krajewska and Ewa Niewiadomska-Szynkiewicz. 2024. Efficient Data Completion and Augmentation. In 2024 IEEE 11th International Conference on Data Science and Advanced Analytics (DSAA). 1–10. https://doi.org/10.1109/DSAA61799.2024.10722793 [32] Antonina Krajewska and Ewa Niewiadomska-Szynkiewicz. 2026. Randomized Approach to Matrix Completion: Applications in Recommendation Systems and Image Inpainting. Machine Learning 115, 3 (2026), 1–32. https://doi.org/10.1007/s10994-026-06995-2 [33] Jiaxuan Li, Xiaoyan Zhu, Hongrui Wang, Yu Zhang, Xin Lai, and Jiayin Wang. 2025. Multi-Label Ranking Loss Minimization for Matrix Completion. Proceedings of the AAAI Conference on Artificial Intelligence 39, 17 (Apr. 2025), 18333–18340. https://doi.org/10.1609/aaai.v39i17.34017 [34] Mu Li, David G. Andersen, Jun Woo Park, Alexander J. Smola, Amr Ahmed, Vanja Josifovski, James Long, Eugene J. Shekita, and Bor-Yiing Su. 2014. Scaling Distributed Machine Learning with the Parameter Server. In Proceedings of the 11th USENIX Symposium on Operating Systems Design and Implementation (OSDI). USENIX. https://doi.org/ 10.1145/2640087.2644155 [35] Gangmuk Lim, Jeongseob Ahn, Wencong Xiao, Youngjin Kwon, and Myeongjae Jeon. 2021. Zico: Efficient GPU Memory Sharing for Concurrent DNN Training. In 2021 USENIX Annual Technical Conference (USENIX ATC 21). USENIX Association, 161–175. https://www.usenix. org/conference/atc21/presentation/lim [36] Canyi Lu, Changbo Zhu, Chunyan Xu, Shuicheng Yan, and Zhouchen Lin. 2015. Generalized Singular Value Thresholding. Proceedings of the AAAI Conference on Artificial Intelligence 29, 1 (Feb. 2015). https: //doi.org/10.1609/aaai.v29i1.9464 [37] Aanchal Mongia, Debarka Sengupta, and Angshul Majumdar. 2019. McImpute: Matrix Completion Based Imputation for Single Cell RNAseq Data. Frontiers in Genetics Volume 10 - 2019 (2019). https://doi. org/10.3389/fgene.2019.00009 [38] Arvind Narayanan and Vitaly Shmatikov. 2008. Robust Deanonymization of Large Sparse Datasets. In 2008 IEEE Symposium on Security and Privacy (sp 2008). 111–125. https://doi.org/10.1109/SP. 2008.33 [39] Deepak Narayanan, Aaron Harlap, Amar Phanishayee, Vivek Seshadri, Nikhil R. Devanur, Gregory R. Ganger, Phillip B. Gibbons, and Matei Zaharia. 2019. PipeDream: generalized pipeline parallelism for DNN training. In Proceedings of the 27th ACM Symposium on Operating Systems Principles (Huntsville, Ontario, Canada) (SOSP ’19). Association for Computing Machinery, New York, NY, USA, 1–15. https: //doi.org/10.1145/3341301.3359646 [40] Deepak Narayanan, Mohammad Shoeybi, Jared Casper, Patrick LeGresley, Mostofa Patwary, Vijay Korthikanti, Dmitri Vainbrand, Prethvi Kashinkunti, Julie Bernauer, Bryan Catanzaro, Amar Phanishayee, and Matei Zaharia. 2021. Efficient large-scale language model training on GPU clusters using megatron-LM. In Proceedings of the International Conference for High Performance Computing, Networking, Stor14

1137/100781894 arXiv:https://doi.org/10.1137/100781894 [54] Christina Teflioudi, Evangelos Makari, and Rainer Gemulla. 2012. Distributed matrix completion. In 2012 IEEE 12th International Conference on Data Mining. IEEE, 655–664. https://doi.org/10.1109/ICDM.2012. 120 [55] Adrien Todeschini, François Caron, and Marie Chavent. 2013. Probabilistic Low-Rank Matrix Completion with Adaptive Spectral Regularization Algorithms. In Advances in Neural Information Processing Systems. 845–853. https://proceedings.neurips.cc/paper_files/paper/ 2013/file/5ef0b4eba35ab2d6180b0bca7e46b6f9-Paper.pdf [56] Firas Trabelsi, David Vilar, Mara Finkelstein, and Markus Freitag. 2024. Efficient Minimum Bayes Risk Decoding using Low-Rank Matrix Completion Algorithms. In Advances in Neural Information Processing Systems, A. Globerson, L. Mackey, D. Belgrave, A. Fan, U. Paquet, J. Tomczak, and C. Zhang (Eds.), Vol. 37. Curran Associates, Inc., 54714– 54733. https://doi.org/10.52202/079017-1734 [57] Hoi-To Wai, Jean Lafond, Anna Scaglione, and Eric Moulines. 2017. Decentralized Frank–Wolfe Algorithm for Convex and Nonconvex Problems. IEEE Trans. Automat. Control 62, 11 (2017), 5522–5537. https://doi.org/10.1109/TAC.2017.2685559 [58] Jialei Wang, Peder A. Olsen, Andrew R. Conn, and Aurelie C. Lozano. 2016. Removing Clouds and Recovering Ground Observations in Satellite Image Sequences via Temporally Contiguous Robust Matrix Completion. In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR). https://doi.org/10.1109/CVPR.2016. 301 [59] Zaiwen Wen, Wotao Yin, and Yin Zhang. 2012. Solving a Low-Rank Factorization Model for Matrix Completion by a Nonlinear Successive Over-Relaxation Algorithm. Mathematical Programming Computation 4, 4 (2012), 333–361. https://doi.org/10.1007/s12532-012-0044-1 [60] Nicholas Wilt. 2013. The cuda handbook: A comprehensive guide to gpu programming. Pearson Education. https://www.cudahandbook.com/ [61] Wencong Xiao, Romil Bhardwaj, Ramachandran Ramjee, Muthian Sivathanu, Nipun Kwatra, Zhenhua Han, Pratyush Patel, Xuan Peng, Hanyu Zhao, Quanlu Zhang, Fan Yang, and Lidong Zhou. 2018. Gandiva: Introspective Cluster Scheduling for Deep Learning. In 13th USENIX Symposium on Operating Systems Design and Implementation (OSDI 18). USENIX Association, Carlsbad, CA, 595–610. https: //www.usenix.org/conference/osdi18/presentation/xiao [62] Kun Xie, Xiaocan Li, Xin Wang, Gaogang Xie, Jigang Wen, and Dafang Zhang. 2019. Active Sparse Mobile Crowd Sensing Based on Matrix Completion. In Proceedings of the 2019 International Conference on Management of Data (Amsterdam, Netherlands) (SIGMOD ’19). Association for Computing Machinery, New York, NY, USA, 195–210. https://doi.org/10.1145/3299869.3319856 [63] Zikai Xiong, Yimin Wei, Renjie Xu, and Yanwei Xu. 2022. Low-rank traffic matrix completion with marginal information. J. Comput. Appl. Math. 410 (2022), 114219. https://doi.org/10.1016/j.cam.2022.114219 [64] Yahoo. 2023. Yahoo Research Outreach Dataset. https://www.yahoo. com/research/outreach-dataset. https://www.yahoo.com/research/ outreach-dataset Accessed: 2025-04-14. [65] Haishan Ye, Ziang Zhou, Luo Luo, and Tong Zhang. 2020. Decentralized Accelerated Proximal Gradient Descent. In Advances in Neural Information Processing Systems, H. Larochelle, M. Ranzato, R. Hadsell, M.F. Balcan, and H. Lin (Eds.), Vol. 33. Curran Associates, Inc., 18308–18317. https://proceedings.neurips.cc/paper_files/paper/2020/ file/d4b5b5c16df28e61124e13181db7774c-Paper.pdf [66] Yuanhang Yu, Dong Wen, Ying Zhang, Xiaoyang Wang, Wenjie Zhang, and Xuemin Lin. 2021. Efficient Matrix Factorization on Heterogeneous CPU-GPU Systems. In 2021 IEEE 37th International Conference on Data Engineering (ICDE). 1871–1876. https://doi.org/10.1109/ ICDE51399.2021.00169 [67] Lianmin Zheng, Zhuohan Li, Hao Zhang, Yonghao Zhuang, Zhifeng Chen, Yanping Huang, Yida Wang, Yuanzhong Xu, Danyang Zhuo,

Eric P. Xing, Joseph E. Gonzalez, and Ion Stoica. 2022. Alpa: Automating Inter- and Intra-Operator Parallelism for Distributed Deep Learning. In 16th USENIX Symposium on Operating Systems Design and Implementation (OSDI 22). USENIX Association, Carlsbad, CA, 559– 578. https://www.usenix.org/conference/osdi22/presentation/zhenglianmin [68] Oliver Zier, Rahul Kannan, Aaron Smith, Mark Vogelsberger, and Erkin Verbeek. 2024. Adapting arepo-rt for exascale computing: GPU acceleration and efficient communication. Monthly Notices of the Royal Astronomical Society 533, 1 (2024), 268–286. https://doi.org/10.1093/ mnras/stae1837

15

Related documents

• Fast Differentiable SVD on GPU via Polar Decomposition#1116836
• Spectral Deflation for Factorization-Free Matrix Filtering in Muon and Semidefinite Programming#1035072
• PRISM: Processing-In-Memory Sparse MTTKRP for Tensor Decomposition Acceleration#238564
• REFINE: Super-efficient 3D Gaussian Splatting Pruning via Rendering-Free Primitive Importance#618102
• REFINE: Super-efficient 3D Gaussian Splatting Pruning via Rendering-Free Primitive Importance#617602