The IMR package implements the Incomplete Matrix
Regression method of Fouda et al. [@fouda2026imr]. It provides a framework for
both matrix completion and regression on response matrices with missing
values. Let \(\boldsymbol{Y} \in \Re^{n\times
m}\) denote the observed incomplete matrix, where missing values
are designated by either NA or zero. The estimator for each
entry of the matrix is specified by any combination of the following
\[\boldsymbol{\hat Y_{ij}} = \boldsymbol{\hat\beta_{oi}} + \boldsymbol{\hat\Gamma_{oj}} + \boldsymbol{X_i\hat\beta_j} + \boldsymbol{\hat\Gamma_i Z^{'}_j} + \boldsymbol{\hat M_{ij}}\]
where \(\boldsymbol{X}\in \Re^{n\times p}\) and \(\boldsymbol{Z}\in \Re^{m\times q}\) are row (\(p\) predictors) and column (\(q\) predictors) covariate matrices, respectively. The vectors \(\boldsymbol{\hat\beta_{o}}\in \Re^{n}\) and \(\boldsymbol{\hat\Gamma_{o}}\in \Re^{m}\) represent the row-level and column-level intercepts. The term \(\boldsymbol{\hat\beta}\) denotes the row covariate coefficients, which may be structured as either an \(p\times m\) matrix (one coefficient of each (predictor,column) pair, or a \(p\)-dimensional vector (one coefficient for each predictor), forcing coefficients to be equal across all columns for each covariate. Similarly, the column covariate coefficients, denoted by \(\boldsymbol{\hat\Gamma}\), can be either an \(n \times q\) matrix or a \(q\)-dimensional vector where all rows share the same coefficient for each covariate. To avoid having too many parameters, we impose Lasso (\(L_1\)) penalties on the row and column covariate coefficients. Finally, \(\boldsymbol{M}\) is a rank-r (r is a hyper-parameter) low-rank matrix subject to a nuclear norm penalty. Together, they yield the following penalty structure:
\[\mathrm{Penalty} = \lambda_\beta {\|\boldsymbol{\beta}\|}_{1} + \lambda_\Gamma {\|\boldsymbol{\Gamma}\|}_{1} + \lambda_m {\|\boldsymbol{S_r^{1/2}MS_c^{1/2}}\|}_{*},\]
where \(\boldsymbol{S_r} \in \Re^{n\times n}\) and \(\boldsymbol{S_c} \in \Re^{m \times m}\) are similarity (or information) matrices that describe the correlation structure among the rows and columns of the response matrix, respectively. In the absence of a known correlation structure, this penalty term reduces to the standard nuclear norm, \({\|\boldsymbol{M}\|}_{*}\). We have 4 penalty parameters: \((\lambda_\beta, \lambda_\Gamma, \lambda_m, r)\). We provide a method to estimate those parameters.
As we said above, we can use any combination of the model components to define our estimator. Examples of these combinations include:
A primary advantage of the algorithm is that it can deal with very large matrices, as long as they are sparse. With the exception of the similarity matrices, the algorithm neither stores nor performs operations on any matrices with dimensions \(\min{(n,m)}\times \min{(n,m)}\) or larger. The matrix \(\boldsymbol{Y}\) is stored as a sparse object, which replaces the dense matrix format with a vector of the observed values and their corresponding coordinate locations. All other matrices and associated operations are restricted to dimensions bounded by \(\max\{m\times r,\, n \times r,\, n\times p,\, m\times q\}\).
To illustrate the standard workflow within the package, assume the matrices \(\boldsymbol{Y}\) and \(\boldsymbol{X}\) are defined as above, and the objective is to fit example 1 above to obtain the complete estimated matrix \(\boldsymbol{\hat Y}\). Then, the implementation is as follows:
# load the library
require(IMR)
# set the hyperparameter value.
lambda_beta <- 0.02
# load the data example (see ?IMR::Bixi_sample for more information)
Bixi <- IMR::Bixi_sample
# create the data object
data <- imr_data(Y = Bixi$Y, X = Bixi$X)
# update the model structure to fit example 1
data <- update(data, row_covariates = TRUE, # turn XBeta on (on by default when X is provided)
shared_beta = TRUE, # make beta a p-dimensional vector (off by default)
low_rank_component = FALSE, # turn M off (on by default)
row_intercept = TRUE) # turn row intercepts on (off by default).
# fit the model
fit <- imr_fit(data, lambda_beta = lambda_beta )
# obtain \hat{Y}
Y_hat <- reconstruct(fit, data)$estimates
#> Constructing XBeta ...
#> Constructing row intercepts ...
#> done.This vignette is structured as follows: we first discuss the data structure and processing steps for the proposed model, using an illustrative example. Second, we show how to define the model structures (i.e., which components to include in the estimator). Third, we show an example of the fit function, followed by a presentation of the results and performance assessment. Next, we describe our procedure for selecting the values of the penalty parameters. We conclude with a discussion on the construction and integration of the similarity matrices, \(\boldsymbol{S_r}\) and \(\boldsymbol{S_c}\).
As an example, we analyze data from Bixi 1, a docked bike-sharing service in Montreal, Canada. We use the data compiled by [@lei2025], which contains normalized daily departure counts for each of 579 stations over 196 days (April 15 to October 27, 2019). The data is accompanied by side information. We select two temporal variables (temperature and precipitation) and one spatial variable (park area). For simplicity, we only consider a sample of 150 stations and 100 days, taken at random. Below is a summary of the data (Note that the data list also contains similarity matrices which we discuss in the final section).
Bixi <- IMR::Bixi_sample
head(Bixi$X,2)
#> x_mean_temp_c x_total_precip_mm
#> [1,] 0.05785124 0.09492635
#> [2,] 0.23140496 0.00000000
head(Bixi$Z,3)
#> [,1]
#> [1,] 0.2629746
#> [2,] 0.0000000
#> [3,] 0.3097933
dim(Bixi$Y)
#> [1] 100 150
Bixi$Y[1:6,11:12]
#> 6423 - Hôpital général juif (de la Côte Ste-Catherine / Légaré)
#> [1,] NA
#> [2,] NA
#> [3,] NA
#> [4,] 0.1881188
#> [5,] NA
#> [6,] 0.1485149
#> 6752 - Hutchison / Beaubien
#> [1,] NA
#> [2,] NA
#> [3,] NA
#> [4,] NA
#> [5,] NA
#> [6,] NA
sprintf("Percentage of observed entries in Y: %.1f%%", 100*sum(Bixi$Y != 0)/length(Bixi$Y))
#> [1] "Percentage of observed entries in Y: NA%"
sprintf("Percentage of entries in the test set: %.1f%%", 100* sum(Bixi$test != 0)/length(Bixi$test))
#> [1] "Percentage of entries in the test set: 14.4%"To facilitate data integration within the proposed framework, we
introduce the imr_data class, which automatically
structures the inputs into a model-compliant format. As in the preceding
output, unobserved entries in the response matrix are denoted by \(0\), though NA is also
supported. Because we treat all zero-values entries as missing data,
structural zeros cannot be accommodated as observed values.
data <- imr_data(Y = Bixi$Y, X = Bixi$X, Z = Bixi$Z, val_prop = 0.2, seed = 2026)
#> 'as(<dgeMatrix>, "dgCMatrix")' is deprecated.
#> Use 'as(., "CsparseMatrix")' instead.
#> See help("Deprecated") and help("Matrix-deprecated").
print(data)
#>
#> == IMR Data Object ==
#> Target Matrix (Y) : 100 x 150 (15000 cells)
#> Observed : 9573 (63.82%)
#> Missing (sparsity) : 5427 (36.18%)
#>
#> -- Data Split --
#> Training : 7659 (80.0% of observed)
#> Validation : 1914 (20.0% of observed)
#>
#> -- Model Configuration --
#> Low-Rank Matrix (M) : [active]
#> Row Intercepts : [inactive]
#> Column Intercepts : [inactive]
#> Row Covariates (X) : [active] 2 vars, per-column
#> Column Covariates (Z) : [active] 1 var, per-column
#> Row Similarity : [not provided]
#> Column Similarity : [not provided]
#> ==========================The response argument Y in the function above accepts
any standard matrix class from base R or the Matrix
package. The framework coerces this input into a dgCMatrix
object using the function as_incomplete(). Run
class?CsparseMatrix 2 for more information about the matrix
arguments.
Bixi$test <- as_incomplete(Bixi$test)
Bixi$test[3:5,3:4]
#> 3 x 2 sparse Matrix of class "dgCMatrix"
#> 6210 - Métro Sauvé (Berri / Sauvé) 6186 - St-Hubert / Laurier
#> [1,] . .
#> [2,] . .
#> [3,] . .Furthermore, the function performs a QR decomposition on both
covariate matrices, \(\boldsymbol{X}\)
and \(\boldsymbol{Z}\), storing the
resultant components as XQ, XR,
ZQ and ZR. This transformation is necessary as
the fit function works on orthonormal matrices (see the main article for
more details on the methodology). However, in our final coefficient
estimates, the estimates are projected back onto the original covariate
spaces.
The latter portion of the preceding data summary shows the model
configuration. The first column contains the supplied auxiliary
information, indicating, for instance, the presence of two row
covariates and one column covariate. The second column specifies the
active components in the estimator. By default, the low-rank matrix and
any provided side information are included in the active model
specification. The designation (Unshared) indicates that
the covariate coefficients are permitted to vary across rows or columns;
that is, the coefficients take the form of matrices rather than vectors.
This default specification can be modified via the update()
function. To illustrate, we introduce row intercepts and constrain the
row covariate coefficients to be shared across all columns, which
enforces a fixed effect of temperature and precipitation across all
stations.
data <- update(data, shared_beta = TRUE, row_intercept = TRUE)
print(data)
#>
#> == IMR Data Object ==
#> Target Matrix (Y) : 100 x 150 (15000 cells)
#> Observed : 9573 (63.82%)
#> Missing (sparsity) : 5427 (36.18%)
#>
#> -- Data Split --
#> Training : 7659 (80.0% of observed)
#> Validation : 1914 (20.0% of observed)
#>
#> -- Model Configuration --
#> Low-Rank Matrix (M) : [active]
#> Row Intercepts : [active]
#> Column Intercepts : [inactive]
#> Row Covariates (X) : [active] 2 vars, shared
#> Column Covariates (Z) : [active] 1 var, per-column
#> Row Similarity : [not provided]
#> Column Similarity : [not provided]
#> ==========================The update function modifies the model architecture
through logical arguments that activate or deactivate specific
components. We proceed to fit this model.
The fit function requires fixed values for the penalty parameters: the three lambda terms and the rank of the low-rank component. For now, let’s set arbitrary values for these parameters. Tuning procedure is detailed in a subsequent section.
The print function yields diagnostic information
regarding the estimation procedure, while the summary
function provides an overview of the resultant estimates. The estimated
coefficients can be extracted using the coef function, or
directly by accessing the fit$coefficients object. It
should be noted that the summary function show the
covariate coefficients with respect to the original covariate matrices
\(\boldsymbol{X}\) and \(\boldsymbol{Z}\), while the raw extracted
coefficients in coef() correspond to the orthonormal
matrices XQ and ZQ.
print(fit)
#>
#> ====================================================
#> === Incomplete Matrix Regression (IMR) Fit ===
#> ====================================================
#> Formula : Y ~ 𝛃₀ + X𝛃 (shared) + 𝚪Z + M
#> Status : Converged (in 194 iterations)
#> Target : 100 x 150 matrix (36.18% missing)
#>
#> -- Dimensions --
#> 𝛃 (row covariates) : 2 x 1 matrix
#> 𝚪 (col covariates) : 100 x 1 matrix
#> 𝛃₀ (row intercepts) : length 100 vector
#> M (latent factors) : U(100 x 10), D(length 10), V(150 x 10)
#>
#> -- Fit (in-sample, on observed entries) --
#> RMSE : 0.0554
#> Pseudo R2 : 0.9496
#>
#> -- Penalties & Hyperparameters --
#> Rank (r) : 10
#> Lambda M : 0.1000
#> Lambda Beta : 0.0000
#> Lambda Gamma : 0.0200
#> Row Similarity : None
#> Column Similarity : None
#> ====================================================
summary(fit)
#>
#> =====================================================
#> === Summary of Incomplete Matrix Regression (IMR) ===
#> =====================================================
#>
#> Goodness of Fit (in-sample, on observed entries)
#> -------------------------------------------------------------
#> RMSE : 0.0554
#> MSE : 0.0031
#> Pseudo R2 (1 - RSS/SST) : 95.0%
#>
#> Variance Decomposition (share of explained variance)
#> -------------------------------------------------------------
#> Latent Matrix (M) : 85.8%
#> Row Covariates (X𝛃) : 59.3%
#> Col Covariates (𝚪Z) : 0.1%
#> Row Intercepts (𝛃₀) : 0.2%
#> Overlap / non-additive : -45.4%
#> (Overlap is due to regularization.)
#> -------------------------------------------------------------
#>
#> Model Estimates
#> -------------------------------------------------------------
#>
#> -- Row Covariates --
#> Mode: Shared across columns (p = 2)
#> Estimate
#> x_mean_temp_c 0.7188
#> x_total_precip_mm -0.0424
#>
#> -- Column Covariates --
#> Mode: Row-specific (100 x 1 matrix)
#> Summary of effects across 100 rows
#> Mean SD Min Max Sparsity L2.Norm
#> 1 -0.0046 0.0527 -0.1699 0.1721 40.00% 0.5265
#>
#> -- Intercepts --
#> Row Intercepts (n=100) | Mean: 0.0609 | SD: 0.1430 | Min: -0.3201 | Max: 0.3159
#>
#> -- Latent Component --
#> Rank (r): 10
#> Singular Values:
#> [1] 28.0966 7.6497 5.6739 3.6104 3.2412 2.9792 2.7055 2.5040 2.3693
#> [10] 2.1823
names(coef(fit))
#> [1] "u" "d" "v" "beta" "gamma" "beta0" "gamma0"Since the fit function runs an iterative estimation
algorithm, we need to specify convergence parameter. There are two
parameters to set: the maximum allowed iterations and the tolerance
threshold. We suggest setting the threshold to \(1e-5\) or less, and set the number of
iteration to be large enough to reach this threshold. These parameters
are specified using the imr_convergence object.
convergence <- imr_convergence(maxit = 15, thresh = 1e-5, trace = TRUE); print(convergence)
#>
#> == IMR Convergence Parameters ==
#> Max Iterations: 15
#> Threshold: 1e-05
#> Initialization: Least Squares
#> Trace Progress: Enabled
#> ================================
fit2 <- imr_fit(data, rank = 10, lambda_m = 1e-1, lambda_beta = 0, lambda_gamma = 0.02,
convergence = convergence)
#> 1 obj= 0.02705 ratio= 0
#> 1 obj= 0.0035 ratio= 0.1810626
#> 2 obj= 0.00299 ratio= 0.03820343
#> 3 obj= 0.00277 ratio= 0.02050178
#> 4 obj= 0.00265 ratio= 0.01296901
#> 5 obj= 0.00258 ratio= 0.01024488
#> 6 obj= 0.00252 ratio= 0.008465095
#> 7 obj= 0.00249 ratio= 0.006684274
#> 8 obj= 0.00246 ratio= 0.006121865
#> 9 obj= 0.00243 ratio= 0.005670869
#> 10 obj= 0.00241 ratio= 0.004994333
#> 11 obj= 0.00239 ratio= 0.00477189
#> 12 obj= 0.00238 ratio= 0.004683454
#> 13 obj= 0.00237 ratio= 0.00424778
#> 14 obj= 0.00235 ratio= 0.0041631
#> 15 obj= 0.00234 ratio= 0.004076078
#> Did not converge in 15 iterations.In the output above, obj denotes the value of the
objective function at the current iteration and ratio
represents the relative change in the parameter estimates between
successive iterations, evaluated using the Frobenius norm. The
optimization algorithm terminates when either the maximal iteration
limit is attained or the Frobenius difference between subsequent
estimates falls below the specified tolerance threshold.
We begin by showing the two imputation methods we have, beginning with the full reconstruction of the response matrix.
data_out <- reconstruct(fit, data)
#> Constructing M ...
#> Constructing XBeta ...
#> Constructing GammaZ ...
#> Constructing row intercepts ...
#> done.
print(names(data_out))
#> [1] "beta" "gamma" "M" "beta0" "gamma0" "xbeta"
#> [7] "gammaz" "estimates"The reconstruct function outputs all estimated model
components alongside the fully imputed response matrix
($estimates). It should be noted that \(\beta\) and \(\Gamma\) estimates are projected back onto
the original covariate matrices, \(X\)
and \(Z\). For illustration, the
explicit computation for the first element of the response matrix is
cat(sprintf("True value of the entry (1,1) int he data matrix is %.4f and the estimated value is %.4f",
data$Y[1,1],
# row intercept +
data_out$beta0[1] +
# row covariates (X beta) +
Bixi$X[1,] %*% data_out$beta +
# column covariates (gamma Z) +
data_out$gamma[1,] %*% t(Bixi$Z)[,1] +
# low-rank matrix (M)
data_out$M[1,1]))
#> True value of the entry (1,1) int he data matrix is 0.4158 and the estimated value is 0.3763This reconstruction function can be computationally prohibitive in high-dimensional settings. Moreover, we may only be interested in imputing a subset of the response matrix. Consider the test matrix we have. Let’s only impute the entries that corresponds to that non-missing entries in this matrix.
preds <- reconstruct_partial(fit, data, Bixi$test@i, Bixi$test@p,
trace=TRUE, return_matrix = TRUE)
#> Constructing M ...
#> Constructing XBeta ...
#> Constructing GammaZ ...
#> Constructing row intercepts ...
print(preds[1:9,1:9])
#> 9 x 9 sparse Matrix of class "dgCMatrix"
#>
#> [1,] . . 0.2560806 . . . .
#> [2,] . . . 0.1648926 . . .
#> [3,] . . . . . . .
#> [4,] . . . . . 0.07572037 .
#> [5,] . . . . 0.4394715 . .
#> [6,] . . . . . . .
#> [7,] . . . . . 0.16555796 0.3490084
#> [8,] 0.6451307 . 0.3124619 . . . .
#> [9,] 0.7952965 0.5528109 . . . . .
#>
#> [1,] 0.7098218 0.1362061
#> [2,] . 0.1227303
#> [3,] . 0.3374960
#> [4,] . .
#> [5,] . 0.3128002
#> [6,] . 0.1808324
#> [7,] . .
#> [8,] . 0.1955837
#> [9,] . .The reconstruct_partial function takes two additional
arguments, the \(i\) and \(p\) vectors, which denote the row and
column indices of the entries to be estimated. If the
return_matrix argument is TRUE, the function
returns a sparse matrix containing the estimates. If it’s
FALSE, then it returns a vector of the estimates
corresponding directly to the @x slot of the sparse data
structure. We also provide an evaluation function that compute standard
error metrics between two numeric vectors.
| RMSE | Rel_RMSE | MAE | MAPE | Spearman_Rho |
|---|---|---|---|---|
| 0.0978 | 0.2182 | 0.0702 | 30.0065 | 0.9285 |
We now address the tuning of the penalty parameters. Our model has
four parameters: three regularization weights \((\lambda)\) and the rank of the low-rank
component, \(r\). We begin by
constructing a parameter grid. For the lambdas (\(\lambda_\beta,\lambda_\Gamma, \lambda_m\)),
the grid is defined by a minimum, a maximum, and a sequence length. For
the rank, the grid requires a minimum, a maximum, and a step size. Both
the nuclear norm penalty and the rank require an early-stopping
tolerance parameter; a value of 2 or 3 is generally sufficient. If the
maximum for any of the lambdas is unknown, it may be initialized to
NA. We later estimate their upper bounds using the
Karush-Kuhn-Tucker (KKT) conditions.
grid <- imr_tune_grid(
beta = c(0, NA, 40), # (min, max, length)
gamma = c(0, NA, 40), #(min, max, length)
nuclear = c(0, NA, 40, 2), #(min, max, length, early-stopping tolerance)
rank = c(2, 15, 1,2) #(min, max, step, early-stopping tolerance)
)
print(grid)
#>
#> == IMR Hyperparameter Configuration ==
#> Beta: Range: 0 -> auto (Grid: 40 points)
#> Gamma: Range: 0 -> auto (Grid: 40 points)
#> Nuclear: Range: 0 -> auto (Length: 40, Streaks: 2)
#> Rank: Range: 2 -> 15 (Step: 1, Streaks: 2)
#> ===========================================================To fix any parameter, it can be supplied as a scalar. Parameters omitted from the specification will automatically default to standard configurations, as illustrate below:
imr_tune_grid(
beta = 3,
gamma = c(0, NA),
nuclear = c(0),
rank = c(0, 10)
)
#>
#> == IMR Hyperparameter Configuration ==
#> Beta: Fixed at 3
#> Gamma: Range: 0 -> auto (Grid: 20 points)
#> Nuclear: Fixed at 0
#> Rank: Range: 0 -> 10 (Step: 2, Streaks: 2)
#> ===========================================================Returning to the initial grid specification, we need to find upper limits for the \(\lambda\) parameters. This is achieved with a two-step procedure for each unspecified maximum:
1- An initial upper bound is derived using the KKT conditions. 2- Because the KKT bound can occasionally be too big, we perform a subsequent bisection search to find the minimum \(\lambda\) required to shrink all corresponding coefficients exactly to zero.
When determining the maximum for either \(\lambda_\beta\) or \(\lambda_\Gamma\), the penalty for the alternate covariate is set to infinity, while the low-rank penalty parameters are held at default values. Conversely, when determining the upper limit for \(\lambda_m\), the two covariate penalties are fixed at default values.
convergence <- imr_convergence(maxit = 1000, thresh = 1e-5)
grid <- imr_set_grid_limits(data, grid,
default_rank = 2, default_lambda_m = 0,
default_lambda_beta = 0, default_lambda_gamma = 0,
bisection_iter = 10, # number of iteration in the second step
verbose = TRUE,convergence = convergence )print(grid)
#>
#> == IMR Hyperparameter Configuration ==
#> Beta: Range: 0 -> 0.005547 (Grid: 40 points)
#> Gamma: Range: 0 -> 0.303631 (Grid: 40 points)
#> Nuclear: Range: 0 -> 15.985542 (Length: 40, Streaks: 2)
#> Rank: Range: 2 -> 15 (Step: 1, Streaks: 2)
#> ===========================================================With the grid fully specified, we proceed to find the optimal penalty parameters. The tuning function internally runs one of the three estimation processes:
One Covariate: if only one of \(\lambda_\beta\) and \(\lambda_\Gamma\) requires tuning, the algorithm runs in parallel across the respective covariate grid. For each candidate value, it executes the first scheme above to determine the optimal low-rank parameters.
Alternating Optimization: If both \(\lambda_\beta\) and \(\lambda_\Gamma\) require tuning, an alternating minimization procedure is employed. The algorithm fixes \(\lambda_\Gamma\) and tunes \(\lambda_\beta\) (along with \(\lambda_m\) and \(r\)) using scheme 2 above; it then fixes \(\lambda_\beta\) and tunes \(\lambda_\Gamma\), again using scheme 2. It iterates until the difference in estimated coefficients is smaller than some threshold.
Given that all parameters in our example require tuning, the algorithm uses the third scheme. We use the fast mode, which generally yields good performance. For an exhaustive final tuning, however, the slow model is recommended, as it evaluates the complete the parameter space of \(\lambda_m\) and \(r\).
cv_out <- imr_tune(data, grid, fast_nuclear = TRUE, convergence = convergence,
n_cores= 1, seed = 2026, verbose=1)summary(cv_out$fit)
#>
#> =====================================================
#> === Summary of Incomplete Matrix Regression (IMR) ===
#> =====================================================
#>
#> Goodness of Fit (in-sample, on observed entries)
#> -------------------------------------------------------------
#> RMSE : 0.0709
#> MSE : 0.0050
#> Pseudo R2 (1 - RSS/SST) : 91.8%
#>
#> Variance Decomposition (share of explained variance)
#> -------------------------------------------------------------
#> Latent Matrix (M) : 78.5%
#> Row Covariates (X𝛃) : 61.4%
#> Col Covariates (𝚪Z) : 0.4%
#> Row Intercepts (𝛃₀) : 0.2%
#> Overlap / non-additive : -40.6%
#> (Overlap is due to regularization.)
#> -------------------------------------------------------------
#>
#> Model Estimates
#> -------------------------------------------------------------
#>
#> -- Row Covariates --
#> Mode: Shared across columns (p = 2)
#> Estimate
#> x_mean_temp_c 0.7191
#> x_total_precip_mm -0.0450
#>
#> -- Column Covariates --
#> Mode: Row-specific (100 x 1 matrix)
#> Summary of effects across 100 rows
#> Mean SD Min Max Sparsity L2.Norm
#> 1 -0.0224 0.1174 -0.3779 0.321 18.00% 1.1889
#>
#> -- Intercepts --
#> Row Intercepts (n=100) | Mean: 0.0613 | SD: 0.1427 | Min: -0.3171 | Max: 0.3191
#>
#> -- Latent Component --
#> Rank (r): 5
#> Singular Values:
#> [1] 27.1736 6.4034 3.8979 2.7565 1.7213
print(cv_out$fit)
#>
#> ====================================================
#> === Incomplete Matrix Regression (IMR) Fit ===
#> ====================================================
#> Formula : Y ~ 𝛃₀ + X𝛃 (shared) + 𝚪Z + M
#> Status : Converged (in 24 iterations)
#> Target : 100 x 150 matrix (36.18% missing)
#>
#> -- Dimensions --
#> 𝛃 (row covariates) : 2 x 1 matrix
#> 𝚪 (col covariates) : 100 x 1 matrix
#> 𝛃₀ (row intercepts) : length 100 vector
#> M (latent factors) : U(100 x 5), D(length 5), V(150 x 5)
#>
#> -- Fit (in-sample, on observed entries) --
#> RMSE : 0.0709
#> Pseudo R2 : 0.9176
#>
#> -- Penalties & Hyperparameters --
#> Rank (r) : 5
#> Lambda M : 0.5322
#> Lambda Beta : 0.0000
#> Lambda Gamma : 0.0156
#> Row Similarity : None
#> Column Similarity : None
#> ====================================================preds <- reconstruct_partial(cv_out$fit, data, Bixi$test@i, Bixi$test@p, return_matrix = FALSE)
knitr::kable(evaluate(preds, Bixi$test@x),format = "simple",digits = 4)| RMSE | Rel_RMSE | MAE | MAPE | Spearman_Rho |
|---|---|---|---|---|
| 0.0869 | 0.1938 | 0.0635 | 27.7056 | 0.9433 |
We now turn to the final component of the proposed methodology: the incorporation of similarity matrices. In many applications, the response matrix exhibits correlation structures across its rows, its columns, or both. These structures can be explicitly described with similarity or covariance matrices. The framework adjusts the estimates of the low-rank component for these correlation structures.
Returning to the Bixi example, the data are supplemented by two
distance matrices: a spatial distance matrix derived from the geographic
coordinates of the stations (columns) and a temporal distance matrix
representing the elapsed days (rows). For the purpose of this
illustration, we impose a Matern 5/2 covariance structure on the columns
and a Gaussian (radial basis function) kernel with a length-scale of 6
on the rows. These kernels are defined using the
imr_similarity function.
similarity_cols <- imr_similarity(x = "matern", d = Bixi$spatial_distance,
invert=TRUE, jitter = .1,
matern_smoothness = 5/2, matern_range = .018); print(similarity_cols)
#>
#> == IMR Similarity Decomposition ==
#> Source: Matern Kernel (Normalized) (Inverted)
#> Dimensions: 150 x 150
#> Jitter value: 0.1
#> Parameters: smoothness=2.5, range=0.018
#> Condition Number: 11.00
#> Top 5 Eigenvalues: 10, 10, 10, 10, 10, ...
#> ==================================
similarity_rows <- imr_similarity(x = "rbf", d = Bixi$temporal_distance,
invert=TRUE, jitter = .1,
rbf_ell = 6); print(similarity_rows)
#>
#> == IMR Similarity Decomposition ==
#> Source: RBF Kernel (Normalized) (Inverted)
#> Dimensions: 100 x 100
#> Jitter value: 0.1
#> Parameters: ell=6
#> Condition Number: 11.00
#> Top 5 Eigenvalues: 10, 10, 10, 10, 10, ...
#> ==================================Key notes regarding the imr_similarity function:
1- The primary argument, x, accepts either a character
string (“Matern” or “RBF”) or a matrix (covariance or information). If a
string is provided, a corresponding distance matrix d must
be supplied. If a matrix is supplied, then d is
omitted.
2- The framework uses an information (precision) matrix.
Consequently, if a covariance kernel is specified or provided, the
argument invert=TRUE must be declared to compute the
inverse.
3- Setting jitter to a small positive value helps
reducing the condition number. The supplied value will be added to the
diagonal elements of the covariance matrix before inverting it.
4- The functions returns the components U and
d, corresponding to the singular vectors and singular
values of the target similarity matrix, respectively. This
representation is the required input for the framework.
We proceed to integrate these matrices into the data object and update the model specification
data <- imr_data(Bixi$Y, Bixi$X, Bixi$Z,
similarity_rows = similarity_rows,
similarity_cols = similarity_cols,
val_prop = 0.2, seed = 2026)
data <- update(data, shared_beta = TRUE, row_intercept = TRUE)
print(data)
#>
#> == IMR Data Object ==
#> Target Matrix (Y) : 100 x 150 (15000 cells)
#> Observed : 9573 (63.82%)
#> Missing (sparsity) : 5427 (36.18%)
#>
#> -- Data Split --
#> Training : 7659 (80.0% of observed)
#> Validation : 1914 (20.0% of observed)
#>
#> -- Model Configuration --
#> Low-Rank Matrix (M) : [active]
#> Row Intercepts : [active]
#> Column Intercepts : [inactive]
#> Row Covariates (X) : [active] 2 vars, shared
#> Column Covariates (Z) : [active] 1 var, per-column
#> Row Similarity : [active]
#> Column Similarity : [active]
#> ==========================We rerun the parameter tuning. However, this time, we only tune \(\lambda_m\) and we fix the values of \(r\), \(\lambda_\beta\), and \(\lambda_\Gamma\) to their previously determined optimal values.
grid <- imr_tune_grid(
beta = cv_out$params$lambda_beta,
gamma = cv_out$params$lambda_gamma,
nuclear = c(0, NA, 40, 2),
rank = cv_out$params$rank_in
)
print(grid)
#>
#> == IMR Hyperparameter Configuration ==
#> Beta: Fixed at 0
#> Gamma: Fixed at 0.0155708384166746
#> Nuclear: Range: 0 -> auto (Length: 40, Streaks: 2)
#> Rank: Fixed at 5
#> ===========================================================grid <- imr_set_grid_limits(data, grid, verbose = TRUE,
default_rank = cv_out$params$rank_in,
default_lambda_beta = cv_out$params$lambda_beta,
default_lambda_gamma = cv_out$params$lambda_gamma,
convergence = convergence)
cv_out <- imr_tune(data, grid,
n_cores= 4, seed = 2026, verbose=1,
nuclear_log_scale = TRUE,
convergence = convergence)summary(cv_out$fit)
#>
#> =====================================================
#> === Summary of Incomplete Matrix Regression (IMR) ===
#> =====================================================
#>
#> Goodness of Fit (in-sample, on observed entries)
#> -------------------------------------------------------------
#> RMSE : 0.0695
#> MSE : 0.0048
#> Pseudo R2 (1 - RSS/SST) : 92.1%
#>
#> Variance Decomposition (share of explained variance)
#> -------------------------------------------------------------
#> Latent Matrix (M) : 88.0%
#> Row Covariates (X𝛃) : 5.4%
#> Col Covariates (𝚪Z) : 0.4%
#> Row Intercepts (𝛃₀) : < 0.1%
#> Overlap / non-additive : 6.1%
#> (Overlap is due to regularization.)
#> -------------------------------------------------------------
#>
#> Model Estimates
#> -------------------------------------------------------------
#>
#> -- Row Covariates --
#> Mode: Shared across columns (p = 2)
#> Estimate
#> x_mean_temp_c 0.2078
#> x_total_precip_mm -0.0802
#>
#> -- Column Covariates --
#> Mode: Row-specific (100 x 1 matrix)
#> Summary of effects across 100 rows
#> Mean SD Min Max Sparsity L2.Norm
#> 1 0.0242 0.1235 -0.3246 0.4156 17.00% 1.2525
#>
#> -- Intercepts --
#> Row Intercepts (n=100) | Mean: 0.0121 | SD: 0.0519 | Min: -0.1048 | Max: 0.1215
#>
#> -- Latent Component --
#> Rank (r): 5
#> Singular Values:
#> [1] 48.7679 6.5452 4.2564 3.5545 1.9540
print(cv_out$fit)
#>
#> ====================================================
#> === Incomplete Matrix Regression (IMR) Fit ===
#> ====================================================
#> Formula : Y ~ 𝛃₀ + X𝛃 (shared) + 𝚪Z + M
#> Status : Converged (in 27 iterations)
#> Target : 100 x 150 matrix (36.18% missing)
#>
#> -- Dimensions --
#> 𝛃 (row covariates) : 2 x 1 matrix
#> 𝚪 (col covariates) : 100 x 1 matrix
#> 𝛃₀ (row intercepts) : length 100 vector
#> M (latent factors) : U(100 x 5), D(length 5), V(150 x 5)
#>
#> -- Fit (in-sample, on observed entries) --
#> RMSE : 0.0695
#> Pseudo R2 : 0.9206
#>
#> -- Penalties & Hyperparameters --
#> Rank (r) : 5
#> Lambda M : 0.0579
#> Lambda Beta : 0.0000
#> Lambda Gamma : 0.0156
#> Row Similarity : Active
#> Column Similarity : Active
#> ====================================================preds <- reconstruct_partial(cv_out$fit, data, Bixi$test@i, Bixi$test@p)
knitr::kable(evaluate(preds, Bixi$test@x),format = "pipe",digits = 4)| RMSE | Rel_RMSE | MAE | MAPE | Spearman_Rho |
|---|---|---|---|---|
| 0.085 | 0.1895 | 0.0619 | 27.05 | 0.9442 |