Stochastic Pattern-Completion Systems

• Stochastic pattern-completion systems are computational frameworks that reconstruct incomplete tensors using nonnegative constraints and stochastic gradient methods.
• They employ alternating optimization to decompose the tensor completion task into manageable nonnegative matrix subproblems, enhancing scalability.
• The approach leverages parallel accelerated stochastic gradients with adaptive step-sizes for rapid convergence and efficient large-scale performance.

Stochastic pattern-completion systems constitute a class of computational methods designed to reconstruct incomplete multidimensional arrays (“tensors”), particularly under nonnegativity constraints and in the presence of stochasticity introduced via sampling, noise, or algorithmic randomness. In the context of tensor completion, these systems are distinguished by their use of stochastic optimization schemes and pattern-exploiting decompositions to efficiently and scalably infer dense structure from sparse, partially observed, or corrupted data. The accelerated stochastic system proposed in "Accelerated Stochastic Gradient for Nonnegative Tensor Completion and Parallel Implementation" exemplifies this approach by integrating alternating optimization over CPD decomposed tensor factors with a parallel and locally adaptive stochastic gradient framework (Siaminou et al., 2021).

1. Mathematical Foundation of Nonnegative Tensor Completion

Let ${\mathbfcal X}^o \in \mathbb{R}_+^{I_1 \times \cdots \times I_N}$ be an unknown nonnegative tensor of canonical polyadic decomposition (CPD) rank $R$, represented as

${\mathbfcal X}^o = \llbracket U^{o(1)}, \dots, U^{o(N)} \rrbracket = \sum_{r=1}^R u_r^{o(1)} \circ \cdots \circ u_r^{o(N)},$

with each factor $U^{o(i)} \in \mathbb{R}_+^{I_i \times R}$. Only a noisy, partially observed version is accessible:

${\mathbfcal X} = {\mathbfcal X}^o + {\mathbfcal E}, \qquad {\mathbfcal M}(i_1, \dots, i_N) = \begin{cases} 1 & (i_1, \dots, i_N) \in \Omega \ 0 & \text{otherwise} \end{cases}$

where $\Omega$ indexes observed entries and ${\mathbfcal M}$ is the binary mask.

The nonnegative tensor completion (NTC) task is then formalized as:

$$\min_{\{ U^{(i)} \geq 0 \}_{i=1}^N} F(U^{(1)}, \dots, U^{(N)}) + \frac{\lambda}{2} \sum_{i=1}^N \|U^{(i)}\|_F^2,$$

where

$F(U^{(1)}, \dots, U^{(N)}) = \frac{1}{2} \| {\mathbfcal M} \circledast ( {\mathbfcal X} - \llbracket U^{(1)}, \dots, U^{(N)} \rrbracket ) \|_F^2.$

Mode-wise unfolding yields an equivalent nonnegative matrix completion (NMC) subproblem for each factor.

2. Alternating Optimization and Subproblem Decomposition

The tensor completion system operates via an outer alternating optimization (AO) loop cycling through all tensor modes $i = 1, \dots, N$. At each step $k$, the system computes the Khatri–Rao product $K_k^{(i)}$ of all modes except $i$, splitting the full problem into a sequence of nonnegative matrix completion subproblems:

$$U_{k+1}^{(i)} = \mathsf{S\_NMC}(X_{(i)}, M_{(i)}, K_k^{(i)}, U_k^{(i)}, \lambda).$$

Here, $X_{(i)}$ and $M_{(i)}$ are the unfolded tensor and mask in mode $i$, and $\mathsf{S\_NMC}$ denotes the stochastic matrix solver.

This approach enables scalability for high-dimensional and large-scale tensor data by leveraging separability across modes and subspace optimizations tailored to the NMC context.

3. Stochastic Accelerated Gradient for NMC

The nonnegative matrix completion subproblem addressed at each AO step is defined as:

$$\min_{A \geq 0} \; f_\Omega(A) = \frac{1}{2} \| M \circledast (X - A B^T) \|_F^2 + \frac{\lambda}{2} \|A\|_F^2,$$

where the optimization variables and data are $A \in \mathbb{R}_+^{P \times R}$, $X \in \mathbb{R}_+^{P \times Q}$, $B \in \mathbb{R}_+^{Q \times R}$.

At each inner iteration $l$, a subset $\widehat{\Omega}_l \subset \Omega$ is generated by subsampling a fraction $c$ of nonzeros within each row. The corresponding stochastic cost and gradient estimates are computed:

$$\nabla f_{\widehat\Omega_l}(Y_l) = -(\widehat M_l \circledast X - \widehat M_l \circledast (Y_l B^T)) B + \lambda Y_l,$$

executed row-wise.

Row-wise local smoothness constants $L_p$ are determined via the spectral properties of each block Hessian:

$$H_{l,p} = B^T \operatorname{diag}(\widehat M_l(p,:)) B + \lambda I_R,$$

defining individual row step-sizes and momentum parameters:

$$A_{l+1}(p,:) = (Y_l(p,:)-\tfrac{1}{L_p} \nabla_p f_{\widehat\Omega_l}(Y_l(p,:)))_+, \ \beta_{l,p} = \tfrac{\sqrt{L_p}-\sqrt{\lambda}}{\sqrt{L_p}+\sqrt{\lambda}},$$

with extrapolation:

$$Y_{l+1}(p,:) = A_{l+1}(p,:) + \beta_{l,p}[A_{l+1}(p,:) - A_l(p,:)].$$

Theoretically, the update exhibits an $O(1/l^2)$ convergence rate in the strongly convex, smooth regime, with rapid empirical convergence in practical settings (Siaminou et al., 2021).

4. Parallelization and Computational Scalability

The independence of row-wise stochastic accelerated gradient updates enables highly efficient parallelization. Within each inner iteration, each thread processes a distinct set of rows, with no overlap in writes to $A_{l+1}$ or $Y_{l+1}$. Synchronization is relegated to an implicit barrier at the end of each inner loop over $l$.

OpenMP is the parallelization framework utilized, supporting nearly linear speedup up to the physical core count. The table below summarizes scaling efficiency observed on benchmark datasets:

# Threads	Dataset	Speedup Factor
1	--	1×
16	Uber pickups (4th order, real)	10-15×
16	Uniform synthetic	>15×

On a $183 \times 24 \times 1140 \times 1717$ Uber pickups tensor ($3.3$M nonzeros), using 1 vs. 16 threads yields a $10$–$15\times$ reduction in per-epoch time, with stronger scaling on load-balanced synthetic data (Siaminou et al., 2021).

5. Empirical Performance and Benchmarks

Extensive empirical evaluation demonstrates the system’s effectiveness in several application regimes:

• Image completion: RGB image tensor $(1063 \times 1599 \times 3)$, $90\%$ missing, rank $R=50$, subsampling $c=0.02$, and $500$ epochs result in visually successful reconstructions even under minimal stochastic sampling.
• Large-scale recommendation: MovieLens 10M tensor $(71567 \times 65133 \times 730)$ and synthetic analogs show that, in high-noise environments, larger stochastic fractions $c$ yield faster reduction in relative reconstruction error (RRE) over epochs, while in low-noise regimes, smaller $c$ suffices:

$\mathrm{RRE}(k) = \frac{ \| {\mathbfcal M} \circledast ( \widehat{\mathbfcal X}_k - {\mathbfcal X}^o ) \|_F }{ \| {\mathbfcal M} \circledast {\mathbfcal X}^o \|_F }.$

• Scalability: Per-epoch computation time diminishes nearly linearly as threads increase from $1$ to the core limit, for all tested ranks $R$.

6. System Architecture Summary and Significance

The stochastic pattern-completion system integrates an AO outer loop centered on CPD-based tensor factorization with an inner, accelerated stochastic per-row solver that leverages local smoothness constants to adaptively set step-sizes. Full row-wise independence is harnessed via OpenMP, leading to parallel efficiency and scaling on extremely large, sparse tensors. This construction enables recovery of missing tensor structure in high-noise, high-dimensional settings with practical computational expenditure (Siaminou et al., 2021).

A plausible implication is that such systems will be critical for domain applications where data is both large and highly incomplete, such as recommendation, remote sensing, and scientific data fusion, given their robust empirical convergence for ill-posed and strongly regularized completion tasks.