Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
53,901 characters · 14 sections · 45 citation commands
Learning Nonlinear Factor Models with Unknown Monotone Links from Incomplete and Noisy Data
Keywords: Non-linear Factor models, Reproducing Kernel Hilbert Space, Block-Coordinate Descant, Incomplete Data
Advances in data collection technologies have led to the increasing availability of high-dimensional data, which are typically characterized by a large number of cross-sectional units and a high-dimensional feature space. Analyzing such data poses significant statistical and computational challenges. Factor models provide an effective dimension-reduction framework to address these challenges. Originally introduced in the field of psychometrics, factor models have since become a fundamental tool for summarizing and modeling complex, high-dimensional data sets, with applications in statistics, economics, and other data science disciplines ding2021high,fan2021recent,barigozzi2023quasi.
At their core, factor models assume that high-dimensional observations are generated from an underlying low-dimensional latent structure, often represented by a low-rank matrix. In high-dimensional settings, this low-rank structure acts as a natural form of regularization and often leads to more interpretable representations, particularly in scientific applications. Beyond traditional statistical modeling, low-rank assumptions play a central role in many machine learning problems. Prominent examples include matrix regression negahban2011estimation,chen2015fast, rank-$r$ principal component analysis johnstone2009consistency,birnbaum2013minimax, low-rank and sparse matrix decomposition candes2011robust,chandrasekaran2011rank,hsu2011robust, matrix completion candes2010matrix,sun2016guaranteed,zheng2016convergence, and recommendation systems alshbanat2025survey,sankagiri2025recommendations. These problems can be unified under the following optimization problem with $\mathcal{L}$ as a loss function and $r\ll n,T$,
Classical (linear) factor models, which have been extensively studied over the past decades, assume the following representation for the observed vector at time $t$:
where $y_t, u_t\in\mathbb R^{n}$, $f_t$ is the $t$-th column of $F$, and $u_t$ is the error term, commonly referred to as the idiosyncratic component, which is assumed to be uncorrelated with both $B$ and $f_t$. The matrices $B$ and $F$ denote the latent loading factor and latent factors, respectively, and are the primary parameters of interest in literature of linear factor models. Within this framework with Gaussian errors, the loss function $\mathcal{L}$ is often a quadratic loss which leads to the following optimization: $\min_{B,F} \sum_{i,t}(y_{i,t}-b_i^\top f_t)^2+\text{Re}(B,F)$, where $\text{Re}(B,F)$ denotes a regularization terms that could impose additional structure on the factors and loading factors.
Linear factor models are widely used for estimating high-dimensional covariance and precision matrices. In particular, suppose that the factors $\{f_t\}$ in (ref) are realizations of a multivariate distribution with unknown stationary covariance $\text{cov}(f_t)$. It is straightforward to see from (ref) that the covariance matrix of the observed data can be expressed as $\hat\Sigma\!:=\!\text{cov}(y_t)\!\!=\!\!B^\top \text{cov}(f_t) B\!+\!\text{cov}(u_t)$, where the first term is low-rank and the second is often sparse, reflecting the weak dependence of idiosyncratic errors. A central objective in the linear factor model literature is to estimate $\text{cov}(y_t)$ accurately, especially in high-dimensional settings where the number of variables $n$ exceeds the number of observations $T$. Consequently, different methods have been proposed to exploit the factor structure of $\hat\Sigma$ in order to produce more accurate estimators. For instance, formulating the covariance matrix estimation under factor models as a low-rank plus sparse matrix decomposition problem barratt2023covariance. Prominent application of this framework include portfolio management and risk assessment in financial economics de2021factor, graphical models meinshausen2006high,zhao2024high, high dimensional classification hastie2009elements, covariance estimation kereta2021estimating,sihag2022covariance, and testing the capital asset pricing model sentana2009econometrics.
A key assumption underlying linear factor models is that the observed data are complete; that is, $y_{i,,t}$ is observed for all $i\in[n]$ and $t\in[T]$. In practice, however, this assumption is often violated due to missing observations. One contribution of this work is to relax this assumption, allowing entries to be observed only at random indices and random time points.
While successful in many applications, linear factor models cannot capture nonlinear dynamics in real-world data. As an illustrative example, consider the problem of learning from comparison data, where the goal is to estimate users’ preferences over items based on observed pairwise comparisons. A widely used approach models such comparisons using probabilistic choice models, such as the Bradley–Terry–Luce model maystre2015fast,negahban2012iterative,sankagiri2025recommendations. In this setting, the probability that user $u$ prefers item $i$ over item $j$ is given by $\sigma(x_{u,i}\!-\!x_{u,j})$, where $\sigma$ is the sigmoid link function. It is further assumed that users and items are embedded in a low-dimensional latent space $\mathbb R^r$. Interpreting $B$ as the matrix of user feature vectors and $F$ as the matrix of item feature vectors, implies $x_{u,i}\!-\!x_{u,j}\!=\!b^\top_u(f_i\!-\!f_j)$. Consequently, the recommendation problem reduces to learning $B$ and $F$ by optimizing a loss function of the form $\mathcal{L}(\sigma(\langle S_k,B^\top F\rangle))$, where $S_k$ is the sampling matrix encoding the random observed comparisons sankagiri2025recommendations. A crucial distinction in this setting is that the observations are no longer linear due to the presence of the link function. To manage this complexity, existing works typically assume that either this function is known a priori saha2021optimal,sankagiri2025recommendations or the latent factors are observed caner2025deep. In contrast, this work strengthens the existing literature by removing both of these assumptions.
Linear Factor Models. As noted earlier, linear factor models have been extensively studied in the context of covariance matrix estimation. For example, fan2013large introduce the Principal Orthogonal Complement Thresholding (POET) estimator, which exploits a low-rank plus sparse decomposition of the covariance matrix. caner2023sharpe propose a node-wise regression approach for covariance matrix estimation and establish consistency results for precision matrix estimation under a linear factor model. An alternative is the shrinkage-based methods. Linear and nonlinear shrinkage estimators, developed by ledoit2012nonlinear, ledoit2022power, improve the sample covariance matrix by shrinking it toward a structured target. In the linear framework, the estimator takes the form $\lambda I + (1-\lambda)\hat\Sigma$, where $\lambda$ is the shrinkage intensity. In contrast, nonlinear shrinkage estimators apply data-driven, nonlinear transformations to each eigenvalue of $\hat\Sigma$ individually.
Recently, fan2024factor propose the Factor-Augmented Sparse Throughput (FAST) model, which combines latent factors with sparse idiosyncratic components for nonparametric regression. Using diversified projections to estimate the latent factors, the authors employ truncated deep ReLU networks for nonparametric factor regression. This approach, however, relies on a linear factor model assumption with known factors. Beyond estimation, there is a growing literature on the statistical properties and applications of linear factor models; see, for example, filipovic2024fundamental,fletcher2024examination. Dynamic linear factor models extend static linear models by allowing factor dependence to operate through time-varying filters rather than fixed loading factors, in a manner analogous to autoregressive models breitung2006dynamic, altissimo2010new. Such dynamic models have been applied to high-dimensional time-series data, including structural economic analysis and macro-financial modeling barigozzi2016non.
Matrix Completion & Factorization. One component of our problem, namely, recovering the factors $F$ and loadings $B$ from their subsampled and masked (due to the link function) product, that is, from a small subset of entries of $\phi(B^\top\! F)$, is closely related to the matrix completion problem. Matrix completion concerns the recovery of a low-rank matrix from a small subset of its entries, potentially corrupted by noise. The low-rank structure is typically enforced either through nuclear-norm regularization candes2009exact, candes2010matrix, negahban2012restricted or via an explicit matrix factorization formulation. The latter is motivated by the observation that any low-rank matrix $X$ can be factorized as $X\!=\! U^\top V$ mnih2007probabilistic. Although this reformulation leads to a nonconvex optimization problem, it is computationally more efficient and has demonstrated strong empirical performance on real-world datasets koren2009matrix. This empirical success has motivated several theoretical studies on the convergence properties of the resulting nonconvex optimization problem keshavan2010matrix,chen2015fast,sun2016guaranteed, zheng2016convergence. A common technical ingredient across this literature is a concentration result by candes2009exact, which relies on two key assumptions: (i) the ground-truth matrix is incoherent, ensuring that no single row or column dominates the matrix, and (ii) the observed entries are sampled uniformly at random. Moreover, these methods require that the iterates remain incoherent throughout the optimization process. To enforce this property, prior works introduce either explicit regularization sun2016guaranteed or projection steps chen2015fast,zheng2016convergence, ma2018implicit. In this work, we adopt a similar strategy to recover the latent factors and loadings and thus rely on the same assumptions with the additional challenge of simultaneously learning the unknown link function.
Nonlinear Factor Models. To overcome the linearity assumption, caner2025deep consider a nonlinear factor model $y_t\!=\!\varphi(f_t)\!+\!u_t$ (slightly different from our model in (ref)) for asset returns in large portfolios. In their model the link function $\varphi$ in unknown but the responses $\{y_t\}$ and the factors $\{f_t\}$ are observed. As a result, their problem reduces to a regression task, which they solve using neural networks and prove the consistency of their solution. In fact, several works study similar formulations. For example, the multi-index model aims to recover both an unknown index space $B$ (loading factors in our formulation) and an unknown link function from observed pairs $(y_t,f_t)$. Multi-index models have been widely used to mitigate the curse of dimensionality in high-dimensional regressions li2022dimension, klock2021estimating, steffen2025pac. Notably, klock2021estimating establish a connection between multi-index models and sufficient dimension reduction methods lee2013general, li2018sufficient, and propose an estimator for the index space along with sharp concentration bounds.
The nonlinear factor model considered in this paper, e.g., Equation (ref) has also been studied by chen2014nonlinear, chen2021nonlinear, with the key distinction that the nonlinear link function is assumed to be known in those works. Specifically, chen2014nonlinear study estimation and inference in semiparametric nonlinear panel single-index models via maximum likelihood. By proposing an iterative two-step likelihood maximization procedure and showing that the objective is concave in each step, they establish convergence to a local optimum for this class of models. In contrast, chen2021nonlinear develop an EM-type algorithm, building on chen2016estimation, to jointly estimate the factors and loading factors. However, neither of these works provides theoretical guarantees for the convergence of their proposed algorithms.
\paragraph{Our Contributions} We study nonlinear factor models with an unknown link function and propose a projected gradient–based method that jointly learns the latent factors, loading factors, and the link function. We consider several observation settings, including complete, random, noisy, and noiseless. The random observation setting, which allows for randomly missing entries, is more realistic than even the classical linear factor model assumption of fully observed data. By assuming a nonparametric link function that lies in an RKHS with smooth and bounded kernels, and imposing mild incoherence conditions on the factors and loadings (which are standard assumptions for handling missing data, borrowed from the matrix completion literature), we provide a theoretical analysis of the proposed algorithm. While the primary contribution of this work is to lay the foundation for learning nonlinear factor models from incomplete and noisy data, we also present experimental results that support our theoretical findings and analyze the sensitivity of the method with respect to its hyperparameters.
There is a finite number of factors, $\{f^*_1,...,f^*_T\}$ and finite number of loading factors $\{b^*_1,...,b^*_n\}$ such that $f_t, b_i\in\mathbb R^r$ for all $t\in[T]$ and $i\in[n]$ and $n,T\gg r$. The date generating process is given by
where $\phi^*:\mathbb R\rightarrow\mathbb R$ is the link function and $u_{i,t}$ denotes the error, also referred to as the idiosyncratic component, and is assumed to be uncorrelated with the factors. Errors are distributed according to a distribution $P$ with variance bounded by $\sigma^2$. Depending on wether $\sigma^2$ is positive or zero, we consider the noisy and noiseless settings, respectively. $y_{i,t}$ is the observed response for the $i$-th individual at time $t$.
We distinguish between two settings depending on the observation pattern: complete observation, where $y_{i,t}$ is observed for all $i\in[n]$ and $t\in[T]$, and random observation, where observations are available only for random subsets of $[n]\times [T]$. The random observation setting is intended to model realistic scenarios with randomly missing data.
To have a unified model for both complete and random observation settings, we use $y_k$ to denote the $k$-th observation,
where $X^*= (B^*)^\top F^*\in\mathbb R^{n\times T}$, $B^* = [b_1,...,b_n]$, and $F^* = [f_1,...,f_T]$. In the above equation, $M (\leq nT)$ denotes the number of observations. Indices $i_k\in[n]$ and $t_k\in[T]$ represent the factor and loading factor, respectively in the $k$-th observation. We also use $\Omega:=\{(i_k,t_k): \ k\in [M]\}$ to denote the set of observed indices. Then, the complete observation setting is when $\Omega=[n]\times [T]$.
In this paper, we focus on recovering $(\phi^*, B^*, F^*)$ from the observed data $\{y_k\}$ and observed indices $\Omega$. To this end, let $U^* \Sigma^*_r (V^*)^\top$ be the rank-$r$ SVD of matrix $X^*$, that is $U^*\in\mathbb R^{n\times r}, V^*\in\mathbb R^{T\times r}$ and $\Sigma^*\in\mathbb R^{r\times r}$ such that $(U^*)^\top U^*=(V^*)^\top V^*=I_r$ and $\Sigma^*$ is a diagonal matrix with entries $\sigma_1^* \geq \ldots \geq \sigma_r^* > 0$.
We define $Z^*$ and $Y^*$ as
Using this notation, we can rewrite the $k$-th observation in (ref) as follows
where $\langle A, B \rangle :=\text{tr}(A^\top B)=\sum_{i,j}a_{i,j}b_{i,j}$ denotes the matrix inner product and $A_k$ is the sampling matrix for the pair $k:=(i_k,t_k)\in\Omega$ defined by
where $e_{i_k}$ and $\tilde{e}_{t_k}$ denote unit vectors in $\mathbb R^n$ and $\mathbb R^{T}$, respectively. In the above equation, $0$ denotes matrices with all entries zero of the appropriate size. In the random observation setting, $\{A_k\}$ are i.i.d. random matrices of form (ref), with the index $i_k$ being chosen uniformly at random from $[n]$, and $t_k$ being chosen uniformly at random from $[T]$.
Given the relation between matrices $X^*$, $Y^*$, and $Z^*$, estimating the ground-truth factors and loading factors is equivalent to estimating $Z^*$. The major advantage of this reformulation is that it reduces the number of parameters from $nT$ in $X^*$ to $(n + T)r$ in $Z^*$.
\paragraph{Condition Number} We define $\kappa:=\sigma^*_1/\sigma^*_r$ as the condition number of the data. Note that $\kappa$ is also the condition number of $Z^*$ as the singular values of $Z^*$ are $\sqrt{2\sigma^*_1}, \ldots \sqrt{2\sigma^*_r}$.
\paragraph{Incoherence} For matrix $Z^*$ with $n+T$ rows, let $\|Z^*\|_{2, \infty}$ denote the maximum of the $\ell_2$ norm of its rows and let $\|Z^*\|_{F}$ denote the Frobenius norm of $Z^*$. Define the incoherence parameter of the matrix as
In principle, $\mu$ can take values from $1$ to $n+T$. However, the sample complexity worsens with $\mu$.
\paragraph{Link Function} The unknown link function $\phi^*(\cdot)$ is assumed to be monotone and to belong to a reproducing kernel Hilbert space, defined below, associated with the kernel $K(\cdot,\cdot)$. For details see scholkopf2018learning.
Given constants $\xi, \Xi>0$, we also define $\mathcal{H}_{\xi,\Xi}\subset\mathcal{H}$ as the subset of functions in $\mathcal{H}$ with bounded derivatives, i.e.,
Since this set is convex, the projection onto it is well-defined $\mathcal{P}_{\mathcal{H}_{\xi,\Xi}}(\phi):=\arg\min_{f\in\mathcal{H}_{\xi,\Xi}}\|f - \phi\|_\mathcal{H}$. As there is no closed form solution to this projection, a practical approach is to discretize the support of $\phi$ and impose inequality constraints to ensure that its derivative remains between $\xi$ and $\Xi$. This reduces the problem to a quadratic program with box constraints. Nevertheless, our empirical study show that Algorithm (ref) converges even without any projections.
Under this assumption, the infinite-dimensional problem of estimating the link function can be reduced to a finite-dimensional one via the representer theorem scholkopf2018learning. This approach is standard in RKHS-based methods and is widely used in the machine learning literature when the parameter of interest is a continuous function.
Recall that our goal is to estimate the pair $(\phi^*,Z^*)$. We do so by maximizing the log-likelihood. Equivalently, we formulate a loss function in terms of the negative log-likelihood, and minimize this function using a gradient-based method. Under the assumption that the errors are independent and identically distributed according to a sub-Gaussian distribution, the negative log-likelihood function is given by
\paragraph{Tikhonov Regularization} Function $\phi$ belongs to $\mathcal H$, which can have many degrees of freedom, possibly infinite. To favor smoother and simpler representations within $\mathcal H$, it is common to regularize the objective function using $\|\phi\|_\mathcal{H}$, also known as Tikhonov regularization scholkopf2018learning. This is because $\|\phi\|_\mathcal{H}$ measures the “roughness” of $\phi$, and the regularization penalizes large coefficients in the RKHS expansion of the solution, thereby encouraging smoother solutions. Another important consequence is that this regularizer enables the use of the representer theorem, which guarantees that the optimal solution is unique and can be expressed as a finite linear combination of kernels evaluated at the observed data points.
\paragraph{Scale Invariance} The generative model, and consequently the log-likelihood function, is invariant to a certain transformation in the $Z$ matrix. Below, we explore such symmetry and introduce a regularizer to favor solutions that have two factors with comparable second-order magnitudes.
Note that the mapping from a matrix $X$ to its factor representation $Z=[U; V]$ is generally non-identifiable. In particular, for any invertible $P\in\mathbb{R}^{r\times r}$, the reparameterization $\tilde Z=[UP^\top ; VP^{-1}]$ yields the same matrix $X$, and therefore the same likelihood, i.e., ${\mathcal{L}}(\phi, Z)={\mathcal{L}}(\phi, \tilde Z)$. To distinguish “imbalanced” factorizations from “balanced” ones, we augment the loss with the regularizer $\|U^\top U - V^\top V\|_F^2$, which encourages solutions in which the two factors have comparable second-order magnitudes. In compact form, this regularizer can be written as
Overall, the loss function that we use in this work to recover the unknown parameters is given by
where $\lambda$ and $\alpha$ are the regularization coefficients, which will be specified later.
\paragraph{Equivalent Solution Set} Beyond scale invariance, the problem exhibits rotational symmetry in the latent dimension. Let $R \in \mathbb{R}^{r \times r}$ be any orthogonal matrix, i.e., $RR^\top= I_r$. Then, $ZR=[UR; VR]$ produces exactly the same loss function, i.e., $\tilde{\mathcal{L}}(\phi, Z)=\tilde{\mathcal{L}}(\phi, ZR)$. Hence, the true features cannot be uniquely recovered; they are determined only up to an orthogonal rotation. Accordingly, we define the equivalence class of ground-truth feature matrices as
Using this equivalent solution set, we can measure the "goodness" of an arbitrary $Z$ via $\|\Delta(Z)\|$, where
As mentioned in the previous section, we apply a gradient-based method to optimize the loss function in (ref) which is presented in Algorithm (ref). However, for proving theoretical guarantees, we need to use projected gradient descent. Notably, two projection steps are designed: first projecting $Z$ onto a set of "incoherent matrices" $\mathcal{C}$, defined below and then projecting $\phi$ onto $\mathcal{H}_{\xi,\Xi}$ (defined in (ref)).
Note that this is a convex and closed set containing matrices that are nearly as incoherent as $Z^*$ when $\|{Z_0}\|_F \approx \|{Z^*}\|_F$. For any $Z$, the projection of $Z$ onto $\mathcal{C}$ is obtained by clipping the rows of $Z$ as follows:
where $\beta = (4/3)\sqrt{({\mu}/{n+T})}\|{Z_0}\|_F$. As noted in sankagiri2025recommendations, the projection onto $\mathcal C$ is often unnecessary in practice; however, it is essential for deriving the theoretical guarantees. The second projection ensures that the learned link function preserve monotonicity in $\mathcal H$.
In this algorithm, the gradients are given by
where $x_k:=\langle A_{k}, ZZ^\top\rangle$ and $g_k:=y_{k} - \phi\big(x_k\big)$.
\paragraph{Initialization} Let $\mathcal{B}(\epsilon):= \{Z: \|{\Delta(Z)}\|_{F}^2 \leq \epsilon\sigma^*_r\}$ denote a neighborhood around the true solution. We assume that $Z_0\in\mathcal{B}(\epsilon)$ for some $\epsilon>0$. Following the initialization schemes in zheng2016convergence, sankagiri2025recommendations, we initialize $Z_0$ via a rank-$r$ SVD of the data matrix, with missing entries filled with zeros.
In this section, we present our theoretical results under different settings. We begin with the noiseless observation setting, which is further divided into complete and random cases. We then present our results for the noisy setting.
Recall that in this setting, the idiosyncratic components have zero variance, i.e., $y_{i,t}=\phi^*((b_i^*)^\top f_t^*)$. We begin by showing the properties of the loss function.
This follows directly from the definition of the loss function.
Proof is in Appendix (ref). Next, we show that both $\nabla_\phi \tilde{\mathcal{L}}$ and $\nabla_Z \tilde{\mathcal{L}}$ are Lipschitz continuous with respect to $Z$ for any fixed $\phi$. To establish this result, we assume that the kernel function itself is Lipschitz continuous and smooth. This assumption is satisfied, for example, by RKHSs induced by Gaussian kernels.
Proof is in Appendix (ref). In the paper, we can simply set $\tilde L_{\max}=\tilde{\mathcal{L}}(\phi_0,Z_0)$ when applying Algorithm (ref).
Proof is in Appendix (ref). As noted earlier, the loss function is non-convex in $Z$. However, when $\phi$ is fixed, the subsequent results characterize its local properties, which in turn facilitate the convergence analysis of the Algorithm (ref).
Proof is in Appendix (ref). In the random observation setting, the required number of observations $M$ given in (ref) matches the order established in matrix completion results candes2009exact, zheng2016convergence.
Proof is in Appendix (ref). Note that both Lemmas (ref) and (ref) assume that $Z\in\mathcal B(\epsilon)$. This is ensured by the initialization assumption, i.e., $Z_0\in\mathcal B(\epsilon)$ together with the result of Proposition (ref) that $\|\Delta_t\|_F$ remains bounded. Using the above results and by introducing a Lyapunov potential function, see below, we establish the convergence of Algorithm (ref). To define the Lyapunov function, we first introduce the $\phi$-subproblem minimizer. Specifically, for any fixed $Z$, let $ \phi^\sharp(Z):= \arg\min_{\phi\in {\mathcal{H}}} \tilde{\mathcal{L}}(\phi, Z). $ Since the objective function is strongly convex with respect to $\phi$, (Lemma (ref)), then $\phi^\sharp(Z)$ is the unique stationary point of $\tilde{\mathcal{L}}$ for a given $Z$. Define the bias radius as
We do not require a closed-form expression for $\chi$, only its existence and finiteness. But if a model-specific bound is available for $\chi$ (e.g., in terms of kernel eigenvalues and $\alpha$), it can be used to tighten the final error bounds.
\paragraph{Lyapunov-based Analysis} Define the Lyapunov potential function as
where $E_t:=\|\phi_t-\phi^\sharp(Z_t)\|_\mathcal{H}^2$ and $D_t:=\|\Delta(Z_t)\|_F^2$. The constant $\gamma>0$ will be chosen explicitly later. The term $D_t$ quantifies the closeness of $Z_t$, at iteration $t$, to the equivalent solution set $\mathcal E$. The term $E_t$ measures the closeness of the $t$-th estimate of the link function to the subproblem minimizer corresponding to $Z_t$. Next result shows that $D_t$ remains bounded and $E_t$ converges to zero.
Proof is in Appendix (ref). The following convergence result establishes a recursion for the potential function, implying that when $\phi_t\to\phi^*$, the potential function, and consequently $D_t$, converges to zero.
A proof including the constant terms are in Appendix (ref).
\paragraph{Regret Analysis} Herein, we define and bound the corresponding regret of Algorithm (ref). Let $\tilde R_{\bar T}:=\tfrac{1}{\bar T}\sum_{t=1}^{\bar T} \big(\tilde{\mathcal{L}}(\phi_t,Z_t)\!-\!\min_{\phi\in \mathcal{H}_{\xi,\Xi},Z\in\mathcal{C}}\tilde{\mathcal{L}}(\phi,Z)\big)$. This quantity measures the cumulative regret of Algorithm (ref) relative to the global minimum of the loss function. However, since the loss function is non-convex with respect to $Z$ for any $\phi$, convergence of the algorithm’s iterates to the global minimum cannot be guaranteed. In particular, when the initialization is closer to a local minimizer than to the global one, the algorithm may converge to such a local solution. As a result, $\tilde R_T$ grows linearly with $T$. Therefore, we introduce the following modified notion of regret with respect to the global minimum of the loss function when $Z$ is fixed,
Using this regret definition, we derive a sublinear regret bound for Algorithm (ref) (proof in Appendix (ref)).
\paragraph{Analyzing the Stationary Points.} Let $(\phi, Z)$ be a stationary point of the $\tilde{\mathcal{L}}(\phi, Z)$, i.e., the gradients of $\tilde{\mathcal{L}}$ with respect to both $\phi$ and $Z$ vanish at $(\phi, Z)$. The following result shows that, at such a point, the average residuals remain bounded and tend to zero as $\alpha$ decreases. Proof is in Appendix (ref).
Recall that in this setting, the idiosyncratic components have non-zero but bounded variance.
As the results from the noiseless setting carry over to this case with minor adjustments to account for the noise, we provide the auxiliary lemmas in Appendix (ref) and present below only the main result for the potential function.
This result implies that when $\alpha$ is small enough, $\lim\limits_{t \to \infty}\!\|\Delta_t\|$ $\!\in\mathcal{O}(\chi(Z_\infty)+\sigma)$, where $Z_\infty$ is the limit point of $Z_t$. See Appendix (ref) for details.
In this section, we evaluate the performance of Algorithm (ref) in both synthetic and real-world settings, study the effect of key parameters, and empirically validate our theoretical results.
We consider two settings: a known link function and an unknown link function. The former isolates the impact of mismatch in the link function and clarifies the role of the remaining parameters.
\paragraph{Learning the latent factors and loadings with a known link function} In this setting, no learning of the link function is required, consequently, Algorithm (ref) does not update $\phi_t$.
Effect of the noise level: We applied the algorithm only on the $Z$-block under varying noise levels, where the idiosyncratic terms are $u_{i,t}\stackrel{\text{i.i.d.}}{\sim} \mathcal{N}(0, \sigma^2)$. For different $\sigma$, Figure (ref)(right) depicts the convergence behavior for varying $\sigma$. In this experiment, the link function is known and given by the sigmoid function $\sigma(x):=1/(1+\exp(-x))$, with $M=500000$ samples. Moreover, to illustrate the dependency between the convergence point and the noise level, the algorithm is run for a fixed number of iterations, and the resulting normalized estimation error $\|\Delta_t\|_F/(\|\Delta_t\|_F|_{\sigma=0})$, where the denominator is the final error for noiseless observation ($\sigma=0$), is shown in Figure (ref)(left). This is consistent with the prediction of Theorem (ref), which implies a linear dependence of the latent factor estimation error on the noise variance when the link function is known. These results validate the tightness of the theoretical bound in the known link function setting.
Effect of the sample size: To investigate this effect, we set $n\!=\!T\!=\!1000$, $r\!=\!10$ and assume the link function is known and equal to the identity, thereby isolating any mismatch in the link. Figure (ref) shows the convergence of $\|\Delta_t\|_F$ for varying sample sizes $M\in\{1,2,..,8\}\times rn$. This results shows that as $M\geq 3nr$, the iterates converges.
\paragraph{Joint learning of latent factors, loadings, and the link function} Herein, we considered the full problem and applied Algorithm (ref). We examined the convergence behavior of the factors estimation error $\|\Delta_t\|_F$, the link function estimation error $\|\phi_t-\phi^\star\|_{\mathcal H}$, and the Lyapunov $\mathcal V_t$ across different noise levels $\sigma$ and regularization parameters $\alpha$. In these experiments, the shared parameters are $n=T=100,r=3, M=5000,\lambda=0.5,\zeta=10^{-5},\eta=10^{-4}$.
Effect of the noise level: Figure (ref) illustrates the estimation errors and potential function over iterations for various noise levels with $\alpha=10^{-3}$. Our results show that, for a fixed $\alpha$, increasing the noise level $\sigma$ leads to a degradation in both $\|\Delta_t\|$ and $\|\phi_t-\phi^\star\|_{\mathcal H}$, while the Lyapunov potential $\mathcal V_t$ converges to a neighborhood whose radius grows with $\sigma$. In the noiseless case ($\sigma=0$), the limiting value of $\mathcal V_t$ provides an empirical estimate of the constant $C_\phi(t)$ appearing in Theorem (ref). For additional plots see Appendix (ref)
Effect of regularization coefficient $\alpha$: To further investigate the role of the $\alpha$, we fix the noise level at $\sigma=0.1$ and vary $\alpha$. Figure (ref) shows the final estimation error of Algorithm (ref) as function of $\alpha$. The result aligns with our discussion in Appendix (ref) and highlights the practical impact of the regularization parameter in nonlinear factor recovery.
We consider two real-world datasets in this section to demonstrate the applicability and scalability of our algorithm. Moreover, our experimental results show that the projection steps in the algorithm do not adversely affect its scalability.
MovieLens\footnote{\href{https://grouplens.org/datasets/movielens/}{MovieLens dataset link}.} uses a $610\times9724$ user-movie matrix with 5-star ratings, yielding a total of 5,931,640 entries. Among these, 100,836 ratings are observed. We use a row-stratified 10% validation split, with 91,104 ratings for training and the remainder for validation. $\beta$ is used for projecting the $Z$ matrix. We report results at the algorithms' best checkpoints in Table (ref); for example, for Id. Link, this occurs at iteration 1000. The Val RMSE and Train RMSE correspond to the validation and training Root Mean Squared Errors (RMSE), respectively.
Id. Link and NL. Link denote the linear factor model and our model, respectively. Z proj. method performs only projection of matrix $Z$ while NL monotone performs both projections, as also indicated in the second and third columns. Clearly, the monotone method (which performs both projection) and Z proj. method (which only does $Z$ projection) have similar runtimes. This further supports our analysis that projecting the link function does not impact scalability. An important observation is the advantage of link function learning that decreases the validation RMSE compared to the identity link function.
Table (ref) shows the range of the inferred factor-loading products, i.e., the matrix $X$. For the linear model (Id. Link), the range is roughly 0 to 5, corresponding to the scale of the ratings. In contrast, for the nonlinear model, the range is much smaller, approximately [-0.2, 0.2], and the link function maps this range to the 5-star ratings through a nonlinear transformation. The inferred link functions are presented in Figure (ref).
Jester\footnote{\href{https://www.kaggle.com/datasets/vikashrajluhaniwal/jester-17m-jokes-ratings-dataset}{Jester dataset link}.} dataset dataset has 24,983 users and 100 jokes, with continuous ratings from -10 to 10. After cleaning the data, the matrix contains 2,498,300 possible user-joke entries, of which 1,810,455 are observed. We use a 10% hold-out for validation, with the remaining entries used for training. Table (ref) presents the results.
We studied the problem of jointly learning the latent factors, loading factors, and the link function under several observation settings, including complete, random, noisy, and noiseless. Our proposed learning method is a projected block coordinate descent. By assuming a nonparametric link function that lies in an RKHS with smooth and bounded kernels, and imposing mild incoherence conditions on the factors and loadings, we provided a theoretical analysis of the proposed algorithm. The analysis in this work relies on several assumptions that may be relaxed, and thus represent limitations of the present study. Moreover, we conjecture that the projection steps in our algorithm are unnecessary.
\nocite{langley00}