Panel Coupled Matrix-Tensor Clustering Model with Applications to Asset Pricing

Panel Coupled Matrix-Tensor Clustering Model with Applications to Asset Pricing

Introduction

Common factor models are central to empirical asset pricing, capturing time-series co-movement and cross-sectional return variation fama_five-factor_2015. However, estimating asset-specific betas (factor loadings) remains statistically challenging; high idiosyncratic volatility in sparse datasets often masks underlying risk signals, a problem exacerbated in segmented markets. Market segmentation leads to significant cross-sectional variation in factor risk premia hou2011factors, limiting the explanatory power of a common factor model across asset classes. To address this, patton2022risk and cong2023sparse propose different clustering methods to group assets by within-group factor loadings, revealing significant cross-sectional heterogeneity. Similarly, GIGLIO2024TestAssetsandWeakFactors show that a factor's explanatory power depends on test asset selection, as factor loadings vary across asset classes. These results highlight the potential of group-specific factor models to better explain variations in asset returns. However, these group-specific factor models often rely solely on asset returns ($Y_{t}$) and factor returns ($f_t$), neglecting asset characteristics (${\cal X}_{ t}$), which provide valuable incremental information kelly2019IPCA.

From a statistical perspective, these group-specific models can be viewed as a matrix clustering task, modeling excess returns of $p_1$ assets using $m_1$ factors with latent groups:

where ${\mbox{\boldmath $ Y$}} \in \mathbb{R}^{p_1 \times T}$ is the return matrix, ${\mbox{\boldmath $ F$}} \in \mathbb{R}^{m_1 \times T}$ are factors, ${\mbox{\boldmath $ B$}} \in \mathbb{R}^{r_1 \times m_1}$ denotes group-level loadings, and ${\mbox{\boldmath $ M$}}_1 \in \{0,1\}^{p_1 \times r_1}$ encodes asset memberships. Clustering methods, including $k$-means and spectral clustering jain2010data, vonluxburg2007tutorial, zhang2024leave, and extensions to structured/high-order data gao2022iterative, han2022exact, provide consistent recovery under suitable conditions but typically rely on a single data source.

Asset pricing researchers observe returns ${\mbox{\boldmath $ Y$}}$, factors ${\mbox{\boldmath $ F$}}$, and asset-specific characteristics that contain information beyond returns alone lettau20243d. These characteristics naturally form a tensor with latent group structure:

where {$\mathcal{X}\in \mathbb{R}^{p_1\times p_2\times T}$ collects $p_2$ asset characteristics for $p_1$ stocks over $T$ periods}, ${\mbox{\boldmath $ M$}}_2 \in \mathbb{R}^{p_2 \times r_2}$ is the membership matrix for characteristics, and ${\cal S} \in \mathbb{R}^{r_1 \times r_2 \times T}$ is a core tensor capturing cluster centroids. The shared first mode, ${\mbox{\boldmath $ M$}}_1$, provides a direct link between the outcome matrix ${\mbox{\boldmath $ Y$}}$ and the characteristics tensor ${\cal X}$. Here the $k$-mode product of ${\cal X}\in\mathbb{R}^{p_1\times p_2\times \cdots \times p_K}$ with a matrix ${\mbox{\boldmath $ U$}}\in\mathbb{R}^{r_k\times d_k}$, denoted as ${\cal X}\times_k {\mbox{\boldmath $ U$}}$, is an order $K$-tensor of size $d_1\times \cdots \times d_{k-1} \times r_k\times d_{k+1}\times \cdots \times d_K$ such that $ ({\cal X}\times_k {\mbox{\boldmath $ U$}})_{i_1,...,i_{k-1},j,i_{k+1},...,i_K}=\sum_{i_k=1}^{d_k} {\cal X}_{i_1,i_2,...,i_K} {\mbox{\boldmath $ U$}}_{j,i_k}.$

Asset characteristics are typically incorporated in two approaches. The first follows the classical portfolio‐sorting scheme Fama1992crosssection, where a small set of characteristics, size and value, for example, is used to sort individual assets into portfolios and evaluate factor models under the implicit assumption that all assets within a portfolio share the same loading vector ${\mbox{\boldmath $ \mathnormal\beta$}}$. This approach, however, relies on only a limited subset of available characteristics, depends on ad hoc sorting breakpoints, and may introduce selection bias. The second approach uses conditional factor models, which specify factor loadings as explicit functions of characteristics, as in kelly2019IPCA and gu2021autoencoder. While flexible, these conditional models use characteristics to define loadings rather than to identify latent asset groupings. While asset characteristics are known to proxy for risk exposures and expected returns, traditional sorting and conditional models fail to fully integrate this high-dimensional information into a unified grouping structure. Our paper proposes a complete integration of the return matrix and asset characteristic tensor, enhancing the accuracy of cross-sectional clustering and factor loading estimation compared to methods that rely on a single return source.

In this paper, we propose a coupled matrix-tensor clustering framework integrates the heterogeneous factor model (ref) with a characteristics tensor (ref) by enforcing a shared membership matrix (${\mbox{\boldmath $ M$}}_1$) across both data sources. Our objective is to jointly exploit the shared clustering structure of the matrix and tensor, thereby improving accuracy over single-source approaches. To achieve this, we develop a two‐stage estimation strategy. The first stage employs Panel Coupled Matrix-Tensor Spectral Clustering (PMTSC) to obtain a warm initialization, while the second stage refines the clustering via the Panel Coupled Matrix-Tensor Lloyd (PMTLloyd) algorithm. PMTSC relies on a coupled low‐rank factorization, implemented through Panel Coupled High‐Order Orthogonal Iteration (PCHOOI), which extends HOOI to settings where a tensor and a matrix are jointly modeled. PMTLloyd then iteratively updates cluster assignments using innovative orthogonal projection-based refinements. Importantly, even when applied to tensor data alone, PMTLloyd improves upon existing refinement procedures han2022exact, indicating that the algorithmic contribution is not limited to the coupled setting. The recovered clustering structure further enables accurate estimation of the group-level factor loading matrix.

Our theoretical analysis reveals that coupling strengthens the signal in the shared mode by increasing relevant matricized singular values, which relaxes signal-to-noise ratio (SNR) requirements under appropriate conditions. A particularly striking result is that PMTLloyd achieves sharp misclustering error bounds that are uniformly superior to those of single-source clustering, regardless of the relative signal strengths in the tensor and matrix components. Additionally, our error bounds feature exact constants in the exponent, significantly improving upon previous tensor co-clustering results han2022exact. For factor loading estimation, we establish convergence rates under both observed and latent factor scenarios, demonstrating that increases in either $p_1$ or $T$ enhance the estimation of the loading matrix. These bounds dominate ungrouped factor analysis results, with the notable property of guaranteeing consistency even with finite sample sizes, improving over existing group panel regression results su2016identifying. Simulation studies confirm that PMTLloyd achieves the lowest clustering error, even in weak SNR regimes, compared with other popular methods.

Empirical results further highlight the practical value of the proposed algorithms when applied to the Panel Tree portfolios of cong2025growing for the period 1980-2024. Our method yields a higher out-of-sample total $R^2$ than return-based clustering and traditional univariate or bivariate characteristic sorting methods. By jointly exploiting information in returns and asset characteristics, the method identifies sharper and more stable latent asset groups, resulting in more accurate factor-loading estimates and improved predictive performance. The clusters reveal economically meaningful patterns: differences in factor exposures align with underlying characteristics, providing interpretable links between the empirical “factor zoo” and characteristic-driven asset behavior. These findings demonstrate that coupling ${\mbox{\boldmath $ Y$}}$ and ${\cal X}$ enhances clustering precision, providing a more robust and interpretable framework for understanding cross-sectional return variation.

Related Literature

Our paper contributes to recent advances in low-rank tensor decomposition. Tucker-type models offer efficient representations for multi-way data, enabling theoretical analysis of estimation and recovery zhang2018tensor,zhang2019optimal. Recent progress in tensor clustering has delivered sharp guarantees for identifying structured groups in high-order data under various noise regimes Luo2021Sharp, Hu2022Multiway, Luo2022Tensor,lyu2023optimal, Lyu2025Optimal. We extend this line of research by introducing a coupled matrix-tensor framework that achieves uniformly superior misclustering error bounds and improved computational efficiency compared to state-of-the-art single-source tensor clustering methods.

In addition to the matrix clustering paradigm, model (ref) can be viewed from the perspective of coefficient homogeneity. Under this view, assets correspond to separate regression tasks, and the grouping structure induces equality constraints on regression coefficients within each group. This connects our framework to the literature on homogeneity pursuit, where pairwise (fused) penalties recover latent group structures among coefficients shen2010grouping, zhu2013simultaneous. Recent work extends these ideas to multi-task and panel settings with adaptive and robust formulations duan2023adaptive, cui2025do.

Our framework also draws on data fusion methodologies that incorporate auxiliary information to improve statistical efficiency, including Covariate-Assisted Sparse Tensor Completion ibriga2023covariate and Covariate-Assisted Spectral Clustering binkiewicz2017covariate. These approaches demonstrate how covariates or side information can enhance clustering and representation learning in multi-way settings. Building on this insight, our method treats the characteristics tensor as auxiliary linked data, developing a unified estimation framework that jointly exploits shared clustering structure across outcomes and characteristics. By explicitly modeling the shared latent group structure across the outcome matrix and characteristics tensor, our Panel Coupled Matrix-Tensor Clustering (PMTC) model provides a statistically efficient mechanism for information fusion, thereby enhancing recovery of the shared mode.

We also build on coupled factorization methods that jointly model multiple structured datasets. Matrix-matrix approaches lock2013joint,fan2019distributed,tang2021integrated,ma2024optimal, matrix-tensor factorizations acar2011all,de2017coupled, and tensor-tensor formulations liu2023joint,chen2025distributed provide flexible tools for capturing shared latent structures across heterogeneous sources. These methodologies illustrate the benefits of coupling information across multiple modes, a principle underlying the proposed PMTC approach.

Notation, Preliminaries and Organization

Let $[n]$ denote the set $\{1, 2, \ldots, n\}$. For a vector $x = (x_1, \ldots, x_p)^{\top}$, we define its $\ell_q$-norm as $\| x \|_q = (\sum_{i=1}^p |x_i|^q)^{1/q} $ for $ q \geq 1 $. For a matrix ${\mbox{\boldmath $ A$}} = (a_{i,j}) \in \mathbb{R}^{m \times n}$, denote its singular values as $\lambda_1({\mbox{\boldmath $ A$}}) \geq \lambda_2({\mbox{\boldmath $ A$}}) \geq...\geq \lambda_{\min\{m,n\}}({\mbox{\boldmath $ A$}}) \geq 0 $. The subspace spanned by the first $r$ left singular vectors is denoted as $ {\mbox{\boldmath $ U$}}_r = {\rm LSVD}_r ({\mbox{\boldmath $ A$}}) $, and the spectral norm is $\|{\mbox{\boldmath $ A$}}\|_2 = \lambda_1({\mbox{\boldmath $ A$}})$. We denote the $i$-th row and $j$-th column of ${\mbox{\boldmath $ A$}}$ as ${\mbox{\boldmath $ A$}}_{i:} $ and ${\mbox{\boldmath $ A$}}_{:j} $, respectively. We also use $a \wedge b = \min\{a, b\}$ and $a \vee b = \max\{a, b\}$.

For any two orthonormal matrices ${\mbox{\boldmath $ U$}}, \widehat{\mbox{\boldmath $ U$}} \in \mathbb{O}^{m \times r}$, the distance between their column spaces is measured by the spectral norm of their projection difference $\ell_2({\mbox{\boldmath $ U$}}, \widehat{\mbox{\boldmath $ U$}}) = \|\widehat{\mbox{\boldmath $ U$}} \widehat{\mbox{\boldmath $ U$}}^\top - {\mbox{\boldmath $ U$}} {\mbox{\boldmath $ U$}}^\top\|_2 = \sqrt{1 - \lambda_r({\mbox{\boldmath $ U$}}^\top \widehat{\mbox{\boldmath $ U$}})^2},$ which equals the sine of the largest principal angle between the subspaces. Let $\mathrm{vec}(\cdot)$ denote the vectorization operator. The mode-$k$ unfolding (matricization) of tensor $\mathcal{A}$ is defined as $\hbox{\rm mat}_k(\mathcal{A})$, mapping $\mathcal{A}$ to a matrix in $\mathbb{R}^{m_k \times m_{-k}}$ where $m_{-k} = \prod_{j \neq k}^K m_j$. For example, if $\mathcal{A} \in \mathbb{R}^{m_1 \times m_2 \times m_3}$, then $(\hbox{\rm mat}_1(\mathcal{A}))_{i,(j+m_2(k-1))} = (\hbox{\rm mat}_2(\mathcal{A}))_{j,(k+m_3(i-1))} = (\hbox{\rm mat}_3(\mathcal{A}))_{k,(i+m_1(j-1))} = \mathcal{A}_{ijk}.$ For a $d$-way tensor $\mathcal{A}$, we define its minimal matricized singular value as $\lambda_{\min}(\mathcal{A}) = \lambda_{\min}(\hbox{\rm mat}_i(\mathcal{A})),i = 1,...,d,$ the smallest singular value across all mode-$i$ matricizations.

The remainder of the paper is organized as follows. Section (ref) develops the PMTC methodology and presents two estimation algorithms: PMTSC and PMTLloyd. Section (ref) establishes theoretical properties and convergence guarantees for the proposed estimators. Section (ref) reports simulation evidence on clustering and factor loading estimation performance. Section (ref) applies the method to empirical asset-pricing data. Section (ref) concludes with a discussion of potential extensions.

Panel Coupled Matrix-Tensor Clustering

The Model

We consider a general Panel Coupled Matrix-Tensor Clustering (PMTC) model for the asset pricing problem. One observes a $(d+1)$-order characteristics tensor $\mathcal{X} \in \mathbb{R}^{p_1 \times \cdots \times p_d \times T}$ and a panel outcome matrix $\mathbf{Y} \in \mathbb{R}^{p_1 \times T}$. The model takes the form

or equivalent, ${\cal X}_t = {\cal S}_t \times_1 {\mbox{\boldmath $ M$}}_1 \times_2 \cdots \times_d {\mbox{\boldmath $ M$}}_d + {\mbox{\boldmath $ E$}}_t, \quad Y_{i,t} ={\mbox{\boldmath $ \mathnormal\beta$}}_{i}^\top f_t + \eta_{it}, \quad {\mbox{\boldmath $ \mathnormal\beta$}}_{i} = \sum_{k=1}^{r_1} {\mbox{\boldmath $ b$}}_{k} \cdot \mathbf{1}\{i \in \mathcal{G}_{1k}\}$, where ${\cal S} \in \mathbb{R}^{r_1 \times \cdots \times r_d \times T}$ is a low rank core tensor capturing latent block centroids, ${\mbox{\boldmath $ B$}} \in \mathbb{R}^{r_1 \times m_1}$ is a group-level factor loading matrix, ${\mbox{\boldmath $ F$}} \in \mathbb{R}^{m_1 \times T}$ contains $m_1$ observed or latent common factor processes, $\mathcal{G}_{1k}$ is the $k$-th group in the mode-1, and ${\mbox{\boldmath $ b$}}_i $ is the $i$-th row of ${\mbox{\boldmath $ B$}}$. The error terms are denoted by ${\cal E}$ and ${\mbox{\boldmath $ \mathnormal\eta$}}$. Each ${\mbox{\boldmath $ M$}}_i \in \{0,1\}^{p_i \times r_i}$ is a membership matrix that maps the $p_i$ objects in mode $i$ into $r_i$ latent clusters such that $({\mbox{\boldmath $ M$}}_i)_{j,a} = \mathbb{I}\{ j \text{-th fiber in mode-}i \text{ belongs to cluster } a\} $, and $r_i \ll p_i$. By construction, every row of ${\mbox{\boldmath $ M$}}_i$ contains exactly one nonzero entry, indicating the unique cluster assignment of each entity in mode $i$. In the characteristics tensor ${\cal X}$, the temporal mode (mode $d+1$) typically lacks cluster structure, though our framework can accommodate such extensions. We further define ${\mbox{\boldmath $ S$}}_Y = {\mbox{\boldmath $ B$}} {\mbox{\boldmath $ F$}}$ to represent the centroids matrix of ${\mbox{\boldmath $ Y$}}$. In this formula, our model not only nests the common factor structure but is also sufficiently general to accommodate other cases where ${\mbox{\boldmath $ Y$}}$ is not strictly driven by factors but still admits a latent group representation.

In the existing literature patton2022risk, han2022exact, the group structure of panel units in ${\mbox{\boldmath $ Y$}}$ is typically derived through one of two approaches: either by employing group factor models or by clustering based on the characteristics tensor $\mathcal{X}$. Our framework integrates these two approaches into a unified model (ref). By leveraging the group dependency between panel outcomes ${\mbox{\boldmath $ Y$}}$ and characteristics ${\cal X}$, we enable bidirectional information sharing that significantly improves group estimation along the first tensor mode. Furthermore, this proposed coupled dependency framework provides more interpretable results in real applications compared to methods that analyze panel outcomes ${\mbox{\boldmath $ Y$}}$ or characteristics ${\cal X}$ in isolation.

For applications in empirical asset pricing, the main objectives are twofold: (i) to recover accurate estimates of the membership matrices ${\mbox{\boldmath $ M$}}_1,{\mbox{\boldmath $ M$}}_2,\cdots, {\mbox{\boldmath $ M$}}_d$, thereby identifying the latent grouping structures in both $\mathcal{X}$ and ${\mbox{\boldmath $ Y$}}$; (ii) to leverage these estimated groupings for improved estimate on factor loadings ${\mbox{\boldmath $ B$}}$.

Methodology

We estimate membership matrices ${\mbox{\boldmath $ M$}}_i,i=1,...,d$ (clustering) and factor loadings ${\mbox{\boldmath $ B$}}$ via a two-stage strategy. The proposed clustering procedure includes two steps: an initialization step using the Panel Coupled Matrix-Tensor Spectral Clustering (PMTSC) procedure, presented in Algorithm (ref), and an iterative refinement step using the Panel Coupled Matrix-Tensor Lloyd algorithm (PMTLloyd), presented in Algorithm (ref). Following the estimation of group memberships, we employ PCA or the least squares method to estimate the group-level factor loading matrices.

Under PMTC model (ref), we aim to jointly estimate the membership matrices ${\mbox{\boldmath $ M$}}_1, {\mbox{\boldmath $ M$}}_2, \ldots, {\mbox{\boldmath $ M$}}_d$, by solving the following optimization problem:

While problem (ref) is nonconvex and computationally intractable when optimizing all parameters simultaneously, the objective function is convex in each individual parameter when the others are held fixed. This multi-convex structure naturally lends itself to the use of an efficient alternating optimization algorithm, combined with a warm initialization procedure.

A key advantage of problem (ref) is its ability to exploit the shared clustering structure in the first mode by integrating information from both ${\cal X}$ and ${\mbox{\boldmath $ Y$}}$. Specifically, for the first mode, we construct the augmented matrix ${\mbox{\boldmath $ Z$}}_1 = [\hbox{\rm mat}_1({\cal X}), {\mbox{\boldmath $ Y$}}] = {\mbox{\boldmath $ M$}}_1 [\hbox{\rm mat}_1({\cal S})(\otimes_{j=2}^d {\mbox{\boldmath $ M$}}_j \otimes {\mbox{\boldmath $ I$}}_T), {\mbox{\boldmath $ S$}}_Y]$. This formulation demonstrates that ${\mbox{\boldmath $ M$}}_1$ can be consistently estimated through the joint analysis of ${\cal X}$ and ${\mbox{\boldmath $ Y$}}$, rather than relying on either source independently. For modes $i \neq 1$, we define $ {\mbox{\boldmath $ Z$}}_i = \hbox{\rm mat}_i({\cal X}) = {\mbox{\boldmath $ M$}}_i [\hbox{\rm mat}_i({\cal S})(\otimes_{j=1,j\neq i}^d {\mbox{\boldmath $ M$}}_j \otimes {\mbox{\boldmath $ I$}}_T)] $. Although estimation of ${\mbox{\boldmath $ M$}}_i$ for $i \neq 1$ depends on ${\cal X}$, the coupling ${\cal X}$ and ${\mbox{\boldmath $ Y$}}$ improves the accuracy of ${\mbox{\boldmath $ M$}}_1$. This improvement propagates to the other tensor modes, resulting in more reliable cluster recovery for the entire model. This propagation effect underscores a primary advantage of the coupled framework: enhanced recovery in the shared mode strengthens recovery across all remaining modes in practice.

Panel coupled matrix-tensor spectral clustering algorithm

We first develop a warm initialization procedure that generalizes high-order spectral clustering with specific coupled-mode innovations. Similar to tensor block models, the core tensor ${\cal S}$ in model (ref) may exhibit degenerate ranks, meaning the rank of the unfolded matrix $\hbox{\rm mat}_i({\cal S})$ is strictly less than the dimension $r_{i}$. This characteristic allows us to reformulate model (ref) as an equivalent low-rank Tucker decomposition,

where ${\mbox{\boldmath $ U$}}_i\in\mathbb{R}^{p_i\times m_i}$ are orthogonal matrices with $m_i\le r_i$. Building on the High Order Orthogonal Iteration (HOOI) framework zhang2018tensor, we introduce the Panel Coupled High Order Orthogonal Iteration (PCHOOI) algorithm, which jointly captures the low rank structures present in both ${\cal X}$ and ${\mbox{\boldmath $ Y$}}$. For the coupled first mode, PCHOOI employs a distinctive approach. During initialization, we estimate the leading $m_1$ left singular matrix as $\widehat{\mbox{\boldmath $ U$}}_1^{(0)}={\rm LSVD}_{m_1} ([\hbox{\rm mat}_1({\cal X}),{\mbox{\boldmath $ Y$}}])$. In subsequent iterations, we refine this estimate using the coupled matrix, $$\widehat{\mbox{\boldmath $ U$}}_1^{(k)} = {\rm LSVD}_{m_1}([\hbox{\rm mat}_1(\mathcal{X}\times_2 \widehat{\mbox{\boldmath $ U$}}_2^{(k-1)\top} \times_3 \cdots \times_d \widehat{\mbox{\boldmath $ U$}}_d^{(k-1)\top} ), {\mbox{\boldmath $ Y$}}]).$$ For the remaining uncoupled modes, the procedure follows the standard HOOI approach. The complete PCHOOI algorithm is detailed in Algorithm (ref).

Following the paradigm of classical spectral clustering, which performs low-rank projection before applying $k$-means zhang2024leave, our initialization procedure combines PCHOOI with a $k$-means clustering step. This integrated approach, termed Panel Coupled Matrix-Tensor Spectral Clustering (PMTSC), operates in two stages:

{\bf Stage 1}: We apply PCHOOI (Algorithm (ref)) to estimate the orthogonal matrices ${\mbox{\boldmath $ U$}}_i$, obtaining estimates $\widehat{\mbox{\boldmath $ U$}}_i$ that span the principal subspaces of ${\cal X}$ and ${\mbox{\boldmath $ Y$}}$.

{\bf Stage 2}: We employ a modified high-order spectral clustering algorithm (Algorithm (ref)) using these $\widehat{\mbox{\boldmath $ U$}}_i$ estimates to recover the latent group structures $\widehat{\mbox{\boldmath $ M$}}_i^{(0)}$. Analogous to PCHOOI, for the coupled first mode, we construct the augmented matrix, $$\widehat{\mbox{\boldmath $ Z$}}_1 = \widehat{\mbox{\boldmath $ U$}}_1\widehat{\mbox{\boldmath $ U$}}_1^{\top}[ \hbox{\rm mat}_1(\mathcal{X} \times_2 \widehat{\mbox{\boldmath $ U$}}_2^{\top} \times_3 \cdots \times_d \widehat{\mbox{\boldmath $ U$}}_d^{\top}), {\mbox{\boldmath $ Y$}}]. $$ Given the computational complexity of exact $k$-means, we implement a relaxed version that can be efficiently solved using approximation algorithms such as $k$-means++ with relaxation factor $\kappa=O(\log \bar r)$. The pseudo-code is presented in Algorithm (ref).

Our algorithm requires the ranks $r_k$ as inputs, which are permitted to scale with tensor dimensions in our theoretical framework. While our simulations assume known true ranks for simplicity, practical applications should employ data-driven selection criteria such as BIC or leverage domain-specific prior knowledge.

Panel coupled matrix-tensor Lloyd algorithm

Following warm initialization via PMTSC (Algorithm (ref)), we employ the Panel Coupled Matrix-Tensor Lloyd algorithm (PMTLloyd, Algorithm (ref)) to refine the membership matrices ${\mbox{\boldmath $ M$}}_i$ and estimate the core tensor ${\cal S}$ and centroid matrix ${\mbox{\boldmath $ S$}}_Y$. PMTLloyd extends the iterative projection framework developed for Tucker and CP factor models han2024tensor, han2024cp to our PMTC model (ref). The key insight motivating our approach involves orthogonal projections. Define ${\mbox{\boldmath $ P$}}_i = {\mbox{\boldmath $ M$}}_i({\mbox{\boldmath $ M$}}_i^{\top} {\mbox{\boldmath $ M$}}_i)^{-1}$ and consider

Since ${\mbox{\boldmath $ M$}}_i^\top{\mbox{\boldmath $ P$}}_i=I_{r_i}$, model (ref) yields

where $\widetilde{\cal X}_{i}$ has dimensions $r_1 \times \cdots \times r_{i-1} \times p_i \times r_{i+1} \times \cdots \times r_d$. While this dimension reduction from $p_j$ to $r_j$ (where $r_j\ll p_j$) improves efficiency, the non-orthogonality of ${\mbox{\boldmath $ P$}}_j,j\neq i$ creates heterogeneous noise in $\widetilde{\cal E}_{i}$. To address this issue, we introduce an orthogonal projection approach. Let ${\mbox{\boldmath $ W$}}_i$ contain the normalized columns of ${\mbox{\boldmath $ M$}}_i$, and ${\mbox{\boldmath $ \mathnormal\Lambda$}}_i^2={\mbox{\boldmath $ M$}}_i^\top {\mbox{\boldmath $ M$}}_i$. Define

Then model (ref) implies that

Unlike (ref), this formulation preserves signal strength while applying homogeneous noise reduction. Given the true core tensor ${\cal S}\times_{\ell\neq i}^d {\mbox{\boldmath $ \mathnormal\Lambda$}}_\ell$ and centroids matrix ${\mbox{\boldmath $ S$}}_Y$, we estimate ${\mbox{\boldmath $ M$}}_i$ via nearest neighbor assignment:

The membership matrix is then reconstructed as $({\mbox{\boldmath $ M$}}_i)_{j,a} = 1$ if and only if $(g_i)_j = a$, for all $i=1,...,d$. The operation in (ref) achieves two critical objectives: it dramatically reduces the dimensionality by projecting onto all modes except the $i$-th, and it effectively averages out noise. Under proper conditions on the combined noise tensor ${\cal E}_{i}^*$, estimation of the membership matrix ${\mbox{\boldmath $ M$}}_i$ based on ${\cal X}_i,{\mbox{\boldmath $ Y$}}$ can be made significantly more accurate, as the statistical error rate now depends on $p_i\prod_{j\neq i}^d r_j$ rather than $p_1p_2\ldots p_d$.

In practice, we do not know ${\cal S}$, ${\mbox{\boldmath $ S$}}_Y$ and ${\mbox{\boldmath $ W$}}_{i}$, $1\le i\le d$. Similar to back-fitting algorithms, we iteratively estimate the membership matrix ${\mbox{\boldmath $ M$}}_{i}$ at iteration $k$ based on

using estimates $\widehat{\mbox{\boldmath $ W$}}_{j}^{(k-1)},~ j\neq i$, from the previous iteration. And the centers $\widehat{\cal S}_i^{(k)},\widehat {\mbox{\boldmath $ S$}}_Y^{(k)}$ are estimated through block-wise averaging and projections $\widehat{\cal S}_i^{(k)}={\cal X} \times_{i} \widehat{\mbox{\boldmath $ P$}}_i^{(k-1)\top}\times_{j\neq i}^d \widehat{\mbox{\boldmath $ W$}}_j^{(k-1)\top}$ and $\widehat{\mbox{\boldmath $ S$}}_Y^{(k)}=\widehat{\mbox{\boldmath $ P$}}_1^{(k-1)\top}{\mbox{\boldmath $ Y$}}$. As we shall show in the next section, such an iterative procedure leads to a much improved statistical rate in the high dimensional panel coupled matrix-tensor clustering scenarios, as if all ${\mbox{\boldmath $ W$}}_{i}$, $1\le i\le d$, and ${\cal S},{\mbox{\boldmath $ S$}}_Y$ are known and we indeed observe ${\cal X}_{i}$ following model (ref).

Estimation of factor loading matrix

Let $\widehat {\mbox{\boldmath $ M$}}_i$ denote the final membership matrix estimates from the PMTLloyd algorithm, and define $\widehat{\mbox{\boldmath $ W$}}_i=\widehat{\mbox{\boldmath $ M$}}_i(\widehat{\mbox{\boldmath $ M$}}_i^\top \widehat{\mbox{\boldmath $ M$}}_i)^{-1}$. For latent factors, we estimate the factor loading matrix and latent factors as

where $\widehat {\mbox{\boldmath $ U$}}_B$ represents an estimate of the left singular matrix of the loading matrix ${\mbox{\boldmath $ B$}}$. As in conventional factor models, this estimate is subject to rotational ambiguity. When factors are observable, we directly estimate the factor loading matrix via least squares method:

Theoretical Analysis

We now establish the statistical consistency and error rates for the proposed estimators and factor loading matrix ${\mbox{\boldmath $ B$}}$ under proper regularity conditions. The membership matrix ${\mbox{\boldmath $ M$}}_i$ encodes cluster assignments through the relationship $({\mbox{\boldmath $ M$}}_i)_{j,a} = 1$ if and only if the cluster label $(g_i)_j = a$, for $i=1,...,d$. Define the SVD ${\mbox{\boldmath $ M$}}_i={\mbox{\boldmath $ W$}}_i{\mbox{\boldmath $ \mathnormal\Lambda$}}_i{\mbox{\boldmath $ Q$}}_i^\top$, where ${\mbox{\boldmath $ W$}}_i$ contains the normalized columns of ${\mbox{\boldmath $ M$}}_i$, and define the rescaled core tensor as

Because rearranging the cluster labels leaves the clustering outcome unchanged, the cluster label vector $g_i\in\mathbb{R}^{p_i}$ for mode-$i$ can only be determined up to a label permutation. Starting with an initial labeling $g_i^{(0)}\in\mathbb{R}^{p_i}$, we denote by $\pi_i^{(0)} : [r_i] \to [r_i]$ the best permutation that minimizes the discrepancies between $g_i^{(0)}$ and $g_i$, specifically:

where $(\pi \circ g_i)_j := \pi((g_i)_j)$ and $\Pi_{r_i}$ is the collection of all permutations on $[r_i]$. Let $k$ denote the iteration step in the PMTLloyd algorithm. We define $h_i^{(k)}$ as the mode-$i$ misclustering error rate (CER) at iteration $k$,

To complement the clustering error rate, we introduce the following misclustering loss,

We also impose a non-degeneracy condition on the distance between block centers (defined by the core tensor and center matrix) to ensure the identifiability of clustering.

for $i=1,...,d$. Then $\Delta_{1}^2 \ge \Delta_{1,x}^2+\Delta_y^2$. Specifically, we define $\Delta_i^2 = \infty$ when $r_i = 1$. Define $\Delta_{\min}=\min_{1\le j\le d} \Delta_j$. Let $p_* = \prod_{i=1}^{d} p_i$, $\bar{p} = \max\{p_1, ...p_d\}$, $\underline{p} = \min\{p_1, ...p_d\} $ and $ p_{-i} = p_*/p_i$. Analogous notation applies for ranks, i.e., $r_*$, $\bar{r}$, $r_{-i}$, and $m_*$, $\bar{m}$, $m_{-i}$.

To present theoretical properties of the proposed procedures, we impose the following assumptions.

Assumption (ref) parallels noise conditions commonly imposed in the clustering literature gao2022iterative, loffler2021optimality, han2022exact. While we assume independent sub-Gaussian entries for technical convenience and to ensure fast statistical error rates, this condition could theoretically be relaxed to accommodate sub-Gaussian noise with mode-wise additive covariance or Gaussian noise with general cross-sectional dependence. However, such generalizations would substantially complicate the mathematical formulations, statistical results, and technical requirements of our framework. Our choice thus strikes a balance between analytical tractability and the preservation of the core insights of our study.

For analytical tractability, Assumption (ref) imposes a “balanced cluster" condition, which ensures that no single group becomes too sparse to allow for consistent recovery of its centroid gao2022iterative, loffler2021optimality, han2022exact.

Assumption (ref) imposes a standard mixing condition that accommodates a broad class of time series models, including causal ARMA processes with continuously distributed innovations fan2008nonlinear, tsay2018nonlinear. This assumption requires the tail probabilities of $f_{t}$ to decay exponentially; notably, when $\gamma_2=2$, the process $f_{t}$ becomes sub-Gaussian.

We begin by analyzing the PCHOOI estimators in Algorithm (ref), which form the foundation of the PMTSC algorithm (Algorithm (ref)). Theorem (ref) establishes performance bounds for the orthogonal loading matrices. This result extends the HOOI framework of zhang2018tensor to coupled matrix-tensor analysis and is of independent interest.

Notably, identifiable cores ${\cal S},{\mbox{\boldmath $ S$}}_Y$ in ${\cal X},{\mbox{\boldmath $ Y$}}$ may exhibit degenerate ranks with $m_i<r_i$. Theorem (ref) accommodates such degeneracy. While the Frobenius norm error bounds for $\widehat{\cal X}$ and $\widehat{\mbox{\boldmath $ Y$}}$ show no improvement over zhang2018tensor through coupled matrix-tensor analysis, the bounds for orthogonal loading matrices $\widehat{\mbox{\boldmath $ U$}}_i$ reveal important distinctions. For the shared mode 1, the error depends on a weighted average of the spectral norms of both the error tensor and error matrix, scaled by the minimum singular value of the coupled signal matrix. In contrast, for non-shared modes ($i\ge 2$), the error depends solely on the error tensor and the minimum singular values of the signal tensor, consistent with zhang2018tensor.

This coupled matrix-tensor analysis provides substantial benefits. In the extreme case where ${\mbox{\boldmath $ Y$}}$ is noiseless ($\sigma_y=0$), the shared component estimation error for mode 1 becomes significantly smaller than in uncoupled analysis for ${\cal X}$, since $\lambda_{m_1}([\hbox{\rm mat}_1({\cal X}),{\mbox{\boldmath $ Y$}}])> \lambda_{m_1}(\hbox{\rm mat}_1({\cal X}))$. Even with noisy ${\mbox{\boldmath $ Y$}}$, improvement persists: since $\lambda_{m_1}([\hbox{\rm mat}_1({\cal X}),{\mbox{\boldmath $ Y$}}])> \lambda_{m_1}(\hbox{\rm mat}_1({\cal X})) + \lambda_{m_1}({\mbox{\boldmath $ Y$}})$ typically holds, the statistical error for mode 1 improves compared to uncoupled analysis, provided $\sigma_y\asymp \sigma_x$.

We now establish theoretical guarantees for our clustering algorithms, PMTSC and PMTLloyd. We first present convergence rates for the PMTSC algorithm (Algorithm (ref)), which serves as initialization for the subsequent PMTLloyd algorithm.

Theorem (ref) provides a starting point for our further theoretical analysis. It establishes rough upper bounds for both the misclustering error $h_i^{(0)}$ and the loss measure $\ell_i^{(0)}$ ($i=1,...,d$), providing the foundation for our subsequent analysis of the PMTLloyd algorithm. While these polynomial rates in Theorem (ref) could be sharpened in the context of the spectral clustering literature, we focus instead on the error bounds of the iterative algorithm, where PMTLloyd substantially improves these initial estimates.

Next, we examine the statistical performance of the iterative PMTLloyd algorithm (Algorithm (ref)) following PMTSC initialization. The dimension reduction operation in (ref) projects ${\cal X}$ in other modes of the tensor from $\mathbb{R}^{p_j}$ to $\mathbb{R}^{r_j}$ for all modes $j\neq i$, preserving cluster centers and signal strength while reducing noise. The following theorem establishes the conditions under which ideal rates, based on population projections, are achieved.

The SNR requirements in (ref) reveal key differences between coupled and uncoupled modes. For the coupled mode 1, the condition incorporates signals from both the characteristics tensor ${\cal X}$ and panel matrix ${\mbox{\boldmath $ Y$}}$. The first term mirrors the requirement for uncoupled modes, while the second resembles conditions from (sub-)Gaussian mixture model clustering gao2022iterative, zhang2024leave. In contrast, uncoupled modes ($i\ge 2$) depend solely on ${\cal X}$. This signal aggregation provides substantial practical advantages. By Remark (ref), when $\sigma_x\asymp \sigma_y$ and $p_i\gg r_i^2$ (as typically holds), the coupled approach achieves weaker SNR requirements than clustering on ${\mbox{\boldmath $ Y$}}$ alone. Similarly, when $\Delta_y$ is large, the requirements are weaker than clustering on ${\cal X}$ alone. Thus, coupling never strengthens SNR requirements and often relaxes them considerably.

Combining Theorems (ref) and (ref) yields the following exact clustering results.

Beyond recovering cluster memberships in the PMTC model (ref), another important task is estimating the factor loading matrix ${\mbox{\boldmath $ B$}}$. We establish guarantees for the estimated loading matrix under both observable and unobservable factor scenarios.

In the latent factors case, the factor loading matrix exhibits rotational ambiguity, requiring that estimation accuracy be measured by the distance between column subspaces, as is standard in the latent factor literature bai2003,lam2012,han2024tensor. Remarkably, consistency of $\widehat{\mbox{\boldmath $ U$}}_B$ holds even with finite or slowly growing $T$, provided weak factor strength $\lambda_B\asymp \sigma_y$ and dimensionality requirement $p_1\gg r_1^2$. A stronger factor strength further relaxes the dimensionality requirement for finite $T$. Similarly, for observed factors, consistency requires only slowly growing $T$ when $p_1\gg m_1r_1$. Intuitively, the grouping structure effectively provides $p_1/r_1$ repeated samples per cluster pattern, reducing noise and weakening requirements on $T$. This is a unique advantage of grouped panel data.

Simulation Studies

We evaluate the clustering recovery and loading matrix estimation accuracy of our methods via simulations under two data-generating processes. By varying parameters such as SNR, we benchmark our coupled approach against methods using only the tensor ${\cal X}$ or the matrix ${\mbox{\boldmath $ Y$}}$. These experiments demonstrate how joint estimation enhances clustering assignments, factor loading estimation, and reconstruction accuracy.

Simulation for PCHOOI Algorithm

In this subsection, we evaluate the proposed PCHOOI algorithm under model (ref) with $d=2$. The orthogonal loading matrices ${\mbox{\boldmath $ U$}}_i$ ($i=1,2$) are obtained as the left singular matrix of $p_i \times m_i$ matrices with i.i.d.\ $\mathcal{N}(0,1)$ entries. The noise tensor $\mathcal{E} \in \mathbb{R}^{p_1 \times p_2 \times T}$ and error matrix ${\mbox{\boldmath $ \mathnormal\eta$}} \in \mathbb{R}^{p_1 \times T}$ have i.i.d.\ entries from $\mathcal{N}(0, \sigma_x^2)$ and $\mathcal{N}(0, \sigma_y^2)$, respectively. We set $p_1 = p_2 = 50$, $T = 40$, $\sigma_x = \sigma_y = 1$, $m_1 = m_2 = 5$, $\lambda_{\min}(\mathcal{S}) = c_x \sqrt{p_1 + m_*T}$, and $\lambda_{\min}({\mbox{\boldmath $ S$}}_Y) = c_y \sqrt{p_1 + T}$. Two scenarios are considered: (i) $\log c_x = 1$ with $\log c_y$ varying from 1 to 3; and (ii) $\log c_y = 2$ with $\log c_x$ varying from 0 to 2. We compare Algorithm (ref) against SVD on ${\mbox{\boldmath $ Y$}}$ and HOOI on ${\cal X}$, reporting average performance over 100 replications. Performance is measured by the subspace distance $\ell_2({\mbox{\boldmath $ U$}}_i, \widehat{\mbox{\boldmath $ U$}}_i) = \|\widehat{\mbox{\boldmath $ U$}}_i \widehat{\mbox{\boldmath $ U$}}_i^\top - {\mbox{\boldmath $ U$}}_i {\mbox{\boldmath $ U$}}_i^\top\|_2$, $i=1,2$.

By leveraging information from both ${\mbox{\boldmath $ Y$}}$ and ${\cal X}$, PCHOOI delivers more accurate loading space estimates, particularly for the shared first mode. For the coupled mode (${\mbox{\boldmath $ U$}}_1$), PCHOOI achieves accuracy exceeding that of the better-performing baseline, with the largest gains occurring when SVD and HOOI perform similarly. In contrast, for the uncoupled mode (${\mbox{\boldmath $ U$}}_2$), improvements are minimal, with $\ell_2({\mbox{\boldmath $ U$}}_2, \widehat{\mbox{\boldmath $ U$}}_2)$ decreasing by less than 0.1%. This results from the theoretical structure of the PCHOOI estimator; coupling directly increases the singular values for the shared mode but only provides secondary benefits to uncoupled modes.

Simulation for Clustering and Factor Loading Estimation

In this subsection, we evaluate the proposed clustering algorithms under model (ref) with $d=2$. The noise tensor $\mathcal{E}$ has independent entries from a zero-mean sub-Gaussian distribution with sub-Gaussian norm $\sigma_x$, and the error matrix $\eta$ has independent entries from a zero-mean sub-Gaussian distribution with sub-Gaussian norm $\sigma_y$. Group assignments are balanced across clusters, where each entity has an equal probability of being assigned to any cluster.

For the core tensor, factor loading matrix, and factors, we set $\mathcal{S}_{j_1j_2t} \sim \mathcal{N}(0, \sigma_s^2)$ with $\mathcal{S} \in \mathbb{R}^{r_1 \times r_2 \times T}$, ${\mbox{\boldmath $ b$}}_{i } \sim \mathcal{N}(\mu_B, \sigma_B^2 {\mbox{\boldmath $ I$}}_{m_1})$ for $i=1,...,r_1$ where ${\mbox{\boldmath $ b$}}_i$ is the $j$-th row of ${\mbox{\boldmath $ B$}}$, and $f_t \sim \mathcal{N}(\mu_f, \sigma_f^2 {\mbox{\boldmath $ I$}}_{m_1})$ with $f_t \in \mathbb{R}^{m_1 }$, where ${\mbox{\boldmath $ I$}}_{m_1}$ denotes the $m_1 \times m_1$ identity matrix. Throughout, we set $r_1 = r_2 = m_1 = 5$, $p_1 = p_2 = 200$, $T = 120$ and $\sigma_x = \sigma_y = \sigma_s = \sigma_B = \sigma_f = 1$. We specify $\mu_B = (1, 1, 1, 0, 0)$ and $\mu_f = 0.03 \cdot \mathbf{1}_5$, where $\mathbf{1}_5$ denotes the 5-dimensional vector of ones. After generating $\mathcal{S}$ and ${\mbox{\boldmath $ B$}}$, we normalize their magnitudes to satisfy the following SNR constraints motivated by (ref):

where $\bar{p} = \max\{p_1, p_2\}$, $p_* = p_1 p_2$, and $\Delta_{x}^2 = \min_{i} \Delta_{i,x}^2 $ with $\Delta_{i,x}$ and $\Delta_{y}$ defined in (ref). We consider two scenarios: (i) fix $\gamma_x = -0.50$ and vary $\gamma_y$ from $-0.45$ to $0.10$; (ii) fix $\gamma_y = -0.10$ and vary $\gamma_x$ from $-0.7$ to $-0.3$.

We compare our coupled algorithms, PMTSC (“X+Y: PMTSC”) and PMTLloyd refinement with PMTSC initialization (“X+Y: PMTSC+PMTLloyd”), against several benchmarks that utilize ${\cal X}$ and/or ${\mbox{\boldmath $ Y$}}$. For convenience, HLloyd denotes the High-order Lloyd algorithm representing the projection idea in (ref), applicable to coupled ${\cal X}, {\mbox{\boldmath $ Y$}}$, or ${\cal X}$ alone; similarly, PMTLloyd represents the orthogonal projection in (ref) and can also be applied solely to ${\cal X}$. The benchmarks are: spectral clustering on ${\mbox{\boldmath $ Y$}}$ zhang2024leave, HLloyd refinement with High-order Spectral Clustering (HSC) initialization on ${\cal X}$ han2022exact, PMTLloyd refinement with HSC initialization on ${\cal X}$ (“X: HSC+PMTLloyd”), and HLloyd refinement with PMTSC initialization on coupled ${\cal X}, {\mbox{\boldmath $ Y$}}$ (“X+Y: PMTSC+HLloyd”). As expected, increasing $\gamma_x$ or $\gamma_y$, which corresponds to higher ${\rm SNR}_x$ and ${\rm SNR}_y$, reduces the Clustering Error Rate (CER) for all methods.

The left two panels of Figure (ref) show clustering accuracy for the first mode (coupled mode). Methods based on coupled ${\cal X}, {\mbox{\boldmath $ Y$}}$ achieve uniformly lower CER than methods using ${\cal X}$ alone, and also outperform methods using only ${\mbox{\boldmath $ Y$}}$, except for PMTSC; PMTSC+PMTLloyd is consistently the best performer. In the second panel with varying $\gamma_x$, on coupled ${\cal X},{\mbox{\boldmath $ Y$}}$, PMTSC+HLloyd improves over PMTSC, while PMTSC+PMTLloyd further improves over PMTSC+HLloyd, demonstrating both the benefit of orthogonal projection and the advantage of iterative refinement. The right two panels show clustering accuracy for the second mode (uncoupled mode). Surprisingly, methods based on coupled ${\cal X}, {\mbox{\boldmath $ Y$}}$ still achieve uniformly lower CER than methods using ${\cal X}$ alone. This likely stems from improved first-mode clustering, which enhances the projection updates in (ref) and propagates benefits throughout the algorithm. In the first, third, and fourth panels, all coupled methods coincide. Although PMTLloyd and HLloyd applied solely to ${\cal X}$ are indistinguishable in these figures, Appendix (ref) demonstrates the advantages of PMTLloyd over HLloyd in tensor co-clustering, particularly with imbalanced clusters. These simulation results demonstrate the consistent superiority of coupled clustering methods and the benefits of PMTLloyd refinement.

We next evaluate the accuracy of estimated factor loadings. For observable factors, we compute $\sqrt{\sum_{i=1}^{p_1}\| \widehat{\mbox{\boldmath $ b$}}_i - {\mbox{\boldmath $ b$}}_i \|_{2}^2}$ as the estimation error, applying a time-series demeaning step to ${\mbox{\boldmath $ Y$}}$ before estimation. For latent factors, we measure accuracy using the subspace distance $\| \widehat {\mbox{\boldmath $ U$}}_B \widehat {\mbox{\boldmath $ U$}}_B^{\top} - {\mbox{\boldmath $ U$}}_B {\mbox{\boldmath $ U$}}_B^{\top} \|_2$, where $\widehat{\mbox{\boldmath $ U$}}_B = \text{LSVD}_{m_1}(\widehat{\mbox{\boldmath $ M$}}_1 \widehat{\mbox{\boldmath $ B$}})$ and ${\mbox{\boldmath $ U$}}_B = \text{LSVD}_{m_1}({\mbox{\boldmath $ M$}}_1 {\mbox{\boldmath $ B$}})$.

Figure (ref) displays the factor loading matrix estimation errors. The superiority of coupled methods over those using ${\cal X}$ or ${\mbox{\boldmath $ Y$}}$ alone is reflected in the clustering results, where PMTSC+PMTLloyd consistently outperforms PMTSC+HLloyd, remaining the best performer. For observable factors, when $\gamma_y$ is very small (low SNR in ${\mbox{\boldmath $ Y$}}$), all clustering algorithms outperform no-clustering estimation, even with high CER. However, as $\gamma_y$ increases, ${\cal X}$-based methods deteriorate due to clustering errors, while coupled methods avoid this degradation. Although spectral clustering on ${\mbox{\boldmath $ Y$}}$ eventually improves as its CER approaches zero, our coupled methods converge faster and maintain the best performance throughout.

Appendix (ref) presents additional SNR settings under balanced clusters, Appendix (ref) extends to imbalanced clusters, and Appendix (ref) addresses smaller dimensions. Across all settings, we observe the same phenomena. Overall, coupling $\mathcal{X}$ with ${\mbox{\boldmath $ Y$}}$ yields significantly lower clustering errors and superior factor loading matrix estimation.

Applications to Empirical Asset Pricing

\paragraph{Data.} We apply our methodology to U.S. equity portfolios constructed via the Panel Tree (P-tree) approach \citep*{cong2025growing}. This flexible method allows for portfolio construction based on multiple asset characteristics, capturing nonlinear and interactive effects. The dataset comprises monthly returns and 61 characteristics for 400 portfolios spanning from January 1990 to December 2024. From this data, we construct a characteristics tensor $\mathcal{X}$ and a return matrix ${\mbox{\boldmath $ Y$}}$. To ensure comparability across firms, all characteristics within $\mathcal{X}$ are rank-normalized on a cross-sectional basis.

For modeling observable factors, we employ the Fama-French five-factor model, given its established relevance in explaining cross-sectional variations in stock returns. The five factors are the excess market return (Mkt-RF), size (SMB), value (HML), profitability (RMW), and investment (CMA). This allows us to benchmark the performance of our proposed methodology against a widely accepted framework in empirical asset pricing.

\paragraph{Empirical Design.} We set the number of clusters in the first mode to $r_1 \in \{2, 5, 10, 25\}$, while the second mode is fixed at $r_2 = 6$. This choice aligns with the widely recognized clustering of stock characteristics into six distinct themes: momentum, value, investment, profitability, frictions related to size, and intangibles.

To assess model performance, we use the total $R^2$, a standard metric for evaluating cross-sectional model fit feng2024deep. This measure captures the proportion of return variation explained by the factor model relative to a market benchmark: \[ \text{total } R^2 = 1 - \frac{\sum_{i=1}^{p_1} \sum_{t=1}^T \left( Y_{i,t} - \sum_{k=1}^{r_1} \hat{{\mbox{\boldmath $ b$}}}_{k}^\top f_t \cdot \mathbf{1}\{i \in \mathcal{G}_{1k}\} \right)^2}{\sum_{i=1}^{p_1} \sum_{t=1}^T \left( Y_{i,t} - R_t^{\text{mkt-rf}} \right)^2}, \] where $Y_{i,t}$ represents observed returns for asset $i$ at time $t$, $\hat{{\mbox{\boldmath $ b$}}}_{k}$ are estimated factor loadings for group $k$, $f_t$ denotes factor realizations, and $R_t^{\text{mkt-rf}}$ represents excess market returns. The denominator benchmarks the variation explained by the market factor. A positive $R^2$ signifies the factor model's capacity to explain additional variation.

We benchmark PMTC model (ref) against a range of alternatives, including univariate sorts (e.g., book-to-market ratio “BM" and market equity value “ME"), $5 \times 5$ bivariate sorts, and economically motivated specifications with pre-defined clusters on the second mode ($\mathcal{G}_2^F$). These comparisons evaluate PMTC relative to return-based baselines, widely used sorting methods, and economically grounded groupings.

For data splitting, we employ two evaluation schemes to assess our methodology. The first approach uses an in-sample (INS) and out-of-sample (OOS) split, where the dataset is divided into a training sample covering the first 35 years (January 1980-December 2014) and an OOS validation period spanning the subsequent 10 years (January 2015-December 2024). The training sample is used to estimate the latent clustering structures ($\widehat{\mbox{\boldmath $ M$}}_1$ and $\widehat{\mbox{\boldmath $ M$}}_2$) and the factor loading matrix ($\widehat{\mbox{\boldmath $ B$}}$). Predictive performance is then evaluated on the OOS period. Under this design, the algorithm inputs are $\mathcal{X}_{\text{train}} \in \mathbb{R}^{400 \times 61 \times 420}$, ${\mbox{\boldmath $ Y$}}_{\text{train}} \in \mathbb{R}^{400 \times 420}$, and ${\mbox{\boldmath $ F$}}_{\text{train}} \in \mathbb{R}^{5 \times 420}$. An alternative scheme is provided in Appendix (ref).

\paragraph{Empirical Performance Comparison.} PMTLloyd consistently outperforms competing methodologies in both in-sample (INS) and out-of-sample (OOS) evaluations. At $r_1 = 10$, a standard benchmark in empirical asset pricing Fama1992crosssection, it enhances the OOS total $R^2$ by 2.9 percentage points compared to the “Only ${\mbox{\boldmath $ Y$}}$" baseline under a simple in-sample and out-of-sample split.

Relative to traditional sorting methods, PMTLloyd demonstrates considerable improvements. At $r_1 = 10$, the OOS total $R^2$ is 10 percentage points higher than BM and 5.8 percentage points higher than ME. This highlights its capacity to outperform standard approaches across both static and dynamic settings. When compared to the $\mathcal{G}_2^F$ specification, which relies on predefined characteristic clusters, PMTLloyd generally achieves comparable or superior performance, further validating its adaptability and effectiveness in diverse settings.

Using a standard INS-OOS split, we classify the $r_1 = 10$ clusters by ranking them in descending order based on within-group average market equity (ME) and profit margin (PM). This methodology underscores significant cross-group heterogeneity. Clusters characterized by high ME and PM exhibit market betas close to 1 and positive RMW loadings, suggesting that larger firms not only tend to co-move with the market but also demonstrate stronger profitability.

The left panel illustrates a monotonic increase in SMB (size) exposure as market equity declines, while Groups 7 and 8 -- characterized by the highest SMB betas -- simultaneously exhibit the most negative RMW (profitability) loadings. The right panel shows a consistent decline in RMW betas as PM decreases, aligning with the definition of the profitability factor. This trend highlights the systematic reduction in RMW beta with lower PM values. Moreover, clusters with similar market, SMB, and RMW betas often exhibit substantial differences in their HML and CMA loadings, highlighting variation in value and investment characteristics across groups. The systematic alignment between latent clusters and fundamental characteristics, such as market equity and profit margins, indicates that PMTC effectively captures the low-rank structure of risk premia, leading to the observed gains in $R^2$.

Final Discussion

While this study focuses on a panel matrix-tensor setting with no time-mode clustering, our framework naturally extends to general cases that relax this restriction. Coupling a matrix and a tensor, even when they share only a single-mode grouping structure, can still be highly informative: the shared mode stabilizes estimation and, in practice, improves recovery of latent clusters in other modes. This flexibility of our approach suggests that the benefits of joint modeling extend beyond the panel setup considered here. More broadly, settings with multiple modes sharing a common grouping structure could further enhance the precision of structure recovery.

Another natural statistical extension of our framework arises when the group structures in $\mathcal{X}$ and ${\mbox{\boldmath $ Y$}}$ are not perfectly aligned. In practice, economic characteristics may cluster firms slightly differently from return dynamics, yet the information in $\mathcal{X}$ can still serve as a powerful auxiliary signal to guide clustering in ${\mbox{\boldmath $ Y$}}$. This setting resembles transfer learning: one could apply a debiasing step to explicitly correct the partial misalignment between $\mathcal{X}$-based and ${\mbox{\boldmath $ Y$}}$-based clusters, while still exploiting the shared latent structure. Such an approach would broaden applicability to heterogeneous but related datasets.

The empirical results reveal substantial heterogeneity in factor exposures across groups. Some groups primarily load on size and value factors, while others exhibit strong exposure to profitability or investment. This heterogeneity suggests that enforcing a uniform factor structure across firms may obscure significant cross-sectional variation. The PMTC framework can be extended to accommodate high-dimensional “factor zoos" by incorporating sparse estimation techniques to select group-specific relevant factors. Sparse estimation techniques such as LASSO can isolate the most relevant factors within each cluster, aligning with the literature on uncommon factors and asset heterogeneity cong2023sparse. Recognizing that clusters may be driven by distinct factors, the coupled matrix-tensor framework provides a natural setting for factor selection.

{\onehalfspacing }