Bivariate Matrix-valued Linear Regression (BMLR): Finite-sample performance under Identifiability and Sparsity Assumptions
Authors: Nayel Bettache
Venue: NeurIPS 2025 poster ()
Keywords: Matrix-valued regression
High-dimensional statistics
Sparse estimation
Multivariate analysis
Structured linear models
Primary area: theory
TL;DR: We study bilinear matrix regression with matrix-valued inputs and outputs, propose optimization-free estimators, derive non-asymptotic error bounds, and validate performance on synthetic and image data.

Abstract: This paper studies a bilinear matrix-valued regression model where both predictors and responses are matrices. For each observation $t$, the response \( Y_t \in \mathbb{R}^{n \times p} \) and predictor \( X_t \in \mathbb{R}^{m \times q} \) satisfy $Y_t = A^* X_t B^* + E_t,$
with \( A^* \in \mathbb{R}_+^{n \times m} \) (row-wise \(\ell_1\)-normalized), \( B^* \in \mathbb{R}^{q \times p} \), and \( E_t \) independent Gaussian noise matrices. The goal is to estimate \( A^* \) and \( B^* \) from the observed pairs \( (X_t, Y_t) \).

We propose explicit, optimization-free estimators and establish non-asymptotic error bounds, including sparse settings. Simulations confirm the theoretical rates and demonstrate strong finite-sample performance. We further illustrate the practical utility of our method through an image denoising application on real data.

Bivariate Matrix-valued Linear Regression (BMLR): Fimite-sample performance under Identifiability and Sparsity Assumptions
Nayel Bettache
Department of Statistics and Data Science Cornell University, Capital Fund Management, 23 rue de l'Université, 75007 Paris, France com
Abstract
This paper studies bilinear matrix-valued regression model where both predic- tors and responses are matrices. For each observation the response YÚ R"xp and predictor X, E R=*x satisfy Y, = X, + E,, with A E RTXm (row- wise (1-normalized), 8 E IR9*P, and E, independent Gaussian noise matrices. The goal is to estimate A and B" from the observed pairs (Xe, Y.). We propose explicit, optimization-free estimators and establish non-asymptotic error bounds, including sparse settings. Simulations confirm the theoretical rates. We illustrate the practical utility of our method through an image denoising application.
Introduction
Supervised learning is core task in modern data analysis, often applied to high-dimensional datasets. With recent advances in data acquisition technologies, many real-world datasets now ex- hibit intrinsic matrix structures. Examples include spatiotemporal measurements, dynamic imaging, and multivariate longitudinal studies, where rows and columns encode distinct dimensions such as time, space, or experimental conditions B (7). In such settings, both the covariates and the re- sponses may be naturally represented as matrices.
Traditional linear regression models are designed to predict a scalar outcome from a vector-valued covariate.. When the response is vector-valued, naive approach consists in fitting separate linear model to each coordinate. However, this strategy ignores the multivariate structure of the response. To address this, several works 66 18 3 low-rank matrix, leading to the multivariate linear regression model. Other works have focused on predicting scalar responses from matrix-values covariates [38 28), leading to the frace regression model where the parameter becomes a matrix, usually assumed to have a low-rank structure. In our setting, both the covariates and responses are matrix-valuee independent trace regression models to each entry of the output matrix. Thus, propose a bilinear model that explicitly couples the row and column structures the input and output matrices.
propose modeling the stacked regression coefficients
It is again not appropriate to fit
In this paper, we extend these lines of work by studying a bivariate matrix-valued linear regres- sion (BMLR) model, where both the predictor and the response are matrix matrix-valued. This framework captures richer structural dependencies and arises naturally in applications such as spatiotemporal forecasting and dynamic network modeling. The BMLR model considers independent observa- tions (X,, Y:)F=,, where each predictor matrix X, E Rmxx and response matrix Y, E Raxx satisfy the relationship
Y X E,, 1 [= (1..,,T)
(1)
The unknown parameter matrices are A E Raxm and B" E Re*r, The non non-negativity con- straint helps ensuring identifiability and interpretability. In many applications, A" represents mixing
39th Conference on Neural Information Processing Systems (NeurIPS 2025).

---

weights, attention coefficients, or interaction strengths that are naturally non non-negative. Note that the model in ① remains unchanged if A is scaled by a positive constant and B is divided by the same constant. To resolve this identifiability issue, we impose the constraint that each row of A has unit fo-rrmm The noise matrices (Ex)²=1 E R=xp have independent and identically distributed entries, each drawn from centered Gaussian distribution with variance o². We dis- cuss more general noise structures in Appendix A The goal is to estimate the parameters A and B from the observed data (X,, Y=) Once estimated, these parameters can be used to pre- dict the out-of-sample response YT-1 given a new covariate matrix XT+1 E Rm=x, This model builds growing literature dedicated to matrix-valued regression with matrix- valuud observations 41 and complements recent advances in matrix autoregression and dynamic 41 50]. These approaches leverage the matrix structure of the data to improve interpretability and predictive performance.
32 23
16
29
7
systems 9
49
24 43 30
Notice that for any í E (T), A captures the link between the rows of Y and the rows of B, while B captures the link between the columns of Y, and the columns of A A This bilinear structure encodes interactions along both matrix dimensions and is inherently richer than what stan- dard multivariate regression models capture, If one ignores the matrix structure and vectorizes both sides, the model becomes vec(Y,) = (B*) @ A" vec(X,) + vec(E,), where vec(-) denotes the column-wist wise vectorization of matrix and @ is the Kronecker product 221. In this formulation, the regression reduces to standard multivariate linear model with T vectorized observations and coefficient matrix M = (B*) @ A", which can be estimated directly using ordinary least squares [. However, this vectorized approach hides the problem structure. Estimating A" and B from an estimate M of M amounts to solving Kronecker factorization problem. It reduces to finding the nearest rank-one matrix in the Kronecker product space [33). This is a non-convex problem that discards the individual matrix roles of A and B². Moreover, the vectorized formulation leads to estimators with high variance when nmpg >> T, as it fails to exploit the low-dimensional structure induced by the bilinear form.
In contrast, our approach preserves and exploits the bilinear structure directly in the matrix domain, allowing for interpretable and efficient estimation of A and B without solving the intractable Kronecker decomposition problem.
1.1 Context
Given the observations (X,, Y.) the linear regression framework models the relationship be- tween the responses and the predictors using an unknown linear map mx + Rnxp, For each t E [T] this relationship is expressed as Y; = ((x)) + E, where F2 represents the noise term. The goal is to learn the unknown function f only using the given observations. In a parametric setting assume that f belongs to predefined class of functions.
A naive approach would be break down the linear relationship by focusing on each entry of the responses individually. This leads to consider for E [m] and j E [p] linear functional fiy : Rmxxy - R. such that
Miliy fx,( [Edaaj t [T).
where [Y:]yy denotes the coefficient on the 7th row and jth column of the matrix Y. Riesz represen- tation theorem ensures the existence of unique matrix Mt₃ E TRmxq such that
[Y]]e- = Tr (x, Mt) + [Br]as,
LE [T),
where X. stands for the transposed of X, and Tr (x, My) denotes the trace of the square matrix x, Mi. In this model, the objective is to estimate the np matrices (MG)) This problem is equivalent to considering np independent trace regression models (38). Hence, this naive approach ignores the multivariate nature of the possibly correlated entries of each response Y,.
To overcome this issue, we consider for E [n] and j E [p] the vectors a; R and B3 E R2 and assume Mf₃- = a;33. For identifiability issues we assume that for all i E [m]. ||a=]] = 1. This model ensures that the matrices (M,)) share a common structure. It also implies that they share the same rank, namely one here. Hence this model now accounts for the multivariate nature of the problem. It has n(m-- - 1) + pq free parameters and rewrites as follows:
[Yly = Tr (x.T (,,3)) + [Felag. t t [T).

---

Consider A the matrix obtained by stacking the n row vectors and B" the matrix obtained by stacking the p column vectors By. This leads to the BMLR model (
1.2 Related works
The BMLR model, first introduced in [23] has gained notable attention powerful framework for examining relationships between matrix-structured responses and predictors. the development of estimation techniques for this model, two principal methods have been explored: alternating mini- mization and spectral aggregation.
Alternating Minimization, usually presented in Maximum Likelihood Estimation [16] or least squares Estimationn(SE)) context [23]. As the objective is non convex in both parameters, à two-step iterative algorithm is usually derived construct the estimators.
Spectral Methods, presented in a factor model framework in 7. offers an alternative that lever- ages the spectral properties of the target matrices. Authors propose new estimation method, the a-PCA, that aggregates the information in both first and second moments and extract it via a spectral method. They show that for specific values of the hyperparameter a, namely method corresponds to the least squares estimator. However, the procedure is non convex and they rely on approximate solution by alternating minimization. More specifically they maximize row and column variances respectively after projection.
-1, the a-PCA
Kronecker Product Factorization: The Kronecker Product Factorization IKRO-PRO-FAC method has recently introduced new possibilities for estimation within matrix-valued linear regres- sion. [1]] present this approach, which leverages Kronecker products to decompose complex ma- trices into simpler components. A key advantage of this method is that it circumvents the need to estimate the covariance between individual entries of the response matrix, significantly reduc- ing computational complexity. The KRO-PRO-FAC algorithm is accompanied by non-asymptotic bounds on the estimation of the Kronecker product between the parameters. However, an important limitation of this approach is the lack of direct control over the estimation accuracy of each param- eter separately. This restricts its upplicability in scenarios where parameter-wise interpretability or precision is critical, leaving a gap that motivates alternative methodologies, such the optimization- free estimators proposed in this work.
Autoregressive Frameworks: In autoregressive settings. where X, Y,-1 in I the pri- mary focus is on capturing temporal dependencies by minimizing the residual sum of squares 48
24 43 These frameworks are partieularly relevant for modeling dynamic systems and matrix-valued time series. However, parameter estimation in autoregressive frameworks typically relies on computationally intensive procedures, such as iterative optimization or matrix decomposi- tions. These methods become increasingly impractical as the dimensions of the data grow, creating significant bottleneck in high-dimensional settings. This challenge is especially pronounced in applied studies, where the scale and complexity of datasets continue to expand, underscoring the need for more efficient and scalable estimation techniques.
49
1.3 Summary of Contributions
This work studies the estimation problem of the matrix parameters A and B in 1 Our contributions are the following:
Noiseless Case Analysis: In Section 2 we study oracle case and establish that, in the absence of noise, the true parameters can be exactly retrieved. This analysis highlights the fundamental identifiability properties of the model.
3 propose explicit, optimization-free estimators B and Á defined in and (4) respectively. They significantly simplify the estimation process and are particularly advantageous in high-dimensional settings where traditional optimization- methods become computationally prohibitive. Theoretical guarantees are provided in Theorems 3.3 and 3.6, We establish non-asymptotic bounds characterizing the dependence of estimation accuracy on the problem dimensions (n. p. m, q) and the sample size T'. For Á E R"xm, the performance improves with larger sample size T and larger values of p and q, the column dimensions of the response and predictor matrices, respectively. However, the performance deteriorates as n and m, which determine the size of A", increase. In contrast, for B E RR4xP, the convergence rate improves with increases in T and n. showcasing a "blessing of dimensionality" effect in the row dimension of the target matrices. Nonetheless, the performance decreases with larger values of p, q, and m.
Optimization-Fret Estimators: In Section
3

---

Numerical Validation: In Section we validate our theoretical findings through extensive numerical simulations. On synthetic datasets, we show that the empirical convergence rates align closely with theoretical predictions, see Figure [ For real-world data from the CIFAR-10 dataset, we demonstrate the practical effectiveness of our procedures by introducing controlled noise, estimating the correction matrices A and B on the training set, and evaluating denoising performance on the test set. The results, presented in Figures 2 and 3 highlight the effectiveness of the proposed estimators in practical scenarios.
Sparse Adaptive Estimators: Extending the framework to sparse settings, introduce in Appendix 3.3 hard-thresholded estimators, BS and ÁS, defined in (5) and (6) respectively. These estimators exploit the sparsity of the true parameters to achieve improved convergence rates, as established in Theorems 3.7 and 3.8 Moreover, we demonstrate that these estimators can recover the exact support of the true parameters with high probability, providing strong practical guarantees. These estimators exhibit the same dependency on the problem dimensions (n. p: m, q) and the sample size T their dense counterpart / and Á.
2 Analysis at the population level
In this oracle case we consider T observations (X,, M,) that satisfy the relationship
= A',,", t E [T),
(2)
where A E " and B" E R are the unknown parameters to be recovered. The matrix A is assumed to have rows with unit L -noom,, ensuring identifiability of the model.
Understanding this oracle case will allow to establish baseline results that will guide the analysis at the sample level (T) involving noise. The primary goal here is to derive conditions under which the matrices A and 8 can be uniquely recovered, as well as to propose efficient algorithms for their recovery. This involves leveraging the structure of the observed predictors (Xe)i= and the constraints on the parameters.
In this section, we assume that the matrices (X1)T= form a generating family of Rm*x, which implies that T is larger than mq. We define two matrices, M == R. RTxxp xnp and X = (vve((X))..... vec(Xr)) E RT'xm.
vec(
Remark 2.1. When the design matrices (Xt)T= are generated randomly, X is invertible under mild conditions [13 36 39 40).
We note F( the canonical basis matrix with 1 at entry (k, and define the unobserved matrix C mm C :=
()), ke(m), (E(g) where each row of C corresponds to the vectorized form of A B". The entry of C located at the k (1- 1)m-th row and i+ (j 1)n-
th column is denoted by [C]((.) By construction, each entry of the matrix C is defined as [C] [A*] [B',y, for all i E [n]]. k E [m]. E [4]. and [p). Moreover, in the model 2 the matrix C can be exactly reconstructed from the observations, as shown in Lemma C.2 in the supplementary material. Corollary 2.4 shows how A and 3* and can be exactly recovered from this quantity
(A.1)
(k.1)
Proposition 2.2. In the model (2). where the design matrices (X,)²= form generating family of
Rmxx, the parameter matrices A Rrxm and B RRxx satisfy the following relationships:
m
v(L.) E lgl
x (p)
BByl
(3J) (R.))
Wi [n],
(1,2)
(2,k) E [n] x
[m]
[A*]
(k.I)
for any (1,5) [a] x [p] such that (B")i # 0.
[B
Remark 2.3. For fixed (1,j) E [g] (p). the entries (191 (i.j) (M.I) (L.b)
share the same sign, as each entry
is the product of [A*] which is non-negative and [B*]y-
The following corollary provides a representation of the entries of A and B as averages. This characterization will be particularly useful at the sample level, leading to plug-in estimators.

---

denote the set of indices (j. 1) such that [B" * 0 and let F: [n] x mm] Then the entries of the matrices A and B* can be expressed as:
Corollary 2.4. Let Do C [p] x [9]
Σ
(C (R,))
[A] Jak
U
1 [B
Σ
[C]
(C)
(kd))
Σ (i.jDDo
BB*)y"
(1,R)ET
Remark 2.5. The magnitude of [B*]i can be expressed as [B>)|s]
Σ
Σ
[
[0
3 Analysis at the sample level
We now consider I' observations (X,, Yi) that satisfy !
Our objective is to estimate A and B
From the observations (Xe, Yi) we construct Y (vee(Y)))..,,vve(YY))) E Txpp and
X (vec(X1),.... vec(XT)) RTxxm
We assume that the design matrices (X+)== form generating family Rmxq, This assumption implies that T mg and ensures that X X mm is invertible. Following Remark 2.1 this is mild assumption. We further define the unobserved noise matrix E : (vec(E)))... (E))..... .ece)..)......)) ()) Txxp and the unobserved signal matrix M := (vec( (.................... vec(MT))
E Tx where
(EAT= are the noise matrices defined in T and for E [T]. M₁ == A" X,
Following Lemma C.2 we define the unobserved matrix C = X x
X E Rmqxxp To
analyze the influence of noise in the estimation process, we define ID ==
X
X pmpxnp
Finally, we derive from the observations C = xx X Y Rmgxxp This leads C C+D.
(6.))
We assume that for all (i. k)) E [m] [p] x [g]. the row sums Σ
are nonzero. It ensures
(k.1)
that the plug-in estimator for A is well-defined. Notably, this assumption holds almost surely when the noise matrices (Ea) are drawn from a continuous distribution, the case in this study.
3.1 Definition of the estimators
Leveraging the results from Corollary 2.4 we define the plug-in estimators / E R2*p of B", defined for all (j.1) E D == [p] [q] and the plug-in estimator À E Rexm of A", defined for all (i, k) [n] x [m]. follows:
Σ
A (h.))
(k..)
[B]y =
Σ (1,R)E
[A]ia =
(k.1)
Σ
[B]
(..)ED
where [B1(]
= n-1 Σ
1
n
m
(rd) (a.1)
In the definition of [Ā]ER, we ensure that the terms in the
numerator do not appear in the denominator to preserve statistical independence of both terms. important to note that the entries of À are defined as the ratio of random variables. While the
behavior of such ratios has been studied in the literature, particularly in Gaussian cases [17 34), the results obtained through this approach would require heavy assumptions and remain challenging to interpret in context.
When the Gaussian variables in the ratio are centered, the distribution of the ratio known as the Cauchy distribution 25 However, for non-centered Gaussian variables, the probability density function of the ratio takes significantly more complex form [22), making the analysis cumbersome. Under certain conditions, it is possible to approximate the ratio with normal distribution [15), but this requires additional assumptions on the model and would still yield results that are difficult to decipher. Consequently, we opt for different approach that avoids these complications while retaining interpretabiliity Specifically, observe that the plug-in estimator À does not fully exploit
5

---

key property of the model: the entries of the matrix A are constrained to lie between and 1. This additional structure could be leveraged improve the estimator's performance. Hence we define the estimator À T of A", defined for all (i, k) E [n] [ml] as:
[A]i == max (0, min
3.2 Theoretical analysis with known variance
We present the matrix normal case with known fixed variance under the ORT assumption. The matrix normal distribution generalizes the multivariate normal distribution to matrix-valued random variables 20 35]. The sparse case is presented in the section 3.3
Lemma 3.1. Under the assumption that X TTx full rank, the matrix C
( X
X x
R.paxp satisfies C ~~ MM mqxxp
C,
X X o TTTp)
Following the vast literature on linear regression and Gaussian sequence models 44 19
37
12
make the ORT assumption to capture a better understanding of the phenomenon at play. This assumption serves primarily to facilitate the theoretical analysis. Notably, the numerical experiments in Section are conducted without relying the ORT condition. We discuss relaxation of the ORT assumption in Appendix(A)
Assumption 3.2 (ORT assumption). We assume that the design matrix X satisfies X T
X = T -mq-
Under Assumption 3.2 Lemma 3.1 ensures that the entries è are independent and normally
distributed. The following theorem establishes non-asymptotic upper bounds the convergence rates of B under various norms.
Theorem 3.3. Under Assumption 3.2 the estimator B introduced in 33 satisfies the following non- asymptotic bounds for any 0:
P P[IB 1-11++++
pqo 22m
exp(
Tne 2ma2 ), PP P >
S 2expl
Tne
+
2p
TvTnn
3mo
3
eVVT
[1IB- - -B'llop > c] 52exp
Vna
20 2mm
Here, Vi.. VF+ V || ||+ is the elementwise maximum norm, || llop the operator norm, and ||= the Frobenius norm.
We observe that the convergence rate of the estimator В exhibits the anticipated dependence on the sample size T, improving as the number of observations increases. Notably, our analysis reveals "blessing of dimensionality" effect, wherein the convergence rate accelerates as the row dimension of the observed target matrices (Y:) grows. Conversely, the convergence rate is negatively affected by the size pq of B* Furthermore, the variance parameter also exerts a detrimental influence on the convergence rate, as intuitively expected.
The following Lemma provides probabilistic control over the event where the plug-in estimator Á. defined in (3), coincides with its modified version Á defined in (4). We assume that the entries of A are strictly bounded from below by 0 and from above by 1. We also assume that the sum of the entries of B* positive.
Assumption 3.4. We assume that for all (i, ()) E [n] x [m]] we have 0 < [A*] < 1. In addition we
1
assume that 3 0 where 3
B*yy-
pq
In model D. A" has non-negative entries with rows summing to one, its entries lie in [0, 1]. Assumption 3.4 strentttss this to entries in j0, 1[. The sparse case is discussed in Appendix 3.3
Lemma 3.5. Under Assumptions 3.2 and 3.4 the estimators À and Å introduced
and
satisfy the following property for any (i, k) E [n]
[m]:
[A]i # [A]1a]
T 2
T ? Vik T
T
exp
exp
Pak T
252
20²

---

where P4k
1 Pq
EL=
E)=1 [B [A*)a)) M4k [a and ā v/11 -
We first note that under Assumption 3.4 the quantities and Pi are positive. We then observe that as the sample size T increases, the plug-in estimator more likely to coincide with its modified version A. The intuition behind this result is straightforward: as the sample size grows, [A]is con- verges to [A*] which inherently lies between and Consequently, the modifications introduced Á become unnecessary as [A]ik naturally satisfies the constraints of the model. similar phe- nomenon occurs
the number of columns of the response matrices, and q, the number of columns of the predictor matrices, increase. Larger and effectively provide more information about the structure of the model, leading to improved accuracy the plug-in estimator and reducing the need for corrections by A.
The following theorem provides sanite-sample analysis of the performance of Å. Theorem 3.6. Under Assumptions
3.2 and 3.4
the estimator À introduced in (4) satisfies for any
<> 0:
2c
nmo /
Tpqc² 20²
Vm
Tpq(n 1)c²
P
A
A'
exp
exp
|+
<V VVpTTT
(n -
2mo2
Tpqu²
Tpqu²
exp
a gexp
nm
26
20²²
where v ming,k Vik- / mingk Pik with paks Vik and ē defined in Lemmal 3.5
We observe that the T. reflecting the benefit of more observations. Additionally, the convergence rate is positively in- fluenced by increases in the column dimension p of the observed target matrices (Yı)²=1 and the column dimension q of the predictor matrices (X,)²=1- This reflects a "blessing of dimensionality" effect, as additional columns provide richer information for the estimation process. Notably, this behavior contrasts with that of B, detailed in Theorem 3.3 where increases in p and q have detrimental impact. Moreover, the magnitude of 3 plays a crucial role in determining the conver- gence rate, with larger values of | leading to faster convergence. Conversely, the performance of À deteriorates as the size nm of A" increases and the variance parameter o negatively affects the convergence rate. Higher noise levels degrade the accuracy. We note that the degradation of both rates with increasing m breaks the symmetry between A" and B". This is because of the assumption we impose on A
performance of the estimator Å improves with the sample size
3.3 Sparse Case with Known Variance
In this part, extend theoretical analysis to incorporate sparsity assumptions the parameters A" and B" By leveraging the sparse structure, we aim to develop estimation strategies tailored for high-dimensional settings, where many entries of the parameters are expected to be zero. This sparse framework addresses practical scenarios where dimensionality reduction is critical.
We propose the whre-hheeshlding estimator [19 37 BS defined as follows for all (4,5) E [g] x [p]:
[B> = [ (11))
(5)
where [B] are the entries of the initial estimator B defined in (3 and 0 is user-defined threshold. The hand-thresholding operation enforces sparsity by setting small entries of ß to zero, aligning the estimator with the sparse structure of B*
The following theorem establishes a non-asymptotic upper bound on the convergence rate of D$ for the Frobenius norm under sparsity assumptions on the parameter 1²
Theorem 3.7. Under Assumption
3.2 for any å (0, the estimator ss introduced in (5)
with the threshold - == V. nT (Vlog(2pg) + vlog (2)) enjoys the following non-asymptotic properties on the same event holding with probability at least 1 1 ó:
7

---

1. If lf denotes the number of nonzero coefficients in B", then
|
B -
16
2. If the entries of B satisfy jElle > 3r, then the support of is perfectly matches that of B3, namely supp( supp( where supp(-) denotes the set of in- dices corresponding to the nonzero entries of matrix.
Theorem 3.7 highlights the performance of the sparse estimator 8 under sparsity assumptions
on the true parameter B" The convergence rate of BS improves the sparsity of B* increases, demonstrating the benefits of leveraging sparse structures. The threshold exhibits favorable scaling with the sample size I and row dimension both of which contribute to reducing т, enhancing the estimator's performance. Conversely, T increases with the dimensions p and q, reflecting the greater challenge of estimation in higher-dimensional settings. Moreover the threshold scales with o,, capturing the adverse impact higher noise levels on the estimation accuracy. This emphasizes the observations from Theorem 3.3 Finally the condition the minimum magnitude of the entries of B ensures that the threshold T enables exact recovery of the support B", with high probability.
We now propose the hard-thresholding estimator A5, defined as follows for all (1, k) E [n]] [m]]
(b.f)
where 7ik
Pq
(N:0))
[A = [A]as (1ie 22r)
The following theorem establishes a non-asymptotic upper bound on the convergence rate of A for the Frobenius norm under sparsity assumptions on the parameter A".
Theorem 3.8. Under Assumption
3.2 for any ô E (0, 1). the estimator A defined in O with the (Vlog( 2mm) + Vlog (3) sutisfies the following non-asymplotic prop-
threshold == VV Tpg
2
erties on an event holding with probability at least 1 26:
1. If ||A" denotes the number of nonzero coefficients in A", then
I - 1 S | (2ts + 37)²
where
3" =
Σラ El BB and ts satisfies a
exp
Tygt 2mo
+
1 Vnpppp
1s 1/p''
ay2
exp
Tpqt2
6.
2. If the entries of A satisfy (ni [e|m)][A"))) [|A))a] > 37 and 3/ 1, then the support of A perfectly matches that of A", namely supp(A (AS) (A*)(A").
Remark 3.9. The parameter to, which determines the concentration properties of the estimator, de- creases as T, and increase. Consequently, the estimator A benefits from improved performance the sample size T' and latent dimensions p and g grow.
Theorem 3.8 characterizes the non-asymptotic properties of the sparse estimator ÁS under sparsity assumptions the true parameter matrix A". The first result provides a Frobenius that scales with the sparsity level || A This bound inversely related to the squared magnitude of 3* (the average entry B*), indicating that stronger signals in the underlying parameter matrix 3 lead to improved estimation accuracy. This is similar to the phenomenon described in Theorem 336 The error bound also depends on both the threshold parameter and the concentration parameter t.s which capture the impact of noise and sample size on the estimation performance. The threshold exhibits several important dependencies. It decreases with the sample size T, reflecting improved estimation with more observations. It similarly decreases with dimensions p and g, showcasing beneficial effect of higher dimensionality. Conversely it scales with the noise level a, capturing the detrimental impact of increased noise levels. Finally it grows logarithmically with the matrix dimensions n and m, indicating mild sensitivity to the size of A". The second result establishes conditions for perfect support recovery of A". Specifically, when the minimum magnitude of the
error bound
8

---

nonzero entries in A exceeds three times the threshold T, and the average effect 3 is sufficiently large (3,5* 1). the sparse estimator exactly recovers the support of A*
As noted in Remark 3.9 the concentration parameter is decreases with larger values of T, p. and q- This property, combined with the similar behavior of т, demonstrates that the sparse estima- tor benefits from both increased sample size and higher latent dimensions, a particularly favorable characteristic for high-dimensional settings.
4 Numerical Simulations
4.1 Synthetic data
Now we evaluate the performance of the proposed estimators through numerical simulations. Simulation Setup: The simulations involve the generation of matrices A", B*, Xx,, Es and Y,. By 2000, and = 1. analyze the effects of n, m, p. 9. T, and
default, the parameters are set as n = 15, m = 13, = 14. q = 12, These default parameter values adjusted performance of the proposed estimators. A is a n x m matrix with random entries sampled from uniform distribution over [0, 1) and rows then normalized random entries sampled from uniform distribution over [0, 1). (X₂) is sequence of I' matrices of size X q with random entries sampled from uniform distribution over [0, 1). (E,)
the
to 1. B" is a q X p matrix with
a sequence of T noise matrices of size n X with entries drawn from a Gaussian distribution with mean 0 and 1. (Y)) is sequence of T observation matrices, where Y, A" X Eg-
standard deviation
Estimation and Evaluation: The estimators Å and / computed using 4 and B respectively. To evaluate their performances, we vary the parameters n, m, p, and 7 to observe their impact on the estimation accuracy. For each parameter setting, we compute the Frobenius norm of the errors
||A-- A|F and Bl F together with the Operator norm of the errors ||Å-Allop ol and ||
Bllap. These errors are plotted as functions of the varying parameters in Figure !
From the plots, we observe that as T increases, Figure Te the errors in À and 13 decrease, indicating improved estimation accuracy with more data. As p and q increase. Figures Ic and ld the errors in A decrease and the errors in B increase. As m increases, Figures lb both the errors in Á and B increase. As n increases, Figure la the errors Å increase and the errors in 13 decrease. These results confirm the theoretical properties of the estimators, detailed in (rcorems{..6aamg 3.3 Additionally, we have performed numerical simulations to support the statement from Corollary
Validation of Theoretical Properties:
2.4 Appendix B.2 provides additional experiments.
(a) EV w.E.
(b) EV W.E. m
(c) EV w.t
(d) EV w.t (c) EV wat T
Figure 1: Evolution (EV) of the Frobenius norm (resp. operator norm) of Å - A (in blue, resp. in green) and of B - B" (in orange, resp. in red) with respect to (w.zt.) different parameters.
We observe that the empirical error rates align closely with the theoretical rates (derived under ORT). It suggests that deviations from orthogonality may lead only to mild degradation in performance. Moreover, the degradation of both rates with increasing is confirmed by the simulations (Fig- ure 1b). Thus, this observed asymmetry is not an artifact of loose analysis, but a consequence of the model's 's structural assumptions.
4.2 Real-world data
We also evaluate our proposed methods on real-world data using the CIFAR-10 dataset. It contains 50,000 training and 10,000 test RGB images, each of size 32 x 32 x 3. The pixel values are nor-

---

malized to [0, 1] for computational consistency. Our goal is to simulate noisy image transformations and assess the effectiveness of the correction techniques.
Noisy Transformations: We simulate noisy image transformations via left and right matrix multi- plications with A and B respectively. First, define A T32 + c E1,, where Fg is a 32 x 32 matrix with i.idd entries from a standard Gaussian distribution, and A is obtained by normalizing each of sum to I. Similarly, matrix 3 is given by / 132 + E22 where Eg is a 32 32 matrix with i.i.d. standard Gaussian entries. The parameter controls the noise level in both trans- formations. For both training and test images, the noisy transformation is applied independently to each color channel. The transformation for given channel c E f1,2, 3} (corresponding to red, green, and blue) is defined as X (A*)-1.Xd) X(r) -(B*)-1 - where X represents the tth channel of the original image and x is the corresponding noisy version.
noisy
Correction Process: The correction process is learnt on the training set. From the noisy transformed images, we estimate the correction matrices A and ((e) using (4) and (3), processing color train
channel independently. Once the correction matrices A and B(o)) are computed, they applied to the noisy test images to reconstruct the corrected test images. For each channel in a test image, the corrected channel is computed X A X noisy B Figure 2 shows example from the test set, illustrating the original image, its noisy version for 0.02, and the corresponding corrected image respectively.
Figure 2: Original, noisy, and corrected versions of the 11th image from the test set for e = 0.02.
Evaluation of Correction Quality: To evaluate the effectiveness of the correction process, we ll.p. the Frobenius distance between the original image and its Xoosy
compute Da, Σ
(x x(cc)
3
noisy version, and Do.- Σ
||
1? the Frobenius distance between the original image x
and its corrected version. We plot in Figure DDo, and Do as functions of the noise factor averaged over the entire test set. Appendix BB. presents additional plots showing how reconstruction accuracy varies with the effective signal-to-norse ratio (SNR) under both Frobenius and max norms.
3000
Ongintiia t03p
1500
3000
500
Figure 3: Dan (blue) and Doc (orange), averaged on the test set, functions of Error bars indicate standard deviations.
Conclusion
The results demonstrate that the proposed correction process effectively mitigates the impact of noise. The corrected images closely approximate the original images, as shown by both qualitative (Figure 2 and quantitative (Figure3) metrics. This methodology generalizes well real-world data, underscoring the applicability of our framework beyond synthetic simulations.
10

---

Acknowledgments
This research was carried out while the author was affiliated with Cornell University, where the main theoretical development and experiments were conducted. The completion of the manuscript and submission process took place after joining Capital Fund Management (CFM). The author gratefully acknowledges the academic environment at Cornell for fostering this work and the support of CFM during the finalization of the paper.
References
[1] Shane Barratt. A matrix gaussian distribution. arXiv preprint, arXiv:1804. 1010, 2018.
[2] Nayel Bettache. Matrix-valued Time Series in High Dimension. PhD thesis, Institut Polytech- nique de Paris, 2024.
[3] Nayel Bettache and Cristina Butucea. Two-sided matrix regression. arXiv:2303. 04694, 2023.
arXiv preprint
[4] Esther Boyle, Luca Regis, and Petar Jevtic. Matrix variate regression tool for insurers. Variance, 17(1), 2024.
[5] Haïm Brezis. Functional Analysis, Sobolev Spaces and Partial Differential Equations. Springer Science Business Media, 2011.
[6] Florentina Bunea, Yiyuan She, and Marten H. Wegkamp. Optimal selection of reduced rank estimators of high-dimensional matrices. The Annals of Statistics, 39(2):1222-139, 2011.
[7] Elynn Y. Chen and Jianqing Fan. Statistical inference for high-dimensional matrix-variate factor models. Journal of the American Statistical Association, 118(542):0038-1055 2023.
[8] Elynn Y Chen, Xin Yun, Rong Chen, and Qiwei Yao. Modeling multivariate spatial-temports data with latent low-dimensional dynamics. arXiv preprint ar arXiv:2002.01.30: 2020.
[9] Rong Chen, Han Xiao, and Dan Yang. Autoregressive models for matrix-valued time series. Journal of Econometrics, 222(1):539955) 2021.
[10] Rong Chen, Dan Yang, and Cun-Hui Zhang. Factor models for high-dimensional tensor time series. Journal of the American Statistical Association, 117(537):94-116, 2022.
[11] Yin-Jen Chen and Minh Tang. Regression for matrix-valued data via kronecker products fac- torization. arXiv preprint, arXiv:2404. 192220 2024.
[12] Julien Chhor, Rajarshi Mukherjee, and Subhabrata Sen. Sparse signal detection in het- croscedastic gaussian sequence models: sharp minimax rates. Bernoulli, 30(3):2127-2153, 2024.
[13] Kenneth Davidson and Stanislaw Szarek. Local operator theory, random matrices and banach spaces. In Handbook of the geometry of Banach spaces, volume pages 317-366. Elsevier, 2001.
[14] A. Philip Dawid. Some matrix-variate distribution theory: notational considerations and bayesian application. Biometrika, 68(1) (8(1))265-2744 1981.
[15] Eloísa the distribution the ratio of two independent normal random variables. Statistical Papers, 54:309-322 2013.
and Francisco Rubio. On the existence of normal approximation to
[16] Shanshan Ding and R. Dennis Cook. Matrix variate regressions and envelope models. Journal of the Royal Statistical Society: Series (Statistical Methodology), 30(2):387-408, 2018.
[17] Robert Charles Geary. The frequency distribution of the quotient of two normal variates. Journal of the Royal Statistical Society, 93(3):442-444 446, 1930.

---

[18] Christophe Giraud. Low rank multivariate regression. Electronic Journal of Statistics, 5:775- 799, 2011.
[19] Christophe Giraud. Introduction to High-Dimensional Statistics. Chapman and Hall/CRC,
2021.
[20] Arjun K. Gupta and Daya K. Nagar. Matrix variate distributions. Chapman and Hal//CRC, 2018.
[21] Harold Henderson and Shayle R. Searle. The vec- rec-permutation matrix, the vec operator and kronecker products: A review. Linear and Multilinear Algebra, 9(4):271-288, 1981.
[22] David V Hinkley. On the ratio of two correlated normal random variables. Biometrika, 56(3):635 639, 1969.
[23] Peter D. Hoff. Multilinear tensor regression for longitudinal relational data. The Annals of Applied Statistics, 169-1199, 2015.
[24] Hangjin Jiang, Baining Shen, Yuzhou Li, and Zhaoxing Gao. Regularized estimation of high- dimensional matrix variate autoregressive models. arXiv preprint, arXiv:2410. 1320, 2024.
[25] Norman L Johnson, Samuel Kotz, and Narayanaswamy Balakrishnan. Continuous univariate distributions, volume 2, volume 289. John wiley & sons, 1995.
[26] lain M. Johnstone. Gaussian estimation: Sequence and wavelet models.
Unpublished
Manuscript, December 2011.
[27] Eukasz Kidziński and Trevor Hastie. Modeling longitudinal data using matrix completion. Journal of Computational and Graphical Statistics, 33(2):551-566, 2024.
[28] Olga Klopp. Rank penalized estimators for high-dimensional matrices. Electronic Journal of Statistics, 5(none): 161 1183, 2011.
(29] Dehan Kong, Baiguo An, Jingwen Zhang, and Hongtu Zhu. L2rm: Low-rank linear regression models for high-dimensional matrix responses. Journal of the American Statistical Associa- fion, 2020.
[30] Clifford Lam and 7etai Cen. Matrix-valued factor model with time-varying main effects. arXiv preprint arXiv:24066.00/28, 2024.
[31] Béatrice Laurent and Pascal Massart. Adaptive estimation of a quadratic functional by model selection. Annals of Statistics, pages 1302 1338, 2000.
[32] Bing Li, Min Kyung Kim, and Naomi Altman. On dimension folding of matrix-or array-valued statistical objects. The Annals of Statistics, 38(2): 2010.
[33] Charles F. Van Loan. The ubiquitous kronecker product. Journal of Computational and Applied Mathematics, 123(1):85-100 2000. Numerical Analysis 2000. Vol. Ill: Linear Algebra.
[34] George Marsaglia. Ratios of normal variables and ratios of sums of uniform variables. Journal of the American Statistical Association, 60(309):193-222) 204, 1965.
[35] Arak Mathai, Serge Provost, and Hans Haubold. Chapter 4: The matrix-variate gaussian distri- bution. In Multivariate Statistical Analysis in the Real and Complex Domains, pages 217-288. Springer, 2022.
[36] Madan Lal Mchta. Random matrices. Elsevier, 2004.
arXiv preprint, [37] Philippe Rigollet and Jan-Christian Hütter. High-dimensional statistics. arXiv:2310.192.. 2023.
[38] Angelika Rohde and Alexandre B. Tsybakov. Estimation of high-dimensional low-rank matri- ces. The Annals of Statistics, 39(2):887-33,, 2011.
[39] Mark Rudelson. Invertibility of random matrices: norm of the inverse. Annals of Mathematics, pages 575-600, 2008
12

---

[40] Mark Rudelson and Roman Vershynin. Invertibility of random matrices: unitary and orthogo-
nal perturbations. Journal of the American Mathematical Society, 27(2):993-338, 2014.
[41] Yaser Samadi and Lynne Billard. On a matrix-valued autoregressive model. Journal of Time Series Analysis, 2024.
[42] Michel Talagrand. Upper and lower bounds for stochastic processes, volume 60. Springer, 2014.
[43] Ruey S. Tsay. Matrix-variate time series analysis: A brief review and some new developments. International Statistical Review, 92(2):246-262, 2024.
[44] Alexandre B. Tsybakov. Nonparametric estimators. Introduction Nouparametric Estinaaiiniiiiiiiittttnnioo pages 1-76. 2009.
[45] Stephen Tu.
Upper and lower tails for gaussian maxima.
https:
//stephentu, wtt.:iii.tthbbb..oooeeeeeeeeeeee thee
upper and aat---nd-loos lower tails watss-ggassiiiaaiiiaiaa html 2017. Accessed: 2024-12-11.
[46] Roman Vershynin. Introduction to the non-asymptotic analysis of random matrices. arXiv preprint 10X.3,101 2010.
[47] D. J. De Waal. Matrix-valued distributions. Encyclopedia of Statistical Sciences, 2004.
[48] Cheng Yu, Dong Li, Xinyu Zhang, and Howell Tong. Two-way threshold matrix autoregres- sion. arXiv preprint, 2XXi:24477.10222 2024.
[49] Ruofan Yu, Rong Chen, Han Xiao, and Yuefeng Han. Dynamic matrix factor models for high dimensional time series. arXiv preprini, arXiv:2407.05624, 2024.
[50] Xu Zhang, Catherine C Liu, Jianhua Guo, KC Yuen, and AH Welsh. Modeling and learning high-dimensional matrix- variate sequences. Journal of the American Statistical Association, pages 1-16, 2024.
13

---

A On the Validity and Relaxation of Assumptions
In the main text, we assume that the noise matrices (Lo)f= have independent entries drawn from centered Gaussian distribution with constant variance o². This corresponds to matrix normal dis- tribution with isotropic row and column covariances, and facilitates clean non-asymprotic analysis. However, our results can be extended to more general settings under mild additional technical effort.
In some proofs, we extend the analysis to a more general setting in which each E, is drawn indepen- dently from matrix normal distribution MNN (0 E,, Σx) [14 47 20 (5], where 2, R=** and Eg E R2XP capture row-wise and column-wise dependencies in the noise. Imposing the nor- malization condition Σg @ E, o² mp corresponds to the setting presented in the core of the paper.
Under this condition, the general matrix normal setting is equivalent to the isotropic case presented in the core paper.
This formulation allows us to present certain proofs, such as Lemma 3.1 within a broader frame- work.
ORT assumption. Consider the setting where each E, follows matrix normal distribution (0, Σ, Ex). where Eg @ Σ, = I but without the ORT assumption. In this relaxed setting, Lemma 3.1 remains valid and continues characterize the distribution of the central quan- tity C, although the matrix now exhibits an anisotropic row covariance structure.
Extending the proofs of Theorems 3.6 and 3.3 to this setting is conceptually straightforward, but technically more involved. One would need rely on concentration bounds for random matrices with independent but non-isotropic columns. These are available, for instance, in Section 5.5 of [46). In this case, collinearities in X X would naturally appear in the resulting bounds, illustrated in results such as Theorem 5.62 of the same reference
Homoskedasticiity Assuming Σ, & Σ, = o or equivalently that (E₂) have independent entries drawn from a centered Gaussian distribution with constant variance o² plays role analogous the homoskedasticity assumption in classical linear regression. It reflects a setting in which the signal is entirely explained by the model, and the residuals are unstructured. While restrictive, this assumption is standard in theoretical analysis and provides a tractable foundation for deriving error bounds.
Relaxing this assumption to allow for general row-wise dependence while keeping column indepen- dence, still under the ORT Assumption, would require concentration results for random matrices with independent but non-identically distributed rows (see Section 5.4 of (46)) In this case, the noise term in our analysis would be governed by the full tensor-product covariance Eg 8 Ex, which would explicitly appear in the resulting bounds (e.g., Theorem 5.44 in 46)).
Relaxing both assumptions simultaneously presents a more significant challenge. To the best of knowledge, current probabilistic tools do not yet offer sharp and tractable results in this fully general setting. However, this remains a promising direction for future work, potentially requiring new matrix concentration inequalities tailored to specific structured settings.
In summary, although the core analysis assumes isotropic Gaussian noise for clarity and tractability, the main techniques extend to more general noise structures under appropriate assumptions and with access to suitable matrix concentration inequalities.
B Additional Numerical Analyses
This appendix presents complementary numerical studies that assess the robustness of the proposed estimators and quantify the alignment between empirical convergence rates and the theoretical pre- dictions.
B.1 Robustness to the Distribution and Normalization of A and B*
To verify that the simulation results Section are not sensitive to the amplitude or normalization of the true parameters, repeated all experiments with entries of both A" and 3² independently
14

---

drawn from a Uniform[0, c) distribution with E f1,2,3,5, followed by row-wise 4 normal- ization of A to ensure identifiability. Across all values of c, the empirical convergence rates and qualitative dependencies with respect to the parameters (n, m, P. 9; 7)) remain unchanged. This con- firms that the results reported in Figure aree not specific the original choice c 1.
Theoretical analysis shows that the model is identifiable once one of the parameter matrices is nor- malized. If instead the normalization were imposed on B*, the algebraic expressions of Proposi- tion 2.2 and Corollary 2.4 would be modified by exchanging the roles of m (the number of columns of A") and q (the number of rows of B*). Consequently, the dependencies on these dimensions in the non-asympiotic bounds of Theorems 3.3 and 3.6 would also be interchanged. This structural asymmetry stems from the identifiability constraint itself and highlights that several normalization choices are possible. The decision to normalize A is primarily motivated by interpretability its nonnegative, row-stochastic structure aligns with the notion of activation or mixing weights- and by analytical tractability of the resulting expressions.
B.2 Quantitative Comparison Between Empirical and Theoretical Rates
This section reports the quantitative comparison between the empirical convergence slopes of the estimators (Á, B) and the theoretical predictions derived from the finite-sample analysis. The ob- jective is to evaluate how the estimation emor scales with each model dimension (n, m. q) and with the sample size T under three norms: Frobenius, operator, and maximum absolute.
Experimental setup. For each parameter d E {n, m. q, T), we generated independent datasets while keeping all other quantities fixed, recomputed (A, B), and measured their reconstruction er-
rors
Erro(d) (d) = |M Mo De {Il ||r, || llap, lmmax The dependence of log (d)) on was then fitted with a linear model
log(Err((d) 0n
(0) ((d) +
where the function (dd corresponds to:
((d flogd. d, d, dE {n, m.p. (linear-log regression).
d= (log log regression),
This distinction reflects that, in the simulations, both axes were represented logarithmic scale for T, while for (n, m. q) the r-axis was linear and only the y-axis (error) was plotted in log scale.
((7) The fitted slope sd quantifies the empirical rate of variation of the estimation error with respect to d. For the sample size T, the theoretical prediction is s
Empirical slopes. Table summarizes the fitted slopes for all parameters and norms. Positive slopes indicate that the error increases with the corresponding dimension, while negative slopes indicate a decrease.
Table 1: Empirical slopes s
(0)
of log(Err) vs. log(d) for Å and ò under the three norms.
A
B
Parameter
Max
Frobenius
Operator
Max
Frobenius
Operator
+0.003
+0.011
+0.008
-0.011
-0.011
-0.012
m
+0.007
+0.015
+0.011
+0.018
+0.016
+0.016
-0.012
-0.011
-0.011
+0.002
+0.011
+0.008
q
-0.005
-0.005
-0.005
+0.008
+0.017
+0.013
T
-0.332
-0.496
-0.474
-0.579
-0.526
-0.526
Findings. Several patterns emerge consistently across all norms:
Sample-size dependence. For both Å and B, the slopes with respect T lie between -0.58 and -0.47, matching the theoretical prediction sT = - This confirms that esti- mation errors decay the expected 112 rate.
15

---

Dimensional dependencies. For A, enrors increase with m and decrease with (p, q), whereas for B the opposite trend holds-enrors increase with (m,p, q) but slightly de- crease with n. This asymmetric pattern is consistent with the theoretical structure of the model, where the identifiability constraint on A leads to mirrored dependencies in A and B.
Norm-specific behavior. The Frobenius norm shows the strongest sensitivity to dimen- sional changes, reflecting its dependence all matrix entries; the operator norm shows a weaker dependence consistent with a (vD v scaling; and the maximum norm lies between these two regimes.
Goodness of fit. The linear relationships between log(Err) and log(d) exhibit high ex- planatory power for most dimensions, with R² values above 0.9 in the majority of cases, confirming the robustness of the observed trends.
Plots. Figure illustrates the empirical error curves for all parameters. Each subplot reports the average reconstruction error (in logarithmic scale) for Å and B under the three norms. The first four panels correspond variations in the structural dimensions (n, m, q). while the bottom panel shows the dependence on the sample size T. The trends confirm the fitted slopes in Table I errors decrease approximately linearly on the log scale as T grows, consistent with the 7-1/2 rate, and vary smoothly with (n, m, p. according to the theoretical sign pattern. Notably, / exhibits larger sensitivity to and smaller sensitivity to n, while Å shows the opposite, reflecting the asymmetric normalization of A Across all panels, the Frobenius norm produces the steepest gradients, the operator norm the weakest, and the maximum norm lies in between, reproducing the hierarchy predicted by the theoretical bounds.
Summary. The empirical analysis confirms that the proposed estimators obey the predicted non- asymptotic scaling laws. The 1/2 convergence rate holds precisely, and the dependence on (n, m. p, 9) follows the signs and magnitudes expected from the theoretical bounds. These quantita- tive findings validate that the theoretical error expressions accurately capture the dominant sources of variation in sample performance.
16

---

Erior VS.
Error VS.
Operaarr
Counter
10a
Frrer log
Einer (tgl:/ggaaaa scaaee)
100
Error VS.
Erro' VS.
Mix Abesuse Ahoter
Imbertat
Frner (lg scale)
Ereer (1g scale))
Focker Cperrsrr
30-1
10.2
20
60
80
100
20
80
100
Paramette
Paremeter
101
Froo dog scale
10-3
107
Parometerr
Figure 4: Empirical reconstruction errors of and B versus model dimensions and sample size T, displayed on logarithmic scale. Each curve corresponds given norm (Frobenius, operator, or max).
B.3 Reconstruction Error as a Function of Noise Level
To further evaluate robustness noise, we examined the reconstruction quality of the bilinear esti- mator under controlled perturbations of A" and B". For each noise level e, we generated perturbed matrices
A. where ZA, Z are Gaussian random matrices with i.i.d. (0, 1) entries, followed by row-wise nor- malization of A, to preserve identifiability. Each perturbed pair (A,, B,) was used to transform the training and test datasets, and the reconstruction was then estimated using the procedure described in Section 4
=
17

---

Effective signal-to-noise ratio (SNR). Rather than plotting the reconstruction error directly against e, we parameterize the x-axis in terms of an effective signal-to ratio (SNR), defined as
|X||
SNRF == 20 log
|X
Xooisylle
where is the original image and noisy its transformed version. This reparametrization provides scale-invariant measure of perturbation strength and directly reflects the degradation of signal energy in Frobenius norm.
Results. Figure reportss the reconstruction error as a function of the effective SNR for the Frobe- nius distance, while Figure O shows the analogous result for the element maximum norm. In both cases, the distance between the original and corrected images (orange) is consistently below that between the original and noisy images (blue), confirming that the estimator effectively compen- sates for multiplicative perturbations in A and B" The monotonic growth of both curves as SNR decreases quantifies the degradation rate, with the Frobenius error emphasizing global reconstruc- tion quality and the max-norm capturing local, glement-wise discrepancies.
Revardrurtios (robesi
Figure 5: Reconstruction error versus effective SNR (Frobenius norm). Lower distances indicate improved correction quality. The x-axis is expressed in dB following the definition SNRr 20 (o) -Xoony l()).
Recendrrction
Norm
Figure 6: Reconstruction error versus effective SNR (element-wise maximum norm). This comple- mentary metric highlights the preservation of fine-scale details under increasing noise.
Overall, the analysis confirms that the proposed correction procedure maintains stable performance across broad range of noise intensities and that the reconstructed images preserve both global and local structure in accordance with the theoretical robustness guarantees.
18

---

C Proofs
C.1
Proofs of Section 2
We first analyze a simplified case of (2) where the design matrices (x,)== are the elements of the canonical basis of TRmx In this setup predictor X, is a matrix with all entries set to zero except for a single entry equal to onc. This setting well studied in the vector case, typically referred
given T mq observations (X,, M, where each
Sequence Model. Such models have been widely explored in the literature [4 26]. By focusing on this simple scenario, we aim to uncover insights that can extend to more complex cases.
Formally, for (k. I) E [m] x [g]. let M(R.I) corresponds to the basis element E(,,), where ((.,)) is the matrix with all entries equal to zero except for one in the k-th row and 1-th column. The model 2 can then be expressed as:
M(a A ((k.)) B"
which corresponds to (2 for L = h (1 - 1)m. The matrix (( acts as a selector, isolating the effect of the k-th row of A and the 1-th row B
Lemma C.I. In the model (7), the parameters A and B can be explicitly recovered from the observed matrices M(K.t) as follows:
v(45) [9] x pl: [B >
Manle Vi e [n].
k=1
((i [n] [m] [A [Maanl
[B*
for any (1.j) E [g] lol such that [B*]i +
Proof of Lemma C.7 From the model 7 we have
M(r A E B",
where E(E.)) is the matrix with a single entry of 1 at the (As, l) position and zeros elsewhere. Expanding this relationship element-wise gives:
[Ma.nlls [A [B",,, Vii,k,,;)) [m]
[m]] [g].
Step 1: Recovery of [B*,
Summing over & E [m] for fixed i, we observe that:
m
Σ
Mannle = ΣIA (B" lyy-
k=1
Since the rows have unit Lorm,, it follows that:
m
[B*) = > Maanle v(I.. E [g] [p] [n].
Step 2: Recovery of [A*]
From the model equation, for each (Is, l), we isolate [A by dividing [Maanls by [B*]y. provided that [B*) bi # 0:
[Mann v(i,k) E [n] x [m], where [B'] + 0. [Ba)
[A
Lemma C.2. In the model
2 where the design matrices X form generating family of RRxx the unobserved matrix lg maqupp satisfies the following equality:
T
Here, M = (ve((M))....e(((AM)) Txxp and X = (vec(X))...,ee(((x))) E Txmq
19

---

Proof of Lemma C.2 We notice that M XC by construction, which concludes the proof.
Proof of Proposition
From Lemma
the unobserved matrix can be expressed as:
T
This means that C be computed directly from the observations (X₁, Me)eeftr.. provided that T invertible, which is ensured by the assumption that the design matrices (Xt) form generating family of mxxx
Step 1: Relating C A and 8
By the definition of C, each entry (C)) corresponds to:
(i.f)
[C] (+.) (k.l)
= [A*]a [B*],
for all i E [n], k [m]. E [g]. and j E [p). This bilinear structure allows us recover A and B separately by exploiting their roles in this product.
Step 2: Recovering B*
To isolate [B*]i,, sum [C] over all k [m]] for a fixed L, d. Specifically;
Σ
Σ ((A")s [B*)
Since the rows of A satisfy the L1 norm constraint (i.e., E [A7] 1), the summation simpli- fies to:
m
Σ
(9 [B")y.
This establishes the second equation in the proposition:
(4)) BBa
(2.)
Step 3: Recovering A
To isolate [A") lik we use the bilinear relationship:
[C [A*]a [B"]---
For fixed i, k, 4,j, solve for [A*] provided that [B*] # 0:
C (k.l)
[A sk
This establishes the first equation in the proposition:
[C] (i.j)
[A*]
(k.0)
B # O. for
B 1g
Proof of Corollary 2.4 The result on 8 follows immediately from Proposition 2.2 the result on A we need prove that for any n E No, if there exist (01, an) E R and (SI..., S.)) E R such that:
To prove
On = On
On
31
B
then it follows that:
n
Σ
1=1
On
Σ
3
1=1
20

---

Specifically, since the entries of A can be expressed as different ratios, all being equal, applying this result to the equations satisfied by A in Proposition 2.2 completes the proof.
We prove this result by induction on n. Initially the result is trivially true.
Assume the statement holds step n, i.e., if
In =
B1
n
then:
n
Σ Oi
7/
Σ
31
f=1
We now prove that the statement holds at step n + 1. Suppose:
an+1
B1
3
By the definition of 7 and the assumption of step n being true, we can write:
aa+1
n++1
and
7n+1-
Bn+1
Σ
B
j=1
This allows write:
Σ
02
Σ
3j
On+1
But1
f=1
f=1
Hence we deduce:
n Σ Bj
Σ + On+1
Gnti + ant1 Bati
Σ By + Bn+1
By + Bm+t
1=1
8, 1=1
+ 1
Pn+1
0n+1
Σ
Bg + Bn+1
J=1
Finally multiplying the numerator by B,+1 and dividing an+1 by the same factor provides:
04 + an+1
On+1
1nt1.
Bnt1
Σ j=1
Bg + Bu+1
Thus, the statement holds at step n + 1.
C.2
Proofs of Section 3
Proof of Lemma
3.1
We first recall that, by construction, we have:
C = [x xx)
x
M
(x
C
21

---

Next we prove the following, which will conclude:
MN mqxnp (0moxpps (x x)
-1 Ee 8 x
First we recall that the noise matrices (Fx) are assumed to be independent and follow the same matrix normal distribution, namely MMN nxp (Onxp: Ex. De). Hence by definition of the matrix normal distribution:
vec(Er)
itd. ~ Nup (Omp; Ee @ ).)
Txx defined in (??) as follows: This leads to derive the distribution of the matrix E E R
E MNTxxp (0Txxps IXT: Le E,).
Finally, the matrix (x "0)00x1x x being of full rank my and (x xxx))= being symmetric, by def- inition of the matrix D myx and by property of affine transformations of matrix gaussian distributions, we conclude that
D ~ MN maxup (OTxmp: X Xx
Ee 2 Σ,,
Proof of Theorem 3.3
We start by deriving the finite-sampled inequalities on B. From Lemma 3.1
Eg X E,
and Assumption 3.2 we deduce that vec
( Nuagmp
(C)
As mentioned in the
introduction D. Σ. o² Lnp which ensures that all entries of vec (C) are independent and follow
the same Gaussian distribution. Hence for all (k, (, E [m]] x 1q] x [n] x [p] we have C(k.))(6.))
(0. o²
((.)). [d] x [p].
1id N
Using the results from Proposition 2.2 and 3 lead to, for all (1,5)
[B]r - [B") Lid
0, mo²
(8)
nT
Mill's inequality, Theorem E. ensures that for all (1,j) E [q] x [p]. for any 0,
P [B} [B") (
2m
nl''²
nTw
exp
2mo
The first inequality is then deduced by using a union bound.
The second inequality is immediately derived from the concentration of extreme singular values of Gaussian matrices with independent entries, see Corollary 5.35 in [46..
[B], [B))s) follows mo²
For the third inequality, we first note that || 1/3-B% B 11
Σ
n7
[=1j=1
(2(g)) distribution. Then, Lemma from ensures that random variable Z following a x² 2pg) distribution satisfies, for any >> O:
(Z
2 pq p/p + 2e) & exp(-e).
Hence we deduce that
m
IP
[B - B" 1
(pg + 2/ppe + +2e) s exp(-c).
n
We notice that, for any c > O, the following inequality holds:
mo²
nT (ag + 2//qc + 26) S mo² (2pg + 3c)
nT
It leads to:
mo² n'r'
(1 BT
(2pg 3e) s exp(-e))
The stated result follows by change of variable.
22

---

Proof of Lemma 3.5 Let us fix (i, E [n] x [m]. Using the definitions of A and A from
4 and
B), respectively, find that [A]xx [A]a if and only if [A]a [0, 1l. Hence [A]u [A] if and only if the event holds, where:
(GJ)
f=
(k.))
A
0
ΣΣ
Σ
1
This event is realized and only if the event A+ or the event A_ holds, where:
(i,5)
p
A+
B
j=1
j=1
1=1
and
(6,5)
A_
- 0
(k.))
1=1
j=1 1=1
Using Fréchet inequalities stated Theorem
wc get:
P [A] 2 max IP
[AA
()-))
Under Assumption 3.2 Lemma
3.1 gives:
o² 7(0.0) ~ N (R.))) pqT
where:
((.,)) = 1
(2,))
pq j=1 1=1
(R.1)
and
Ma.) =
(k..0..6.)))
pq
j=1 1=1
Under the same assumptions, as detailed in the proof of Theorem
3.3
o2m
B ~ N
8", (n 1) p)pT
(10)
where:
B.
p
B
and
8
1
P
9
BByg
pq
pq
By construction, 9(..,) and B, are independent. In addition, T(n.)) ensures that their expected values satisfy:
[A 11* Assumption 3.4
((K.)) < g
Using the symmetry of the Gaussian distribution, we deduce:
and P[A2P((A)).
The event A+ holds if B, and B, hold simultaneously, where:
B = f0 < ) B3 = (AUM S 3..}
Using Fréchet inequalities, stated in Theorem E.2
P [44+] >
P
+ PP 1.
23

---

We now bound from above the probability of the events B3 and B2 separately. For B1, we note:
B1 = { [T(ke)) 2 7(.00) - }(..0)
The complementary event Bi is:
B1 (3,.0)) S S(()))) S(m.)))
As we have (kk,) O, using ( and Mill's incquality stated in Theorem E.1 provides:
Tpq (J(
P [1 <
exp
Mn,i) T
202
For Bg we note that:
B2 -7 2 (8" - -((k)) - (3: 9((ms))
The complementary event B2
B2 8 T(k,)) S (B* - -(ks))
- (3: ((k))
As we have 3 (((.)) 0, using independence of ((a,)) and B,, results from 9 inequality ensures:
10 and Mill's
Tyq B ((Me))
exp
n-
20²
[2
(3+ M(. 12pq/
Finally, using:
P[A-]21-P [B.]] [B2]
we deduce the result.
Let us fix ki,k) [n]] x [m]. We work on the event A = {LA)ue IAlmn
Proof of Theorem 3.6
Using we get:
p
ik
pg
pd
(k.1)
j=1
1=1
Proposition 2.2 guarantees:
A
Hk
B
C (1.0.115)))
pq
Pq
j=1
1=1
1=
By definition of 3 and using the reverse triangle's 's inequality we get:
188
[A]as - [A"]sx
S 9(k.)) - Rime)) +
[
B,
3r
where:
1
9
(1,5)
3 =
P4
pq
(k.I)
1=1
and
1
M(..)
[ (l
pq
j=1==
24

---

Under Assumption 3.2 99 and TOP hold. Using get for any 0:
and Mill's inequality stated in Theorem E.1 we
[0c.0 7(s0)) > > S
ov2
Epull" 20²
cvpqTT exp
Using 10 and Mill's inequality, we find for any 0:
1 3+ >e - o
@ 2m e/pqTin
EpeI'n
cv
exp
2mo²
where à = (n 1).
Finally, on the event A. we have [[]] s 1. We then use the result from Lemma 3.5 inequality stated in Theorem E.2 and we conclude with a union bound.
Fréchet
Proof of Theorem 3.7 For all (1,5) E [gl
[p]. g ensures that:
[B [B0] +
Elji
where:
ay itd. N (0. nT
mo²
From the results in Appendix where A the event
D obtain that the event A holds with probability at least
2m nT
Viog((20)) + log
max
Etyl
S
V
jEP
We recall the definition of the threshold:
2m
VTog(2pp)+
(2)
log
nT
On the event A, we observe the following:
If |B]]ا > 27, then:
[B"] 2 [B]bs] -leyl >
If [[ااز] 2r, then:
[B*uy| s [[]bu] + |eyyl 3T.
From these observations, we deduce:
1 - 1B if [ > 2T,
[B9 - [B"1
[1*)bb if [11B)yy] S 2T.
Rewriting this equality with indicator functions, we have:
| [B")y| =ley (1(8)c) 27)
|1B")yA (1(31c) s 27)
This leads to the bound:
| [B")
Sr1 $r1 (B |ul -)
3r1 (l 1u - 37).
25

---

Finally, we conclude:
[ - [B"ly S min (r, (r,1(1")s))
By definition of the Frobenius norm, obtain:
8 - B
S 16 16
This proves the first point.
For the second point, on the event A and under the stated assumption, we observe the following:
If [B*]| # 0, then |[B*| 3r. This provides the bound:
[Byl = |[B"), 15 + 4ny > 2r.
Hence [BB))1 0.
If [[[[B)] 1y # 0, then [B]|s]
> 2r. This provides the bound:
1ul 2 [B]yy] |a,l 2 T.
Hence [B*]s] [isl * O.
Proof of Theorem 3.8 For all (i. ki) E [n] x [m]. 99) ensures that:
[+] ik 17 "iik + Esk,
where:
700
(l(a0))00
pd
and
lid N
(0,
o²
esk
PqT
Proposition 2.2 ensures that
2ik 9A")"
From the results in Appendix where A is the event
we obtain that the event holds with probability least
(Viog(2mm) + log 8
max iE ke(m)
Cakl I sov Tpg
We recall the definition of the threshold:
2 5
Mogg(2m)) + log
VV
Tpg
On the event A, we observe the following:
If Hiak > 2r, then:
1r6el 2 Ffaa |eax 7.
Thus we have [A"]i # 0 and because [A*]i is bounded from above by we also have
3" ZT.
If [inl S 2T, then:
hnde S 1reel + lecal S 3T.
26

---

From these observations, we deduce:
[A]ix [AA
if liskl > 2T,
|A9 sk
[A*]a
if liskl S 2T.
Rewriting this with indicator functions and using the previously stated implications, we have: |[A) Jak - [A*]ax
S [A]ss
(bos > T), [[A"]ER| (hrxal $ 37)
[A*]
Multiplying all terms by 1381 leads the bound:
18-1 |[A"] - [A =|S [A]e [A7 1 (17Ek) > T),
+ hadl (htal - 3r).
Using the equality
(00 -1) + -B1AAe | [A] - [A4"|s
and the triangle inequality leads
|8*| |[A" [A S(M + -M (hrtal +)
+ 1 (hrul 33)),
where
[ = =18- - 31 [AAla
and
M ||AAss (
Moreover we notice that
M (2)
Pisk 7áll
Equations and
10
together with Mill's inequality, stated in Theorem E.
provide for any t > 0:
ov2m
Tpt
P 8r
> t S
cxp
T
202m
and
0/2
Tppt?
PP [1 Fskl $
exp
t/pql'T
20
Finally, we get for any 5 E (0, 1) with probability at least 5:
|881- [[AS] [A]]ul - 1l (Aj, + 0) 1 (188) 2 7) + (hill < 37).
Dividing by 18-1-1 and using that 7ie aa ensures that with probability least 1 6:
| lac
[A*]**
<2t 18 (4) 0) 1 (18*)~- -T (18)
+ [ (IA*)i - 331881111
Finally we get with probability least 1-6: 8:
||AS] [A*] --3-1-11 (A + 0) (2t, + 3r).
By definition of the Frobenius norm, we obtain with probability at least 1 - 8:
A - 1 < 137122 |l (2t, + 37)
For the second point, on the event A and under the stated assumptions, observe the following:
27

---

If [[<)) & O, then [A*| > 3T. Using the equality 7/x = [5* and the assumption
on the bound satisfied by 8 provide:
IRad = his + eaul 3,3 7 T > 0.
Hence |[AS]al o O.
If [[A4 * 0, then Iaal > 2r. This provides the bound:
Indp 2 Kuil |eyl 2
Hence hyy * 0 and thus [A* 7 O.
D Tails of the Maximum of Absolute Gaussian Variables
In this appendix, adapted from 45 we analyze the upper and lower tails of the random variable maxigien [X&el. where x.... X are independent and identically distributed (i.i.d)) random vari- ables following standard normal distribution ((0,, 1). Specifically, we aim establish bounds for the upper and lower tails of this random variable with high probability, which is crucial for understanding the behavior of maxima in Gaussian settings.
We begin by stating the main result.
Theorem D.1. Fix (0, 1), and let X1, Xn be i..dd N/O, 1). With probability least 6,
π
(og(n/2) log log 2/0) < max [Xil,
2
15i5n
XX S v2 (1og((2) + vlog(2/8))
max ISiSn
The asymmetry in the tails arises from the fact that maxi X; is bounded below by zero almost surely, is not bounded above by any fixed constant. To establish the theorem, separately analyze the upper and lower tails and combine these results with a union bound.
D.1 Upper Tail
The upper tail is analyzed using concentration results for the suprema of Gaussian processes. The key tool is Talagrand's 's concentration inequality, Lemma 2.10.6 in [2), stated in a simpler version follows:
Lemma D.2. Consider C R and let ~ I). Then, for any u > 0,
u S exp
u²
P s) E sup(t, 9)
u
LET
lET
2s
where suPter (L( 2//
To apply this inequality, we consider T = (El, (eI...., Cn. -l,..... J., where e; i-th canoni- cal basis vector in R". For g NO, (0,1), I). this gives
sup g) = max 19.1, and (t(tt = Elg?1 = seT
1515n
Applying Lemma D.2 obtain, for any > O:
P
max
Xx||
> E
max
[x||
/2lg)(1/)) < 8.
Isi5n
15:5n
Thus, with probability at least
max X
S
E
max X
+ lolo((//))
(11)
Snn
The next step is to bound E [max
X from above.
28

---

D.2 Expected Maximum of Gaussian Variables
We now bound mmaxigign X;]]- [X,]] Proposition D.3. Let Zi,..., Zn be random variables (not necessarily independent) with marginal distribution (0,,1) 1). Then.
E
max Zj
v log n.
Isisn
Proof. Fix any O. Observe that
E
e Sig
2, 5
E [e>2a) = =nc
1=1
Taking the logarithm of both sides,
A2
log E | maxisisa 21| S log n +
2
Applying Jensen's inequality,
12
AE
max 1
Z
- log EE' [ @amaxisisn Zi S log n1 +
2
Dividing through by X,
log
E
max
Z
IStS
Optimizing by setting À = /
log ni we conclude
E
max Z,
1 / logn. n.
IiSn
For XI,..., Xn i.i.d. following standard normal distribution (N 1), considering the 2n variables (X1,..., X,; X1 XI...., Xn), we have:
maxx
Xx|
S
222 log(2n).
ISiS
Substituting this into
completes the proof of the upper tail.
D.3
Lower Bound
We now analyze the lower tail of maxiçisal Xz). Fix a positive 7 > 0. Then,
Pr
max Isisn
X
ST = Pr ( S ..... Xx s <T).
Using the independence XI,..., Xn, we get:
Pr max
X1 S
II
Pr (IX, ST)
The Gauss error function, defined as
2
orf
:
c dt,
V J
allows then write:
T
Pr max ISisn
[Xil S 7
erf
2
29

---

The inequality exf(2)² < - 4x
which holds for all x 2 0, provides:
(1- a/2
Pr
max isiS
X
Finally, the inequality 1 S e
which holds for all R ensures:
Pr
max 18Sa
[x
ST S exp
2 -
To derive the lower bound, set
--- =ō, 6,
exp
2
and solve for We find that with probability at least 1 - 8,
max
[
2
log(n/2)
log log(1/6).
SSSA
E Probability bounds and inequalities
E.1 Mill's inequality
Theorem E.1 (Mill's Inequality). Let X be a Gaussian random variable with mean p and variance a². Then, for any > O, the following inequality holds:
P(X pa p> 1)
e
20
I 2m
By symmetry, we also have:
P(X p -))
tv 2m
Proof.
A proof of this theorem can be found in 37
E.2
Fréchet inequalities
Theorem E.2 (Fréchet inequalities). Let A and B be two events. The probability of their intersection satisfies:
(O + ( -) (PAN min(P((A),P(
The probability of their union satisfies:
max(P(A). & P(AUB) & min (1,P((((()) + +(B)).
30

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
The main claims in the abstract and introduction are carefully aligned with the paper's actual contributions, as detailed in Section 1.3 These claims are neither overstated nor misleading. They accurately describe the scope of the work, including the introduc- tion of optimization-free estimators, non-asymptotic error bounds, and the applicability to regression problems. The theoretical and empirical results presented in the core sections fully support these contributions.
Justification:
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made the paper and important assumptions and limitations. No or NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
fine to include aspirational goals motivation as long it is clear that these goals are not attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: The limitations of the work are discussed in Appendix A address the simplifying assumptions made in the theoretical analysis, such as the use of isotropic Gaussian noise, and explain how the results could be extended under more general settings. We also acknowledge the technical challenges involved in such extensions and provide references to the appropriate concentration inequalities required for handling non- isotropic noise or more complex covariance structures.
where we explicitly
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create separate "Limitations" section in their paper.
The paper should point any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The au- thors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested few datasets with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the ap- proach. For example, a facial recognition algorithm may perform poorly when image resolution is low images taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to ad- dress problems of privacy and faimess.
31

---

While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]
Justification: Each theoretical result in the paper is stated with complete and explicit set of assumptions. The main assumptions are clearly outlined, and each theorem or lemma is accompanied by full proof provided in Appendix C The proofs are rigorous, and where simplifications are made (e.g., assuming isotropic noise), we clarify this explicitly and discuss how the results can be extended under milder assumptions in Appendix
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should numbered and cross- referenced.
All assumptions should be clearly stated or referenced in the statement of any theo- rems.
The proofs can either appear in the main paper or the supplemental material, but they appear in the supplemental material, the authors are encouraged to provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be comple- mented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main experimental results of the paper to the extent that it affects the main claims and/or conclu- sions of the paper (regardless of whether the code and data provided or not)?
Answer: [Ys]]
The paper provides detailed description of the experimental setup, including the data generation process, model parameters, and evaluation metrics used in all simula- tions. Each figure in the results section is directly tied clearly defined experimental protocol. The code is included in the supplementary material and all information necessary reproduce the results relevant to the main claims and conclusions is fully disclosed. The experiments serve to validate the theoretical findings, and no critical step is omitted in their description.
Justification:
Guidelines:
The answer NA means that the paper does not include experiments. If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is a dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data often one good way to accomplish this, reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
32

---

While NeurIPS does not require releasing code, the conference does require all sub- missions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily new algorithm, the paper should make it clear how to reproduce that algorithm.
(b) the contribution is primarily a new model architecture, the paper should describe the architecture clearly and fully.
(c) If the contribution is new model (e.g., a large language model), then there should either be way to access this model for reproducing the results or a way re- produce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case au- thors are welcome to describe the particular way they provide for reproducibility. the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers have some path reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [Ys]]
The full implementation is provided as a Jupyter notebook in the supple- mentary material, including all code and data needed reproduce the main experimental results. The notebook contains clear documentation and step-by-step instructions, ensuring that the simulations and figures presented in the paper can be faithfully reproduced. This supports the transparency and reproducibility of the empirical claims made in the work. See /ttpps://gii //githbb.. com http//ayybbeeacc//BBI...
Justification:
Guidelines:
The answer NA means that paper does not include experiments requiring code. Please see the NeurIPS code and data submission guidelines (https://nips) cc/
onPolicy) for more details.
While we encourage the release of code and data, understand that this might not be possible, so "No" an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https:
//nips. cttccc/ppuuco
Sume itsSSsmmssiiooooiiiii.ii. onPolicy) for more details.
The authors should provide instructions on data access and preparation, including how access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only subset of experiments reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code is permitted.
6. Experimental setting/details
Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Question:
Answer: [Yes]
Justification: For the synthetic experiments, all parameters, such as sample size, matrix dimensions, noise levels, and distributional assumptions, are fully specified. The estimators
33

---

are closed-fone and do not involve optimization or hyperparameter tuning. no additional training procedure required. For the real-data experiments on CIFAR-10, the paper uses the standard train/test split provided by the dataset's official loader, which is explicitly mentioned in the text. No additional tuning or fine-tuning is performed, and all relevant details are disclosed to ensure full understanding of the results.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to level of detail that is necessary appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, or as supplemental material.
7. Experiment statistical significance
Does the paper report error bars suitably and correctly defined or other appropri- ate information about the statistical significance of the experiments?
Question:
Answer: [Yes]
Justification: For the simulated data, the paper evaluates performance across a wide range of parameter settings (e.g., sample size, matrix dimensions, noise levels), which helps mit- igate the effects of randomness and ensures that the empirical results robustly reflect the theoretical predictions. This systematic variation serves form of sensitivity analysis. For the real-world CIFAR-10 experiments, the paper explicitly reports error bars to convey the variability of the estimators and support the statistical significance of the results. Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper.
The factors of variability that the error bars capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall
run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).
It should be clear whether the error bar is the standard deviation the standard error of the mean.
It is OK to report I-sigma error bars, but one should state it. The authors should prefer- ably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not show in tables figures symmetric error bars that would yield results that are out of range (e.g. negative error rates).
If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
The experiments consist primarily of synthetic simulations and closed-fome estimators, which are computationally lightweight and do not require specialized hardware. No GPU is needed, and all experiments can be run standard laptop CPU-based ma- chine. While exact runtimes are not reported, the simplicity and efficiency the estimators ensure that the results reproducible without significant computational resources.
Justification:
Guidelines:
34

---

The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs as well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g.. preliminary failed experiments that didn't 't make it into the paper).
9. Code ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics htt.iiggtttte//rrrp.pp.pp///aaaa inne
Answer: [Yes]
Justification: The research adheres fully to the NeurIPS Code of Ethics, The work is the- oretical and empirical nature, with experiments conducted on synthetic data and the publicly available CIFAR-10 10 dataset. No personal, sensitive, proprietary data is used. All methods are described transparently, reproducibility supported through open supple- mentary materials, and no foreseeable negative societal impacts or misuse of the research have been identified.
Guidelines:
The answer NA means that the authors have not reviewed the NeurlPS Code of Ethics.
the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., if there is special consid- eration due to laws or regulations their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
Justification: No dedicated section is included in the paper to discuss societal impacts. However, the work is theoretical in nature and focuses on a methodological contribution to matrix-valued regression. As such, it does not directly target any high-risk applica- tion domains. While the methods may have positive downstream impact in areas such as spatiotemporal modeling image analysis, no immediate or foreseeable negative societal impacts are associated with the research.
Guidelines:
The answer NA means that there no societal impact of the work performed.
If the authors answer NA or No. they should explain why their work has no socictal impact why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious unintended (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (c.g., deployment of technologies that could make decisions that unfairly impact spe- cific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
35

---

If there are negative societal impacts, the authors could also discuss possible mitiga- tion strategies (e.g., gated release of models, providing defenses in addition attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility of MI.).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification:
The paper does not involve the release of high-risk models or datasets. It fo- cuses theoretical analysis and simulation-based evaluation, along with experiments on the publicly available CIFAR-10 dataset. No pretrained models, sensitive data, or genera- tive systems used.
Guidelines:
The answer NA means that the paper
such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by re- quiring that users adhere usage guidelines or restrictions access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [Yes]
Justification: All external assets used in the paper are publicly available and properly cred- ited. The CIFAR-10 dataset cited appropriately and used in accordance with terms of Any third-party libraries or tools used for simulations TensorFlow, scikit-learn) follow open-source licenses and are acknowledged either directly in the text within the supplementary notebook. No proprietary restricted-use resources are employed.
experiments (e.g., NumPy,
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package or dataset.
The authors should state which version of the asset is used and, if possible, include URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in the pack- age should be provided. For popular datasets, tattasets.ppperswiiit snperswwthcode com/ datasets has curated licenses for some datasets. Their licensing guide can help determine the li- cense of dataset.
For existing datasets that are re- re-packaged, both the original license and the license of the derived asset (if has changed) should be provided.
If this information not available online, the authors are encouraged to reach out to the asset's s creators.
13. New assets
36

---

Question: Are new assets introduced in the paper well documented and is the documenta- tion provided alongside the assets?
Answer: [Yes]
Justification: The paper introduces new simulation code and experiments, which are pro- vided in a well-documented Jupyter notebook included the supplementary material. The code includes clear instructions, parameter settings, and explanations necessary to repro- duce the results. No new datasets or models involving human subjects or requiring consent are introduced. All assets anonymized to comply with the double-blind review policy.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details the dataset/code/mooell as part their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and research with human subjects
For crowdsourcing experiments and research with human subjects, does the pa- per include the full text of instructions given participants and screenshots, if applicable, as well as details about compensation (if any)?
Question:
Answer: [NA]
The paper does not involve crowdsourcing or research with human subjects. All experiments conducted on synthetic data and the publicly available CIFAR-10 dataset, which does not contain personally identifiable information or require participant interaction.
Justification:
Guidelines:
The answer NA means that the paper does not involve crowdsourcing research with human subjects.
Including this information in the supplemental material is fine, but if the main contri- bution of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According to the NeurlPS Code of Ethics, workers involved in data collection, cura- tion, or other labor should paid at least the minimum wage the country of data collector.
15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: The paper does not involve human subjects or study participants. All ex- periments are conducted using synthetic data or the CIFAR-10 dataset, which is publicly available and does not involve any human interaction or identifiable information. There- fore, no IRB approval was required.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing research with human subjects.
Depending on the country in which research is conducted, IRB approval (or equiva- lent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
37

---

We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
16. Declaration of LLM usage
Does the paper describe the usage of LLMs if it is an important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing, editing, formatting purposes and does not impact the core methodology. scientific rigorousness, or originality of the research, declaration is not required.
Question:
Answer: [NA]
The core methods, theoretical analysis, and experimental components of the paper do not involve the use of large language models (LLMs) in any important, original, or non-standard way. Any use of LLMs was limited to editing support and does not affect the scientific contributions of the work.
Justification:
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, non-standard components.
Please refer to our LLM policy (attps://neurips.cc//oo http:/pps.c///oofeeenneeeee//e/s 2005/11 for what should or should not be described.
38


Reviews:
[{"review_id": "WJ6s9F6E40", "reviewer": "Reviewer_ivfu", "summary": "This paper studies the problem of bivariate matrix-values linear regression (BMLR). The authors first studies the explicit solution in the noiseless case, then subsequently proposed an explicit, optimization free estimators and established finite sample convergence analysis. The paper also extends the estimator to exploit the sparsity structure and achieve better convergence rate. Empirical evaluation demonstrates the practical value of the proposed estimator.", "questions": "I'm mainly curious about the importance of various assumptions. See \"Weaknesses\".", "limitations": "Yes.", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strength\n\n1. This paper addresses a significant gap in multi-dimensional regression by tackling the fully bilinear matrix-to-matrix regression problem (BMLR) rather than vectorized or low-rank approximations. In particular, the authors proposed optimization-free, closed from estimators, which is a very different approach from previous iterative algorithms. This provides great practical value as people don't need to tweak the algorithm for any specific problem.\n\n2. Strong theoretical guarantee. This paper presents thorough theoretical analysis for the proposed estimator, and provide finite sample convergence results. This is the first work that provide both non-asymptotic analysis and closed-form solution for this non-convex problem. It is also interesting that the convergence rates can be further improved by using hard-thresholding to exploit sparse structure.\n\n3. The presentation of this paper is pretty clear. The problem is clearly formulated, and the intuition are well presented for the noiseless case solution, the general estimators construction, and for the convergence results.\n\nWeaknesses\n\n1. More discussion on assumptions. While the overall results are reasonable, I think readers will benefit from more discussion on various assumptions. For instance (a) the problem assumes the matrix A^* to be non-negative. This assumption doesn't seem very natural to me, and I didn't find thorough discussion on this assumption. It would be good if the authors could clarify more on this non-negative assumption; (b) the noise is assumed to be Gaussian, but does the same analysis also cover other sub-gaussian noise distribution? (c) the assumption relies on the ORT assumption which also seems strong -- can we relax the assumption to be just full-rank? Is there more explicit dependency on the condition number of the X^T X matrix? I think it would be interesting to see more discussions along this line.\n\n2. Limited empirical comparison. The empirical results provide evidence that the proposed estimator solves the problem, but it doesn't provide insights on the relative performance to other methods proposed for the similar problem. The CIFAR-10 corruption & recovery experiment also seems a bit artificial. Is there other real-world dataset that the BMLR applies more naturally? I think that would also help to clarify the practical significance.", "was_revised": true, "final_justification": "The author's response addressed my questions. I don't have further concerns. I think this paper makes good theoretical contribution and I've raised my rating to 5.", "rebuttal": "We thank the reviewer for the positive assessment of our paper's clarity, novelty, and theoretical contributions. We appreciate the constructive feedback and provide detailed responses below to the concerns regarding assumptions and empirical validation.\n\n***\"The problem assumes the matrix A* to be non-negative. This assumption doesn't seem very natural... It would be good if the authors could clarify more on this.\"***\n\nThe non-negativity constraint on $A^\\*$, combined with row-wise $\\ell_1$-normalization, serves a crucial role in ensuring **identifiability** of the BMLR model. Without such constraints, the bilinear form $A^* X_t B^*$ exhibits fundamental scaling ambiguities: any transformation $(A^\\*, B^\\*) \\to (\\alpha A^\\*, \\alpha^{-1} B^\\*)$ for $\\alpha > 0$ produces identical observations, making individual recovery of $A^\\*$ and $B^\\*$ impossible (see line 35).\n\n**Why non-negativity specifically?**\n\n1. **Interpretability**: In many applications, $A^\\*$ represents mixing weights, attention coefficients, or interaction strengths that are naturally non-negative (e.g., in image processing, spatial modeling, and neural attention mechanisms).\n\n2. **Theoretical tractability**: Non-negativity combined with row-normalization naturally bounds entries to $[0,1]$, enabling our projection-based estimator $\\hat{A}$ in Equation (7) and simplifying the finite-sample analysis.\n\n3. **Closed-form feasibility**: Alternative identifiability constraints (e.g., orthogonality or spectral normalization) would likely require iterative procedures, defeating our main contribution of optimization-free estimation.\n\n**Important note**: As we show in **Weakness 1 response to Reviewer 1juu, our empirical results remain robust when $A^\\*$ is generated from distributions extending beyond Uniform on $[0,1]$ before normalization, suggesting the constraint is not artificially restrictive in practice.\n\n***\"The noise is assumed to be Gaussian, but does the same analysis also cover other sub-Gaussian distributions?\"***\n\n**Yes, our analysis extends to general sub-Gaussian noise.** The key insight is that our concentration inequalities (Theorems 3.3 and 3.6) rely primarily on **tail bounds** rather than specific distributional properties.\n\n**Technical extension**: For sub-Gaussian noise with parameter $\\sigma^2$, we can replace Gaussian tail bounds (Mill's inequality, Theorem E.1) with standard sub-Gaussian concentration inequalities. The $T^{-1/2}$ scaling component and the qualitative dimensional dependencies remain the same under sub-Gaussian noise extensions, with the overall bounds differing by factors involving the sub-Gaussian parameters that may depend on problem dimensions.\n\n**Practical significance**: This extension covers many realistic noise models including bounded noise, Laplace distributions, and mixtures with sub-Gaussian tails, significantly broadening the applicability of our results.\n\nAs noted in **Appendix B**, we focus on Gaussian analysis for clarity, but the techniques naturally generalize to this broader class.\n\n***\"The assumption relies on the ORT assumption which also seems strong — can we relax it to full-rank? Is there explicit dependence on the condition number?\"***\n\nThe ORT assumption ($X^{\\top}X/T = I_{mq}$) is indeed strong and primarily serves to **simplify exposition** and obtain clean rate expressions. \n\n**Relaxation to full-rank case**: As detailed in **Appendix B (lines 514-522)**, our analysis extends to the general full-rank setting $X^{\\top}X \\succ 0$ using concentration inequalities from **Vershynin (2010, Section 5.5)**. In this setting:\n\n1. **Condition number dependence**: The bounds would explicitly involve $\\kappa(X^{\\top}X)$, the condition number, appearing as multiplicative factors in the concentration inequalities.\n\n2. **Rate preservation**: The $T^{-1/2}$ scaling component of the error bounds and the overall dimensional structure remain unchanged, but with multiplicative factors involving the condition number $\\kappa(X^{\\top}X)$.\n\n3. **Technical complexity**: The analysis becomes significantly more involved, requiring matrix concentration bounds for non-isotropic random matrices, which is why we present the ORT case in the main text.\n\n**Practical relevance**: In many controlled experimental settings and when design matrices are approximately well-conditioned, the ORT assumption provides valuable theoretical insight while remaining reasonably realistic.\n\n***\"The CIFAR-10 experiment seems a bit artificial. Are there other real-world datasets where BMLR applies more naturally?\"***\n\nWe appreciate this observation and acknowledge that our CIFAR-10 setup involves controlled perturbations. However, this design choice was **intentional**:\n\n**Why controlled experiments?**\n1. **Theory validation**: Our primary contribution is **theoretical** — providing the first closed-form estimators with non-asymptotic guarantees for a fundamentally non-convex problem. Controlled experiments allow precise validation of theoretical predictions (Figure 1 shows excellent agreement between empirical and theoretical convergence rates).\n\n2. **Ground truth availability**: Unlike purely observational datasets, our setup enables quantitative assessment of estimation accuracy, directly demonstrating that our estimators recover true parameters as predicted by theory.\n\n3. **Proof of concept**: The image correction application showcases practical feasibility while maintaining experimental rigor.\n\n**Regarding more natural applications**: \nWe agree that additional real-world applications would strengthen practical relevance. However, such evaluation would require:\n- Comparative analysis with existing iterative methods [16, 23]\n- Application-specific performance metrics beyond parameter recovery\n- Domain expertise for proper evaluation\n\nSince our contribution focuses on **methodological innovation** (first optimization-free approach) rather than empirical benchmarking, we designed experiments to validate theoretical claims rather than demonstrate superior performance across applications.\n\n**Future work**: The reviewer's suggestion motivates valuable future research directions, including applications to spatiotemporal data, multi-view learning, and tensor completion where BMLR structure arises naturally.\n\n***Conclusion***\n\nWe thank the reviewer for the thoughtful feedback. Our responses demonstrate that: (1) the modeling assumptions, while restrictive, are well-motivated and can be substantially relaxed, (2) our theoretical framework extends beyond the presented special cases, and (3) our experimental design appropriately validates the core theoretical contributions. \n\nGiven the acknowledged strengths — particularly the novelty of optimization-free estimation and rigorous theoretical analysis — we believe this work represents a meaningful methodological advance for the matrix-valued regression community."}, {"review_id": "7WcWEgQ6uO", "reviewer": "Reviewer_1juu", "summary": "The paper deals with bivariate matrix valued linear regression problems ($Y_t = A^\\ast X_t B^\\ast + E_t$). There is a growing literature for this non-convex optimization problem and the paper makes new theoretical contributions on analyzing the non-asymptotic estimation error rates under ORT assumption and $\\ell_1$ normalized rows for $A^\\ast$. The paper develops theorems and proofs for optimization free estimators for $A^\\ast$ and $B^\\ast$ and analyses both sparse and non-sparse cases. In noiseless case, exact recovery is established implying identifiability. Some numerical simulations are also provided in support of the theoretical claims.", "questions": "(Q1) Lines 274-275: Does $\\beta^\\ast$ hide some dimensional dependency? I find it surprising that increasing $q$ does not degrade rate for $\\hat A$, but increasing $m$ does degrade rate for $\\hat B$.\n\n(Q2) Addressing (W3), (W4), (W5), (W6) would improve clarity.\n\n(Q3) Addressing (W1), (W2) would improve quality.", "limitations": "yes", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 3, "significance": 4, "originality": 4, "strengths_and_weaknesses": "Strengths\n\n(S1) The theoretical parts of the paper are of high quality and clarity. Proofs and assumptions are clear and correct. Intuitions are provided for several results which is helpful to better interpret the theretical results. The ORT assumption simplifies the interpretation and technical work needed for the proof, but hints on relaxing it are provided in the Appendices. It also recognizes that current probabilistic tools are inadequate for simultaneously relaxing both ORT and homoskedasticity.\n\n(S2) Highly significant work. Optimization-free estimators are computationally efficient in high-dimensions and theoretical analysis of the corresponding estimation error rates is valuable. Error rates are developed not just for Frobenius norm, but also for operator norm and element-wise max norm on $\\hat B - B^\\ast$ and $\\hat A - A^\\ast$. Some error rates are shown to exhibit blessing of dimensionality effect.\n\n(S3) Highly original work. This seems to be the first work to report non-asymptotic error rates for optimization free estimators for the BMLR problem under the assumption that rows of $A^\\ast$ are $\\ell_1$ normalized. Relevant literature with theoretical results has been referenced and the paper's contributions have been positioned very well against prior art.\n\nWeaknesses\n\nAlthough the paper's contribution is theoretical analysis, I still find the numerical simulations and the commentary around it inadequate for a theory paper.\n\n(W1) Lines 285-286: Generating entries of $B^\\ast$ from Uniform[0,1) could become an implicit confounder to performance assessment, since Eq (7) utilizes this aspect of $A^\\ast$ (i.e. shrinkage to [0,1)) to improve $\\hat A$ estimation. It would be more reflective of the theory to simulate for a scenario where entries of $B^\\ast$ are not in [0,1) with at least constant probability. I suspect that if normalization criteria was imposed on $B^\\ast$ instead of $A^\\ast$, the error rate expressions both in theory and in simulation would look very different.\n\n(W2) Line 305: This comment is not adequately supported in the text. Since plots in Fig 1 are approximately linear in the depicted representation (i.e. logarithmic y-axis and linear x-axis), it should be straightforward to fit a linear model to the data in each plot and determine effective the slope/coefficients of the estimation error rates as informed by the simulation. These should then be compared with the theoretical predictions and the results should be reported. This would more explicitly quantify the sharpness of the analysis, expose where large gaps remain, and would support the intended claim.\n\n(W3) The plot in Fig 3 is a bit strange. I am not quite sure why $\\epsilon$ should be on the x-axis. Since BMLR model is for $Y_t = A^\\ast X_t B^\\ast + E_t$, I would expect some analog of signal-to-noise ratio (which should depend on the Frobenius norm of $E_t$) to be on the x-axis.\n\n(W4) In Fig 3, it would be instructive to add the analogous plot for the element-wise max norm $||\\cdot||_{+}$.\n\nThe paper could benefit from a little bit of reorganization.\n\n(W5) The theorem for the sparse Gaussian case in Appendix A merits to part of the main text in my opinion. Weaknesses with the simulation results presentation should be addressed. Although the writing of the overall paper is clear, it could be made a bit crisp to accommodate these additional suggestions for the main text. For example, statement of contributions is repeated near line 63 and section 1.3 and could be reduced in verbosity.\n\n(W6) From proof of Lemma C.1, I feel that $\\ell_1$ normalization of rows of $A^\\ast$ is critical for optimization-free estimators to admit simple expressions. If true, it would be good to highlight this in the main text to provide additional justification for the $\\ell_1$ normalized rows assumption.", "was_revised": true, "final_justification": "The authors have provided a clearly written rebuttal that addresses almost all questions raised in my review. (W1) and (W2) have been adequately addressed and this significantly improves the quality of the paper. (Q1) has been explained well and the authors have agreed to the reorganization suggestions in (W5) and (W6) for improved clarity. I have raised the quality, clarity and overall recommendation scores.", "rebuttal": "We thank the reviewer for their thorough and thoughtful assessment of our paper. We are encouraged by the recognition of the theoretical quality and clarity (S1), the significance (S2), and the originality (S3) of our work. We also appreciate the constructive feedback on improving the simulation design and presentation. Below, we address each of the reviewer’s comments and suggestions in detail.\n\n***(W1) On the use of Uniform[0,1) for generating $A^\\*$ and $B^\\*$***\n\nWe appreciate the reviewer’s observation. To test the robustness of our simulation results, we conducted additional experiments where entries of both $A^\\*$ and $B^\\*$ were independently drawn from a Uniform[0, c) distribution, with $c \\in$ {1, 2, 3, 5}, followed by row-wise $\\ell_1$-normalization of $A^\\*$ to ensure identifiability. The empirical convergence rates remained consistent across these values of $c$, confirming that the behavior observed with $c=1$ was not specific to the original choice. These plots will be added to the appendix of the revised version.\n\nRegarding the reviewer's suggestion to normalize $B^\\*$ instead of $A^\\*$, we confirm that this would indeed lead to different expressions in Proposition 2.2 and Corollary 2.4. In this alternate setup, the roles of $m$ (the number of columns in $A^\\*$) and $q$ (the number of rows in $B^\\*$) would be swapped in the analytical expressions. As a result, the dependence on these dimensions in the non-asymptotic error bounds of Theorems 3.3 and 3.6 would also change accordingly.\n\nThis structural asymmetry reflects the fact that identifiability in the BMLR model requires normalizing one of the two parameter matrices. Our choice to normalize $A^\\*$ is motivated by both interpretability (e.g., non-negative activation weights) and analytical tractability, but we acknowledge that alternative normalizations are possible. We will clarify this point in the revised version.\n\n***(W2) Quantifying empirical slopes vs. theory***\n\nWe thank the reviewer for this excellent suggestion. Following their recommendation, we have fitted linear models to the log-linear data in each plot to extract the empirical slopes and compare them quantitatively with our theoretical predictions.\n\n**Methodology**: For parameter T, we used log-log regression (since both axes are effectively on log scale), while for parameters n, m, p, q, we used linear-log regression (y-axis on log scale). The empirical analysis is made for the three norms: Frobenius, maximum Element-wise, Operator. From our concentration inequalities, we derived the expected theoretical slopes by analyzing the dominant scaling terms in each bound.\n\n**T-dependence validation**: All bounds contain $T$ in the exponent terms (e.g., $\\exp⁡(−T× \\epsilon^2)$) and balancing these exponents to obtain meaningful bounds yields the characteristic $T^{−1/2}$ scaling. This reflects the standard statistical principle that estimation error decreases as $1/\\sqrt{T}$ with $T$ independent observations. Our empirical results confirm this almost exactly, with slopes ranging from -0.538 to -0.569 compared to the theoretical -0.5, demonstrating that our estimators achieve the expected rate of convergence.\n\n**Parameter dependencies**: The empirical slopes match qualitatively the theoretical predictions. For example, for matrix $B$ Frobenius error norm with typical parameter values of 55 (parameter range is [10, 100]):  \n- n: theoretical -0.009, empirical -0.011  \n- m: theoretical +0.009, empirical +0.016  \n- p: theoretical +0.009, empirical +0.011  \n- q: theoretical +0.009, empirical +0.017  \n\n**Key findings**:\n1. **Signs match perfectly** across all parameters and norms\n2. **R² values > 0.85** for all fits\n\n**Norm-specific behavior**: The numerical analysis reveals important distinctions between norms, as theoretically expected:\n- Frobenius norm shows strong $p,q$ dependencies as predicted by the $pq$ scaling for $B$.  \n- Operator norm shows weaker $p,q$ dependencies due to the predicted ($\\sqrt p + \\sqrt q$) structure for $B$.  \n- Max norm shows intermediate behavior.\n\n**Asymmetric effects**: While theory often predicts similar behavior for $p$ and $q$, empirical results consistently show $q$ having larger slopes than $p$ for $B$ and the opposite for $A$, suggesting structural asymmetries not fully captured by the main bound terms.\n\n**Conclusion**: This analysis confirms that our theoretical rates align closely with the observed trends. In all tested regimes, the signs of the empirical slopes match theoretical expectations, and the rates scale correctly with respect to dimensional parameters. The max-norm, Frobenius, and operator-norm behaviors differ in expected ways, validating the distinction made in our analysis. These findings demonstrate that our theoretical bounds, though conservative, meaningfully capture the dominant trends in estimator performance.\n\n***(W3) Plotting $\\epsilon$ on the x-axis in Figure 3***\n\nWe acknowledge this could be confusing. Here, $\\epsilon$ controls the perturbation strength applied to the identity in defining $A^\\*$ and $B\\^*$ in the CIFAR experiment, and thus indirectly controls noise severity. That said, we agree that rephrasing the x-axis as an effective signal-to-noise ratio (SNR) or plotting performance versus $\\| E_t \\|_F / \\|A^\\ast X_t B^\\ast \\|_F $ would better reflect model structure. We will revise this plot in the final version.\n\n***(W4) Adding max-norm plots in Figure 3***\n\nWe will include the element-wise max norm in Figure 3 to better align with the norms analyzed in our theory.\n\n***(W5) Placement of sparse case (Appendix A)***\n\nWe agree that the results in Appendix A on sparse recovery are important. We will promote them to the main paper (Section 3.2) in the camera-ready version. This will streamline the narrative and clarify the contribution on sparsity-aware estimators.    \n⁡\n\n***(W6) On the critical role of row-wise $\\ell_1$-normalization***\n\nIndeed, the row-wise $\\ell_1$-normalization of $A^\\*$  is instrumental in deriving the closed-form expressions (e.g., Proposition 2.2) and in ensuring identifiability without iterative scaling. We will explicitly emphasize this connection in Section 2 and the introduction to further justify this modeling choice.\n\n***(Q1) Dimensional dependencies in error rates***\n\nThe parameter $\\beta^\\*$ is defined as:\n$$\\beta^\\* := \\frac{1}{pq} \\sum_{j=1}^p \\sum_{l=1}^q [B^\\*]_{lj}.$$\n\nThis means $\\beta^\\*$ represents the average entry value of matrix $B^\\*$. Crucially, $\\beta^\\*$ does contain implicit dimensional dependencies through the averaging process over the $pq$ entries of $B^\\*$.\n\n**Why increasing $q$ helps $\\hat{A}$ but hurts $\\hat{B}$:**\n\n1. For $\\hat{A}$: The bound shows that larger $|\\beta^\\*|$ improves the convergence rate (appears in denominator as $2\\epsilon/|\\beta^\\*|$). When $q$ increases:\n  - More columns are added to $B^\\*$\n  - Increasing $q$ thus provides more \"constraints\" or \"equations\" to solve for each entry of $A^\\*$\n  - The key benefit comes from having more observations (effectively #$\\mathcal{D}_0$ terms) to estimate $A^\\*$, creating a \"blessing of dimensionality\" effect (see the expression of $A^\\*$ in Corollary 2.4). \n\n2. For $\\hat{B}$: The Frobenius bound shows degradation with increasing $q$ through the $+\\frac{2pq}{3}$ term in the exponent, which directly reflects that:\n\n  - $B^*$ has more entries to be estimated\n  - More parameters generally require more data to estimate accurately\n  - This is the classical \"curse of dimensionality\" in parameter estimation\n\n**Why increasing $m$ hurts both $\\hat{A}$ and $\\hat{B}$:**\n\n\n1. For $\\hat{A}$: \n- Matrix $A^\\*$ has dimensions $n \\times m$, so increasing $m$ directly increases the number of parameters in $A^\\*$\n- Each row of $A^\\*$ must sum to 1 (identifiability constraint), creating coupling across all $m$ columns\n- This constraint means estimating any entry $[A^\\*]_{ik}$ requires information about all other entries in that row, making the problem harder as $m$ grows\n\n2. For $\\hat{B}$:\n - When $m$ increases, the entries $[\\hat{B}]_{lj}$ involve a sum over more indices (see the expression of $B^\\*$ in Corollary 2.4). \n- Each $[\\hat{C}]_{(k,l),(i,j)}$ has variance $\\sigma^2/T$, thus the sum is over $nm$ terms, but only $n$ of them are independent for each $(l,j)$ pair and the normalization only involves $n$.\n- Larger $m$ effectively dilutes the signal relative to noise in the estimation of each $[\\hat{B}]_{lj}$\n\n**Clarification for the paper**:\n\nWe should emphasize that while $\\beta^\\*$ is normalized by $pq$, it still captures the \"signal strength\" of matrix $B^\\*$, and in practice, this signal strength interacts with dimensionality in ways that can benefit $A^\\*$ estimation while simultaneously making $B^\\*$ estimation more challenging due to the increased parameter space.\n\nThis explains the seemingly counterintuitive asymmetric dimensional effects: the blessing/curse of dimensionality manifests differently for each matrix due to their distinct roles in the bilinear model structure and the identifiability constraints imposed.\n\n***Conclusion***\n\nWe thank the reviewer again for their constructive feedback, which has helped us significantly improve both the clarity and depth of our empirical and theoretical presentation. In particular, the slope-based analysis strongly supports our theoretical predictions, and the requested presentation updates (e.g., improved figures, appendix material, and reorganization) will be reflected in the revised version. We hope these clarifications further strengthen the case for this work's contribution to the theoretical understanding of BMLR."}, {"review_id": "75vYtZA0iD", "reviewer": "Reviewer_4ujm", "summary": "The paper presents closed-form estimators for the BMLR model that requires no iterative optimization. The paper also provides the associated non-asymptotic error bounds of the estimator. This work claims it is the first to offer a complete finite-sample analysis for estimating both A* and B* individually. Empirical results on synthetic and real-image datasets confirm the method can somehow estimate the A and B up to some accuracy.", "questions": "I reiterate the weakness here, please address the questions below. I appreciate it if the questions may be clarified here. Thanks. \n\nWeakness 1: Presentation. The connection between that question on page 1 and the actual problem solved in Equation 2 is unclear. More clarification here is needed. \n\nWeakness 2: Clarity. In line 162, what is matrix E? Is it the noise? Since the noise is unknown. What is the implication of having C constructed using unknown A, B, E? Why the estimator here of A and B in Equation 6 has unknown parameters C? Is it solved in iterative manner? Why entries in matrix in A are constrained to lie between 0 and 1? These derivations and the reasoning behind them are not clearly motivated and explained in the paper. \n\n\nWeakness 3: Experiments. The paper is lacking experiment comparisons with other denoising methods in the literature. Actually, it is not compared with any other method. Also, in Figure 2, it seems the denoising is not very successful where the reconstruction is still noisy. \n\nWeakness 4: (Minor) Instead of “distance” in the y axis, usually, paper would describe it as mean squared error (MSE). \n\nWeakness 5: The sparse settings is in the appendix. I personally recommend that this part should be included in the main paper, while many of the derivation details can be instead put in the appendix, leaving the main conclusive theories and results in the main paper.", "limitations": "Yes.", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strength:\n\nThe paper provides the associated non-asymptotic error bounds of the estimator. This work offers a finite-sample analysis for estimating both A* and B* individually. Empirical results on synthetic and real-image datasets confirm the method can somehow estimate the unknown A and B up to some accuracy.\n\nWeakness: I appreciate it if the weakness below may be clarified during rebuttal. Thanks. \n\nWeakness 1: Presentation. The presentation of the paper can be improved. The current version is hard to follow in terms of the story flow. The paper opens with the equation on page 1. But the connection between that question on page 1 and the actual problem solved in Equation 2 is unclear. More clarification here is needed. \n\nWeakness 2: Clarity. In line 162, what is matrix E? Is it the noise? Since the noise is unknown, what is the implication of having C constructed using unknown A, B, E? Why the estimator here of A and B in Equation 6 has unknown parameters C? Is it solved in iterative manner? Why entries in matrix in A are constrained to lie between 0 and 1? These derivations and the reasoning behind them are not clearly motivated and explained in the paper. \n\n\nWeakness 3: Experiments. The paper is lacking experiment comparisons with other denoising methods in the literature. Actually, it is not compared with any other method. Also, in Figure 2, it seems the denoising is not very successful where the reconstruction is still noisy. \n\nWeakness 4 (minor): Instead of “distance” in the y axis, the error here is described as mean squared error (MSE). \n\nWeakness 5: The sparse settings is in the appendix. I personally recommend that this part should be included in the main paper, while many of the derivation details can be instead put in the appendix, leaving the main conclusive theories and results in the main paper.", "was_revised": true, "final_justification": "I read the rebuttal from the authors, and the discussion has clarified most of my concerns. I am looking forward to more experiments in the camera ready version.", "rebuttal": "We thank the reviewer for their careful assessment and detailed feedback. We are grateful for the opportunity to clarify the presentation and address the technical questions raised. Below, we respond point by point to each concern.\n\n**Weakness 1: Presentation and Story Flow**\n\nWe agree that the connection between the introductory equation and the core estimation problem can be made clearer. Here is the intended narrative.\n\n**The progression from general to specific:**\n1. **Page 1 motivation**: We begin with the general BMLR model $Y_t = A^\\* X_t B^\\* + E_t$ to establish the fundamental problem of estimating matrix parameters $A^\\*$ and $B^\\*$ from matrix-valued observations.\n\n2. **Section 2 (Equation 2)**: We first analyze the **noiseless oracle case** $M_t = A^\\* X_t B^\\*$ to establish fundamental identifiability properties and derive the theoretical foundation for our estimators.\n\n3. **Section 3**: We then return to the **practical noisy case** from page 1, using insights from the oracle analysis to construct optimization-free estimators.\n\n**Why this progression matters**: The oracle analysis (Section 2) is crucial because it reveals the exact mathematical relationships between $A^\\*$, $B^\\*$, and the observations, which we then exploit in Section 3 to design plug-in estimators that simply replace unknown quantities with empirical counterparts.\n\nWe will revise the introduction to more clearly explain this progression and include a roadmap paragraph outlining how the oracle case enables our estimation procedure in the presence of noise.\n\n**Weakness 2: Clarity of Technical Details**\n\nWe address each specific question:\n\n**Q: What is matrix $E$ in line 162?**\n\nThe matrix $E_{(k,l)}$ denotes a canonical basis matrix: it has a 1 at entry $(k,l)$ and zeros elsewhere. We agree that this notation was not clearly explained and will make it precise in the revised manuscript.\n\n**Q: What is the implication of $C$ being constructed using unknown $A$, $B$, $E$?**\n\nIn the oracle setting (Section 2), matrix $C$ is defined as $C = (X^{\\top}X)^{-1}X^{\\top}M$ and is observable under the noiseless model $M = A^\\* X B^\\*$. Corollary 2.4 shows how $A^\\*$ and $B^\\*$  can be exactly recovered from this quantity. In the noisy setting, we compute the empirical counterpart $\\hat{C} = (X^{\\top}X)^{-1}X^{\\top}Y$ which differs from $C$ by a noise term. Our estimators in Equation (6) are built directly from $\\hat{C}$ requiring no knowledge of $A^\\*$, $B^\\*$ or $E_t$. \n\n**Q: Why does the estimator in Equation (6) have unknown parameters?**\n\nThis is a misunderstanding. In Equation (6), $\\hat{B}$ and $\\hat{A}$ are defined using $\\hat{C}$, which is which is computed directly from the data $(X_t, Y_t)$ via Equation (5). The estimators are fully explicit and require no iterative procedures.\n\n**Q: Why are entries of $A^\\*$ constrained to $[0,1]$?**\n\nThis follows from two assumptions stated on lines 34-36:\n1. $A^\\*$ has non-negative entries (modeling assumption)\n2. Each row sums to 1 (identifiability constraint).\n\nThese two constraints naturally bound entries to $[0,1]$. Our estimator $\\hat{A}$ defined in Equation (7) respects this structure via row-wise projection, improving finite-sample performance, as stated on lines 211-214.\n\n**Weakness 3: Experimental Comparisons and Denoising Quality**\n\nWe acknowledge that our empirical evaluation does not benchmark against classical denoising methods. This is intentional: our primary goal is to demonstrate the theoretical soundness of our closed-form estimators for the BMLR model, not to propose a general-purpose image denoising technique.\n\nFurthermore, our method targets a specific matrix regression setting, where the corruption is modeled as bilinear noise over structured matrix inputs — quite distinct from the assumptions behind traditional denoisers (e.g., wavelet shrinkage, CNN-based filters). Thus, direct comparisons are not straightforward or meaningful.\n\nRegarding Figure 2, the reconstructed images visually appear noisy because the transformations are applied on compressed representations. The apparent \"noise\" in the corrected image reflects two factors:\n1. We use a relatively high noise level ($\\epsilon = 0.02$) to demonstrate robustness\n2. Some residual error is expected - our method provides **statistical estimation** with finite-sample guarantees, not perfect reconstruction\n\nDespite visible noise, Figure 3 shows a clear and consistent reduction in MSE, demonstrating denoising effectiveness in the intended setting. To address this concern further, we will:\n\n- Add visual results at different noise levels\n- Add a short discussion on the statistical limits of recovery\n\n**Weakness 4: Terminology in Figure Axis**\n\nAgreed. We will change \"distance\" to \"Mean Squared Error (MSE)\" for clarity and consistency with standard terminology.\n\n**Weakness 5: Sparse Results Placement**\n\nWe agree. The sparse recovery results (Appendix A) represent significant theoretical contributions and merit inclusion in the main text. We will:\n1. **Move Theorems A.1 and A.2** to a new Section 3.3 in the main paper\n2. **Relocate detailed proofs** in the appendix\n3. **Streamline other sections** to accommodate this addition without exceeding page limits\n\n**Additional Improvements Based on Your Feedback:**\n\n1. **Enhanced Introduction**: Clear problem motivation and solution roadmap\n2. **Methodology Overview**: Step-by-step explanation of our estimation approach\n3. **Improved Figures**: Clearer labeling and additional examples\n4. **Reorganized Structure**: Main theoretical results in the body, technical details in appendices\n\n**Conclusion:**\n\nWe believe these clarifications and improvements directly address your concerns about presentation clarity, technical exposition, and experimental evaluation. The theoretical contributions remain novel and significant, and with improved presentation, we hope the paper will be accessible to a broader audience while maintaining its theoretical rigor.\n\nWe thank you for the constructive feedback and look forward to submitting a substantially improved revision."}, {"review_id": "gonk2rhQOz", "reviewer": "Reviewer_Bwup", "summary": "This paper studies a bilinear matrix-valued regression model. \nThey study an oracle case and establish that, in the absence of noise, the true parameters can be exactly retrieved. This analysis highlights the fundamental identifiability properties of the model.\nThey propose explicit, optimization-free estimators. \nIn theoretical analysis, non-asymptotic bounds are presented to characterize the dependence of estimation accuracy on the problem dimensions (n, p,m, q) and the sample size T.", "questions": "I was surprised that a similar result hasn't appeared in the literature. I choose to believe that the authors have done a thorough literature search and it's the first appearance of these results.", "limitations": "No separated paragraph in the paper that discussed limitations of their approach.", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "The paper is well written and easy to follow. It's a pleasure to read. \n\nThe novelty of the result is not very clear. The main estimate is based on a typical least square solution in (5) with additional linear transformation in (6) and (7). These expression is so similar to the ones in the existing linear regression setting, the related results seem to be easily reachable by applying standard concentration inequalities.", "was_revised": true, "final_justification": "See my responses to authors. No additional comments.", "rebuttal": "We thank the reviewer for their positive evaluation of our paper’s clarity and technical soundness. We appreciate the comment that it is “well written and easy to follow,” and we respond below to the concerns regarding the novelty and limitations of our contributions.\n\n**On the Novelty and Originality of Our Results**\n\nWhile the final estimators may resemble classical least-squares formulas, their application to the bilinear matrix-valued regression (BMLR) setting is non-trivial, and we respectfully highlight several important points of novelty:\n\n**1. Structural Complexity Beyond Standard Linear Regression:**\n\nThe BMLR model $Y_t = A^* X_t B^* + E_t$ introduces unique challenges not present in classical regression:\n- **Non-convexity**: The bilinear dependence on $A^\\*$ and $B^\\*$ yields a non-convex estimation problem.\n- **Identifiability constraints**: The need for row-wise $\\ell_1$-normalization of $A^*$, necessary to resolve scaling ambiguities, creates coupling across parameters.\n- **Matrix-valued observations**: Our approach maintains and leverages the full matrix nature of the inputs and outputs, unlike vectorized or rank-constrained alternatives.\n\n**2. Technical Innovation in Estimation:**\n\nWhile Equation (5) resembles a least-squares solution, the downstream estimators (Equations (6) and (7)) involve several non-standard techniques:\n- **Ratio estimator** in Equation (6) for $\\hat{A}$, where we avoid the complex distributional analysis of ratios of Gaussian variables by using the projection approach in Equation (7).\n- **Sample-splitting strategy** in the definition of $[\\tilde{B}]^{(i)}_{lj}$ (line 198) to maintain statistical independence (detailed on line 199).\n- **Matrix structure preservation** that enables exploitation of dimensional effects (blessing of dimensionality for some parameters, curse for others). The estimators are explicit and optimization-free, yet tailored to the bilinear structure — a setting where most prior work requires iterative procedures.\n\n**3. Non-Asymptotic Analysis with Dimension-Dependent Rates**\n\nThe theoretical analysis goes well beyond standard concentration inequalities:\n- **Multi-norm bounds**: We provide simultaneous analysis for elementwise max, Frobenius, and operator norms, each requiring different technical approaches.\n- **Asymmetric dimensional dependencies**: The analysis reveals counter-intuitive effects where increasing $q$ helps $\\hat{A}$ but hurts $\\hat{B}$, while $m$ degrades both estimators. Thus the analysis shows both blessing and curse of dimensionality effects.\n- **Sparse recovery guarantees**: Our results in Appendix A provide exact support recovery, using tools adapted to the matrix bilinear setting.\n\n**4. Literature Gap Confirmation:**\n\nWe understand the reviewer's surprise that these results have not appeared before. To the best of our knowledge:\n\n- **Existing BMLR work** focuses on iterative optimization methods [23, 16] without non-asymptotic guarantees\n- **Kronecker factorization approaches** [11] provide bounds on the Kronecker product $M^* = (B^*)^{\\top} \\otimes A^*$ but not on individual matrices $A^*$ and $B^*$.\n- **Matrix regression literature** [38, 28] addresses different settings (scalar responses or matrix predictors, but not both matrix-valued).\n\n**5. Methodological Contributions:**\n\n-  The first closed-form, optimization-free estimator for BMLR with finite-sample guarantees,\n- The first individual matrix recovery analysis in this setting,\n- And the first to show how dimension interactions affect estimator performance — e.g., how increasing $q$ helps $\\hat{A}$ but degrades $\\hat{B}$.\n\n**Addressing the Limitations Concern:**\n\nWhile Appendix B discusses the validity and potential relaxation of our key assumptions, we acknowledge that a more prominent and comprehensive limitations discussion would benefit readers. We will add a dedicated \"Limitations\" section covering:\n\n- **ORT assumption**: While this assumption facilitates clean theoretical analysis, it may not hold in practice. As we detail in Appendix B, our theoretical framework can be extended to more general settings where $X^{\\top}X/T \\neq I_{mq}$, though this requires more sophisticated concentration inequalities for matrices with non-isotropic columns and introduces additional complexity in the bounds.\n\n- **Noise structure assumptions**: Our main analysis assumes isotropic Gaussian noise ($\\sigma^2 I_{np}$). However, as shown in Appendix B, we can accommodate more general matrix normal distributions $\\mathcal{MN}_{n \\times p}(0, \\Sigma_r, \\Sigma_c)$ with structured row and column dependencies, albeit with increased technical complexity.\n\n- **Joint assumption relaxation**: As we acknowledge in Appendix B, simultaneously relaxing both the ORT assumption and homoskedasticity represents a significant theoretical challenge that current probabilistic tools cannot adequately address with sharp, tractable bounds.\n\n- **Identifiability constraints**: Our choice to impose $\\ell_1$-normalization on $A^*$ rather than $B^*$ is somewhat arbitrary and affects the dimensional dependencies in our bounds. Alternative normalization schemes would yield different (though related) theoretical guarantees.\n\n- **Computational scalability**: Despite being optimization-free, our method requires inverting $(X^{\\top}X)$ of size $mq \\times mq$, which becomes computationally expensive when the product $mq$ is large relative to available computational resources.\n\n- **Model specification**: The bilinear structure $Y_t = A^* X_t B^* + E_t$ may not capture all relevant dependencies in real applications, where higher-order interactions or non-linear relationships might be present.\n\nWe will incorporate this detailed limitations discussion to provide readers with a comprehensive understanding of our method's scope and applicability.\n\n**Why These Results Matter:**\n\nThe combination of:\n1. **Computational efficiency** (no iterative optimization)\n2. **Theoretical guarantees** (first non-asymptotic bounds for individual matrices)\n3. **Structural exploitation** (blessing of dimensionality effects)\n4. **Practical applicability** (demonstrated on real data)\n\nrepresents a significant advance over existing BMLR methodology. For instance, the counter-intuitive finding that increasing $q$ helps $\\hat{A}$ estimation but hurts $\\hat{B}$ estimation (with slopes of opposite signs in our empirical validation, see rebuttal for reviewer 1juu) reveals fundamental asymmetries in the bilinear structure that would not be apparent from standard regression analysis.\n\n**Conclusion:**\n\nWhile our estimators build on fundamental least squares principles, the technical challenges of the bilinear matrix setting, the sophisticated theoretical analysis required, and the novel insights about dimensional effects constitute substantial contributions to the matrix-valued regression literature. The apparent simplicity of our final estimators reflects successful mathematical engineering that distills complex bilinear relationships into tractable closed-form expressions.\nWe will enhance the paper with a deeper limitations section. We believe these clarifications and the planned improvements directly address the concerns raised and will strengthen the contribution to the matrix-valued regression literature. We thank the reviewer for the constructive feedback."}]

Author rebuttal:

Meta-review: The paper proposes an optimization-free approach for Bivariate Matrix-valued Linear Regression (BMLR), where both predictors and responses are matrices. The authors introduce explicit estimators and provide the first non-asymptotic error bounds for the individual recovery of both parameter matrices under identifiability and sparsity assumptions. The reviewers uniformly praised the paper's theoretical rigor, originality, and clarity. The authors' rebuttal successfully addressed all concerns, including the interpretation of assumptions, the lack of empirical comparisons, and the presentation of key results. Reviewers confirmed that the authors' responses were highly informative and resolved all initial questions, leading to a strong consensus for acceptance.
Decision: Accept (poster)