SA-MJLR: A sparse multivariate joint latent regression statistical framework for ultra-high dimensional prediction: Monte Carlo evaluation under (p, q)>n settings

Introduction

The rapid expansion of high-throughput technologies has transformed many scientific disciplines by enabling the simultaneous measurement of thousands of variables from relatively small numbers of observations, and its applications in biomedical research, genomics, neuroscience, metabolomics, environmental sciences, and human phenotyping increasingly involve datasets where the number of measured characteristics is significantly larger than the available sample size. Datasets of this category create a challenging statistical environment commonly described as ultra-high dimensional modelling, in which traditional regression assumptions become inadequate because the information contained in the observations is not sufficient for direct estimation of the large number of unknown parameters (Li et al., 2024; Fan et al., 2020) [6,13] . Expressing a general multivariate regression framework can be:

𝑊 ∈ ℝ 𝑛×𝑝 𝑎𝑛𝑑 𝑍 ∈ ℝ 𝑛×𝑞 Represent the multivariate predictor and response matrices respectively,

𝛽 ∈ ℝ 𝑝×𝑞

Represents the coefficient matrix being unknown, and the random errors is denoted by  .

Estimation is reliable when the number of observations is sufficiently larger than the number of variables in conventional regression settings, and modern scientific datasets increasingly exhibit the simultaneous highdimensional condition; however, p>n and q > n, that is, where the predictor and response dimensions are greater than the number of observations. This situation introduces severe statistical difficulties since the coefficient matrix contains pq unknown parameters; and when both predictor and response dimensions equal 1000, for instance, the model contains one million potential predictor-response relationships, while the available observations may be limited to only a few tens or hundreds. Developments in high-dimensional regression recently have shown that these challenges cannot be addressed adequately through conventional multivariate regression since covariance matrices become unstable, parameter estimates become nonidentifiable, and prediction performance deteriorates due to excessive model complexity (Ham et al., 2024; Li et al., 2024) [9,13] . These problems are particularly important in biomedical and biological applications, where thousands of molecular, physiological, or imaging measurements are collected together with multiple correlated outcomes. A major characteristic of such datasets is the presence of strong dependence structures, where predictor variables often represent groups of related biological mechanisms, while response variables frequently reflect interconnected physiological processes. Hence, modelling these systems requires approaches capable of simultaneously exploiting predictor correlation, response correlation, dimensional reduction, and variable selection. It is clear that, several statistical strategies have been proposed to address high-dimensional regression problems, and one major direction is regularized estimation, in which penalty functions are introduced to stabilize parameter estimation, and the recent advances in sparse multivariate regression have focused on estimating large coefficient matrices while accounting for dependence among multiple response variables, and these approaches demonstrate that incorporating response correlation can improve estimation efficiency compared with treating responses independently (Ham et al., 2024) [9] . Regularized coefficient estimation, however, remains challenging when both dimensions increase simultaneously. In the setting where: p>n and q > n, the coefficient matrix becomes extremely large, and direct estimation requires strong assumptions regarding sparsity, covariance structure, or low-rank representation. Recent studies on high-dimensional output regression have emphasized that scalable methods must reduce the complexity of the output space while maintaining predictive information (Li et al., 2024) [13] . Sparse regression provides an attractive solution because many high-dimensional systems are believed to contain only a limited number of relevant predictors, and recent sparse multivariate regression frameworks have incorporated L1based penalties to identify informative variables while controlling model complexity. Sparse regression methods that directly estimate regression coefficients may struggle, however, when predictors are strongly correlated because biologically meaningful variables often occur in groups rather than independently (Ham et al., 2024) [9] . Another important challenge is the modelling of multiple correlated responses, and when the response variables are associated, ignoring their dependence structure can lead to inefficient estimation. Studies recently have therefore investigated methods that jointly estimate regression coefficients and response covariance structures to improve prediction accuracy in highdimensional multivariate regression (Ham et al., 2024) [9] ; nevertheless, covariance-based approaches become computationally demanding when the number of responses is extremely large. An alternative approach is to use the latent variable model, which assumes that high dimensional observations are generated from a lower dimensional latent space, methods that are appealing in particular are the ones which many complex scientific observations have hidden structures that cannot be directly observed. Factor-based regression techniques have been gaining more attention in recent years, as they are able to simultaneously reduce dimensionality and predict, and the recent advances in factor-augmented regression extend the applicability of dimensional reduction to prediction with both common variation and sparse effects present in the response structure (Fan et al., 2023) [7] . The reduced-rank regression has also become an important framework for multivariate prediction because it assumes that the coefficient matrix can be represented using a smaller number of latent dimensions. The sparse reduced-rank regression approaches recently have demonstrated improved predictive performance and interpretability by combining low-rank structure with variable selection, particularly in biological applications involving thousands of molecular predictors and multiple outcomes (Kobak et al., 2021) [12] . Even though latent techniques provide important advantages, existing techniques often focus primarily on either reducing predictor dimensionality or modelling response dependence, and they do not fully address the more difficult situation where both predictor and response spaces are simultaneously ultra-high dimensional. The machine learning approaches have become increasingly important for high-dimensional prediction because of their ability to model nonlinear relationships and complex interactions, and among these methods, Extreme Gradient Boosting (XGBoost) has achieved considerable success due to its computational efficiency, regularization mechanisms, and ability to capture nonlinear patterns in large datasets (Chen & Guestrin, 2016) [3] . Applications of machine learning recently in biomedical and scientific prediction demonstrate that ensemble tree methods can provide strong predictive performance when traditional statistical assumptions are violated, and therefore machine learning models generally focus on prediction rather than uncovering interpretable latent relationships among correlated predictors and responses. As for multivariate problems involving hundreds or thousands of responses, predictive accuracy alone may not be sufficient, but researchers often require knowledge about which variables contribute to prediction and how groups of predictors relate to groups of outcomes. Hence, there remains a need for statistical approaches that combine predictive capability with interpretable latent structures.

Material and Methods

Problem Formulation

This study considers an ultra-high-dimensional multivariate regression problem in which the predictor and response dimensions are both greater than the number of observations. Let:

Represent the multivariate predictor and response matrices respectively, where , n p  n q  Then, the general multivariate regression association is represented as stated in Equation ( 1). In ultra-high dimensional settings, the direct estimation of  becomes infeasible since the number of unknown parameters is: pq Which may greatly be much greater than the available observations, and additionally, if n p  , the predictor covariance matrix:

W W T

is singular, preventing the use of ordinary least squares, and similarly, if n q  , the response covariance structure becomes rank deficient, limiting the effectiveness of conventional multivariate regression approaches. To overcome these challenges, this study develops a Sparse Adaptive Multivariate Joint Latent Regression (SA-MJLR) framework that assumes that the coefficient matrix has two complementary structures, which are: a sparse direct-effect component and a low-rank latent component. Thus, the proposed coefficient structure is mathematically defined as follows:

Where A denotes the latent regression association. Thus, the response model is therefore:

This decomposition enables SA-MJLR to simultaneously identify individual predictor-response associations and shared latent dependencies.

Mechanism for Data Generation

The predictor and response matrices were generated according to the proposed sparse-latent structure for simulation experiments, and the predictor matrix is generated as:

Where  𝑇 ∈ ℝ 𝑛×𝑟 denotes latent scores;  w Q ∈ ℝ 𝑝×𝑟 denotes sparse predictor loadings;  w E denotes the predictor noise.

On the other hand, the response matrix is generated as follows:

Where;  𝐴 ∈ ℝ 𝑟×𝑟 denotes latent regression coefficients;

 s B ∈ ℝ 𝑝×𝑞 denotes sparse direct coefficients;   controls the magnitude of direct sparse effects;

 z E denotes the response noise.

Thus, the true coefficient matrix employed for recovery evaluation is:

Sparse Direct-Effect Estimation

The first component of SA-MJLR estimates sparse direct predictor-response associations, and for each response variable

Where B  controls sparsity,

and the resulting sparse coefficient matrix is:

After the sparse direct component estimation, its contribution is removed, that is:

Where r Z denotes the remaining response variation that is explained by latent structures.

Sparse Joint Latent Loading Estimation

The latent component of SA-MJLR is being estimated from the residual response matrix r Z , whereas the initial predictorresponse association matrix and singular value decomposition applied are respectively given in Equations ( 10) and (11) as:

The initial loading matrices are obtained using Equations ( 12) and ( 13):

Where r denotes the selected latent rank, and to introduce sparsity, SA-MJLR applies alternating soft-thresholding updates, whereby the soft-thresholding operator is defined as:

The predictor loading update and the response loading update are respectively represented in Equations ( 15) and ( 16) as:

After each update, the loading vectors are normalized as:

The alternating procedure will continue for a fixed number of iterations until a stable latent loading estimate is achieved.

Latent Regression Estimation

Latent predictor and response scores are computed, following latent loading estimation, as in Equations ( 18) and ( 19):

Then, the latent regression association becomes:

Using ridge regularization, the latent coefficient matrix is estimated as:

 controls numerical stability and I denotes the identity matrix.

Coefficient Reconstruction and Prediction

The latent coefficient component is reconstructed using:

The final SA-MJLR coefficient matrix is given as:

or written equivalently as

Predictions are obtained for new predictor observations, new W as:

This formulation, however allows SA-MJLR in order to add both sparse direct effects and latent multivariate associations.

Benchmark Regression Methods

The proposed SA-MJLR technique is compared with four existing high-dimensional regression methods.

Ridge Regression

The Ridge estimates coefficients as stated in Equation (26):

in which λ controls shrinkage.

Elastic Net Regression

Elastic Net combines 1 L and 2

Where 5 . 0   is the implementation employed.

Reduced Rank Regression

The Reduced Rank Regression (RRR) estimates a low-rank approximation as:

Where r denotes the retained rank.

Sparse Partial Least Squares

The SPLS estimates latent predictor and response directions through singular value decomposition of Z W T and reconstructs:

Computational Algorithm of SA-MJLR

The SA-MJLR algorithm provides a unified framework for modelling high-dimensional multivariate systems which contain both sparse direct associations and latent dependency structures, and is summarized as follows:

Step eight: Then produce/generate predictions: B W Z  2.9 Monte Carlo Simulation Design 2.9.1 Objective of Simulation Study: The Monte Carlo simulation study is designed to evaluate the finite-sample performance of the proposed Sparse Adaptive Multivariate Joint Latent Regression (SA-MJLR) method under controlled high-dimensional settings, and the objective is to examine whether SA-MJLR can effectively recover regression structures and provide accurate predictions whenever the predictor and response dimensions exceed the available sample size. Thus, the simulation framework is specifically designed to reflect situations where the true coefficient matrix contains two simultaneous components, which are: 1. A sparse direct-effect component representing a small number of important predictor-response relationships, and 2. A low-rank latent component representing hidden shared variation between predictor and response variables.

This structure reflects many modern high-dimensional applications in which observed variables are influenced by both individual effects and underlying latent processes.

Data-Generating Mechanism

In each simulation replicate; the predictor and the response matrices 𝑊 ∈ ℝ 𝑛×𝑝 and 𝑍 ∈ ℝ 𝑛×𝑞 respectively were generated based on the model as follows:

Where W denotes the high-dimensional predictor matrix, s B denotes a sparse direct-effect coefficient matrix, w Q and z Q denote, respectively, sparse latent loading matrices for predictor and response spaces, A denote the latent regression relationship, whereas  denotes random response noise; and the predictor matrix is produced (generated) using:

Where T represents latent scores and w  denotes predictor noise, and the true coefficient matrix used for evaluation is defined as:

Where  controls the strength of the sparse direct contribution.

Simulation Parameters and Scenarios

In the experiment, five simulation scenarios were considered to evaluate model robustness under increasing dimensional complexity, and the scenarios were designed such that both predictor and response dimensions exceeded sample size, representing genuine high-dimensional problems. Thus, the simulation parameters are summarized in Table 1. 1 shows the Monte Carlo simulation settings in which the scenarios progressively increased dimensionality while reducing signal-to-noise ratio and sparsity levels, and SC1 represents a relatively favourable high-dimensional condition, whereas SC5 represents an extreme setting with weak signals and highly sparse underlying relationships.

Competing Methods

The proposed SA-MJLR method is compared against four established multivariate regression approaches, known as Ridge Regression (representing dense shrinkage-based regression), Elastic Net Regression(representing sparse regularized regression), Reduced Rank Regression(representing low-rank multivariate modeling), and Sparse Partial Least Squares (representing sparse latent projection modeling). These methods were selected because they represent major modelling strategies commonly used for high-dimensional regression problems.

Monte Carlo Procedure

For each scenario, the simulation experiment was repeated 10 times using independent randomly generated datasets, and a fixed random seed was applied to ensure reproducibility. Then, for each replication, the following were executed, but the final reported results represent the average performance across all replications:  Synthetic predictor and response matrices were generated.  Each competing method was fitted using 80% of observations for training.  The remaining 20% of observations were used for prediction evaluation.  Estimated coefficient matrices were compared with the known true coefficient matrix.  The prediction and recovery metrics were computed.

Evaluation Criteria

Two categories of evaluation metrics were used, which are prediction performance and Coefficient Recovery Performance. The prediction accuracy was assessed using Root Mean Square Error (RMSE), Mean Absolute Error (MAE), and Coefficient of Determination (R 2 ), whereas coefficient recovery performance, which is the ability to recover the true coefficient structure, was evaluated using Frobenius Error (FE) and Correlation (ρ), and their computational formulas are stated in Equations ( 33) to (37).

Here, a lower value indicates a more accurate coefficient recovery.

Here, a higher value indicates a better structural recovery

Monte Carlo Simulation RStudio Codes for the Five Methods:

The Monte Carlo Simulation for the proposed method, as well as the four competing methods, is committed to a GitHub repository accessible from the link https://github.com/judend/Monte-Carlo-Simulation-RStudio-Codes-for-SA-MJLR-and-Four-Competing-Methods/blob/main/SAVE%20EZE-GOLDEN.docx The result in Table 3 shows the coefficient recovery analysis, which provides stronger evidence of the structural advantage of SA-MJLR; while some methods achieved competitive prediction in simple settings, accurate recovery of the underlying coefficient matrix requires successful identification of both sparse and latent components. From SC2 to SC5, SA-MJLR consistently achieved the lowest Frobenius Error and highest coefficient correlation, and this indicates that the proposed method provides a more accurate reconstruction of the true regression structure. The improvement becomes particularly evident under SC4 and SC5, and in SC5, SA-MJLR achieved a Frobenius Error of 0.0171 and correlation of 0.7617, which outperforms Elastic Net, RRR, Ridge, and SPLS, confirming that separating sparse direct effects from latent shared variation improves parameter estimation in ultra-high-dimensional settings. The results in Table 4 shows the aggregated ranking of the four existing methods against the proposed method (SA-MJLR), which confirms the overall superiority of SA-MJLR, and the proposed method achieved the highest performance score (6.9888), substantially exceeding all competing approaches. The SA-MJLR obtained the best overall predictive accuracy, structural recovery, and coefficient correlation, and the results demonstrate that combining sparse estimation with latent joint modelling provides a more complete representation of complex multivariate relationships than methods relying exclusively on shrinkage or low-rank assumptions. In overall, the simulation study provides strong empirical support that SA-MJLR is a robust framework for high-dimensional multivariate regression problems characterized by simultaneous sparsity and latent dependency structures.

Results

Conclusion

In this study, the Sparse Adaptive Multivariate Joint Latent Regression (SA-MJLR) framework was proposed specifically for high-dimensional multivariate regression problems characterized by the simultaneous presence of sparse direct effects and latent dependency structures. But the motivation for developing SA-MJLR is based on the lacunas or limitations of existing approaches, which often focus exclusively on either sparse variable selection or latent structure extraction, while many real-world high-dimensional systems exhibit both characteristics simultaneously. An extensive Monte Carlo simulation study was conducted to examine the performance of SA-MJLR under progressively challenging high-dimensional scenarios, and the simulation framework considered situations where both predictor and response dimensions were more than the sample size, with increasing dimensionality, decreasing sparsity levels, and reduced signal-to-noise ratios. The comparison was based on the developed technique, as against the four existing techniques known as Ridge Regression, Elastic Net Regression, Reduced Rank Regression (RRR), and Sparse Partial Least Squares (SPLS), and the findings revealed that SA-MJLR provided superior overall performance, particularly under complex high-dimensional conditions. Although SPLS showed competitive performance in the simplest simulation scenario, SA-MJLR consistently achieved better predictive accuracy and coefficient recovery as the dimensional complexity increased, but the proposed method obtained lower RMSE and MAE values, higher R 2 values, lower Frobenius reconstruction errors, and stronger correlations between estimated and true coefficient matrices across most challenging scenarios. Hence, the improved performance of SA-MJLR could be attributed to its hybrid modelling strategy, which separates sparse direct relationships from shared latent structures; unlike conventional regularization methods that assume either dense or purely sparse relationships, or latent projection methods that primarily capture hidden variation, SA-MJLR simultaneously exploits both sources of information, and this enables more accurate estimation of complex predictorresponse relationships while maintaining interpretability through sparse coefficient recovery. Thus, the overall ranking analysis further confirmed the robustness of SA-MJLR, with the proposed method achieving the highest aggregate performance score among all competing approaches, and these results provide evidence that SA-MJLR represents a promising framework for high-dimensional multivariate regression applications where traditional approaches may suffer from information loss, poor prediction, or inaccurate structural recovery.

References

1. Bunea F, She Y, Wegkamp MH. Optimal selection of reduced rank estimators of high-dimensional matrices. The Annals of Statistics. 2011;39(2):1282-1309. DOI:10.1214/10-AOS859.
2. Chen J, Huang JZ. Sparse reduced-rank regression for simultaneous dimension reduction and variable selection. Journal of the American Statistical Association. 2012;107(500):1533-1545. DOI:10.1080/01621459.2012.716613.
3. Chen T, Guestrin C. XGBoost: a scalable tree boosting system. In: Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. Association for Computing Machinery; 2016. p. 785-794. DOI:10.1145/2939672.2939785.
4. Chun H, Keleş S. Sparse partial least squares regression for simultaneous dimension reduction and variable selection. Journal of the Royal Statistical Society: Series B (Statistical Methodology). 2010;72(1):3-25. DOI:10.1111/j.1467-9868.2009.00723.x.
5. Fan J, Han F, Liu H. Challenges of big data analysis. National Science Review. 2014;1(2):293-314. DOI:10.1093/nsr/nwt032.
6. Fan J, Li R, Zhang CH, Zou H. Statistical foundations of data science. Annual Review of Statistics and Its Application. 2020;7:1-29. DOI:10.1146/annurev-statistics-031219-041205.
7. Fan J, Li R, Liao H, Wang Y. Learning high-dimensional factor models with structured sparsity. Journal of the American Statistical Association. 2023;118(542):10021016.
8. Hastie T, Tibshirani R, Friedman J. The elements of statistical learning: data mining, inference, and prediction. 2nd ed. New York: Springer; 2017.
9. Ham Y, Price BS, Rothman AJ. Sparse multivariate regression with structured covariance estimation for highdimensional responses. Journal of Computational and Graphical Statistics. 2024;33(2):451-465.
10. Hoerl AE, Kennard RW. Ridge regression: biased estimation for nonorthogonal problems. Technometrics. 1970;12(1):55-67. DOI:10.1080/00401706.1970.10488634.
11. Izenman AJ. Reduced-rank regression for the multivariate linear model. Journal of Multivariate Analysis. 1975;5(2):248-264. DOI:10.1016/0047-259X(75)90058-4.
12. Kobak D, Lomond R, Sanchez A. Sparse reduced-rank regression for high-dimensional multivariate data. Statistics and Computing. 2021;31:76.
13. Li Y, Li L, Wang H, Zhang X. High-dimensional multivariate regression with large-scale response variables: sparse estimation and prediction. Journal of Multivariate Analysis. 2024;204:105341.
14. Rothman AJ, Levina E, Zhu J. A new approach to Cholesky-based covariance regularization in high dimensions. Biometrika. 2010;97(3):539-550. DOI:10.1093/biomet/asq022.
15. Tibshirani R. Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society: Series B (Methodological). 1996;58(1):267-288.
16. Wainberg M, Alipanahi B, Frey BJ. Are random forests truly the best classifiers? Journal of Machine Learning Research. 2016;17(1):3837-3841.
17. Wold S, Sjöström M, Eriksson L. PLS-regression: a basic tool of chemometrics. Chemometrics and Intelligent Laboratory Systems. 2001;58(2):109-130. DOI:10.1016/S0169-7439(01)00155-1.
18. Yuan M, Zhang CH. On tensor completion via nuclear norm minimization. Foundations of Computational Mathematics. 2012;12:279-297.
19. Zhao P, Yu B. On model selection consistency of Lasso. Journal of Machine Learning Research. 2006;7:25412563.
20. Zou H, Hastie T. Regularization and variable selection via the elastic net. Journal of the Royal Statistical Society: Series B (Statistical Methodology). 2005;67(2):301-320. DOI:10.1111/j.1467-9868.2005.00503.x.