Solving Sparse \& High-Dimensional-Output Regression via Compression
Authors: Renyuan Li
Zhehui Chen
Guanyi Wang
Venue: NeurIPS 2024 poster (2024)
Keywords: multi-output regression
sparsity
compression
non-convex optimization
Primary area: learning_theory

Abstract: Multi-Output Regression (MOR) has been widely used in scientific data analysis for decision-making. Unlike traditional regression models, MOR aims to simultaneously predict multiple real-valued outputs given an input. However, the increasing dimensionality of the outputs poses significant challenges regarding interpretability and computational scalability for modern MOR applications. As a first step to address these challenges, this paper proposes a Sparse \& High-dimensional-Output REgression (SHORE) model by incorporating additional sparsity requirements to resolve the output interpretability, and then designs a computationally efficient two-stage optimization framework capable of solving SHORE with provable accuracy via compression on outputs. Theoretically, we show that the proposed framework is computationally scalable while maintaining the same order of training loss and prediction loss before-and-after compression under arbitrary or relatively weak sample set conditions. Empirically, numerical results further validate the theoretical findings, showcasing the efficiency and accuracy of the proposed framework.

Solving Sparse & High--Dimensionll-Ouuuttt Regression via Compression
Renyuan Li
Zhehui Chen Google
Department of Industrial Systems Engineering & Management National University of Singapore renyuan.. 1i li@u. mus.edu
zhehui chongugece@g
Guanyi Wang
Department Industrial Systems Engincering & Management National University of Singapore guanyi www.uuu.uuu...... edu.sg
Abstract
Multi-Outpur Regression (MOR) has been widely used in scientific data analysis for decision-making. Unlike traditional regression models, MOR aims to simultane- ously predict multiple real-valued outputs given input. However, the increasing dimensionality of the outputs poses significant challenges regarding interpretability and computational scalability for modern MOR applications. As a first step address these challenges, this paper proposes a Sparse & High-dimennsioal-Ouupe REgression (SHORE) model by incorporating additional sparsity requirements resolve the output interpretability, and then designs computationally efficient two- stage optimization framework capable of solving SHORE with provable accuracy via compression on outputs. Theoretically, we show that the proposed framework is computationally scalable while maintaining the same order of training loss and prediction loss before-and-after compression under arbitrary relatively weak sample set conditions. Empirically, numerical results further validate the theoretical findings, showcasing the efficiency and accuracy of the proposed framework.
1 Introduction
Multi-Q Multi-Output Regression (MOR) problem is a preponderant tool for factor prediction and decision decision-making making modern data analysis. Compared with traditional regression models that focus a scalar output for each sample, MOR aims to predict multiple outputs y e E simultaneously based given input E RRd, i.e.,
arg min dist(u,ū(z)) with 9 :- arg min wey
1
Σ
l(v -((x))
gEg i=1
where we use {xx }') to denote its given sample set with x E Rd i-th input feature vector and y E K corresponding output vector, define f: g x R + R as the loss function, dist R kg + E as some prediction/distance metric, 3 as the structure/constraint set for multiple outputs, and & as the candidate set for predicting model g gd k Hence, MOR and its variants have been used for numerous regression tasks with structure requirements outputs arising from real applications, such as simultaneous estimation of biophysical parameters from remote sensing images [40], channel estimation through the prediction of several received signals [35], the grounding (e.g., factuality check [66] in the Large Language Model (LLM, [34 era, name but few.
multi-dimennionall
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

In this paper, we are interested in the interpretability issue of high high-dimensional outputs obtained from modern MOR tasks. One typical example is raised from algorithmic trading. In particular. in algorithmic trading, MOR helps to construct the portfolio [31 from a large number of financial
instruments (e.g. different stocks, futures, options, equities, etc [99) based on given historical market and alternative data [22). To be concise, high-dimensional output in this example could be viewed as a "decision". where every component denotes the investment for the corresponding financial instruments. Thus, other than accuracy, quantitative researchers prefer outputs with only few instruments to enhance interpretability for the underlying decision-making reasons, which naturally introduces sparse output condition. Similar scenarios apply to other applications, including offline reinforcement learning in robotics 331, discovering genetic variations based on genetic markers (25).
As a result, the dramatic growth in output dimensions gives rise to two significant challenges: 1. gegh-dimensional-ourrae impedes human interpretation for decision-making 2. Approaches with better computational scalability are desired for training & predicting MOR. Upon these challenges, conceptual question that motivates this research is:
How design framework that predicts output with enhanced interpretability, better computational scalability, and provable accuracy under a modern gegh--imenniooal-utppe setting?
Generally speaking, this paper provides an affirmative answer first step the above question. Before presenting the main contributions, let us first introduce the model that will be studied in this paper. Unlike the classical MOR model, further assume that given outputs are of high-dimensional (i.e., d < K), and to address the interpretability issue, these outputs have at most 8 non-zero components, i.e., 1/y S 8, for all i [n]] with some pre-determined sparsity-level s<< K). Based on such given samples, this paper proposes the (uncompressedd Sparse & Hiigh-dimensional-nur REgression (SHORE) model that aims to predict an interpretable high- -ige-dimensional output y (i.e.. s-sparse in this paper) via any input feature vector æ. In particular, to be concise and still capture the essential relationship, the proposed (uncompressed) SHORE model predicts y from æ under a linear model, i.e., y = arg fyllos dist(y. Zx) for some distance metric (see Section 3.1 prediction stage) and the linear regression Ž is obtained by solving the following linear regression problem:
Ż = :arg min (Z) := Y ZXI
ZERRxé
where X := ===xx1111 I xx)) E dd is the input matrix and Y == (y² ... (1n) E Kxn xn is the corresponding column-sparse output matrix.
1.1 Contributions and Paper Organization
This paper makes the following three main contributions:
1. We propose a two-stage computationally efficient framework for solving SHORE model. Specif- ically, the first training stage offers computationally scalable reformulation on solving SHORE through compression in the output space. The second prediction stage then predicts high- high-dimensional outputs from a given input by solving a specific sparsity-conee constrained minimization problem via an efficient iterative algorithm.
2. We show that for arbitrarily given samples, the training loss in the first stage with compression is bounded by a + 5 multiplicative ratio of the training loss for the original one ( with some positive constant 6. Additionally, the proposed iterative algorithm in the second stage exhibits global geometric convergence within a neighborhood of the ground-truth output, with a radius proportional to the given sample's optimal training loss. Furthermore, if all samples are drawn from a light-tailed distribution, the generalization error bound and sample complexity remain in the same order for SHORE with output compression. This finding indicates that the proposed framework achieves improved computational efficiency while maintaining the same order of generalization error bounds statistically.
3. We conduct rich numerical experiments that validate the theoretical findings and demonstrate the efficiency and accuracy of the proposed framework on both synthetic and real-world datasets.
In summary, this paper studies the SHORE model through computational and statistical lenses and provides a computationally scalable framework with provable accuracy.

---

The paper is organized as follows: Section 2 reviews related literature; Section B presents our proposed framework and provides theoretical results on sample complexity and generalization error bounds; Section I compares the proposed method with existing baselines in a suite of numerical experiments on both synthetic and real instances. Concluding remarks are given in Section 5
Notation. Given a positive integer n, we denote [n] := (l......1)) nt... We use lowercase letters a as scalars and bold lowercase letters a as vectors, where a; its i-th component with i E [d], and bold upper case letters A as matrices. Without specific description, for a m-by-n matrix A, we denote A., as its (i. j-th component, A,,. as its i-th row, A., as its j-th column. For symmetric square matrix A, we denote Amax(A), Amin (A) and A₂ as its maximum, minimum and i-th largest eigenvalue, respectively. We denote ||a|]], |al|22 ||a se ||Al|
IF || as the L,, l2, -norm of vector a, the Frobenius norm and the operator norm of a matrix A, respectively. We denote I(-) as the indicator function, ||a||o Σ
-i I(a; F 0) as the fo-norm (i.e., the total number of nonzero components), supp(a) == {i E id] a, #0 the support set. We denote V.K := fyy E R ||y|llo I s} set of s-sparse vectors, Bp(c:p) := {y E x |yy - cl2 I p} as a closed ballll with center c and radius p. (u,, as Gaussian distribution with mean μr and covariance o².
For two sequences of non-negative reals (fn) 1221 and {gn 122,, we use Jn S In to indicate that there is a universal constant C' > 0 such that fn S Cgn for all 2 1. We use standard order notation fn = Ognn) indicate that fn S In and fn (g.)) ) to indicate that fr, < In In (1//T) for some universal constants T and c. Throughout, we use e, 6, T, c, C1 C2.... and C, C1 Cz... to denote universal positive constants, and their values may change from line to line without specific comments.
2 Literature Review
Multi- Multi-output regression (MOR) and its variants have been studied extensively over the past decades. In this section, we focus on existing works related to our computational and statistical results.
Computational part. Existing computational methods for solving MOR can be, in general, clas- sified into two categories [8), known as problem transformation methods and algorithm adaptation methods. Problem transformation methods (c.g.. Binary Relevance (BR), multi-target regressor stacking (MTRS) method (37], regression chains method (37) aim to transform MOR into multiple single-output regression problems. Thus, any ofate-of-thh-arr single-output regression algorithm can be applied, such as ridge regression [5]. regression trees [9]. and etc. However, these transfor- mation methods ignore the underlying structures/relations between outputs, which leads to higher computational complexities. In contrast, algorithm adaptation methods focus more on the underlying structures/relations between outputs. For instance, (36] investigates input component selection and shrinkage in multioutput linear regression; [ later couples linear regressions and quantile mapping and thus captures joint relationships among variables. However, the output dimension considered in these works is relatively small compared with modern applications, and their assumptions concerning low-dimensional structure of outputs are hard to verify. To overcome these shortages, we consider high-dimenssionll-owe output regression with only additional sparsity requirement on outputs.
Statistical part. There are numerous works concerning statistical properties of traditional or multi- output regressions. [18] gives sharp results on "out-of pout-of-sample" (random design) prediction error for the ordinary least squares estimator of traditional linear regression. 4 proposes an empirical risk minimization framework for large-scale multi-label learning with missing outputs and provides excess risk generalization error bounds with additional bounded constraints. (28) investigates the generalization performance of structured prediction learning and provides generalization error bounds on three different scenarios, i.c., Lipschitz continuity, smoothness, and space capacity condition. [27] designs efficient feature selection procedure for multiclass sparse linear classifiers (a special case for SHORE with sparsity-level s = 1), and proves that the proposed classifiers guarantee the minimax generalization error bounds in theory. A recent paper 42] studies transfer learning via multi-task representation learning, a special case in MOR, which proves statistically optimistic rates on the excess risk with regularity assumptions on the loss function and task diversity. In contrast with these works, our contributions concentrate on how generalization error bounds change before and after the compression under relatively weak conditions on the loss function and underlying distributions.
Specific results in MLC. MLC is an important and special case for MOR with {0, 1 output per dimension, i.c., y E 10, and thus, in this paragraph, we use labels to replace
3

---

outputs. Here, we focus on dimensionality reduction techniques on outputs, in particular, the compressed sensing and low-rank conditions on the output matrix Y. The idea of compressed sensing rises from signal processing, which maps the original iigh-dimensional output space into a smaller one while ensuring the restricted isometry property (RIP). To the best of our loowledge, the compressed sensing technique is first used in [17] to handle sparse expected output El y [r]. Later. [39 [2] propose Principle Label Space Transformation (PLST) and conditional PLST through singular value decomposition and canonical component analysis respectively. More recently, many new compression approaches have been proposed, such as robust bloom filter 13,. log time log space extreme classification [23], merged averaged classifiers via hashing [32,, etc. Additionally, computational efficiency and statistical generalization bounds can be further improved when the output matrix Y ensures a low-rank condition. Under such condition, [45] provides a general empirical risk minimization framework for solving MLC with missing labels. Compared with the above works, this paper studies MOR under sparse & high-dimeensioall-oupee additional correlation assumptions or low-rank assumptions for output space, and then provides a complete story through computational and statistical lens.
setting without
3
Main Results
3.1 Two-Stage Framework
This subsection presents a general framework for solving SHORE and then the computational complexity for the proposed framework with/vithout compression. Given set of training samples {(x', y y') as described in Section [ the framework can separated into two stages: (compressed) training stage & (compressed) prediction stage.
Training stage. In the first training stage, the framework finds compressed regressor by solving a linear regression problem with compressed outputs. In particular, the framework compresses the original large output space (K dim) to a smaller "latent" output space (m-dimm by left-multilyinn a so-called "compressed" matrix o E m to outputs. Thus, the compressed version of training stage in SHORE can be represented follows,
W
== arg min We
(W) = oy WX 1.
(2)
n
We would like to point out that the idea of compressing the output space into some smaller intrinsic dimension has been used in many existing works, e.g., [17 39 [2] mentioned in Section 2
Prediction stage. In the second prediction stage, given any input x E R the framework predicts a sparse output ÿ by solving the following prediction problem based on the learned regressor W in the training stage,
y (W) == arg min by - W Wxll2 s.t. y DKK nF,
(3)
y
where D* is the set of s-sparse vectors in R and F is some feasible set to describe additional requirements of For example, by letting F be RK RK {0, 13 the intersection V* n_F denotes the set of s-sparse output, non-negative s-sparse output, {0, 1}-valued s-sparse output, respectively.
We use D( W) (shorthanded in T to specify that the predicted output is based on the regressor W. To solve the proposed prediction problem (3), utilize the following projected gradient descent method (Algorithm T). which could be viewed as a variant/generalization of existing iterative thresholding methods 6 21 for nonconvex constrained minimization. In particular, step 4 incorporates additional constraints from F other than sparsity into consideration, which leads to non-trivial modifications in designing efficient projection oracles and convergence analysis. Later, we show that the proposed Algorithm I ensures a near-optimal convergence (Theorem 2 and Theorem 4 while greatly reduces the computational complexity (Remark() (2) of the prediction stage for solving compressed SHORE.
Before diving into theoretical analysis, we first highlight the differences between the proposed prediction stage (3), general cpansity optimization (SCO), and sparse regression in the following remark.
Remark 1. Proposed prediction stage v.s. General SCO: To be clear, the SCO here denotes the following mininization problem minjjollosk | B11 Thus, the prediction stage is special case

---

of general SCO problem. In particular. the predicted stage takes a random projection matrix c with restricted isometry property (RIP) be its A and uses Wa with W obtained from the compressed -ttage to be its B.. As a result (Theorem 2 and Theorem A. the proposed Algorithm 7 for prediction stage ensures globally linear convergence to ball with center y (optimal solution of the prediction-stage and radius @ī Wxll2), which might not hold general SCO problems.
Proposed prediction stage v.s. Sparse regression: Although the proposed prediction stage and sparse high- high--dmmensional regression share a similar optimization formulation minjjallo ||Y X B11..
the proposed prediction stage (3 is distinct from the sparse regression in the following parts:
(I) Underlying Model:
Most existing works about sparse high-dimensional regression assume that samples are i..dd generated from the linear relationship = 1TxTB 3 + with underlying sparse ground truth 3 In the proposed prediction stage, do not assume additional underlying models on samples if there is no further specific assumption. The problem we studied in the predicted stage takes the random projection matrix & with restricted isometry property (RIP) as its X (whereas X in sparse regression does not ensure RIP), and uses W x with W obtained from the compressed training stage as its Y.
=x-
(2) Problem Task:
The sparse regression aims recover the sparse ground truth 3 given sample set ((x', yy) "²= with n i..dd samples. In contrast, the task of the proposed prediction stage is to predict a sparse high-dimensional output given a random projection matrix T and single imput x. As a quick summary, some typical and widely used iterative algorithms [38 29 for sparse regression cannot be directly applied to the proposed prediction stage.
Then, we provide the computational complexity with and without the compression for the proposed two-stage framework to complete this subsection.
Remark 2. Training stage: Conditioned XX is invertible, the compressed regressor W has closed form solution W x (XxT)-1 with overall computational complexity
O mmKmm + mnd + nd? + d + md?) zo O (KK mm)).
Compared with the computational complexity of finding 2 from the uncompressed SHORE ① O(K nd + nad2 d KdP) Rs O(K(n + djd),
solving W enjoys a smaller computational complexity on the training stage if m In later analysis (see Section 3.2). m = 1088-2, s log( K ))) with some predetermined constants 8, and sparsity-level d, thus many applications with large output space, the condition m < holds.
Prediction stage: The computational complexity of each step-3 3 in Algorithm 7 is
O(Km - K + Km + K) = O(m)).
The projection in step is polymomially solvable with computational complexity O K min (s, log K}) (see proof in Appendix A.S.I. Thus, the overall computational complexity of Algorithm is
O(K(m min(s, log K))T).
Compared with the complexity O(K(d (d + min {s. log KJ)) of predicting ÿ from the uncompressed SHORE (1. the compressed version enjoys smaller complexity on the prediction stage if
(m + min(s, log K})T < + min(s, log K).
(4)
In later analysis (see Theorem since ((8-2_ log( ))) with predetermined constants d, sparsity-level d, and T Ologl - 122 from inequality 5 we have condition 4 holds.
(D--Wall
Whole computational complexity: Based on the analysis of computational complexity above, we conclude that when the parameters (K, d, m, satisfies
> k1/3 d > 0(8 log(K/r) =TT, K>
the compressed SHORE enjoys a better computational complexity with respect to the original one 1
5

---

Algorithm/ 1 Projected Gradient Descent (for Second Stage)
Input: Regressor W, input sample æ, stepsize n) total iterations T
1: Initialize point ( E V n.F.
2: for t = 0, ..... T 1: do
3: Update ((++1) (tt) 7 T ( W x)).
4: Project ((++)) I(((++)) arg min 1/0 11
5: end for
Output: ((T)
3.2 Worst-Case Analysis for Arbitrary Samples
We begin this subsection by introducing the generalization method of the compressed matrix db Assumption 1. Given an m-by-K K compressed matrix $ all components $i, for S i S m and 1 J & K, are i.i.d. generated from Gaussian distribution (0, 1/m).
Before presenting the main theoretical results, let us first introduce the definition of restricted isometry property (RIP, [1O), which is ensured by the generalization method (Assumption] D].
Definition 1. (V, P).P A m-by-K K matrix + is said to be (V, RRPP over given set of vectors V RRK for every o E V.
(1 (1- 1- < 11 S (1 (1+0)||u11.
In the rest the paper, we use (s, RRIP to denote 5 Recall VKK = for RR 1101|o S 8) is the set of s-sparse vectors.
Remark 3. From Johnson-Lindenstruuss Lemma 43 for any E (0. 1). any T E (0, 1). and any
finite vector set [V] < co. if the number of rows m 2 88-2 log( ) then the compressed
matrix & generated by Assumption 7 satissies (V,8)-RIP with probability at least 1 T.
Now, we are poised to present the first result on training loss defined in 2
Theorem 1. For any 5 E (0, 1) and t E (0, 1), suppose compressed matrix ob follows Assumption 7 with m 2 O log( K )). We have the following inequality for training loss
l W S (1+ Y 2
holds with probability at least 1 - T, where 2, W are optimal solutions for the uncompressed
and compressed SHORE (2). respectively.
The proof of Theorem D is presented in Appendix A In short, Theorem D shows that the optimal training loss for the compressed version is upper bounded within (1 + 5) multiplicative ratio with respect to the optimal training loss for the uncompressed version. Intuitively, Theorem [ implies that SHORE remains similar performances for both compressed and compressed versions, while the compressed version saves roughly O(Kn - m) + Kd) computational complexity in the training stage from Remark Moreover, the lower bound condition on m 2 O log(
K )) ensures that
the generated compressed matrix $ is (1, 8)-RIP with probability at least 1 T. For people of independent interest, Theorem only needs (1, RIP (independent with the sparsity level) due to the unitary invariant property of & from Assumption I (details in Appendix (A.I). Additionally, due to the inverse proportionality between m and 8², for fixed K and т, the result can be written as
x W < (1 + 01)/m))))))mmmmmmmmmmmm))))))) ||Y 2x
which is verified in our experiments
We then present the convergence result of the proposed Algorithm for solving prediction problem 3 Theorem 2. For any (0, 1) and E (0, 1), suppose the compressed matrix & follows Assump-
tion with 2 0 fz log( K) )) With a fixed stepsize 1) E (3-22- 1), the following inequality
| - o 11 s 1 (00) ||2
C2
4 Wxl12
c1
holds for all i E [T] simultaneously with probubility least 1 where := 2 2-21 21p < 1 is
some positive constant strictly smaller than 1, and Co 2n + O is some constant.

---

The proof of Theorem 2 is given in Appendix A.2 Here, the lower bound condition on the number of rows m ensures that the generated compressed matrix (3s, 6)-RIP with probability at least 1 by considering a 6/ 922-net cover of set V = V3. n B2(0; 1) from Johnson-LLindessraane Lemma 4 Moreover, since the number of rows m required in Theorem 2
is greater than the one required in Theorem 1 term | Wællz can be further upper bounded using the uncompressed version (1 + ū Z la with probability at least 1 - T. Then we obtain a direct corollary of Theorem
suppose l
(0) 2 - W if
1 4. = O log |
(0) 1/ |2/ W log (1/ci)).
(5)
the proposed Algorithm guarantees globally linear convergence to a ball (g: [1 W W
In contrast with OMP used in [17] for multi-label predictions, Theorem Saaeemmm@hor holds for arbitrary sample set without the so-called bounded coherence guarantee on Φ. Moreover, as reported in Section 4 the proposed prediction method (Algorithm T has better computational efficiency than OMP.
3.3
Generalization Error Bounds for IID Samples
This subsection studies specific scenario when every sample (x', sr'.y") is i..d. drawn from some underlying subGaussian distribution 'D over sample space IRd x yk Specifically, we use
ED
P fly
and Varp
20 (
E ae Eya
Eay
=: O(d+K)x(d+K)
Eyy
to denote its mean and variance, respectively. Let &a = x pw: Ey == My. (Ee Egy be centered subGaussian random variables of æ. y, (x
respectively. Let a. b E ( {æ, y). we Mab ED jab Eab paal4, Mab := E a' (b')T to denote the population second (oooss))mmmetss and empirical second (cross-)moments, respectively. Then, the population training loss is defined as
min C(Z) == (( [lly ZERK xd
Z
with its optimal solution Z. = My Ma Similarly, given a 4, the compressed training loss is given by
- (w) = B [| W
with optimal solution W, = $ MMM M We then define the following assumption:
Assumption 2.
Let D be o2-subGaussian for some positive constant o² O, i.e., the inequality Ep (exp(Xu E)] S exp (X²o2/2) holds for any X > O and unitary vector u E RdtKK Moreover, the covariance matrix Exn is positive definite (i.e., its minimum eigenvalue Amin (x 0).
Remark 4. Assumption 2 Jenurees the light tail property of distribution DD Note that in some real applicatioons e.g.. factuality check [3TT, algorithmic trading [16), one can normalize input and output vector ensure bounded E2-norm. Under such situation, Assumption 2
naturally satisfied.
Our first result in this subsection gives the generalization error bounds.
Theorem 3. For any 8 E (0, 1) and r E (0, )), suppose compressed mairix & follows Assumption ! with 0 &a log( K )), and Assumption 2 holds, for any constant O, the following results hold:
(Matrix Error), The inequality for matrix error M
1/2 M
M
1/2
< 4 holds with probability at
least 1 the number of samples ni with
24
322 Ps 12202
(2v4+ 1og(1/7) 21
01 max
9X2 mini (Mra) (d+ log(2/T)),
min (Ma))
where some fixed positive constant used in matrix concentration inequality of operator norm.
(Uncompressed). The generalization error bound for uncompressed SHORE satisfies ((2) s [(Z.) + 4e with probability at least 1 - 3T, as the number of samples n 2 max(ni, with
d 2 d logg//// K/T) + +2 2og (K/T) Alliy Z. Z.
d
= max 4(1 4(1 + K)

---

(Compressed). The generalization error bound for the compressed SHORE satisfies L ( - _ (W.) de with probubility least 3r, as the number of sample n 2 max (11, R22) with
A[I] + [+)) d+ 2+22/0 22dddlogggg llg(m/// + 2lo(mm/T)
d c
M2 = max
11 W.pa
The proof of Theorem Biss presented in Appendix A.3 The proof sketch mainly contains three steps: In Step-1, I, we represent the difference ( L ZZ.)) - (W) L. (W.) as product between matrix error (in Theorem B and rescaled approximation enror (see Appendix A.3 for definition), i.c.
( ((z)) or L (W) - - (W.)= I (matrix error) x (rescaled approximation error);
Step-2 controls the upper bounds for matrix error and rescaled approximation error separately, using concentration for subGuassian variables; for Step-3, we combine the upper bounds obtained in Step-2 and complete the proof. Based the result of Theorem ignoring the logarithm term for т, the proposed generalization error bounds can be bounded by
d ((2) S ((Z.) Ö [max(]] 1. 1lly Z. K} m
n1
C (W) (W))& (W.) + (max(] W. 1|
W.p 112: 11 71
d
Remark 5. To make a direct comparison between the generalization error bounds of the uncom- pressed and the compressed version, we further control the norms || WWW.le | W W.Pa 112. + based on additional conditions on the compressed matrix Q. Recall the generalization method of the compressed matrix & mentioned in Assumption we have the following event
7
W I| = 2Z. 5 (1+ 8) 11?
& == E BR
E 1 1F 4 ( Z. Z.pz (1+ pm l|?
holds with probability at least 1 - T due the RIP property for fixed matrix. Moreover since every component Pig is i.idd drawnfromt a Gaussian distribution N (0, 1/m). using the concentration tail bound for chi-squared variables (See Lemma I in (261), we have the following event
K log(1/r) 2 log((1/T)
8 = 1 E RR
| | | SKSK22
m
171
holds with probability at least T. Conditioned on these two events & and &g, the generalization ervor bound of the compressed version achieves the same order (ignoring the logarithm term of r) as the gemeralization error bound of the uncompressed version. That is to say,
d 12
c (*) (1+ C(Z.) ((,)))
(rax(12.1133 I, Z. K}
holds with probability at least 5T.
Comparing with existing results generalization error bounds mentioned in Section 2 would like to emphasize that Theorem guaranttess that the generalization error bounds maintain the order before and after compression. This result establishes on i.i.d. subGaussian samples for the SHORE model without additional regularity conditions on loss function and feasible set as required in [45). Additionally, we obtained O ((d//n) generalization error bound for squared Frobenius norm loss function L or L which is smaller than (K²d/n) as presented in [Theorem 4, [55
We then give results on prediction error bounds.
Theorem 4. For any d E (0, 1) and any T E (0,1/3), suppose the compressed matrix & follows O log( K)), and Assumption holds. Given any learned regressor W from training problem (2), let (ar, y) be new sample drawn from the underlying distribution D, we have the following inequality holds with probability at least -Ti
Assumption 1 with
Ep 1
vl
- 1 1
W
where ÿ is the optimal solution from prediction problem with input vector x.
8

---

A.4 Theorem gives an upper bound of é2- S v v v yll2. combined with 2 when T 2 t. defined in (5) (see Appendix A.5pe)) where the final inequality holds due to the optimality of ÿ. Hence, we achieve an upper bound of (z-norm distance between v(T") and as presented in Theorem see Remark 6
The proof of Theorem is presented in Appendix
norm distance between and Since ||
y
Theorem 2 we have ( ylll O(E
((7)
I1 W
Remark 6. For any (0. and E (0, 1/3), suppose compressed matrix & follows Assump- tion with m 2 O log( )), and Assumption holds, for any constant e > 0, the following inequality holds with probability at least 1 - 3r:
EDlle víl S O(EDIl by w SO(C 4e).
4 Numerical Experiments
In this section, we conduct numerical experiments on two types of instances (i.e., synthetic data sets and real data sets) validate the theoretical results and illustrate both efficiency and accuracy of the proposed prediction method compared with typical existing prediction baselines, i.e., Orthogo- nal Matching Pursuit (OMP, (46)), Fast Iterative Shrinkage-Treeshoddnng Thresholding Algorithm (FISTA, (2) and Elastic Net (EN, 47) Due to the space limit, we put the implemented prediction method (Algorithm 2 in Appendix A.6.1, aforementioned existing prediction baselines in Appendix A.66.
experiment setting details and results for real data in Appendix A.7
Performance measures. Given sample (z, y) with input æ and corresponding true output y. we use to denote the predicted output obtained from any prediction method, and measure the numerical performances based on the following three metrics:
1. For a ground truth y with sparsity-level s, the metric precision selected supports, i.e. Precision@e == //sppp( n supp(y)) measures the percentage of correctly identified supports in the predicted output;
2. The metric output difference, i.e., Output- - diff == |/ v measures the 2n-mrrm distance between the predicted output and the ground-truth,
3. For any given MOR W and compressed matrix ④ the metric prediction loss, i.e., Prediction Loss cvu- Wxll3 computes the prediction loss with respect to Wa a.
Synthetic data generation procedure. The synthetic data set is generated as follows: Every input x for E [n] is i.i.d. drawn from a Gaussian distribution N (Me, Eax), ). where its mean vector p and covariance matrix Eaw are selected based on the procedures given in Appendix A.6.3 For any given sparsity-level. s, underlying true regressor Z. R xd and Signal-to-Nol -ooiise Ratio (SNR), the ground- truth y' (corresponding with its given input x') is generated by y = Iygn Z + +²). where e E K a.i.d. random noise drawn from the Gaussian distribution ((x SNR '1 l Ik).
Parameter setting. For synthetic data, we set input dimension d = 10*, output dimension K = 2 x 10t, and sparsity-level above, i.c., Szyn := {(x', , )² 3x300 with SNR 1 E [1, 0.32, 0.032] to ensure the signal-to-novl decibels (dB, [4] takes values on dB == 10 (SNR) E {0. 10, 30). We select the number of rows for compressed matrix & by E {100, 300, 500, 700, 1000, 2000). For computing the empirical regressor W E Rmmxd we first split the whole sample set Saya into two non-overlap subsets, i.e., a training set Stza with 80% and a testing set Stest with rest 20%. The regressor W is therefore obtained by solving compressed SHORE (2) based on the training set Stra with a randomly generated compressed matrix &. For evaluating the proposed prediction method, Algorithm 22, we pick fixed stepsize 77 0.9, JF K and set the maximum iteration number as 60, and run prediction methods over the set Steat
3. We generate in total
3 x 10*, i.i.d. samples described
Hardware & Software. All experiments are conducted in Dell workstation Precision 7920 with 3GHz 48Cores Intel Xeon CPU and 128GB 2934MHz DDR4 Memory. The proposed method and other methods are solved using PyTorch version 2.3.0 and seikit-leam version 1.4.2 in Python 12.3.
Numerical Results & Discussions. The results are demonstrated in Figure I which does not include the results from the Elastic Net and OMP due to relatively much longer running time.

---

637
030
70
120
150
110
108
100
136
70
130
Pasce115 mande Ziee USt b852
Produc ioseade whar w.55
Protuss sader 33tt
9000 1236 2008
2700
1000
Figure 1: Numerical results on synthetic data. In short, each dot in the figure represents the average value of 10 independent trials (i.e., experiments) of compressed matrices ((1), (100) on a given tuple of parameters (K, d, n, SNR, m)) The shaded parts represent the empirical standard deviations over 10 trials. In the first row, we plot the ratio of training loss after and before compression, i.e.. 2xI versus the number of rows m. is obvious that the ratio converges to one increases, which validates the result presented in Theorem [ In the second row, plot percision @ versus the number rows. As we can observe, the proposed algorithm outperforms CD and FISTA.
I W X / |
Based on Figure | we observe that the proposed algorithm enjoys better computational cost and accuracy on most metries. The running time for the proposed algorithm and baselines are reported in Table[[(see in Appendix (A7)), which further demonstrates the efficiency of the proposed algorithm. The implemented code could be found on Github https: ftttps://gitull com frr//rroo-ooo coom/trom-rran/Sllvvng SHORE_vise Lire__vi__co__ www.aaccopeession.
5 Conclusion and Future Directions
In conclusion, we propose a two-stage framework to solve Sparse & figh-diinensionl--uure REgression (SHORE) problem, the computational and statistical results indicate that the proposed framework is computationally scalable, maintaining the same order of both the training loss and prediction loss before and after compression under relatively weak sample conditions, espe- cially in the sparse and gigh-dimeesiioal-ouppe setting where the input dimension is polynomially smaller compared to the output dimension. In numerical experiments, SHORE provides improved optimization performance over existing MOR methods, for both synthetic data and real data.
We close with some potential questions for future investigation. The first is to extend our theoretical results to nonlincarhnonconvex SHORE frameworks [24). The second direction is to improve existing prediction variable reduction methods for better scalability while maintaining small sacrificing accuracy, e.g., new design and analysis on randomized projection matrices. The third direction is to explore general scenarios when high dimensional outputs enjoys additional geometric structures 30 from real applications in machine learning or operations management other than s-sparsity and its variants as discussed in the paper. Taking result for SHORE as initial start, expect stronger follow-up up work that applies to MOR with additional structures, which eventually benefits the learning community in both practice and theory.
10

---

Acknowledgement
Renyuan Li and Guanyi Wang were supported by the National University of Singapore under AcRF Tier-1 grant (A-8000607-00-00) 22-5539-A0001. Zhehui Chen would like to thank Google for its support providing the research environment and supportive community that made this work possible.
References
[1] Z. Abraham, P.-N. Tan, Perdinan, J. Winkler, S. Zhong, and M. Liszewska. Position preserv- ing multi-output prediction. In Machine Learning and Knowledge Discovery in Databases: European Conference, ECML. PKDD 2013, Prague, Czech Republic, September 23-27, 2013, Proceedings, Purt II 13, pages 320-335. Springer, 2013.
[2] A. Beck and M. Teboulle. A fast iterative shrinkage- thankkgggg-hreshldingg algorithm for linear inverse problems. SIAM journal on imaging sciences, 2(1):183-202, 2009.
[3] D. Bertsimas and B. Van Parys. Sparse high- ditheeinnniionll regression: Exact scalable algorithms and phase transitions (2017). arXiv preprint arXiv:1709. 10029, 2019.
[4] D. Bertsimas and B. Van Parys. Sparse high-dimensional regression. The Annals of Statistics, 48(1):300-3233 2020.
[5] K. Bhatia, K. Dahiya, H. Jain, P. Kar, A. Mittal, Y. Prabhu, and M. Varma. The extreme classification repository: Multi-label datasets and code, 2016.
[6] T. Blumensath and M. E. Davies. Iterative hard thresholding for compressed sensing. Applied and computational harmonic analysis, 27(3):265-277 2009.
[7] H. Boche, R. Calderbank, G. Kutyniok, and J. Vybíral. Compressed sensing and its applications: MATHEON workshop 2013. Springer, 2015.
[8] H. Borchani, G. Varando, C. Bielza, and P. Larranaga. A survey on multi-output regression. Wiley Interdisciplinary Reviews: Data Mining and Knowledge Discovery, 5(5):216-2)) 221166233,2222 2015.
[9] L. Breiman. Bagging predictors. Machine learning, 24:123-140, 1996.
[10] E. J. Candes and T. Tao. Decoding by linear programming, IEEE transactions on information theory, 51(12):4203-4215, 2005.
[11] Y. Chang, X. Wang. J. Wang. Y. Wu, L. Yang, K. Zhu, H. Chen, X. Yi, C. Wang. Y. Wang, et al. survey on evaluation large language models. ACM Transactions on Intelligent Systems and Technology, 15(3):1 15,5, 2024.
[12] Y.-N. Chen and H.-T. Lin. Feature-aware label space dimension reduction for multi-label classification. Advances in neural information processing systems, 25, 2012.
[13] M. M. Cisse, N. Usunier, T. Artieres, and P. Gallinari. Robust bloom filters for large multilabel classification tasks. Advances in neural information processing systems, 26, 2013.
[14] S. S. Dey, G. Wang, and Y. Xie. Approximation algorithms for training one-node relu neural networks. IEEE Transactions on Signal Processing, 68:6696-6706, 2020.
[15] A. E. Hoerl and R. W. Kennard. Ridge regression: Biased estimation for nonorthogonal problems. Technometrics, 12(1):55-67, 1970.
[16] O. Honovich, R. Aharoni, J. Herzig, H. Taitelbaum, D. Kukliansy, V. Cohen, T. Scialom, Szpektor, A. Hassidim, and Y. Matias. True: Re-evaluating factual consistency evaluation. arXiv preprint arXiv:2204.04990
2022.
[17] D. Hsu, S. M. Kakade, Langford, and T. Zhang. Multi-label prediction via compressed sensing, 2009.

---

[18] D. Hsu, S. M. Kakade, and T. Zhang. An analysis of random design linear regression. arXiv preprint arXiv: 1eXiv:100.22666 6, 2011.
[19] C. Hull and S. Basu. Options, futures, and other derivatives. Pearson Education India, 2016.
[20] I. M. Jacobs and J. Wozencraft. Principles of communication engineering. 1965.
[21] P. Jain, A. Tewari, and P. Kar. On iterative hard thresholding methods for high-dimensional m-estimation. Advances in neural information processing systems, 27, 2014.
[22] S. Jansen. Machine Learning for Algorithmic Trading: Predictive models to extract signals from market and alternative data for systematic trading strategies with Python. Packt Publishing Ltd, 2020.
[23] K. Jasinska and N. Karampatziakis. Log-time and log-space extreme classification. arXiv preprint arXiv: 161101964 2016.
[24] A. Kapoor, R. Viswanathan, and P. Jain. Multilabel classification using bayesian compressed sensing. Advances neural information processing systems, 25, 2012.
[25] S. Kim and E. Xing. Tree-guided group lasso for multi-response regression with structured sparsity, with an application to eQTL mapping. The Annals of Applied Statistics, 6(3): 6(3):0955- 1117, 2012.
[26] B. Laurent and P. Massart. Adaptive estimation of quadratic functional by model selection. Annals of statistics, pages 1302-1338, 2000.
[27] T. Levy and F. Abramovich. Generalization error bounds for multiclass sparse linear classifiers. Journal of Machine Learning Research, 24(151):1-35, 2023.
[28) S. Li and Y. Liu. Towards sharper generalization bounds for structured prediction. Advances in Neural Information Processing Systems, 01-26844-26857 2021.
[29] L. Liu, Y. Shen, T. Li, and C. Caramanis. High dimensional robust sparse regression. In International Conference on Artificial Intelligence and Statistics, pages 411-421. PMLR, 2020.
[30] S. Ludwig et al. Algorithms above the noise floor. PhD thesis, Massachusetts Institute of Technology, 2018
[31] Y. Ma, R. Han, and Wang. Portfolio optimization with return prediction using deep learning and machine learning. Expert Systems with Applications, 165:1 13973, 2021.
[32] T. K. R. Medini, Q. Huang, Y. Wang, V. Mohan, and A. Shrivastava. Extreme classification in log memory using count-min sketch: A case study of amazon search with 50m products. Advances in Neural Information Processing Systems, 32, 2019.
[33] M. Riedmiller, R. Hafner, T. Lampe, M. Neunert, J. Degrave, T. Wiele, V. Mnih, N. Heess, and J. T. Springenberg. Learning by playing solving sparse reward tasks from scratch. In International conference on machine learning, pages 4344 4353. PMLR, 2018.
[34] A. Roberts, C. Raffel, and N. Shazeer. How much knowledge can you pack into the parameters of language model? arXiv preprint 100i::20022089100 2020.
[35] M. Sänchez-FFemándee M. de Prado-Cumplido J. Arenas-Gras García, and F. Pérez-Cruz. Svm multiregression for nonlinear channel estimation in multiple-inpur multiple-output systems. IEEE transactions on signal processing, 52(8):2298-2307, 2004.
[36] T. Similä and J. Tikka. Input selection and shrinkage in multiresponse linear regression. Computational Statistics & Data Analysis, 52(1):406-475 2007.
[37] E. Spyromitros-Xioon methods for multi-target regression. arXiv preprint arXiv:1211.6581, pages 1159-1168, 2012.
G. Tsoumakas, W. Groves, and I. Vlahavas. Multi-label classification
[38] Q Sun, H. Zhu, Y. Liu, and J. G. Ibrahim. Sprem: sparse projection regression model for high- dimensional linear regression. Journal of the American Statistical Association, 110(509):289- 302, 2015.
12

---

[39] F. Tai and H.-T. Lin. Multilabel classification with principal label space transformation. Neural Computation, 24(9):2508 2::-2222228222,, 2012.
[40] D. Tuia, Verrelst, L. Alonso, Pérez-Cruz, and G. Camps- Valls. Multioutput support vector regression for remote sensing biophysical parameter estimation. IEEE Geoscience and Remote Sensing Letters, 8(4):804 808, 2011.
[41] M. Wainwright. High-dimensional statistics: A non-asymptotic viewpoint, volume 48. Cambridge university press, 2019.
[42] A. Watkins, Ullah, Nguyen-Tang and R. Arora. Optimistic rates for multi-task representa- tion learning. Advances in Neural Information Processing Systems, 36, 2024.
[43] B. J. William and J. Lindenstrauss, Extensions of lipschitz mapping into hilbert space. Contem- porary mathematics, (6(189-206)3323 1984.
[44] D. Xu, Y. Shi, I. W. Tsang, Y.-. Ong. C. Gong, and X. Shen. Survey on multi-output learning. IEEE transactions neural networks and learning systems, 31(7):2409-244 2429, 2019.
[45] H.-F. Yu, P. Jain, P. Kar, and I. Dhillon. Large-scale multi-label learning with missing labels. In International conference on machine learning, pages 993-601. PMLR, 2014.
[46] Z. Zhang. Matching pursuits with time-frequeng dictionaries, IEEE Transactions on signal processing, 41(12):3397-3415 1993.
[47] H. Zou and T. Hastie. Regularization and variable selection via the elastic Journal of the Royal Statistical Society Series B: Statistical Methodology, 67(2):301-322) 2005.
13

---

A Appendix supplemental material
A.1 Proof of Theorem 1
Recall the following theorem in the main text.
Theorem 1
For any E (0, 1) and E (0, 1), suppose the compressed matrix & follows Assump-
tion with m 2 O O log( K) )) We have the following inequality for training loss
+ W (1 + 6) Y 2
holds with probability least 7, where Z, W optimal solutions for the uncompressed and compressed SHORE (2). respectively.
Proof. Given set of n samples {(x', y )')) and matrix & generated from Assumption | have
bY W ΦY xãx 1 || I
where the inequality holds due to the optimality of W. Let Y' = - 2x. Consider the singular value decomposition: Y Uy Ey Vy,,, and then we have
( xx) 1 4 2 DUUr
= bUyr Er
Set +- u Since the generalization method for & ensures its 61,6)-RIP with probability Lemma), and Uy is a real unitary matrix, then Lemma [ for unitary invariant shows that is also (1, 5)-RIP with probability - r. Now, using T (1, RIIP and all columns in Σy has at most one non-zero component, have
1 - T (from the
UU Ey || g 11 s (1+ + o) 13
(1 + 8)1 Y'I.
Combining the above inequalities together implies
x W 1| (1 8) Y (1 5) |Y 2
which completes the proof.
A.2 Proof of Theorem
Recall the following theorem in the main text.
Theorem
2
For any 8 (0, 1) and (0, 1), suppose the compressed matrix & follows Assump-
tion I with
O fz log( K )). With fixed stepsize 17 E (2-28 1), the following incquality
| v 19 ll2 S 1 ( 1|2
C2
pÿ W W
ci
holds for all t e [T'] simultaneously with probability at least 1 where := 2 2n + 2nj < 1 is some positive constant strictly smaller than 1, and e2 = 2n/1 + d is some constant.
Proof. Suppose Assumption [ holds, then the randomized compressed matrix & ensures (3s,8)-RI) RIP with probability at least 1 - (see Remark 3). Thus, to complete the proof, it sufficient to show the above inequality holds for all [T] under such compressed matrix & with (3s, 6)-RIP. We conclude that the proof is, in general, separated into three steps.
Step-1. We establish upper bound on the fg-norm distance between the current point and the optimal solution. Due to the optimality of the projection (step in Algorithm I]. we have the following inequality
| n (1+1)
(1((+)]] =012
holds for all v VK nF, which further implies that
[|](++|]
(t+ 1| < || 2/12
# ||l (1+1) < 2(v x((+1) (((+1))
holds for all VK n.F.
14

---

Step-2. We show one-iteration improvement based on the above inequality. Since ÿ E V n.F, we can replace v by ÿ. which still ensures the above inequality. Set A (t) := - - v vl for all t E [7]. Based on the updating rule (step in Algorithm the above inequality can be written as
v((+1)
v " + T ( Wa a). Thus,
A (++1) 12 S ( (()
+1 (Do
Wx), A (112))
= 2( (4) A
( + 2n (
Wa x. D (+11)
2
(++1))
$ (
A (++1)) + (D
Wæ, $ (11)).
(6)
Using Lemma 2 with
A (1)
A (++1)
A (t)
A (++1)
A (4)
A (++1)
1122
A 11
| (
A 11
A (0) 2
1 2
all 3s-sparse vectors, and + a (3s, 8)-RIP matrix, we have
(E)
A (1-1)
(e)
(++1)
27
S 2ōn - 2n
(1+1)
A 12
(111)
which implies
2 (4)
A
()+
- 20nl A (8) |2 (++1))
112 - 2m A (4) A ((+!))).
Inserting the above result into inequality gives
(++1) 1 s (2 2n) A (4)
A
((+1)) + 25nll A (4) 12 A (++ 11)
1l2 + ΦD
Wo, T (1+1))
(4)
< (2 2n + 2no)
2 A (++1)
1|2
2ny ワ
W x |||2| 4
2
< (2 - 2n + 2nd)
(0)
|||| (++1)
1| + 2 + 6 Φ6 Wx
|2
[1+1)
where the above inequality (i) requests 7) < 1 to ensure the inequality (2 2n) A (+)
(+))) <
1-22n)ll (+) ||2|| A (c+1)|12 holds. Therefore, dividing ||
(t-1)]|2 on both side implies
a (1+1) 12 S (2 2n 2nj5 (0 1/22 + 2n 1 Φ W
which gives the one-step improvement.
Step-3. Combine everything together. To ensure contractions in every iteration, we pick stepsize such that 2 - 2n + 2nj6 E 10, with < 1, which gives 1) E 211- 1). Using the above inequality ( for one-step improvement, we have
21 1 +
1 n 00 112 S (2 22 2n + 2n5)². 1 -0(0) 122 + 2n 2np 1
lT -Wa Wx 22|2,
which completes the proof.
A.3 Proof of Theorem 3
Recall the following theorem in the main text.
Theorem 3. For any 8 E (0, 1) and (0. 4), suppose compressed matrix & follows Assumption
|
with O log( )), and Assumption 2 holds, for any constant 0, the following results hold:
(Matrix Error). The inequality for matrix error
|M 1/2 M ze
M
holds with probability at <
up
least 1 - 27 as the number of samples 2 n1 with
24
322 |z 12002
(2v4+ + vlog(I/r) ")
n1 == max
912 1 (Maw) (d + log(2/7)),
X2
(Mau
15

---

where C' is some fixed positive constant used in matrix concentration inequality of operator norm.
(Uncompressed). The generalization error bound for uncompressed SHORE satisfies ((2) s (Z.)) + with probability least -3,, as the number of samples n 2 max(m1, n2) with
== max 4(1l + K) d+ 2 log K/T) + (log(K/T) All Z. 11 $
d
(Compressed). The generalization error bound for the compressed SHORE satisfies < s _ (W.) 4e with probability least 3r, as the number of sample n Z max max(n1, 722) } with
M2 = max 4(1]W..13 + ( d 2 dlog(mmT)) + (2 llg(m//T)
1 W. 12 d
Proof. Let us start with the uncompressed version.
Step-1. Note that the optimal solutions for population loss and empirical loss are
Z. M yau M and 2 YX
xX
-1
Myn M 20
respectively. Thus, the generalization error bound is
((2) L(Z.) 1 Z 21
11 Z.) 2
Note that
(Ø - M.)M
M Myx MEi Myx
M
Myn MynMain Man
M M M22
Elyæ ZZ.xx
M 1/2 IM wu wi
El(y Z. M
1
M 1/2
where we use E[ to denote the empirical distribution. Then, the above generalization error bound can be upper- bounded as follows
| 12 M. 1/2
[i - Z. ) M M M-112 M 1/2 x0
E((y - Z. T Madl1 |
M Mau' 12 M 1/2
op
rescaled approximation errur
matrix eTTOr
Step-2. Next, we provide upper bounds on these two terms
Bi(y - Z.
Mau 1/22//1 //////
and
M 1/2 1/2 in the right-hand -ide separately. have
op
For matrix error term
Mar M
1/2
M- 1/2 1
M
11/2,-1 MM/2M-1MM M1/2 ma
Due to Assumption 2 the centralized feature vector Enr == per ensures the following inequality
Ep [exp (Av [e)] < exp (A3112
(Ab11e|Bo2
2
for all v E Rd, Consider the empirical second moment of x,
T
x (x)) T
M Σ
Su
+
See
pr +
m=1
i=1
16

---

Since e -. are ii.d o² subGaussian random vector with zero mean and covariance matrix Exn.. then based on Lemma there exists positive constant C such that for any T E (0, 1),
d log(2/T) d+ log(2/T) })
(5x)
PD
2ew
I Co max
2 21-7,
op
and based on Lemma
for any T E (0,
a
40 4o 20 log (1/)
PD
21-T.
n
vn
Let A =
Ei=1 E T
and E
then M can be represented by Ei=l E.
M
mm un Euz Enn Mrlll
=MMan
and thus we have minimum eigenvalue of M 1/2 M M
1/2
M on M
+ M
Then the can be lower bounded follows
1/2 (A MEE
E M 1/2
1/22
MM 1/2 M M=1/2)
Aminn
1/2 M
M
ax + MxE
+ Eu
M M=1/2 1/2
lp
2
M 1/2
ma 2000
1/2
1/2 (MoE
+ Elr
M 1/2
sea
Co²
d
log(2/T)
211M22112 40 4o 20 log(1/T)
2 1
Amin (M
Amin (Max)
vn
where the final inequality holds with probability at least 2r by inserting the above non-asymptotic bounds. Then, have
1/2 nw
M-
1/2 va
-Amln (M
Maw M M/1
Co²
d + log(2/T)
4ov + 20 V log(1/T)
Amin
Amin (M
Vn
holds with probability at least
2T. It is easy to observe that as
240
32a M 1 (d log(2/T)),
m max
(2Vd + vlog(1/T)
9X2
(Mxx)
X2 min (Max)
have
M m
M- 1
1/2
holds with probability 1 27.
For rescaled approximation enror term E((y Z.z)x M 1/2
F we first compute variance
proxy for the subGaussian vector Z, Note that the j-th component of the subGaussian vector Z.æ can be written
y Z,x,, = (-(2.) e
where [Z.]]. is the j-th row of Z., e. is a K K-dimensional vector with j-th component equals and rest components equal zero. Thus, it is easy to observe that the 4--oorm square of
(-1 -Z. e, is I1,,,,, 113 + 1, and therefore, based on the Assumption 2
we have
Ep exp A22
Z. æ - Alp Z.Huml)) < exp (0l. + +100212).
17

---

i.e., subGaussian with variance proxy o3:= == -1112,1;:11 1)o². error can be upper-bounded bounded by
Thus the rescaled approximation
[ly Z
M M 1/2
Ely Z Maxe 1/2 2
Σ
==1
2
B[(ly Z,zily
Myy Z )=
M M=1/2 1/2
+2
Bly ZuNaoaaa
M M 1/2
TT
=:T
We control term T3 for all j E [K] separately using Lemma as follows: For all E (0, 1), we have
of(d oj (og((K/T) + llg((//))) 2 21-7/K
PD Tj
For the term T7. we have
73
Uy
Z. (xli((')
M M-112 1/2
71
i=1
My
Zupelli
xc
M
M-112 12
M 22
1/2
1n
VT
My
Z.MLali
12
x
M-1/22 22
1n
Vn
Vn
Now let
x² vi
V. be the singular value decomposition of the matrix
x
the above f2-norm be further written
v
v
[y
Z.Mzl
xx_
20
M- 1/2
1
vn
n
2
[y Z. mm 1 Vo
12
Ia
O (n-d) O ((=-))
V 1n
0 (n
n
(ii)
My Z.M d,
n
where the equality (i) holds due to the definition empirical matrix Mes - x (a+)) the equality (ii) holds due to the unitary property of matrix V.. Combining the above two parts implies
B[(y Z_ T M
1/2
Σ
PB Z,arl yz M 1/2
f=1
2
T2 + Σ T
3=1
(111)
K
od 2 +d (olooo(//)) + 22lgg(K/T))
K [yy Z.Hom? d
< 2
+2
n
/=l
1=1
d+ 2 d log dlog(K//)) (/T) + log log(K///) K/T)
= 2([ Z. + K)
+ 2| ply - Z, Z. 12
18

---

with inequality (iii) holds with probability at least 1 - T. Still, it is easy to observe that for any positive constant c, as
d+ =+2/dlog(K/T) 2log(K/T)) Allpa - Z. 2.
n 2 n22 = max 4(1 + K)
Max 11
have |B Z, M 9 1| holds with probability at least
Step-3. Combining two upper bounds together, if n 2 max(n1, n2). the following inequality for generalization error bound
[ ((Z.)
[ Z. M= M M MLMMMMMM S 4e
up
holds with probability at least
3r.
We then study the compressed version.
Step-1'. Similarly, its optimal solutions for population loss and empirical loss are W. = 4 M ya M-1 and W oYX xx eлror bound is
+ M Mz1 respectively. Thus, the generalization
n
n
2 (W.)
W
W. W.ll.a.. |[( W., M. 2
Still, we have W. M Bl by generalization error bound can be upper-bounder by
M W.r) M |
1/2
and therefore, the
( W. W.) 1/2
M 1/2 M M M El W.æ)
op
Step-2'. Next, provide upper bounds these terms [ W. W.
M
and
M M
1/2
M 1/2
the right-hand-side separately. Note that for the matrix error term could use the same upper bounded mentioned in the proof of uncom- op
M M
1/2
M
1/2 mm
pressed version.
Now, to give the upper bound on the rescaled approximation error term E[( by W. M vector by W. which is
1/2
F we first compute the variance proxy for the subGaussian
J3
== 1| 11))0
for all j [m]. Thus, the rescaled approximation error for the compressed version be upper bounded by
m
E Ty W. M 1/2
Σ
E [ W.xljæ 1
1/2
F
=
re
j=1
< 2 Σ
E[ W.xl,
W. W.pww (
M 1/2
==1
+ 2
WaMtee M 1/2
==1
Still using Lemma 5 for all 7 E (0, 1), we have
(d 2 log(m//) llo((m/T))
Pn T <
21 1 T/m
19

---

For the term T7. following the same proof procedures of the uncompressed version implies
T3 fu W M
M
1/2
x:
1
P Buy W. Mm
Therefore, the rescaled approximation error for the compressed version is upper-bounded by
B(( - W.w) M
d 2 lloggm//) + log(m/T)) 2 Bp W. d
$ 2 113 + I|::)
n with probability at least Similarly, it is easy to get that for any positive constant
n
4( W. + 11=) d 2 /l/o(m///m + + llg(mm/T) A11 W. W. max
d
n n2 =
we have | W. M /3 S holds with probability least
T.
Step-3'. Combining two upper bounds together, if n 2 max A2). the following inequality for generalization error bound
ch C (W.) < [ - W.x)x M M=d/1
|M M M
S 4c
up
holds with probability at least - 3T.
A.4 Proof of Theorem
Recall the following theorem in the main text.
Theorem For any E (0, 1) and any E 103), suppose the compressed matrix follows
Assumption with > O log( )), and Assumption 2 holds. Given any learned regressor W
from training problem 2 let (æx, y) be new sample drawn from the underlying distribution D, we
have the following inequality holds with probability at least 1 -
4
EDlli v11
EDll By
Wa
1-8
where ÿ is the optimal solution from prediction problem with input vector
Proof. Due to the optimality of y. we have
Tÿ W x:
py
Wx
Φÿ by + by
W <
oy W W
<>
ΦΦ 4yl < 2 byy by
W x
b6 myl < Φ6 by 2
by Wx x
4
4yll, < Φy W.
2v ? < 4
py Wx
6)
1 v < by
W x
where the final = holds due to the (3x, 5) RIP property of the compressed matrix & with probability
least 1 т. Taking expectations both sides implies
(1 () v S 4F by
W xz
which completes the story.
20

---

A.5
Technical Lemma
A.5.1 Proof of Claim Proposed in Remark 2
Proof. Let us discuss the computational complexity for to be RK RR {0, 1] separately. Given fixed Ü,
F RK the projection method arg ney |/0 012 can be reformulate using the following mixed-ime mixed-integer programming (MIP),
(n., 2.) == arg min
(pp 0p E Ep=1 2p S
k
S.t.
with v, the output of the projection method. Sorting the absolute values (1 in decreasing order such that
أررقا 2 ... |80g)
the output v, of the proposed projection is
[v.l,
D
i E ((I)....,ss) (s)}
o.W.
with computational complexity O(K min log K}). F = R the projection method arg min DEVRRE | 812 can be reformulate using the following mixed-integer programming (MIP),
(v., z.) = arg min s.t.
ΣK k
p=1 Zr 0 Ep=1 2p S
Up 2 0 4 E [
Sorting { j.K in decreasing order such that
1(1) 2 ... 2...297 2((K):
the output u, of the proposed projection is
[v.li =
8g I(0r 0)
if {(1),..., (s))
O.W.
with computation complexity O(K min(s, log KJ). F = {0, 1)K, the projection method min || - ú12 presented in step of Algorithm can be represented
ming
Σ p= (1 - 22 + (ū 1)2
min,
Σ i - (2 - 1)
s.t.
Σ p= Zp S s
s.t. Σ 2p S s
k
2p 1) Ap [K]
2p E (0, V
[K]
Sort (20p - 1) K in decreasing order such that
20(1) > 28(2) 1 2......2.88 20(K) - 1,
then, the optimal z" can be set by
1(20p - -10 0) ifp (())...., (s)}
0
O.w.
For computational complexity, computing the sequence {2e, - 1 needs O(K), picking the top-s elements of the above sequence requires O(K min(s, log K}). setting the optimal solution z"
needs (s)), and thus the total computational complexity O(K) + +OK min(s, log K))+
O(K min(s, log K}).
21

---

A.5.2 Lemma for the Proof of Theorem
Lemma 1. (Unitary invariant). Let o E R be a randomized compressed matrix described in Assumption | and U E Rddd be a real unitary matrix. Then we have = &U is (1,8)-RIP with probability least 1-
Proof. Note that (i, component of & can be represented as P..j - U. - Ei=1 Pu.Ues Since every component P..j in & is i.i.d. drawn from (0, 1/m),, we have
P Σ Uus N 0.
Σ E=1
1 Uē. -
Now, we need to show that any two distinct components Tin.j Ji and P in + are independent. It is easy to observe that Pin.h and & Piz.da independent when i * ig since Di..: and $ ig.: independent. i1 = 12 = i, then the following random vector satisfies
4 Pich
Dr. U... ~ ~N (())))() (1/m 1/m))
Pinh
D.
That is to say, ii and + J2 ja are jointly Gaussian distributed and uncorrelated, which shows that +; in.h and Pig.a are independent. Combining the above together, we have & is randomized matrix with component i.i.d. from (0,, 1/m). Based on the existing result ((7), Theorem 1.5), when 2 Ci- (inee) + ln(2/r] for any d > 0 and (0, 1), we have & ensures (1, (1P8)- with probability least 1 - T.
A.5.3 Lemma for the Proof of Theorem 2
Lemma 2. For any integer parameter s(S d) and positive parameter 5 E (0, 1), let + Rmxd be a (s, R)-PIP matrix. For u1, u2 if u1, U2 u1 u2 U1 U2 are all s-sparse, then the following inequality holds
-25 un 12
|u2 (12) + 4(u1, u2) S 4 bu,,, Tu2) 25( |u212 + 4(u 42))
Proof. Since u1 U2 11 + u2 U2 are s-sparse, we have
(1 - () + uzll3) v (w + u2). 0(u1 + u2)) - (1+ [[ + +22133)
(8)
(1 - ( =22l²) < (D(u u2), D(ui u2)) < (1 + 0)(luu -uall3
(9)
Subtracting (9) from gives
-2 [ u |u2 (12) + 4(u,, u2) < 4(4pu,, 4u2) 25([|uu 11 ||u2)2) 4(un u2),
which completes the proof.
A.5.4 Lemma for the Proof of Theorem 3
Lemma 3. Let .'....,"
n be n i..dd o² oa-subGannnann random vectors with a zero mean and covariance matrix Σ. Then, there exists a positive constant C such that for all r E (0, 1),
(T))) n
{ d log(2/7) d+ log((/7))
)) 2 2-TT
P
Σ op
Co² max
Proof. Lemma is direct corollary from [Theorem 6.5, 41 Is easy observe that the proposed Lemma 3 holds by setting the parameter ő listed in [Theorem 6.5, [1 min(1,cy/nng]) (n(2/r)/n) with some positive constant.
Lemma Let ,!...., Th be n i..dd oussusaassia random vectors with a zero mean and covariance matrix E. Then for any E (0, 1), we have
n g'
40 40 20 log(1/7)
P
s
21-7.
n
n
2
22

---

Proof. We show this lemma by discretizing the unit 22 nomm ball (0; 1). Le12 be a 2 minimum cover of (0; 1) with its cardinality |N1/22 - 5 Since for any vector & Rd we always have
1
||El| = max (v, I max (v',E) + max (o",E) ||
max (v' + 2 v'EN1/2
max (v", .5),
o' N'EN1/2
|
then | S 2 maxx'eNi (v', E). Therefore, for any o'-subGaussian random vector,
max (0',)) S |W1//2 exp (- d'EN1/2
t²
F(IE 2 PP
So S 5d exp
S
which implies that || 40yd + 20 (lg((1/) with probability least 1 Now, since & = Σラ= E is a o² / /n-subGaussian random vector, inserting this variance proxy into the above inequality gives the desired result.
Lemma 5. Let n)(x)) be a zero-mean, sassussaussiann random variable. Let x ....., x. be n ii.d o'-subGaussian random vectors (may not zero-mean) as described in Assumption 2
Conditioned on
Mæ 1 -.=1x x²x))) > Odxd: for any T E (0, 1), we have
PP (1mma))
of + 2 d (1/(( + log(1/7))
2 -1-
n
Proof.
The proof of Lemma can be found in [Lemma 5, [8]].
A.5.5 Discussion after Theorem
Since [| - ulll S (| 5112 + 1 y/l2, combined with Theorem
we have
EDlllu(T) - - s 2F (T) ع + 2ED (1 y
(T)
8
< 00E Φÿ w + 1-8 ED By
Wa Walb
O(Epl D w
A.6
Additional Numerical Experiments on Synthetic Data
A.6.1
Implemented Prediction Method
The implemented prediction method is presented follows, see Algorithm 2 Algorithm proposed in the main content, it adds an additional stopping criteria
Comparing with
(1) 1((-21)|2 < 10
0.01 + 10
to ensure earlier stop than Algorithm
Note, in later numerical experiments, we use the terminology 'early stopping' denote that the iteration generated by the prediction algorithm satisfies the above additional stopping criteria within the maximum iteration number, i.e. T 60 (as listed in Section
A.6.2 Discussions on Baselines
Baselines. We compare our proposed prediction method with the following basclines.
Orthogonal Matching Pursuit.
Orthogonal Matching Pursuit(OMP) is greedy prediction algo- rithm. It iteratively chooses the most relevant output and then performs least-squares and updates the residuals. The built-in function OrthogonalMMtthhinnnuraaa package 'Sklearn. Leteearnn.inerr... Linear Linear_modell "_rdeemmo""" is used in the experiment.
from Python
Correlation Decoding. Correlation decoding multiplication of the transpose of compression matrix & and the learned regressor W. For any test point the algorithm predicts the top labels in & W orded by magnitude.
standard decoding algorithm. computes the
23

---

Algorithm 2 Implemented Projected Gradient Descent (for Second Stage)
Input: Regressor W, imput sample æ, stepsize 7), total iterations T
1: Initialize point (0)) E V. n F and = 0.
2: while and [|0(9) ( 3) (0.01 + ||v(0 > 10 30:3;0m do
3: Update ((++]) (()) Φ (Dr(l) Wx).
4: Project ((++)) [[(((+))) arg mignin 1
5: Update t == ++1.
6: end while
Output: v(T).
Elastic Net.
The elastic is regression method that combines both the é -norm penalty and the f2-norm penalty to guarantee the sparsity and the stability of the prediction. The parameters for the f-mnrrm penalty and --ror penalty in Elastic Net are be 0.1 through the experiments. The built-in function 'ElasticNet' from Python package Linerrrr.niier.....mm Lwwearnnennnne _model' is used in the experiment.
'Sklearn
Fast Iterative TirehshodiTn Algorithm. The fast iterative barinkage-thresholle al- gorithm is advanced optimization method designed to efficiently solve certain classes of unconstrained convex optimization problems. It utilizes momentum-lit like strategies to speed up convergence rates, which particularly effective for minimizing functions that may be non-smooth.
A.6.3 Procedures on Generating Mean & Covariance
In this paper, we give exact procedures on selecting Me and Enz mentioned in Section 4
The mean vector Pe is selected follows. We first generate d-dimensional vector μa from the standard d-dimensional Gaussian distribution (O,, Ia). Then we set
رلسي وامs for all [d]
by taking absolute values over all components.
The covariance matrix Drn is selected as follows. We first generate a d-by-d matrix A, where every component of A i.i.d. generated from the standard Gaussian distribution N/00, 1). Then the covariance matrix Eng is set to be
Men = d
A Ia.
2
A.7 Additional Numerical Experiments on Real Data
In this subsection, we do experiments real data and compare the performance of the proposed pre- diction method (see Algorithm 2 with four baselines, i.c., Orthogonal Matching Pursuit (OMP, (46)), Correlation Decoding (CD, 20 Elastic Net (EN, (47), and Fast Iterative Seiikkage-Thrssooldin
Algorithm (FISTA, 2 see Appendix A.6.2 for detailed explanations.
Real data. We select two benchmark datasets in multi-label classification, Wikil0-31K and EURLex-4KK due to their sparsity property. Table shows the details for the datasets.
Input Dim
Output Dim
Training set
Test set
Dataset EUKLex-4K
K
n
K
n1
K
5,000
3,993
17.413
236.69
5.30
1,935
240.96 659.65
5.32 19.03
Wiki10-31K
101,938
30,938
14.146
673.45
18.64
6,616
Table 1: Statistics and details for training and test sets, where d, K denote their averaged non-zero components for input and output, respectively.
24

---

Prediction time number of rows when = 3. 1/SMR
1.2
1.0
0.9
0.8
230
500
230
1200
1250
1520
1730
2005
Numteer of mos
The figure reports the prediction running time (measured in seconds) on synthetic data with early stopping by the proposed algorithm under different compressed output dimensions. As we can observe, the running time first decreases dramatically, then increases almost linearly with respect to m Such phenomenon has occurred since the max number of iterations is 60 in the implemented prediction method with early stopping, which relatively large; As m increases but is still less than 500, the actual number of iterations drops dramatically due to early stopping criteria; After passes 500, the actual number of iterations stays around 10, and then the running time grows linearly as dimension increases.
Figure 2:
Parameter setting. In prediction stage, we choose E 3) for EURLex-4K and 8 E {3,5) 5} for Wiki10-31K. We choose the number of rows m {100, 200, 300, 400, 500, 700, 1000, 2000) on both EURLex-4K and Wiki Vi-110-31. Ten independent trials of compressed matrices ((1) (10) are implemented for each tuple of parameters (s, m) on both datasets.
Empirical running time. Here, report the running time of the proposed algorithm and baselines on both synthetic and real datasets, see Table
Dataset Synthetic Data
Prop.Algo.
OMP
CD
Elastic Net
FISTA
~ =1 second
200-400 seconds
<| second
700-900 seconds
<3 seconds
EURLex-4K Wiki10-31K
<1 second
20-80 seconds
<1 second
2 1 second
seconds
500-700 seconds
seconds
5-10 seconds
Table 2: Time Complexity Comparison for each prediction
Numerical Results & Discussions. Figure 2 further illustrates that the computational complexity increases linearly with respect to the growth of compressed output dimension synthetic data, when m is greater than 500 to ensure the convergence of the prediction algorithm (see Remark 2
For real data, Figure 3 and 4 the results of their accuracy performances. In particular, the accuracy grows relatively stable with respect to m when the compression matrix satisfies the RIP-property with high probability. Besides, based on the results presented in Figure 3 Figure and Table 2 we observe that the proposed algorithm slightly outperforms the baselines on precision as increases while enjoys significant better computational efficiency, especially on large instances, which demonstrate the stability of the proposed algorithm.
25

---

Culpul dislance namber of alen 1. Eur les
Outguk distance mamber rows when 3, Eat
53
7.5
54
2.5
5.2
0.1
Ortpot
4.0
2000
2750
2000
15:00
3750
Precision number when Eurlex
Preciiion number rows whene Eurlex
0.773
0.65
0.5
0.55
0.30
Algo
CVP
F.STA
0.40
1000
3750
2000
15:00
Figure 3: This figure reports the numerical results on real data EURLex-4K. Each dot in the figure represents 10 independent trials (i.e.. experiments) of compressed matrices ((1), ...., $(10) on given tuple of parameters (s, m.. The curves in each panel correspond to the averaged values for the proposed Algorithm and baselines over 10 trials; the shaded parts represent the empirical standard deviations over 10 trials. In the first row, we plot the output distance versus the number of rows. In the second row, we plot the precision versus the number of rows, and we cannot observe significant differences between these prediction methods.
26

---

wild
wiki
Proc.Algo
0.40
ProcAAge
Serison 00 Ninter of maN
0.3
0.2
0.35
0.30
0.1
0.05
0.0
0.00
250
500
250
1000
1350
1500
1750
2002
2.50
500
750
1000
1750
1500
1750
2000
Number ctrons
Number CTrDRS
Output
wiki
22.0
Pop Aleo
28.0
Pop Aleo
Dauat berdroon
233
23.0
20.5
224
21.5
20.0
250
500
250
1000
1250
1500
1150
2000
230
500
250
1000
1250
15:00
17500
2000
Number ctrons
number
Preditiion
wiki
Pop Alge
Phop Algo
215
005
000
230
500
250
2000
1250
1300
37500
2000
2.0
300
750
1000
1250
1300
1150
2000
Mumber corows
Number c1rc05
Figure 4: This figure reports the numerical results on real data Wiki10-31K. Similar to the plot reporting EURLex-4K above, each dot in the figure represents 10 independent trials (i.e., experiments) of compressed matrices ((1), (10) curves in each panel correspond to the averaged values for the proposed algorithm and baselines over 10 trials; the shaded parts represent the empirical standard deviations over 10 trials. Similarly, in the first row, the precision of the proposed algorithm outperforms the FISTA especially when small. In the second & third rows for output difference and prediction loss, there only slight improvement on the proposed algorithm than CD of output difference.
given tuple of parameters (s, m.. The
27

---

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 paper's abstract is the short highlight of our introduction, covering the problems we study, the theoretical results establish, and the numerical experiments conduct.
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 in the paper and important assumptions and limitations. A 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 be expected to generalize other settings.
is fine include aspirational goals as motivation long 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: We mention assumptions discuss the limitation the proposed method, and compare method with other existing methods.
Guidelines:
The answer NA means that the paper has 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 out 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 authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect the scope of the claims made, e.g., if the approach was only tested on a few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low or images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play 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 of assumptions and complete (and correct) proof?
28

---

Answer: [Yes]
Justification: Assumptions and complete proofs are presented in the Appendix with short proof sketch in the main paper. All the theorems, used formulas, and proofs in the paper should be numbered and cross-referenced.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced the statement of any theorems.
The proofs can either appear 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 complemented by formal proof's provided in appendix 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 ex- perimental results of the paper to the extent that it affects the main claims and/or conclusions the paper (regardless of whether the code and data are provided or not)?
Answer: [Yes]
The numerical experiments are controlled by the default random seed, and it could be reproduced. Moreover, we not only provide the code but also include the detailed steps of the experiment in the main paper.
Justification:
Guidelines:
The answer NA means that the paper does not include experiments. If the paper includes experiments, 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 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 specific model and empirical evaluation, it may be necessary to either make possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way to accomplish this, but 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 large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurlPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make it clear how to reproduce that algorithm.
(b) the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) the contribution is a new model (e.g., a large language model), then there should either be way access this model for reproducing the results or way to reproduce 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 authors are welcome describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (c.g., to registered users), but it should be possible for other researchers to have some path to reproducing verifying the results.
29

---

reproducible, they 5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, described in supplemental material?
Answer: [Yes]
Justification:
We follow the guidance to provide the code and data.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurlPS code and data submission guidelines (https://nips.c cc/
public/guides
for more details.
While encourage the release of code and data, we understand that this might not be possible, so "No" is 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
7/nips cc/ wtc/publicc/ cooollccggiiesss gwt.///ieee//eee// Codeeoomiiiieeee /esd/guuuude//oeuubbiiiiii SanSSuubissss iolbbmmsssioo ssisnnoiccyyy for more details.
The authors should provide instructions on data access and preparation, including how to 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. If only subset of experiments 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 the paper) is recommended, but including URLs to data and code is permitted.
6. Experimental Setting/Details
Question: 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?
Answer: [Yes]
Justification:
We add the details in the numerical experiment section.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, as supplemental material.
7. Experiment Statistical Significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [Yes]
Justification: We provide the confidence interval for the results, which shows the significant improvement.
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, least for the experiments that support the main claims of the paper.
The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall with given experimental conditions).
30

---

The method for calculating the error bars should be explained (closed form formula, call 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 or the standard error of the mean.
is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma enror bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are range (e.g. negative enror rates).
If error bars reported in tables plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables 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 reproduce the experiments?
Question:
Answer: [Yes]
Justification:
The compute resources are mentioned the numerical section.
Guidelines:
The answer NA means that the paper does include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs 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 or failed experiments that didn't make into the paper).
9. Code Of Ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics
Answer: [Yes]
Justification: is anonymous.
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
If 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., there special consid- eration due to laws or regulations in their jurisdiction).
10. Broader Impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [Yes]
Justification: Our method aims solving challenges in the recent hype in large language model, for which hallucination is big issue. Besides, our paper is also can be used in the fairness considerations.
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 societal impact why the paper does not address societal impact.
31

---

Examples of negative societal impacts include potential malicious unintended uses (e.g., disinformation, generating fake profiles, surveillance), faireess considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific 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 direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that 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 out 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 intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms monitor how system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper deseribe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification: Our paper does not have this issue.
Guidelines:
The answer NA means that the paper poses no 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 requiring that users adhere to usage guidelines or restrictions to 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 not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
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?
Question:
Answer: [Yes]
Justification: We cite all papers that are used in our numerical experiments.
Guidelines:
The answer NA means that the paper does use existing assets.
The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset used and, 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 package should be provided. For popular datasets, paperswithcode com/ datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
32

---

For existing datasets that are re-packaged, both the original license and the license of the derived asset (if has changed) should be provided.
If this information is not available online, the authors are encouraged to reach out to the asset's s creators.
13. New Assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Yes]
Justification:
We provide anonymized zip file.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details the dataset/code/moodl 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 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. Crowdsoureing and Research with Human Subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification:
The paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included the main paper.
According the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the 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 equivalent approval/review based the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: The paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Depending the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to 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.
33


Reviews:
[{"review_id": "Lcf3ma9V23", "reviewer": "Reviewer_hQ2t", "summary": "This manuscript proposes a new approach the authors refer to as Sparse & High-dimensional-Output REgression (SHORE) to tackle linear regression problems promoting sparse predictions. A major component of the author's approach is to reduce the computational cost with a random normal matrix compressing the signal first and then producing sparse predictions in the second phase. The manuscript provides a theoretical analysis of the accuracy of the sparse prediction. The authors include further some numerical numerical investigations supporting their theoretical findings.", "questions": "See weaknesses.", "limitations": "yes.", "rating": 6, "confidence": 2, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "The manuscript is well-written and tackles an interesting problem.", "weaknesses": "Comparison to state-of-the-art methods such as FISTA (Fast Iterative Shrinkage-Thresholding Algorithm) or ADMM (Alternating Direction Method of Multipliers) would strengthen the contribution of this work. A clear comparison highlighting the advantages of the proposed method in terms of convergence speed, memory efficiency, or specific problem suitability would be valuable. Additionally, the paper could benefit from exploring connections to randomized linear algebra, particularly regarding the selection of the matrix Φ. The paper Randomized Numerical Linear Algebra: Foundations & Algorithms by Martinsson and Tropp provides a comprehensive overview of how random projections can be leveraged to achieve efficient solutions in high-dimensional problems. Further, Equation (1) is a simple linear regression or least squares problem  ||Y-ZX||_F^2 =  ||\\vec(Y) - (X'  \\kron  I) vec(Z)||_2^2 which does not require a new renaming as SHORE.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "## Response \nThank you for your positive feedback! Figure1 mentioned below can be found in the uploaded pdf file in 'global' author rebuttal.\n\n**[Response to FISTA or ADMM]** Thank you for your comments. We have conducted some preliminary numerical experiments concerning comparison with FISTA, as you suggested. Please check the first panel in Figure 1 in the global response. Besides, it is difficult for FISTA to achieve the stopping criteria in Appendix 6.1 on page 21. We suspect that this is due to the special structure of our proposed model versus the general SCO.\n\n\n**[Response to other randomized linear algebra]** Thank you for providing these references. Actually, as you mentioned above, we realized that how to select a projection matrix is an important open problem that may significantly influence the accuracy of the final predicted output in both theory and practice (also see Blum et al. (2005), Zhang et al. (2015), Pillai et al. (2011)). We will add this part to Section 5 as one possible future direction.\n\n\n**[Response to SHORE name issue]** Thank you for raising this confusion. We will add a discussion to the camera-ready version based on the following explanations. We agree that Equation (1) for the training stage can be reformulated as a linear regression, as you presented. However, the SHORE framework we proposed here refers to a two-stage framework, which includes both the training stage (Equation (2), compressed version of Equation (1)) and the prediction stage (Equation (3)). Moreover, we would like to mention that vectorizing the output matrix $\\boldsymbol{Y}$ and linear regressor $\\boldsymbol{Z}$ might impede its ability to extend to other structured settings. For example, in some specific applications, the linear regressor $\\boldsymbol{Z}$ is assumed to be low-rank, which is more natural to formulate $\\boldsymbol{Z}$ as a matrix rather than a vector.\n\n## Reference\n**[Blum et al. (2005)]** \nAvrim Blum. Random projection, margins, kernels, and feature-selection. In International Statistical and Optimization Perspectives Workshop” Subspace, Latent Structure and Feature Selection”, pages 52–68. Springer, 2005.\n\n**[Zhang et al. (2015)]**\nShengping Zhang, Huiyu Zhou, Feng Jiang, and Xuelong Li. Robust visual tracking using\nstructurally random projection and weighted least squares. IEEE Transactions on Circuits\nand Systems for Video Technology, 25(11):1749–1760, 2015.\n\n**[Pillai et al. (2011)]**\nJaishanker K Pillai, Vishal M Patel, Rama Chellappa, and Nalini K Ratha. Secure and robust\niris recognition using random projections and sparse representations. IEEE transactions on\npattern analysis and machine intelligence, 33(9):1877–1893, 2011.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "bJivVU0MMP", "reviewer": "Reviewer_vYuH", "summary": "In this paper, the author has proposed a new method for solving the high-dimensional output regression (MOR) problem through compression. The authors introduce a two-stage framework that incorporates output compression to achieve computational efficiency while maintaining accuracy. Theoretical results demonstrate the framework's scalability and error bounds, while empirical results validate its performance.", "questions": "- [The influence on the compressed hyperparameter m] With the increase of the hyperparameter $m$ as the compressed dimension. The model performance of the methods has been significantly improved. Therefore, I am wondering if the model training time is depending on m, so it would be interesting to see how the training time increase in real datasets as the hyperparameter m change since the current training time is just a range, while the computational complexity O(KM) in the second stage, so the training time wshould definitely increase with the improvement of m.\n\n- [The selection of m] I can understand there is a trade-off between model performance and compressed output dimension. However, if the dataset is totally unknown, the training time would increase under the trial. Therefore, I am wondering if that is also an issue when this method is finally applied for real application.", "limitations": "N/A", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "- The author clearly defined the question about the MOR, and the proposed methods that trying to achieve the good performance.\n- The proof and the structure of the paper is well-written and well-stated, containing comprehensive theoretical analyses, including training loss bounds and convergence guarantees.", "weaknesses": "- [Comparison Methods] In this paper, the authors has compare the model performance with some existing prediction baselines like Orthogonal Matching Pursuit and Elastic Net. Both of those two methods seems not to be the latest sparse regression methods that would be interesting to compare with. Therefore, I am wondering if some methods related to sparse high dimension regression like [1,2,3,4] can be discussed within the paper, and tested on the model performances. Also, there are plenty of the methods that are mentioned within the literature review parts, those methods are also not taken into considerations instead of elasticnet, which make the conclusion less convincing from the numerical results perspective.\n\n- [Missing of the Sparsity levels] with the synthetic dataset. The paper doesn't mentioned the sparsity level while the sparsity level is pre-defined within the numerical experimental settings since the set up sparsity level tends to be the most correct baseline for evaluating the sparsity level of the model, so it would be really important to be shown within the paper.\n\n[1] Bertsimas, D., and B. Van Parys. \"Sparse high-dimensional regression: Exact scalable algorithms and phase transitions (2017).\" arXiv preprint arXiv:1709.10029 (2019).\n[2] Bertsimas, Dimitris, and Bart Van Parys. \"SPARSE HIGH-DIMENSIONAL REGRESSION.\" The Annals of Statistics 48.1 (2020): 300-323.\n[3] Liu, Liu, et al. \"High dimensional robust sparse regression.\" International Conference on Artificial Intelligence and Statistics. PMLR, 2020.\n[4] Sun, Qiang, et al. \"Sprem: sparse projection regression model for high-dimensional linear regression.\" Journal of the American Statistical Association 110.509 (2015): 289-302.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "## Response\nWe sincerely thank the reviewer for their positive assessment of the paper's clarity, methodology, and theoretical contributions. Figure1 mentioned below can be found in the uploaded pdf file in 'global' author rebuttal.\n\n**[Response to Comparison with Sparse Regression and related references]**\n\nThank you for these suggestions on related literature and possible baselines. We will add the following comparisons and literature review to the revised version. Additionally, some numerical experiments will be conducted as suggested.\n\n**Proposed prediction stage vs. Sparse regression:**  Please refer to the 'global' author rebuttal part \"Proposed prediction stage vs. Sparse regression\" for detailed discussions. \n\nWe will also add the following literature reviews (Bertsimas et al. 2019, Bertsimas et al. 2020, Liu et al. 2020, Sun et al. 2015) to Section 3 in our camera-ready version:\n\n- **[Bertsimas et al. 2019, Bertsimas et al. 2020]:** This work proposes a binary convex reformulation (BCR) of the sparse regression with an additional $\\ell_2$ squared regularizer and devises a cutting plane method for solving the proposed BCR. The authors of this work demonstrate its scalability under their sample-generating model. Compared with our result, the proposed cutting-plane method lacks an iterative convergence guarantee, even though such a method ensures the finding of a globally optimal solution with exponential time complexity.\n\n- **[Liu et al. 2020]:** This paper proposes an algorithm for high-dimensional sparse regression with a constant fraction of corruptions, which has a different problem setting compared with ours.\n\n- **[Sun et al. 2015]:** This paper develops a sparse projection regression modeling (SPReM) framework to perform multivariate regression, focusing on addressing the low statistical power issue of many standard statistical approaches. However, the authors of Sun et al. 2015 do not provide any iterative convergence result for their proposed coordinate descent algorithm (Algorithm 1).\n\n\n\n**[Response to Missing of Sparsity levels]**\n\nThank you for your feedback regarding the importance of including the sparsity levels in the discussion of our synthetic dataset experiments. We would like to point out that the sparsity levels are indeed mentioned on line 298 on page 8, where $s = 3$. However, we understand that this might not have been sufficiently emphasized. To address this, we will make the description more explicit and highlight it in figure captions and subplot titles.\n\n\n\n**[Response to the influence on the compressed output dimension $m$: model performance]**\n\nHere, the performances are relatively stable as $m$ is greater than some worst-case lower bound. Typically, the worst-case lower bound is a fixed number much smaller than the label number $K$ given the training dataset (see response on \"The selection of $m$\" below). In our numerical experiment, due to resource limitations, we only conducted several experiments on relatively small real datasets where $K \\approx 31,000$. Thus, the ratio of the worst-case lower bound on compression number divided by label numbers $K$ is high compared to large datasets, leading to the result that the accuracy of experiments significantly improves as $m$ increases. Now, we have tested on relatively larger $m$; see the third panel in Figure 1 and the Real-world dataset part in the global author rebuttal for detailed discussions.\n\n**[Response to the influence on the compressed output dimension $m$: computational complexity]**\n\nThank you for pointing this out. Yes, the computational complexity will increase with the improvement of $m$. As noted above, the accuracy performance is relatively stable as $m$ is greater than some worst-case lower bound. Thus, there is no strong motivation to increase hyperparameter $m$ beyond such a worst-case lower bound, which leads to an \"upper bound\" on the computational complexity. See the second panel in Figure 1 in the global rebuttal for details.\n\n**[Response to the selection of $m$]**\n\nThank you for your feedback regarding the selection of $m$ for real applications without prior knowledge. As noted in Remark 3 (lines 182-184), we provide a theoretical worst-case lower bound on $m$ with high probability ensuring that the compressed matrix $\\Phi$ satisfies the RIP-property. More specifically, the worst-case bound is given by $m \\ge \\dfrac{12}{\\delta^2}(3s \\ln{9} + 3s\\ln{\\dfrac{e|\\mathcal{V}|}{3s}}+\\ln{\\dfrac{2}{\\tau}})$ where $s$ represents the sparsity level. As implemented in our numerical experiments, we use the worst-case bound above and repeat the experiments ten times to reduce the impact of randomness. In real applications, the dimension of the dataset is typically known, and the sparsity level $s$ can be chosen close to the average of non-zero labels in the dataset, allowing us to compute the worst-case lower bound of $m$​ accordingly.\n\n\n\n## Reference\n\n**[Bertsimas et al. 2019]**\nD Bertsimas and B Van Parys. Sparse high-dimensional regression: Exact scalable algorithms\nand phase transitions (2017). arXiv preprint arXiv:1709.10029, 2019.\n\n**[Bertsimas et al. 2020]**\nDimitris Bertsimas and Bart Van Parys. Sparse high-dimensional regression. The Annals of\nStatistics, 48(1):300–323, 2020.\n\n**[Liu et al. 2020]**\nLiu Liu, Yanyao Shen, Tianyang Li, and Constantine Caramanis. High dimensional robust\nsparse regression. In International Conference on Artificial Intelligence and Statistics, pages\n411–421. PMLR, 2020.\n\n**[Sun et al. 2015]**\nQiang Sun, Hongtu Zhu, Yufeng Liu, and Joseph G Ibrahim. Sprem: sparse projection re-\ngression model for high-dimensional linear regression. Journal of the American Statistical\nAssociation, 110(509):289–302, 2015.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "qsYVB9G6YN", "reviewer": "Reviewer_PtcT", "summary": "The paper proposes the Sparse & High-dimensional-Output REgression (SHORE) model to address challenges in Multi-Output Regression (MOR). SHORE introduces a two-stage framework to handle high-dimensional outputs efficiently. Theoretical analysis shows that the framework maintains training and prediction loss while being computationally scalable. Empirical results validate these findings, demonstrating the framework’s efficiency and accuracy in handling SHORE for modern MOR applications.", "questions": "- Theorem 4: $\\hat{y}$ in the theorem is the optimal solution of (3). How about the algorithmic solution?\n\n- The selection $\\eta$ is much easier than in [2]. However, it would be good to provide some default values for practical usage. Furthermore, is it possible to select $\\eta$ in a more data-driven way, for example, by leveraging the value of the loss function like in [4]?", "limitations": "The authors have discussed the limitations in this paper.\n\n### Reference\n\n- [1] Blumensath, Thomas, and Mike E. Davies. \"Iterative hard thresholding for compressed sensing.\" Applied and computational harmonic analysis 27.3 (2009): 265-274.\n\n- [2] Jain, Prateek, Ambuj Tewari, and Purushottam Kar. \"On iterative hard thresholding methods for high-dimensional m-estimation.\" Advances in neural information processing systems 27 (2014).\n\n- [3] Hastie, Trevor, Robert Tibshirani, and Ryan Tibshirani. \"Best subset, forward stepwise or lasso? Analysis and recommendations based on extensive comparisons.\" Statistical Science 35.4 (2020): 579-592.\n\n- [4] Wang, Zezhi, et al. \"Sparsity-Constraint Optimization via Splicing Iteration.\" arXiv preprint arXiv:2406.12017 (2024).", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- Overall, this paper is well-written and well-organized.\n\n- The two-stage algorithm proposed in the paper is easy to grasp and implement. Notably, by introducing a random projection matrix, it reduces the dimension of $Y$ and largely relaxes the RIP-type condition in sparsity-constrained optimization.\n\n- The theoretical analysis is solid and includes a high-level interpretation of the proof. The imposed assumptions are mild and likely to (approximately) hold in practice.", "weaknesses": "- The prediction stage is significantly different from the classical manner (i.e., $Z X$). It involves solving sparsity-constrained optimization (SCO) problems. It would be good to highlight this difference. Additionally, solving SCO is more computationally intensive than the classical method, so more discussion (in both theory and experiments) is necessary.\n\n- Algorithm 1: Algorithm 1: Steps 3-4 resemble iterative hard thresholding [1, 2]. It would be better to discuss this in the appropriate place. \n\n- Lines 217-218: The setting for $y$ seems to imply that only $s$ columns have non-zero values while the remaining columns always take zero values. Under this assumption, why not filter out some columns according to the magnitude of the columns of $y$?\n\n- It would be better to explicitly describe how OMP is compared. Is it used for solving (3) or (1)?\n\n- Line 310-311: Elastic net may have better performance in a low signal regime [3].\n\n- Considering these assumptions are mild in practice, why not test the algorithm on real-world datasets and present the results in the main text?\n\n- Lines 437, 467, 542: Incomplete sentences. Perhaps they can be changed to \"Recall the following theorem in the main text.\"", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "## Response\n\nWe sincerely thank the reviewer for their positive and encouraging feedback! Table 1 mentioned below can be found in the uploaded pdf file in 'global' author rebuttal.\n\n**[Response to Comparison with SCO]**\n\nThank you for highlighting this aspect! If we understand correctly, the classical manner you mentioned denotes predicting SHORE *without compression*. Please refer to the 'global' author rebuttal part ''Proposed prediction stage v.s. General SCO'' and ''Proposed prediction stage v.s. Classical manner'' for detailed discussions. We will add these discussions in the introduction and problem setup sections.\n\n**[Response to Comparison with iterative hard thresholding]**\n\nThank you for your suggestion. We will add the following discussion to Section 3.1 after introducing the proposed two-stage framework SHORE. Note that both the iterative hard thresholding methods (Blumensath et al. 2009, Jain et al. 2014) and the proposed Algorithm 1 share the same insight from the vanilla projected gradient method for nonconvex constrained minimization.\nCompared with existing iterative hard thresholding, Step 4 (projection step) projects the resulting vector to a more general feasible set $\\mathcal{V}_s^K \\cap \\mathcal{F}$ with additional constraints other than sparsity, which might require non-trivial modifications in designing projection oracles and convex analysis. In particular, we discuss the different projection methods under different $\\mathcal{F}$ through the perspective of optimization (see Appendix A.5.1 on page 19). We also discuss the bound for stepsize $\\eta$ to guarantee the algorithm's convergence.\n\n**[Response to nonzero of y]**\n\nThank you for raising this ambiguity! We only assume that each column of $\\boldsymbol{Y}$, i.e., $\\boldsymbol{y}^{i}$ for $i = 1, \\ldots, n$, has $s$ non-zero components. However, the support sets of two distinct samples, e.g., $\\boldsymbol{y}^{i}$ and $\\boldsymbol{y}^{i'}$ with $i \\neq i'$, might be different, which does not imply the row-sparsity for matrix $\\boldsymbol{Y}$ (i.e., only $s$ rows/labels are non-zero). \n\n**[Response to Comparison with OMP]**\n\nThank you for pointing this out! We will add an explicit description of both the OMP algorithm and how we compare the OMP algorithm with our algorithm in the camera-ready version.\nHere, we would like to point out that OMP is an algorithm to solve the sparse regression problem, which has some differences from the proposed problem in the prediction stage. See the global rebuttal part for details.\n\n**[Response to Comparison with Elastic Net]**\n\nThank you for your insightful observation! We conducted several small experiments and found that Elastic Net works better in a low SNR regime (Refer to Table 1). We will add relevant comparisons in the revision.\n\n**[Response to testing on real-world dataset]**\n\nThank you for identifying this point. Because of resource limitations, we only conducted several experiments on relatively small real datasets where  $\\text{K} \\approx 31,000$ , which are shown in Appendix 7. Currently, we are preparing more computational resources to conduct experiments on larger real-world datasets.\n\n**[Response to typos]**\n\nThank you for pointing this out! We will correct them in the revision and do more proofreading.\n\n**[Response to $\\hat{y}$ and algorithmic solution ]**\n\nThank you for asking this! We actually discussed the algorithmic solution, denoted by $\\boldsymbol{v}^{(T)}$​, a bit below Theorem 4, see lines 273 - 275. We will highlight the results in a remark in the revision.\n\n**[Response to selection of $\\eta $]**\n\nIn our numerical experiments, the default setting for $\\eta$ is $\\eta = 0.9$, which satisfies the convergence condition proposed in Theorem 2 and Theorem 4 for all $\\delta \\in (0, \\frac{4}{9})$.\n\nFor a data-driven setting on $\\eta$, we realize that (Wang et al. 2024) was recently accessible on arXiv after the NeurIPS 2024's submission deadline, which introduces a novel iterative algorithm with a tuning-free property and linear convergence guarantee. We believe that this is a good related future direction and will be discussed in Section 5.\n\n\n\n## Reference\n\n**[Blumensath et al. 2009]**\nThomas Blumensath and Mike E Davies. Iterative hard thresholding for compressed sensing.\nApplied and computational harmonic analysis, 27(3):265–274, 2009.\n\n**[Jain et al. 2014]**\nPrateek Jain, Ambuj Tewari, and Purushottam Kar. On iterative hard thresholding methods\nfor high-dimensional m-estimation. Advances in neural information processing systems, 27,2014.\n\n**[Wang et al. 2024]**\nZezhi Wang, Jin Zhu, Junxian Zhu, Borui Tang, Hongmei Lin, and Xueqin Wang. Sparsity-\nconstraint optimization via splicing iteration. arXiv preprint arXiv:2406.12017, 2024.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:
## Responses to All Reviewers & Area Chairs

We would like to thank the reviewers for their constructive and high-quality feedback. The manuscript has been revised based on the comments given in three reviewers' reports.

### Comparisons & Differences

We first highlight the differences between our proposed prediction stage, general sparsity-constrained optimization (SCO), sparse regression, and the classical manner (i.e., prediction stage of SHORE without compression). The following discussions will be added to the revised version.

- **Proposed prediction stage vs. General SCO:** The prediction stage requests solving a sparsity-constrained optimization (SCO) problem with its general form: $\min_{||\boldsymbol{\alpha}||_0 \leq k} ||\boldsymbol{\beta} - \boldsymbol{A}\boldsymbol{\alpha}||_2^2$. Finding the global optimal solution of SCO, in general (worst-case scenario), is computationally hard. In contrast with the general SCO, the problem we studied in the predicted stage takes the random projection matrix $\boldsymbol{\Phi}$ with restricted isometry property (RIP) as its $\boldsymbol{A}$ and uses $\widehat{\boldsymbol{W}}\boldsymbol{x}$ with $\widehat{\boldsymbol{W}}$ obtained from the compressed training stage as its $\boldsymbol{\beta}$. As a result, and as presented in Theorem 2 (page 5) and Theorem 4 (page 7), the proposed efficient algorithm (Algorithm 1 -- PGD) ensures a globally linear convergence to a ball with center $\widehat{\boldsymbol{y}}$ (optimal solution of the prediction stage) and radius $O(||\boldsymbol{\Phi}\widehat{\boldsymbol{y}} - \widehat{\boldsymbol{W}}\boldsymbol{x}||_2)$, which might not be true for general SCO problems.

- **Proposed prediction stage vs. Sparse regression:** Although both the proposed prediction stage and sparse high-dimensional regression share a similar optimization formulation: 
  $\min_{||\boldsymbol{\beta}||_0 \leq k} ~~ ||\boldsymbol{y} - \boldsymbol{X}^{\top}\boldsymbol{\beta}||_2^2,$
  the proposed prediction stage (equation (3) under line 149 on page 4) is distinct from sparse regression in the following parts:
  
  - **Underlying Model:** Most existing works on sparse high-dimensional regression assume that samples are i.i.d. generated from the linear relationship $\boldsymbol{y} = \boldsymbol{X}^{\top}\boldsymbol{\beta}^* + \boldsymbol{\epsilon}$ with underlying sparse ground truth $\boldsymbol{\beta}^*$. In the proposed prediction stage, we do not assume additional underlying models on samples if there is no further specific assumption. The problem we studied in the prediction stage takes the random projection matrix $\boldsymbol{\Phi}$ with RIP as its $\boldsymbol{X}^\top$ (whereas $\boldsymbol{X}^{\top}$ in sparse regression does not ensure RIP) and uses $\widehat{\boldsymbol{W}}\boldsymbol{x}$ with $\widehat{\boldsymbol{W}}$ obtained from the compressed training stage as its $\boldsymbol{Y}$.
  - **Problem Task:** Sparse regression aims to recover the sparse ground truth $\boldsymbol{\beta}^*$ given a sample set $\{(\boldsymbol{x}^i, \boldsymbol{y}^i)\}_{i = 1}^n$ with $n$ i.i.d. samples. In contrast, the task of the proposed prediction stage is to predict a sparse high-dimensional output $\widehat{\boldsymbol{y}}$ given a random projection matrix $\boldsymbol{\Phi}$ and a single input $\boldsymbol{x}$.
 
  Therefore, due to distinct underlying models and problem tasks, some existing methods/results for sparse regression cannot be directly applied to the proposed prediction stage.
- **Proposed prediction stage vs. Classical manner:** As mentioned in Remark 1 & Remark 2 on page 4, the compressed SHORE (both training stage & prediction stage) enjoys better computational complexity with respect to the classical manner, especially under the setting with high-dimensional outputs.

### Numerical Experiments
We then report some preliminary results of numerical experiments as suggested by reviewers (see Figure 1 and Table 1). Other numerical experiments, including relatively sophisticated optimization methods or large real instances, will be reported later in the camera-ready version due to resource and time limitations.
- **Comparing with FISTA:** The first panel in Figure 1 compares the proposed algorithm (Algorithm 1 -- PGD) with FISTA. As we can observe, the precision of the proposed algorithm outperforms the precision of FISTA, especially when $m$ is small (heavily compressed).
- **Prediction time:** The second panel in Figure 1 reports the prediction running time (measured in seconds) on solving one implemented prediction method with early stopping (see lines 618-620, Section A.6.1 on page 21) by PGD under different compressed output dimensions $m$. As we can observe, the running time first decreases dramatically, then increases almost linearly with respect to $m$. Such a phenomenon has occurred since the max number of iterations $T$ is 30 in the implemented prediction method with early stopping, which is relatively large; As $m$ increases but is still less than 1000, the actual number of iterations drops dramatically due to early stopping criteria; After $m$ passes 1000, the actual number of iterations stays at 5, and then the running time grows linearly as dimension $m$ increases.
- **Real-world dataset:** The third panel in Figure 1 shows that for a real-world dataset with 31,000 labels, the model performance of our proposed algorithm becomes stable as $m$ is greater than some worst-case lower bound.
- **Comparing with elastic net under low SNR regime:** As observed from Table 1, the elastic net has similar performances on precision compared with the performances of the proposed PGD under a low SNR regime. It's true that the elastic net performs much better in numerical experiments in low SNR regime than in high SNR regime. We would like to point out that the elastic net method performs much slower than PGD. The prediction time for the elastic net is around 1,000 seconds, while the time for the PGD is around 10 seconds.

Meta-review: The authors propose in this article a two-stage approach to tackle multi-output regression under a pre-determined sparsity level on the output vector and provide performance guarantees for this computational efficient method. Several points raised by the reviewers regarding the comparison to other methods, the experimental study and the clarity were effectively addressed by the authors' rebuttal, leading to a consensus to accept the paper, with recommendation to incorporate the discussion and the additional results into the final version.
Decision: Accept (poster)