1 Introduction↩︎

In many application domains, such as recommendation systems [1], weather forecasting [2], and network analysis [3], [4], the available data is often incomplete due to missing measurements, transmission failures, or other technical limitations [5]. These gaps in the data make it challenging to analyze and utilize the information effectively, requiring robust approaches to estimate the missing values.

Low-rank matrix completion has emerged as a widely used method to address this issue. It is based on the assumption that the observed data can be represented as a low-rank matrix, enabling the reconstruction of missing entries from a partially observed subset.

Throughout this paper, we denote by \(M \in \mathbb{R}^{m \times n}\) the complete target matrix that we aim to recover. In practice, we only observe a subset of its entries at indices \(\Omega \subseteq \{1,\ldots,m\} \times \{1,\ldots,n\}\).

1.1 Related work↩︎

Low-rank matrix completion is classically formulated as the problem of finding a low-rank matrix that satisfies constraints on the observed entries: \[\min_{X\in\mathbb{R}^{m\times n}} \quad \mathrm{rank}(X) \qquad \textrm{s.t.} \quad \mathcal{P}_\Omega(X) = \mathcal{P}_\Omega(M). \label{eq:RM}\tag{1}\] Here, \(\mathcal{P}_\Omega\) is the projection operator that restricts a matrix to its observed entries: \[[\mathcal{P}_\Omega(X)]_{ij} = \begin{cases} X_{ij} & \text{if } (i,j) \in \Omega \\ 0 & \text{otherwise}. \end{cases}\] Despite its intuitive formulation, this problem is NP-hard [6] due to the non-convexity of the rank function, making it computationally infeasible for large-scale matrices.

To circumvent the challenges posed by rank minimization, a convex relaxation replaces the rank function with the nuclear norm, defined as the sum of the singular values of the matrix. The relaxed problem is expressed as [6]: \[\min_{X\in\mathbb{R}^{m\times n}} \quad \|X\|_\ast \qquad \textrm{s.t.} \quad \mathcal{P}_\Omega(X) = \mathcal{P}_\Omega(M),\] where \(\|X\|_\ast\) serves as a convex surrogate for the rank. While effective for small matrices, nuclear norm minimization scales poorly with the matrix dimensions due to its reliance on repeated singular value decomposition (SVD) during optimization.

To improve computational efficiency, many methods adopt a factorized representation. These methods approximate the matrix as \(X = UW\), where \(U \in \mathbb{R}^{m \times r}\) and \(W \in \mathbb{R}^{r \times n}\). This reformulation enforces a rank constraint implicitly, transforming the problem into [7]: \[\min_{\substack{U\in\mathbb{R}^{m\times r}\\W\in\mathbb{R}^{r\times n}}} \quad \|\mathcal{P}_\Omega(UW) - \mathcal{P}_\Omega(M)\|_F^2 + \lambda(\|U\|_F^2 + \|W\|_F^2), \label{eq:MC2}\tag{2}\] where \(\lambda\) is a regularization parameter that controls the trade-off between fitting the observed entries and penalizing the complexity of \(U\) and \(W\). This factorized formulation significantly reduces the computational cost, making it suitable for large-scale problems [8].

However, the standard factorized approach treats all rows and columns independently, ignoring potential structural relationships in the data. In many real-world applications, additional side information about the relationships between rows and columns is available. This information is often naturally encoded as graphs that describe relationships, such as user similarities in recommendation systems or spatial proximity in weather data.

To leverage this structural information, graph-based regularization can be incorporated into matrix completion models to enforce smoothness along graph structures [9]–[12]. This means that connected nodes in the graph (e.g., nearby weather stations or similar users) should have similar latent representations. This leads to the following formulation [13]: \[\min_{\substack{U\in\mathbb{R}^{m\times r}\\W\in\mathbb{R}^{r\times n}}} \Big\{ \frac{1}{2}\left\|{\mathcal{P}_\Omega(UW)-\mathcal{P}_\Omega(M)}\right\|_\mathrm{2}^2 + \frac{\lambda_\mathrm{w}}{2}\left\|{U}\right\|_\mathrm{F}^2 + \frac{\lambda_\mathrm{h}}{2}\left\|{W}\right\|_\mathrm{F}^2 + \frac{\lambda_\mathrm{L}}{2}\left\|{U}\right\|_\mathrm{\mathcal{D},u}^2 + \frac{\lambda_\mathrm{L}}{2}\left\|{W}\right\|_\mathrm{\mathcal{D},w}^2 \Big\}.\] Here, the Dirichlet semi-norm \(\|\cdot\|_\mathcal{D}\) encodes graph smoothness constraints. For the row graph with Laplacian matrix \(L_\mathrm{u} \in \mathbb{R}^{m \times m}\), we have: \[\|U\|_{\mathcal{D},\mathrm{u}}^2 = \text{Tr}(U^T L_\mathrm{u} U) = \sum_{k=1}^r u_k^T L_\mathrm{u} u_k,\] where \(u_k\) is the \(k\)-th column of \(U\). Since \(L_\mathrm{u} = D_\mathrm{u} - A_\mathrm{u}\) with adjacency matrix \(A_\mathrm{u}\) and degree matrix \(D_\mathrm{u}\), this expands to: \[\|U\|_{\mathcal{D},\mathrm{u}}^2 = \sum_{k=1}^r \sum_{(i,j) \in E_\mathrm{u}} A_{\mathrm{u},ij} (U_{ik} - U_{jk})^2.\] This formulation penalizes large differences between connected nodes, encouraging similar rows (those connected in the graph) to have similar latent representations. The column regularization \(\|W\|_{\mathcal{D},\mathrm{w}}^2\) follows analogously with the column graph Laplacian \(L_\mathrm{w}\). Methods such as Graph-Regularized Alternating Least Squares (GRALS) [13] solve this formulation using the alternating minimization scheme.

Another promising line of work explores manifold-based optimization. The Riemannian Trust-Region Matrix Completion (RTRMC) [14], [15] framework operates on the Grassmann manifold \(\mathrm{Gr}(m, r)\), i.e., the set of \(r\)-dimensional subspaces in \(\mathbb{R}^m\). To handle the observed entries constraint, RTRMC employs a weighted formulation using a confidence matrix \(C\) and the Hadamard product \(\odot\). This approach differs from the projection operator \(\mathcal{P}_\Omega\) used in previous formulations. Specifically, let us define the \(\Omega\)-norm of a matrix \(M\) as: \[\|M\|_\Omega^2 := \sum_{(i,j) \in \Omega} M_{ij}^2,\] where \(\Omega\) denotes the set of observed entries. The confidence matrix \(C\) assigns weights to observed entries, with \(C_{ij} > 0\) for \((i,j) \in \Omega\) and \(C_{ij} = 0\) otherwise. RTRMC solves the following optimization problem: \[\min_{\mathcal{U}\in\mathrm{Gr}(m,r)}\Big(\min_{W\in\mathbb{R}^{r\times n}} \frac{1}{2}\|C\odot (UW-M)\|_\Omega^2 + \frac{\lambda^2}{2} \|\mathcal{P}_{\bar{\Omega}}(UW)\|_F^2\Big). \label{eqn:rtrmc}\tag{3}\]

1.2 Contribution↩︎

Building on RTRMC, we propose a graph-regularized extension that combines the advantages of Riemannian optimization with graph-based regularization. Our key contribution is extending the RTRMC framework by incorporating Dirichlet semi-norms that enforce smoothness along graph-encoded relationships. The proposed method solves the following optimization problem: \[\min_{\mathcal{U}\in\mathrm{Gr}(m,r)}\Big(\min_{W\in\mathbb{R}^{r\times n}} \frac{1}{2}\|C\odot (UW-M)\|_\Omega^2 + \frac{\lambda^2}{2} \|\mathcal{P}_{\bar{\Omega}}(UW)\|_F^2 + \frac{\lambda_\mathrm{u}}{2}\|U\|_{\mathcal{D},\mathrm{u}}^2 + \frac{\lambda_\mathrm{w}}{2}\|W\|_{\mathcal{D},\mathrm{w}}^2\Big), \label{eq:proposed}\tag{4}\] where the graph regularization terms \(\|U\|_{\mathcal{D},\mathrm{u}}^2\) and \(\|W\|_{\mathcal{D},\mathrm{w}}^2\) penalize deviations along graph structures associated with rows and columns respectively.

To demonstrate the effectiveness of our proposed method, we conduct experiments on meteorological datasets from Belgium and France, as well as the MovieLens collaborative filtering dataset. Our evaluation framework uses k-fold cross-validation to ensure robustness and fair comparisons across methods. The results highlight the superior performance of the graph-regularized RTRMC model, showing improved accuracy in recovering missing entries compared to baseline methods. These findings highlight the value of using graph-based structure when filling in missing data in matrices, especially when there are relationships between rows and columns.

A brief preview of this work, centered on the Belgian weather data application, was presented at the BNAIC/BeNeLearn 2024 conference [16]. The present paper extends the conference version with: (i) a complete theoretical and algorithmic development of the proposed model, (ii) experiments on larger-scale French meteorological data and on the MovieLens collaborative filtering dataset, and (iii) detailed case studies analyzing failure modes during extreme weather events.

The proposed method is presented in Section 2. Numerical experiments are reported in Section 3. Conclusions are drawn in Section 4.

2 Graph-Regularized Riemannian Trust-Region Matrix Completion↩︎

We present the Graph-Regularized Riemannian Trust-Region Matrix Completion
(GR-RTRMC) method, which enhances the classical RTRMC approach by incorporating graph-based regularization. This extension leverages structural relationships between rows and columns of the target matrix to improve recovery performance, particularly in applications where such relationships naturally arise.

Building upon the RTRMC framework introduced in 3 , we incorporate the graph regularization terms from 4 to formulate our optimization problem over matrix factors \(U \in \mathbb{R}_*^{m \times r}\) and \(W \in \mathbb{R}^{r \times n}\), where \(\mathbb{R}_*^{m \times r}\) denotes the set of full-rank \(m \times r\) matrices. We use the confidence matrix \(C\) (with \(C_{ij} > 0\) for \((i,j) \in \Omega\) and \(C_{ij} = 0\) otherwise), graph Laplacians \(L_\mathrm{u} \in \mathbb{R}^{m \times m}\) and \(L_\mathrm{w} \in \mathbb{R}^{n \times n}\) for rows and columns respectively, and \(\overline{\Omega} = \{1,\ldots,m\} \times \{1,\ldots,n\} \setminus \Omega\) to denote the unobserved entries:

\[\begin{align} \min_{U,W} h(U,W) &= \frac{1}{2}\|C \odot (UW - M)\|^2_\Omega + \frac{\lambda^2}{2}\|UW\|^2_{\overline{\Omega}} \nonumber \\ &\quad + \frac{\lambda_\mathrm{u}}{2}\text{Trace}(U^TL_\mathrm{u}U) + \frac{\lambda_\mathrm{w}}{2}\text{Trace}(WL_\mathrm{w}W^T). \label{eq:gr-objective} \end{align}\tag{5}\]

The term weighted by \(\lambda > 0\) regularizes unobserved entries, while \(\lambda_\mathrm{u}, \lambda_\mathrm{w} \geq 0\) control graph regularization strength. The graph regularization terms encourage similar rows and columns (as defined by the respective graph structures) to have similar representations in the factorization.

2.1 Optimization on the Grassmann Manifold↩︎

The key insight [14], [15] enabling efficient optimization of 5 lies in recognizing that the factorization \(UW\) exhibits a fundamental invariance property. This observation allows us to reformulate the high-dimensional optimization problem over matrix factors into a more structured optimization problem on a Riemannian manifold.

To understand this reformulation, we begin by noting that for any fixed left factor \(U \in \mathbb{R}_{*}^{m \times r}\), finding the optimal right factor \(W\) reduces to solving a regularized least-squares problem: \[W_{U} = \underset{W \in \mathbb{R}^{r \times n}}{\operatorname{argmin}} h(U,W). \label{eq:wu-definition}\tag{6}\]

By substituting this optimal \(W_{U}\) back into the original objective 5 , we obtain the reduced cost function: \[f(U) = h(U, W_{U}), \label{eq:reduced-cost}\tag{7}\] a technique known as variable projection [17]. The crucial observation is that \(W_U\) depends only on the column space of \(U\), not on the specific choice of basis for that subspace. To see this, consider any invertible transformation \(M \in \mathbb{R}_*^{r \times r}\) and define \(\tilde{U} = UM\). Then the optimal factor for \(\tilde{U}\) satisfies: \[W_{\tilde{U}} = W_{UM} = M^{-1}W_U.\]

This relationship ensures that the matrix product remains invariant under basis changes: \[\tilde{U}W_{\tilde{U}} = (UM)(M^{-1}W_U) = UW_U.\]

This invariance property is fundamental: it shows that the quality of any factorization \(UW\) depends only on the column space spanned by \(U\), not on the particular orthonormal basis chosen to represent it. This naturally motivates a formulation on the Grassmann manifold \(\text{Gr}(m,r)\), which is precisely the set of \(r\)-dimensional subspaces of \(\mathbb{R}^m\).

For computational efficiency, we work with orthonormal representations on the Stiefel manifold \(\text{St}(m,r) = \{U \in \mathbb{R}^{m \times r} : U^TU = I_r\}\). When \(U\) has orthonormal columns, we can exploit the identity \(\|UW_{U}\|^2_{\overline{\Omega}} = \|W_{U}\|^2_F - \|UW_{U}\|^2_\Omega\) (valid due to \(\|UW_{U}\|^2_F = \|W_{U}\|^2_F\) when \(U^TU = I_r\)) to rewrite our objective as:

\[\begin{align} f(U) = \frac{1}{2}\|C \odot (UW_{U} - M)\|^2_\Omega &+ \frac{\lambda^2}{2}(\|W_{U}\|^2_F - \|UW_{U}\|^2_\Omega) \nonumber \\ &+ \frac{\lambda_\mathrm{u}}{2}\text{Trace}(U^TL_\mathrm{u}U) + \frac{\lambda_\mathrm{w}}{2}\text{Trace}(W_{U}L_\mathrm{w}W_{U}^T). \label{eq:f-grassmann} \end{align}\tag{8}\]

Since \(f(U)\) depends only on the column space of \(U\) due to the invariance property demonstrated above, this function is well-defined on the Grassmann manifold \(\text{Gr}(m,r)\). The optimization problem now becomes, with an abuse of notation discussed in [15]: \[\min_{U \in \text{Gr}(m,r)} f(U), \label{eq:gr-problem}\tag{9}\] where we optimize over \(r\)-dimensional subspaces in \(\mathbb{R}^m\), understanding that in practice we represent these subspaces using orthonormal matrices \(U \in \text{St}(m,r)\).

In comparison with the original optimization problem 5 , problem 9 can be thought of as more intricate in the sense that the linear feasible set has been replaced by a nonlinear manifold. However, this is mitigated by the fact that formulation 9 enables the application of Riemannian optimization algorithms [18]–[20], which can efficiently navigate the geometry of the Grassmann manifold while respecting its intrinsic structure. Furthermore, 9 has two advantages over 5 : the dimension of the feasible set is reduced—drastically when \(m\ll n\)—and the Grassmann manifold is compact.

2.2 Gradient Computation↩︎

Since \(W_U\) minimizes \(h(U,\cdot)\), the optimality condition gives us \(\frac{\partial h}{\partial W}(U,W_U) = 0\). By the envelope theorem, the directional derivative simplifies to: \[Df(U)[H] = \frac{\partial h}{\partial U}(U,W_U)[H].\]

To facilitate differentiation, we rewrite \(h(U,W)\) using only Frobenius norms. We define the mask matrix: \[\Lambda_{ij} = \begin{cases} \lambda & \text{if } (i,j) \in \Omega \\ 0 & \text{otherwise}. \end{cases}\]

Using the identity \(\|UW\|^2_{\overline{\Omega}} = \|W\|^2_F - \|UW\|^2_\Omega\) (valid when \(U\) has orthonormal columns), we can rewrite: \[\begin{align} h(U,W) = & \frac{1}{2}\|C \odot (UW - M)\|^2_F + \frac{\lambda^2}{2}\|W\|^2_F \nonumber \\ & - \frac{1}{2}\|\Lambda \odot UW\|^2_F + \frac{\lambda_\mathrm{u}}{2}\text{Trace}(U^TL_\mathrm{u}U) + \frac{\lambda_\mathrm{w}}{2}\text{Trace}(WL_\mathrm{w}W^T). \end{align}\]

This form simplifies differentiation, as for differentiable mappings \(g\): \[D\left(X \mapsto \frac{1}{2}\|g(X)\|^2_F\right)(X)[H] = \langle Dg(X)[H], g(X) \rangle.\]

Defining \(\hat{C} = C^{(2)} - \boldsymbol{\Lambda}^{(2)}\) (where \((\cdot)^{(2)}\) denotes element-wise squaring), we obtain: \[Df(U)[H] = \langle H, [\hat{C} \odot (UW_U - M_\Omega) - \lambda^2 M_\Omega]W_U^T + \lambda_\mathrm{u} L_\mathrm{u} U \rangle.\]

Define the residue matrix: \[R_U = \hat{C} \odot (UW_U - M_\Omega) - \lambda^2 M_\Omega. \label{eq:residue-matrix}\tag{10}\]

This gives us the Euclidean gradient: \[\nabla f(U) = R_U W_U^T + \lambda_\mathrm{u} L_\mathrm{u} U.\]

In view of [15], and with the same abuse of notation and language, the Riemannian gradient on the Grassmannian is obtained by projecting onto the tangent space: \[\text{grad } f(U) = (I - UU^T)\nabla f(U) = (I - UU^T)(R_U W_U^T + \lambda_\mathrm{u} L_\mathrm{u} U).\]

2.3 Hessian Computation↩︎

We now differentiate the gradient expression to obtain the Hessian of \(f\). We define the vector field \(\bar{F}: \mathbb{R}_*^{m \times r} \rightarrow \mathbb{R}^{m \times r}\): \[\bar{F}(U) = (I - UU^T)(R_U W_U^T + \lambda_\mathrm{u} L_\mathrm{u} U).\]

Distributing the projector, we can rewrite this as: \[\bar{F}(U) = R_U W_U^T + \lambda_\mathrm{u} L_\mathrm{u} U - UU^T R_U W_U^T - \lambda_\mathrm{u} UU^T L_\mathrm{u} U.\]

The Riemannian Hessian of \(f\) on the Grassmannian is represented by: \[\text{Hess } f(U)[H] = (I - UU^T)D\bar{F}(U)[H],\] \(\forall H \in \mathbb{R}^{m\times r}\) with \(U^T H = 0\); see [15].

Computing the differential of \(\bar{F}\) term by term yields: \[\begin{align} D\bar{F}(U)[H] &= [\hat{C} \odot (HW_U + UW_{U,H})]W_U^T + R_U W_{U,H}^T + \lambda_\mathrm{u} L_\mathrm{u} H \nonumber \\ &\quad- (HU^T + UH^T) R_U W_U^T - UU^T [\hat{C} \odot (HW_U + UW_{U,H})]W_U^T\nonumber \\ &\quad- UU^T R_U W_{U,H}^T - \lambda_\mathrm{u} (HU^T + UH^T) L_\mathrm{u} U - \lambda_\mathrm{u} UU^T L_\mathrm{u} H. \end{align}\]

Applying the projector \((I - UU^T)\) and using \((I - UU^T)U = 0\) to eliminate terms, the final expression for the Hessian becomes: \[\begin{align} \text{Hess } f(U)[H] = &(I - UU^T)([\hat{C} \odot (HW_U + UW_{U,H})]W_U^T + R_U W_{U,H}^T + \lambda_\mathrm{u} L_\mathrm{u} H) \nonumber \\ &- HU^T(R_U W_U^T + \lambda_\mathrm{u} L_\mathrm{u} U). \end{align}\]

The computation requires \(W_{U,H}\), the differential of the mapping \(U \mapsto W_U\) along \(H\), which is provided in the following subsection.

2.4 Computing \(W_U\) and Its Derivative↩︎

We now derive a clean and explicit formulation for computing \(W_U\), the optimal right factor given a fixed left factor \(U\), as well as its directional derivative \(W_{U,H}\). These computations are central to evaluating the cost function \(f(U)\) and its Riemannian Hessian. Following the approach of [15], we use vectorization and Kronecker product properties to derive closed-form expressions.

For a fixed \(U \in \text{St}(m,r)\), computing \(W_U\) amounts to solving the least-squares problem 6 : \[h(U, W) = \frac{1}{2} \|C \odot (UW - M_\Omega)\|_F^2 + \frac{\lambda^2}{2} \|W\|^2_F - \frac{\lambda^2}{2} \|UW\|^2_\Omega + \frac{\lambda_\mathrm{w}}{2} \text{Tr}(W L_\mathrm{w} W^T).\]

To facilitate vectorization, we use the following notation. Let \(w = \text{vec}(W) \in \mathbb{R}^{rn}\) be the vectorization of \(W\).

Following [15], we define the matrices \(S \in \mathbb{R}^{k \times mn}\) and operators that will allow us to express the objective in standard quadratic form: \[S = I_\Omega \text{diag}(\text{vec}(C)),\] where \(I_\Omega \in \mathbb{R}^{k \times mn}\) is the sampling operator such that \(I_\Omega \text{vec}(M)\) extracts the entries of \(M\) corresponding to indices in \(\Omega\).

We can now express the objective function \(h(U,W)\) in terms of \(w\):

\[\begin{align} h(U,W) &= \frac{1}{2}\left\|{C\odot (UW-M_\Omega)}\right\|_\mathrm{\Omega}^2 + \frac{\lambda^2}{2}\left\|{W}\right\|_\mathrm{F}^2 - \frac{\lambda^2}{2}\left\|{UW}\right\|_\mathrm{\Omega}^2 + \frac{\lambda_u^2}{2}\left\|{U}\right\|_\mathrm{\mathcal{D},u}^2 \\ &\quad+ \frac{\lambda_w^2}{2}\left\|{W^\top\! }\right\|_\mathrm{\mathcal{D},w}^2\\ &= \frac{1}{2}\left\|{S \mathrm{vec}_\mathrm{}\left({UW}\right) - \mathrm{vec}_\mathrm{\Omega}\left({C\odot M_\Omega}\right)}\right\|_\mathrm{2}^2\\ &\quad + \frac{\lambda^2}{2}\left\|{\mathrm{vec}_\mathrm{}\left({W}\right)}\right\|_\mathrm{2}^2 - \frac{\lambda^2}{2}\left\|{\mathrm{vec}_\mathrm{\Omega}\left({UW}\right)}\right\|_\mathrm{2}^2 + \frac{\lambda_u^2}{2}\left\|{U}\right\|_\mathrm{\mathcal{D},u}^2 + \frac{\lambda_w^2}{2}\left\|{W^\top\! }\right\|_\mathrm{\mathcal{D},w}^2\\ &=\frac{1}{2}\left\|{S(I_n \otimes U)\mathrm{vec}_\mathrm{}\left({W}\right)-\mathrm{vec}_\mathrm{\Omega}\left({C\odot M_\Omega}\right)}\right\|_\mathrm{2}^2 + \frac{1}{2}\left\|{\lambda I_{rn}\mathrm{vec}_\mathrm{}\left({W}\right)}\right\|_\mathrm{2}^2 \\ &\quad- \frac{1}{2}\left\|{\lambda I_\Omega (I_n \otimes U)\mathrm{vec}_\mathrm{}\left({W}\right)}\right\|_\mathrm{2}^2 + \frac{\lambda_u^2}{2}\left\|{U}\right\|_\mathrm{\mathcal{D},u}^2 + \frac{\lambda_w^2}{2}\left\|{W^\top\! }\right\|_\mathrm{\mathcal{D},w}^2\\ &= \frac{1}{2}\left\|{\begin{bmatrix} S(I_n \otimes U) \\ \lambda I_{rn} \end{bmatrix} \mathrm{vec}_\mathrm{}\left({W}\right) - \begin{bmatrix} \mathrm{vec}_\mathrm{\Omega}\left({C\odot M_\Omega}\right)\\0_{rn} \end{bmatrix}}\right\|_\mathrm{2}^2 \\ &\quad- \frac{1}{2}\left\|{\lambda I_\Omega (I_n \otimes U)\mathrm{vec}_\mathrm{}\left({W}\right)}\right\|_\mathrm{2}^2+\frac{\lambda_w^2}{2}\mathrm{Tr}(WL_wW^\top\! ) +\frac{\lambda_u^2}{2}\left\|{U}\right\|_\mathrm{\mathcal{D},u}^2\\ &= \frac{1}{2}\left\|{\begin{bmatrix} S(I_n \otimes U) \\ \lambda I_{rn} \end{bmatrix} \mathrm{vec}_\mathrm{}\left({W}\right) - \begin{bmatrix} \mathrm{vec}_\mathrm{\Omega}\left({C\odot M_\Omega}\right)\\0_{rn} \end{bmatrix}}\right\|_\mathrm{2}^2 +\frac{\lambda_u^2}{2}\left\|{U}\right\|_\mathrm{\mathcal{D},u}^2 \\ &\quad- \frac{1}{2}\left\|{\lambda I_\Omega (I_n \otimes U)\mathrm{vec}_\mathrm{}\left({W}\right)}\right\|_\mathrm{2}^2+ \frac{\lambda_w^2}{2}\mathrm{vec}_\mathrm{}\left({W}\right)^\top\! (L_w\otimes I_r)\mathrm{vec}_\mathrm{}\left({W}\right)\\ &=\frac{1}{2}\left\|{A_1w - b_1}\right\|_\mathrm{2}^2-\frac{1}{2}\left\|{A_2w}\right\|_\mathrm{2}^2+\frac{\lambda_w^2}{2}w^\top\! (L_w\otimes I_r)w+\frac{\lambda_u^2}{2}\left\|{U}\right\|_\mathrm{\mathcal{D},u}^2\\ &=\frac{1}{2}w^\top\! (A_1^\top\! A_1 - A_2^\top\! A_2 + \lambda_w^2(L_w\otimes I_r))w - b_1^\top\! A_1 w + \frac{1}{2}\left\|{b_1}\right\|_\mathrm{2}^2 + \frac{\lambda_u^2}{2}\left\|{U}\right\|_\mathrm{\mathcal{D},u}^2\\ &=\frac{1}{2}w^\top\! A w - b^\top\! w + \text{const}, \end{align}\] where: \[\begin{align} A &= (I_n \otimes U^T) S^T S (I_n \otimes U) + \lambda^2 I_{rn} - \lambda^2 (I_n \otimes U^T) I_\Omega^T I_\Omega (I_n \otimes U) + \lambda_\mathrm{w}^2 (L_\mathrm{w} \otimes I_r). \\ b &= (I_n \otimes U^T) S^T \text{vec}_\Omega(C \odot M_\Omega). \end{align}\]

Define \(\hat{C} = C^{(2)} - \Lambda^{(2)}\) where \((\cdot)^{(2)}\) denotes element-wise squaring. Then: \[\begin{align} A &= (I_n \otimes U^T) \text{diag}(\text{vec}(\hat{C})) (I_n \otimes U) + \lambda^2 I_{rn} + \lambda_\mathrm{w}^2 (L_\mathrm{w} \otimes I_r). \\ b &= (I_n \otimes U^T) \text{vec}(C^{(2)} \odot M_\Omega). \end{align}\]

The solution is given by: \[\text{vec}(W_U) = A^{-1} b. \label{eq:wu-solution}\tag{11}\]

The key difference from the original RTRMC formulation is the additional term \(\lambda_\mathrm{w} (L_\mathrm{w} \otimes I_r)\) in the matrix \(A\), which introduces coupling between all columns of \(W\) via the graph Laplacian \(L_\mathrm{w}\). This represents a significant departure from the block-diagonal structure in the original RTRMC, where each column of \(W\) could be computed independently.

Let us compute the directional derivative \(W_{U,H}\) of the mapping \(U \mapsto W_U\) along a direction \(H\): \[\text{vec}(W_{U,H}) = D(U \mapsto \text{vec}(W_U))(U)[H].\]

Using the formula for the differential of the inverse of a matrix and the product rule: \[\begin{align} \text{vec}(W_{U,H}) &= -A^{-1} \cdot D(U \mapsto A)(U)[H] \cdot \text{vec}(W_U) + A^{-1} \cdot D(U \mapsto b)(U)[H]. \end{align}\]

Computing the differentials: \[D(U \mapsto b)(U)[H] = (I_n \otimes H^T) \text{vec}(\hat{C} \odot M_\Omega) = \text{vec}(H^T (\hat{C} \odot M_\Omega)).\] \[\begin{align} D(U \mapsto A)(U)[H] &= (I_n \otimes H^T) \text{diag}(\text{vec}(\hat{C})) (I_n \otimes U) \\&\quad+ (I_n \otimes U^T) \text{diag}(\text{vec}(\hat{C})) (I_n \otimes H). \nonumber \end{align}\]

Let \(B = \text{diag}(\text{vec}(\hat{C}))\). Then: \[\begin{align} \text{vec}(W_{U,H}) &= -A^{-1} [(I_n \otimes H^T) B (I_n \otimes U) + (I_n \otimes U^T) B (I_n \otimes H)] \text{vec}(W_U)\\ &\quad + A^{-1} \text{vec}(H^T (\hat{C} \odot M_\Omega)) \nonumber\\ &= -A^{-1} [(I_n \otimes H^T) \text{vec}(\hat{C} \odot UW_U) + (I_n \otimes U^T) \text{vec}(\hat{C} \odot HW_U))] \nonumber\\ &\quad + A^{-1} \text{vec}(H^T (\hat{C} \odot M_\Omega)) \nonumber\\ &= -A^{-1} [\text{vec}(H^T (\hat{C} \odot UW_U)) + \text{vec}(U^T (\hat{C} \odot HW_U))] \nonumber\\ &\quad + A^{-1} \text{vec}(H^T (\hat{C} \odot M_\Omega)). \nonumber \end{align}\]

Recalling the definition of the residue matrix \(R_U = \hat{C} \odot (UW_U - M_\Omega)\) from equation 10 , we obtain the compact expression: \[\text{vec}(W_{U,H}) = -A^{-1} \left[ \text{vec}(H^T R_U) + \text{vec}(U^T (\hat{C} \odot HW_U)) \right]. \label{eq:wu-h-final}\tag{12}\]

2.5 Initialization Strategy↩︎

The GR-RTRMC algorithm requires an initial point on the Grassmann manifold \(\text{Gr}(m,r)\), represented by an orthonormal matrix \(U_0 \in \text{St}(m,r)\). Unlike factorization-based approaches that require initialization of both matrices, our formulation only requires initializing the left matrix \(U\), since the right matrix \(W\) is computed optimally via the least-squares solution 11 .

We initialize \(U_0\) using the \(r\) dominant left singular vectors of the observed data matrix \(M_\Omega\): \[U_0 = \text{SVD}_r(M_\Omega)_{\text{left}},\] where \(\text{SVD}_r(\cdot)_{\text{left}}\) denotes the matrix of the \(r\) leading left singular vectors. Since \(M_\Omega\) is sparse with \(k = |\Omega|\) nonzero entries, this computation can be performed efficiently using the PROPACK sparse SVD algorithm [21], with computational cost \(O(kr + r^3)\).

3 Numerical Experiments↩︎

All experiments are performed on a desktop with 4.5 GHz AMD Ryzen 9 7950X with 64GB of RAM running Julia 1.10.4. Reported computational times exclude the initial compilation overhead by discarding the first execution.

We conduct experiments on two types of datasets: weather data (Belgium and France) and MovieLens.

For the meteorological datasets, we follow the experimental protocol (data splits, masking scenarios, and evaluation metrics) and station-graph construction described in [22]. The data matrix has rows corresponding to weather stations and columns to time steps, so that the row graph encodes spatial proximity between stations. The spatial graph connecting weather stations is built using inverse distance weighting with altitude corrections, among other approaches explored in that work; see also [9] for graph construction in matrix completion. We consider two missing data scenarios: Block, where contiguous time intervals are missing at selected stations, and Spread, where missing entries are distributed more uniformly across time and stations. For all methods, hyperparameters (including the rank \(r\) and regularization parameters) are selected via cross-validation on a validation set, and the reported results are computed on a separate test set.

3.1 Methods used for comparison↩︎

For the weather datasets, the following methods were employed for matrix completion:

1. IDW (Inverse Distance Weighting) [23], [24]
   A simple interpolation method that estimates missing values by calculating a weighted average of nearby observations. The weights decrease as the distance from the missing value increases.

2. PCA (Principal Component Analysis) [24], [25]
   A dimensionality reduction technique that imputes missing values by projecting the data onto a lower-dimensional space defined by its principal components.

3. MissForest [26]
   A machine learning-based imputation method that uses random forests. Missing values are predicted iteratively based on the relationships between observed variables.

4. LMaFiT (Low-rank Matrix Fitting) [27]
   A matrix factorization approach that models the observed data as a low-rank matrix. Missing values are estimated through iterative optimization.

5. SoftImpute [28]
   A matrix completion algorithm that employs soft-thresholded singular value decomposition (SVD) to approximate missing entries by minimizing the reconstruction error.

6. RTRMC (Riemannian Trust-Region Matrix Completion) [14]
   A matrix completion technique that operates on the Riemannian manifold of low-rank matrices, employing trust-region optimization to estimate missing data.

7. GRALS (Graph Regularized Alternating Least Squares) [13]
   An imputation method that incorporates graph-based regularization with alternating least squares to leverage structural information within the dataset.

3.2 Weather Data↩︎

Missing data in time series from weather stations can result from a variety of reasons like sensor malfunctions, power failures, communication disruptions, maintenance or removal of misleading values. Addressing these gaps is however essential for climate and weather analyses, which mostly rely on complete datasets.

Before reporting results, we make explicit how the weather imputation task instantiates the graph-regularized problem 5 . The data matrix \(M\in\mathbb{R}^{m\times n}\) collects the temperature measurements: its \(m\) rows index the weather stations and its \(n\) columns index the regularly sampled time steps of the considered time window (for instance \(n=1008\) for the Belgian dataset). The confidence matrix \(C\) is set to one on observed entries and to zero elsewhere, i.e., uniform confidence in the recorded values; the only exception is the storm case study (8), where \(C_{i,j}\) is deliberately lowered over an anomalous time window. Since the rows are the stations, the row Laplacian \(L_\mathrm{u}\in\mathbb{R}^{m\times m}\) encodes the spatial graph between stations, following [22]. This graph is a \(k\)-nearest-neighbour graph on the inter-station haversine distance, that is, the surface distance between two points given their latitude and longitude, with inverse-distance edge weights and an altitude-difference threshold. The column Laplacian \(L_\mathrm{w}\in\mathbb{R}^{n\times n}\) encodes the temporal graph, which links nearby time steps so as to promote temporal smoothness. The model therefore exposes parameters to be chosen: the rank \(r\), the three regularization weights \(\lambda\), \(\lambda_\mathrm{u}\), \(\lambda_\mathrm{w}\) of 5 , and the graph construction method and its parameters.

In Belgium↩︎

For our experiments, we used a dataset provided by the Royal Meteorological Institute of Belgium (RMI), which includes temperature measurements recorded every ten minutes across 96 weather stations in Belgium as represented in 1. The analysis covers the period from January 1, 2020 to December 31, 2023. Results are summarized in ¿tbl:table:irm95results95all?.

\begin{table}[h!]

\vskip 0.10in
\begin{small}
\begin{sc}
\begin{tabular}{lcccl}\toprule
    & \multicolumn{2}{c}{Block} & \multicolumn{2}{c}{Spread}
    \\\cmidrule(lr){2-3}\cmidrule(lr){4-5}
        Methods & RMSE  & Time  & RMSE  & Time\\\midrule
    IDW         &0.64   &0.14   &0.57   &0.19   \\
    PCA         &0.61   &0.25   &0.68   &0.35   \\
    MissForest & 0.59 & 10.1 & 0.56 & 9.6 \\\midrule
    LMaFiT      &0.81   &4.9    &0.76   &4.6    \\
    SoftImpute  &0.53   &3.1    &\textbf{0.41}&2.9  \\
    RTRMC       &0.49   &3.4    &0.46   &3.7    \\\midrule
    GRALS       &0.53   &8.0    &0.63   &12.5   \\
    GR-RTRMC    &\textbf{0.45}  & 8.6   &0.43   &   9.6 \\\bottomrule
\end{tabular}
\end{sc}
\end{small}
\vskip\caption{Average RMSE (in $^\circ$C) and average computational time (in seconds) on the test set for every method considered. Results for the two types of missing data generation on the Belgian weather dataset (RMI)}\label{table:irm95results95all} -0.15in

\end{table}

GR-RTRMC achieves the best overall performance on the Belgian dataset, with the lowest RMSE of 0.45°C in the Block scenario. In the Spread scenario, SoftImpute achieves the lowest RMSE (0.41°C), while GR-RTRMC remains competitive (0.43°C). 2 shows the completion quality across multiple test stations, revealing that while GR-RTRMC has the lowest average RMSE, it is not systematically superior at every individual station. The violin plot (3) illustrates the distribution of completion errors across the test set, showing the variability in method performance.

In France↩︎

For our experiments, we used a dataset provided by Météo-France,¹ which includes temperature measurements recorded every six minutes across 272 weather stations in the northwest (NW) region and 545 stations in the southeast (SE) region of France. The data covers the three-year period from 2016 to 2018. The two regions are represented as rectangular areas in 4.

\begin{table}[h!]

\vskip 0.10in
\begin{small}
\begin{sc}
\begin{tabular}{lcccl}\toprule
        & \multicolumn{2}{c}{Block} & \multicolumn{2}{c}{Spread}
        \\\cmidrule(lr){2-3}\cmidrule(lr){4-5}
            Methods & RMSE  & Time  & RMSE  & Time\\\midrule
        IDW         &0.80   &0.21   &0.74   &0.19   \\
        PCA         &0.77   &0.24   &0.67   &1.1    \\
        MissForest  &0.70   &14.3   &0.56  &13.8   \\\midrule
        LMaFiT      &0.94   &9.7    &1.01   &9.3    \\
        SoftImpute  &0.59   &7.3    &0.61   &6.7    \\
        RTRMC       &0.67   &8.2    &0.58   &4.9    \\\midrule
        GRALS       &0.62   &20.1   &0.68   &19.3   \\
        GR-RTRMC        &\textbf{0.57}  &10.1   &\textbf{0.54}  &13.1
        \\\bottomrule
    \end{tabular}
\end{sc}
\end{small}
\vskip\caption{Average RMSE (in $^\circ$C) and average computational time (in seconds) on the test set for every method considered. Results for the two types of missing data generation on the French weather dataset (M\'{e}t\'{e}o-France)}\label{table:meteofrance95results95all} -0.15in

\end{table}

On the France dataset, GR-RTRMC demonstrates consistent performance with an RMSE of 0.57°C (Block) and 0.54°C (Spread), outperforming most baselines including GRALS (¿tbl:table:meteofrance95results95all?). While SoftImpute shows interesting results, particularly in the Block scenario, GR-RTRMC achieves the best overall performance in the Spread scenario, confirming the effectiveness of graph regularization. The advantage of graph-based methods becomes more pronounced at larger scales, as the French dataset with its higher number of stations provides richer spatial relationships for the graph regularization to exploit.

3.3 MovieLens↩︎

The MovieLens 100K dataset [29] is a well-established resource in the realm of recommender systems, containing 100,000 ratings from 943 users for 1682 movies. The data was collected during the seven months from September 19, 1997 to April 22, 1998. Each user in the dataset has rated at least 20 movies, ensuring sufficient user-item interactions for meaningful analysis. The MovieLens 1M dataset will also be used, which contains 1,000,029 ratings from 6040 users for 3952 movies. In both cases, the data matrix has rows corresponding to users and columns to movies, with the row graph encoding user similarity and the column graph encoding movie similarity.

For the MovieLens experiments, we employ the same matrix completion methods as for the weather datasets (LMaFiT, SoftImpute, RTRMC, GRALS, GR-RTRMC), but replace IDW and PCA with User-Based KNN, which is more appropriate for collaborative filtering.

We construct a user–user similarity graph via \(k\)-nearest neighbors (k-NN) on user rating vectors. Both the neighborhood size \(k\) and the similarity function (e.g., Euclidean, cosine, Jaccard, Manhattan) are treated as hyperparameters and selected by validation. In practice, moderate \(k\) values and standard similarity metrics yield comparable performance, indicating robustness to these choices.

We apply the same procedure to build a movie–movie (item–item) similarity graph from item rating vectors, with \(k\) and the similarity function treated as hyperparameters.

Results↩︎

\begin{table}[t!]

\vskip 0.10in
\begin{small}
\begin{sc}
\begin{tabular}{lcccl}\toprule
    & \multicolumn{2}{c}{ML 100K} & \multicolumn{2}{c}{ML 1M}
    \\\cmidrule(lr){2-3}\cmidrule(lr){4-5}
        Methods & RMSE  & Time  & RMSE  & Time\\\midrule
    User-Based KNN& 0.969&  23.4    &   0.961   &  209.4    \\\midrule
    LMaFiT      &0.993  &  5.5  &   0.965   &   68.9    \\
    SoftImpute  &0.973  &  4.6  &   0.932   &   30.1    \\
    MissForest & 0.987 & 3.8 & 0.956 & 24.1\\
    RTRMC       &0.966  &  4.9  &   0.945   &   39.2    \\\midrule
    GRALS       &0.951  &  10.4 &   0.981   &   71.2    \\
    GR-RTRMC    &\textbf{0.942}&  17.3  &   \textbf{0.913}  &   79.5    \\\bottomrule
\end{tabular}
\end{sc}
\end{small}
\vskip\caption{Average RMSE and average computational time (in seconds) on the test set for every method considered. Results for the MovieLens datasets.}\label{table:collaborative95results95all} -0.15in

\end{table}

The MovieLens experiments demonstrate the effectiveness of our approach in collaborative filtering scenarios. As shown in ¿tbl:table:collaborative95results95all?, GR-RTRMC achieves the best performance on both datasets. The method obtains RMSE values of 0.942 for MovieLens 100K and 0.913 for MovieLens 1M, outperforming all baseline approaches including GRALS and traditional collaborative filtering methods like User-Based KNN.

3.4 Computational Complexity↩︎

While GR-RTRMC achieves superior accuracy across nearly all datasets, this performance comes at the cost of increased computational time. As shown in the timing results, GR-RTRMC consistently requires 2-3 times longer computation time than RTRMC due to the additional graph regularization.

For instance, on the Belgian dataset, GR-RTRMC requires 8.6-9.6 seconds compared to RTRMC’s 3.4-3.7 seconds. Similar patterns are observed across the French and MovieLens datasets. Despite this computational overhead, the method remains practical for real-world applications where accuracy improvements justify the additional computation time.

3.5 Case analysis↩︎

When conducting the numerical experiments with the Belgian weather data, we observed isolated cases of abnormally high RMSE at a few stations that appeared only in exceptional runs. We analysed in more detail two of such cases.

Case 1: Storm issue↩︎

In one particular imputation week, the reconstructed temperature at some stations contained implausible, narrow oscillations. An illustrative example of this artifact is shown in 5.

To understand why this week produced oscillations at multiple stations, we analyzed hourly radar overlays and identified an exceptional convective front crossing Belgium (6). Its passage time was spatially heterogeneous, with stations tens to hundreds of kilometers apart affected at different hours (see [fig:selected_stations_map,fig:temperature_analysis_subplots]).

The station map and time series clearly show this heterogeneity: four stations across Belgium all record a sharp temperature drop, but this occurs at different times, staggered by roughly 1–2 hours ([fig:selected_stations_map,fig:temperature_analysis_subplots]). A low-rank model that assumes time-synchronous behavior, especially when further constrained by spatial graph regularization, cannot capture these temporary phase shifts. When forced to reconcile misaligned signals, the method produces the oscillation observed during the storm.

Reducing the confidence assigned to known data during the storm window, by lowering \(C_{i,j}\) for all \(i\) and for \(j\) within the global storm period, removed this behavior and restored a plausible imputation, as shown in 8.

Case 2: Rain issue↩︎

Rain from scattered thunderstorms creates trouble for accurate imputation of air temperature data. Nearby weather stations often see very different amounts of rainfall, even where they are close together, as illustrated in 9 which shows the temperature, precipitation, solar radiation and predicted temperature using GR-RTRMC for the station 6414. In this case, station 6414 observed heavy rain from 3pm to 4pm on June 23, 2023. The rain caused temperature to drop quickly. However, a weather station about 10 km away (10115001) did not experience any rain during that time. This focused influence makes it tricky to fill in missing data using only temperature information. Nearby locations do not share similar weather conditions at that moment.

We filled in gaps in data from station 6414 during the rainfall window, but the completion disregarded the local cooling and produced estimates more consistent with neighboring stations that did not experience any precipitation event.

This situation shows that relying on temperature alone does not fully explain what happens when localized weather creates disturbances. If two weather stations are close together, they might be on either side of a small area with heavy rain. This is illustrated in 10 which shows hourly precipitation accumulation observed by the weather radars on June 23, 2023 at 14:00, 15:00, 16:00 and 17:00. These radar images show that these precipitation events spread eastward, passing through station 6414 but not reaching station 10115001.

Addressing these situations would require incorporating additional meteorological variables or mechanisms designed to detect and preserve such local anomalies, and to treat them separately during the imputation.

4 Conclusions↩︎

We have presented a novel approach called Graph-Regularized Riemannian Trust-Region Matrix Completion (GR-RTRMC). The method extends RTRMC by incorporating graph regularization to exploit structural relationships in the data. Experiments on meteorological and collaborative filtering datasets demonstrate the effectiveness of our approach.

In the meteorological case, highly localized phenomena such as scattered precipitation events complicate the imputation for matrix completion-based methods. These local disturbances violate the spatial coherence assumptions underlying low-rank models, highlighting the inherent limitations of approaches that rely solely on matrix completion techniques.

Several directions for future work emerge from this study. Incorporating additional meteorological variables (e.g., precipitation, humidity, solar radiation) into the model could help detect and handle localized weather events more effectively. Developing mechanisms to automatically detect anomalous periods and reduce the confidence assigned to observations during such events could improve robustness without manual intervention. Finally, alternative graph construction strategies could be explored to further improve completion accuracy.

Computational resources have been provided by the supercomputing facilities of the Université catholique de Louvain (CISM/UCL) and the Consortium des Équipements de Calcul Intensif en Fédération Wallonie Bruxelles (CÉCI) funded by the Fond de la Recherche Scientifique de Belgique (F.R.S.-FNRS) under convention 2.5020.11 and by the Walloon Region.

Statements and Declarations↩︎

This work was funded by the Fonds de la Recherche Scientifique–FNRS under Grant no T.0001.23.

The authors have no competing interests to declare that are relevant to the content of this article.

The code that produced the experimental results is available at https://github.com/bloucheur/CAM-2026.

References↩︎

---

1. https://meteonet.umr-cnrm.fr/↩︎