ReHLine: Regularized Composite ReLU-ReHU Loss Minimization with Linear Computation and Linear Convergence
Authors: Ben Dai
Yixuan Qiu
Venue: NeurIPS 2023 poster (2023)
Keywords: coordinate descent
linear convergence
primal-dual methods
empirical risk minimization
linear constraints
quantile regression
TL;DR: A novel algorithm for empirical risk minimization that significantly surpasses generic solvers and also outperforms specialized solvers including Liblinear, hqreg, and Lightning.

Abstract: Empirical risk minimization (ERM) is a crucial framework that offers a general approach to handling a broad range of machine learning tasks. In this paper, we propose a novel algorithm, called ReHLine, for minimizing a set of regularized ERMs with convex piecewise linear-quadratic loss functions and optional linear constraints. The proposed algorithm can effectively handle diverse combinations of loss functions, regularization, and constraints, making it particularly well-suited for complex domain-specific problems. Examples of such problems include FairSVM, elastic net regularized quantile regression, Huber minimization, etc. In addition, ReHLine enjoys a provable linear convergence rate and exhibits a per-iteration computational complexity that scales linearly with the sample size. The algorithm is implemented with both Python and R interfaces, and its performance is benchmarked on various tasks and datasets. Our experimental results demonstrate that ReHLine significantly surpasses generic optimization solvers in terms of computational efficiency on large-scale datasets. Moreover, it also outperforms specialized solvers such as Liblinear in SVMs, hqreg in Huber minimization, and Lightning (SAGA, SAG, SDCA, SVRG) in smoothed SVMs, exhibiting exceptional flexibility and efficiency. The source code, project page, accompanying software, and the Python/R interface can be accessed through the link: https://github.com/softmin/ReHLine.

ReHLine: Regularized Composite ReLU-ReHU Loss Minimization with Linear Computation and Linear Convergence
Ben Dai
Yixuan Qiu
Department of Statistics The Chinese University of Hong Kong bendai @cuhk. edu. hk
School of Statistics and Management Shanghai University of Finance and Economics qiuyixuan@sufe edu.
Abstract
Empirical risk minimization (ERM) a crucial framework that offers general approach handling broad range of machine learning tasks. In this paper, we propose a novel algorithm, called ReHLine, for minimizing a set of regular- ized ERMs with convex piecewise linear-quadratic loss functions and optional linear constraints. The proposed algorithm can effectively handle diverse combi- nations of loss functions, regularization, and constraints, making it particularly well-suited for complex domain domain-specific problems. Examples such problems include Va, elastic net regularized quantile regression, Huber minimization, etc. In addition, ReHL Line enjoys a provable linear convergence rate and exhibits a per-iteration computational complexity that scales linearly with the sample size. The algorithm implemented with both Python and interfaces, and its perfor- mance benchmarked on various tasks and datasets. Our experimental results demonstrate that ReHLine significantly surpasses generic optimization solvers in terms of computational efficiency on large-scale datasets. Moreover, also outperforms specialized solvers such as LIBLINEAR in SVMs, IIQREG in Huber minimization and LIGHTNING (SAGA, SAG, SDCA, SVRG) smoothed SVMs, exhibiting exceptional flexibility and efficiency. The source code, project page, accompanying software, and the Python/R interface can be accessed through the link: https: //gi thub. cow.soocccooooooooo tmi n/ ReHL. ine
Introduction
Empirical risk minimization (ERM) with different losses is a fundamental framework in statistical machine learning [43) which has demonstrated significant impact on many application areas, such as artificial intelligence, computational biology, computer vision, natural language processing, electronic engineering, and finance. In fact, many common practical ERMs are associated with a piecewise linear- quadrattc (PLQ; (34)) loss function, including support vector machines (SVMs; (7]. the least absolute deviation regression, quantile regression (26.. Huber regression (22). borrowing cost function in portfolio selection (28), and pricing/planning problems in electricity and operational research 16 1. among many others. In addition, ERMs are frequently blended with linear constraints, for instance, the fairness constraints in fair classification [48] monotonicity constraints in classification or regression, etc. More linear constraints according to the prior domain knowledge, such as the sum-to-one constraints in the portfolio variance minimization problem [24] negative positive constraints in por-orrfonnnng loan (36), are considered various real applications. Therefore, has become of utmost importance to develop fast and scalable numerical algorithms to solve ER Ms with PLQ losses based on large-scale datasets and flcxible linear constraints.
"Co-first authorship, and the order is determined by coin toss. Corresponniigg author.
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

In this paper, we consider general regularized ERM based on a convex PLQ loss with linear constraints:
Σ (x(x - 11 s.t. - 2200200 0,
min BER
where x; E / is the feature vector for the i-th observation, 3 E F is an unknown coefficient vector, A E IR 1gx and RK are defined as linear inequality constraints for B, and Li(-) 2 0 a convex PLQ loss function. The PLQ function class greatly generalizes existing popular loss functions, and the convexity L,(-) guarantees the global convergence of optimization. Here, we focus on working with a large-scale dataset, where the dimension of the coefficient vector and the total number of constraints are comparatively much smaller than the sample size, that is, d n and K < n.
Although (T) is a strongly convex problem with affine constraints, thus admitting a unique global optimum, directly solving ( by generic solver may encounter various challenges. For example, L.(-) is in general non-smooth, so gradient-bass based methods may fail, and hence slow-convergent subgradient methods need be used. In addition, the constraint is polyhedron, whose projection operator is non-trivial to compute, making various projection-basse methods difficult to apply.
Of course, for specific forms of Lz(-). highly efficient solvers for ①I indeed exist. One remarkable example the LIBLINEAR solver [12), which has gained great success in solving large- scale SVMs. However, LIBLINEAR highly relies the hinge loss function, so its success does not directly transfer to more general loss functions.
Outline. To this end, in this article we develop a novel algorithm to solve ( with the greatest generality. The proposed algorithm, named ReHLine, has three key ingredients. First, show that any convex PLQ loss function can be decomposed as the sum of a finite number of rectified linear units (ReLU, (15)) and rectified Huber units (ReHU, defined in Section ( Second, based on this decomposition, the dual problem of (T) is shown to be a box-constrained quadratic programming (hox-QP) problem. Finally, special paired primal-dual coordinate descent algorithm developed to
simultancously solve the primal and dual problems, with a computational cost linear in n.
We emphasize that for the second point, the dual box-QP problem is highly structured. Therefore, general-pe purpaal-ppppeeeeeeee solver may not fully exploit the structure, thus leading to slower convergence rates higher computational costs. As direct comparison between the proposed ReHLine solver and existing algorithms, including the projected (sub- )orodderadnnttt descent (P-GD; ( the coordinate descent (CD; (4))), the interior-poii methods (IPMs; [I3
14 18),, the stochastic dual coordinate ascent (SDCA: 39 40 related methods, and the alternating direction method of multipliers (ADMM; (5)), Table summarizes the required number of iterations and per-iteration cost of each method. More details and discussions for this comparison are given in Section
Table Overview existing algorithms for solving ① Column COMPLEXITY (PER ITER ATION
shows the computational complexity of the algorithm per iteration. Here, we focus only the order of n since d & n assumed in our setting. Column #ITERATIONS shows the number of iterations needed to achieve accuracy of to the minimizer.
ALGORITHM
COMPLEXITY
#ITERATIONS
COMPLEXITY (TOTAL)
(PER ITERATION)
P-GD
((n)
(
6
O(ne-1) O O log(s-1)) O (n2 log(e
CD
((n²)
Olog(
31
IPM
(nn))
Ollog(
18
log(e-1)
ADMM
((n2)
o 20
o ((---1)
SDCA
((n)
( 39
Onn
ReHLine (ours)
O(n)
(log(E
On (!))
Contribution.
Compared with other specialized ERM solvers or general-purposs box-QP solvers, the proposed ReHLine solver has four appealing "linear properties".

---

1. It applies to any convex piecewise linear linear-quadradn quadratic loss function (potential for non -on-smoothness included), including the hinge loss, the check loss, the Huber loss, etc.
2. In addition, supports linear equality and inequality constraints on the parameter vector.
3. The optimization algorithm has a provable linear convergence rate.
4. The per-iteration computational complexity is linear in the sample size.
Moreover, ReHLine is designed to be a computationally efficient and practically useful software package for large-scale ERM problems. As an example, experiments in Section have shown that ReHL significantly surpasses generic solvers in terms of computational efficiency on large-scale datasets, and provides reasonable improvements
specialized solvers and problems, such as LIBLINEAR, HQREG and LIGHTNING, given that LIBLINEAR, HQREG and LIGHTNING are one of the fastest and most widely-used solvers for SVM, Huber minimization, and smoothed SVM, respectively, Finally, ReHLine is available as user-friendly software in both Python and R, which facilitates ease of use for researchers, practitioners, and wider audience.
2 The ReLU-ReHU decomposition
As introduced in Section | we attempt to solve ERM with a general loss function, and one sensible choice is the class of convex PLQ functions. A univariate loss function f R R>o is PLQ if there exist a finite number of knots 1 < 12 < .... < tx such that f is equal to a linear or quadratic function on each of the intervals (-00, f1), [2]. tgl.-....... ..... [tx-] 11. tk], [tx. zo). In what follows merely consider nonnegative convex PLQ functions, which also imply the continuity of the loss function. It is well known that PLQ functions
universal approximators of continuous functions 38) but the canonical definition of PI.Q is not convenient for optimization algorithms. Instead, in this section prove that convex PL.Q functions actually be decomposed into series of simple and basic units.
Before proceeding, we introduce generalized quadratic version ReLU, the so-called rectified Huber unit, ReHU:
12 ReLU
ReHU-1
ReHU-2
10
ReHU-3
0, 22/2,
2 2 0
ReQU/ReHU-s
(2)
0 0<2ST, E>T
4-
where E: [0, oo] is a hyperparameter of ReHU switching from quadratic linear form, see ex- amples in Figure [ ReHU also incorporates the rectified quadratic unit (ReQU) as special case, that is, ReHU (2) ReQU(z)/2
(2+) / when 7 = xc, where 24 ma 0}. Then we give the characterization of a composite ReLU- ReHU function.
Figure 1: The ReLU function and the ReHU function with r == 1,2,3, 2,2,3,00. co.
Definition (Composite ReLU-ReHU function), A ffnntion L(z) is composite ReHU, if there exist u, V L and T,s,t E H such that
L(2) Σ ReLU(ux v1)
Σ
ReHU,, (sh2 + tt),
where ReLU(z)= = and ReHU,,, T. (=) is defined in (2).
Clearly, the composite ReLU-ReHU function class is closed under affine transformations, as is indicated by Proposition[]
Proposition (Closure under affine transformation). If L(z) a composite ReLU-ReHU function as in (3, then for any e > 0, E R, and E I, cL(pz + q) is also composite ReLU-ReHU, that is,
H
cL(pz + + = Σ (efu + +f)
ReHlU o3 (5/, + tty),
1=1
h=1
where ui cpu, 1/2 cuq + cua, Th VoTh. s Vepsn. and t velsng + tn).
3

---

More importantly, we show that the composite ReLU-ReHU function class is equivalent to the class of convex PLQ functions.
Theorem 1. loss function
> Rz0 is convex PLQ if and only if is composite ReLU-ReHU
Based on the fact that the loss functions of various SVMs (the hinge loss, the squared hinge loss, and the smoothed hinge loss), quantile regression (the check loss), and least absolute deviation (the absolute loss) all convex PLQs, Theorem I suggests that Definition incorporates numerous widely-used losses in statistics and machine learning applications. Table provides a convenient translation from some commonly loss functions to their ReLU-ReHU representation.
Table 2: Some widely-used composite ReLU-ReHU losses as in 3 Here SVM weighted SVMs
based on the hinge loss 77. sSVM is smoothed SVMs based the smoothed hinge loss [33. SVM weighted squared SVMs based on the squared hinge loss 771. LAD the least absolute deviation regression, SVR is support vector regression with the g-insensitive loss [4 and QR is quantile regression with the check loss [26].
PROBLEM
Loss (Le(ze))
COMPOSITE RELU-REHU PARAMETERS
SVM
c(1 1:21)+
"ii = Cille Uli
sSVM
cRRHUU1( 1-(1i2 1))
$1 -iiii U1i VG. r VG
SVM
ci((1 +)²
81i 22gii 41i /2,,
LAD
cily zil
W1: = c. th -43li, u2 = -G,, U2i = C33li
SVR
c((ly: - 24| - E)
W1 c, U1 = -(y: + E), -c. Un = -
QR
Cipk(yi 21) A.3
01: CsM, U1: NCIlh, U2: a1 n), 02: = -c,(1 (1 -n)
Table 3: Some widely-used linear constraints of the form A b 0, where U-bidiag(1-11 1]) is an upper bidiagonal matrix with main diagonal being 1 and upper diagonal being 1.
CONSTRAINT
FORM
REHLINE PARAMETERS
Nonnegative
320 0
I, 0
Box
1 Su
1:-1]. b [-1:
Monotonicity
3: ... Ba
A B B = U (-biddia((-1111)) 0a-1
Faimess A.I
|2'X3/ml Sp
= Z' /I/n, b= pp:pl
Taken together, our algorithm is designed to address the empirical regularized ReLU-ReHU mini- mization problem, named ReHLine optimization, of the following form:
min BERd
ReLU(unx]B vn)) +
ReHU (Snx + bNs) +
2 11812
1=1 1=1
1=1
s.t. A 222220 0,
(Ths) Hxn the ReLU-ReHU loss parameters, as illustrated in Table 2 and (A, b) are the constraint parameters, as illustrated in Table This formulation has a wide range of applications, including statistics, machine learning, computational biology, and social studies. Some popular examples include SVMs with fairness constraints (FairSVM), elastic net regularized quantile regression (ElasticQR), ridge Huber minimization (RidgeHuber), and smoothed SVMs (sSVM). See Appendix A for details.
where U = (wx), V (un) L and S (shi), T
Despite appearing impose limitations the l2 regularization in 4 our ReHLine formulation can actually be extended to more general regularization forms, including the elastic net penalty 49
Proposition 2. If L is a composite ReLU-ReHU loss function with parameters (un ww)l and (Sha, the T then its elastic net regularized minimization with \ 2 and X2 > 0,
A2 1|31
min BER
Σ
Lxxx +
s.t. A,3 A/+b20 0,

---

can be rewritten as the form of ReHLine optimization with
U
OLxd
1
V
OLxd
)
U-
Oī
11
V
Aa
O
X
X I
q+d
OT
11
Aa
81
0
n+d
S ( V1
X2
0
T
X2
T
0
T
(vir 03
where +- indicates assigning a value to a parameter.
3 Optimization algorithm
In this section, we present the ReHLine algorithm tackle the optimization in 4 The proposed algorithm is inspired by LIBLINEAR [12. which has gained great success in solving standard S SVMs. The key idea of ReHLine leverage the linear property of Karush Kuhi-KKhnnn--aa Tucker (KKT) conditions, thus simplifying the iterative complexity of CD and achieving rapid convergence.
3.1 Primal and dual formulations
To proceed, first demonstrate equivalent formulations for ReLU and ReHU, respectively, where the equivalence means that the value on the left is equal to the minimum value of the right-hand side:
min
min 0.0 s.t. 0 0
0 To
ReLU(z)
ReHU(z)
(5)
0.
0 2.0
is worth noting that when t oc, the optimal solution for and hence the last constraint (ar 2 0) can be removed. Then according to the definition of the composite ReLU-ReHU loss function and the equivalent forms in (5). the ReHLine optimization problem (4) can be rewritten as:
T
H
H
min
TK
The + 7 =113113
811
1=1
h=1
s.t.
A,
0, Th x G + Ups
Oni + ohi 2 sn. B + this
2 0,
0, for all 1i,
(6)
where II (TT) R A (0ns) R and Σ (ani) RR slack variables.
Theorem 2 shows the dual form of the primal formulation (6) by using the Lagrangian duality theory, and the relation between primal and dual variables is also provided.
Theorem 2. The Lagrangian dual problem of (6) is:
A. argmin A,I) E.A.T
s.t.
E O, E2A20, T2T20,
(7)
C(g,A,T) 2 STTATTE (e (3) U (3)
vec vec(TT) (S)))) + I)ve((T)
A0 + + Tr(AVT) THTTT),
(3 + ver ( S (e)
(8)
where & E R A= (Atc) E Bixn and r (mh) ll are dual variables, Ú ((3) dxnt. and S E RdxnH are the mode-3 unfolding the tensors Ú (wuj) RLxaxd and S (ska)) Rlxxxd, respectively, Uli UgZij.. Shij = ShyJiy, I is the identity matrix, and all inequalities are elementwise.
Moreover, the optimal poins B of 66) can be recovered as:
K
n
1.
ERaF
Σg
x;
Σ
Auau + Σ
Thosks
A (()eee((A) -
(9)
h=1
5

---

3.2 ReHLine update rules
Indeed, 9 is KKT condition of 6 and the main idea of the proposed ReHLine solver is to update the dual variables (5. A. by CD on (7) with an analytic solution, and simultaneously update the primal variable 3 by the KKT condition (9. Most importantly, the computational complexity of the CD updates on 7 can be significantly reduced by using the information of B3.
Canonical CD updates. As first step, we consider the canonical CD update rule that directly optimizes the dual problem 7 with respect to single variable. For brevity, in this section we only illustrate the result for \ variables, and the full details are given in Appendix B
By excluding the terms unrelated to Ati we have A = argminosssi Lu where
(A)
u
(x]xx)12
E
(x
Σ Eous (a]x)
(1,8)((1)))
Wil x AGA
Therefore, by simple calculations we obtain
new = Pl..1l !xx1 EE Exac Evreyaa App
K
Ex Ncs un
[x
(10)
where (2) = max(a, min(b, x)) means projecting a real number x to the interval [a, 6].
Clearly, assuming the values x] as and ||x; 112 are cached, updating one A value requires O + +nd nL+ n/) computation, and updating the whole matrix requires (nL + +nd + +n nl)). Adding all variables together, the canonical CD update rule for one full cycle has a computational complexity of O((K- ndd + nL + nH)(K + nL + nH)).
ReHLine updates. The proposed ReHLine algorithm, on the other hand, significantly reduces the computational complexity of canonical CD by updating (3 according to the KKT condition (9) after each update of dual variable. To see this, let := (E, I) denote all the dual variables, and define
K
B(u)
ELaR
Σx x,
AuaN Σ ShiShi
Then it can be proved that (VxLs) = ( + un). Therefore, when μ is fixed at oold (guld A old rold) and let gold =3( (10) can be rewritten as
(Vx, (V,, (2old) u 1x
uix gold +
Xno Po. Aold
= P(o.1]
Apd +
Ut
u |x,
Accordingly, the primal variable 3 is updated as
nnew = gold
(Nyow old Jun:i
which can then be used for the next dual variable update. Simple calculations show that this scheme only costs (dd) of computation for one A variable.
The update rules for other dual variables are similar, with the complete algorithm shown in Algorithm I and the full derivation details given in Appendix B Overall, the ReHLine update rule has computational complexity of O(K' + nL + nd) for one full cycle. Given that (K, L, H) are all small numbers, the complexity is linear with respect to both the sample size n and the dimension d.
We emphasize that the linear relationship ( between the primal and dual variables, as well as its consequence of reducing the computational cost, is highly non-trivial. This technique has been proposed for specialized problems such as SVMs via the LIBLINEAR solver (21) and here we show that its success be greatly generalized to any convex PLQ loss functions.

---

Algorithm 1: The ReHLine algorithm that solves (4).
Input :X E Ruud U, V E L S.T. T RHxn A E RXxd AERRxd bERK
Initialize 0, E A 0,
r 0: 0:
2 Compute 3 A Li=1
x, Σi=1 im + Eb Mhsha);
3 Compute and store = (r;) E R" r; ||x 12 and p = (px) E RK = =laa/
while maximum iteration do for k K [K] do
E" + max ( Sk. Pk 1 (a + bx)); & + Ea + B B + t
7
for (1, i) E f = [I. x [n] do
6" f- max
Au. min (1
Ais u (Mis t + uix! +,xx B)));
Au - Au ; B s
Ei
for (h, i) E H = [H] X [n] do
10
f max ( This min (The This (s 1) 1tni + shx -hh)));
This t Mhi + E": B - B E² chiniii
Output: :f. A, g
3.3 Global convergence rate
Finally, we show that the ReHLing algorithm has a fast global linear convergence rate.
Theorem 3. Let n = (g(1), A r(s) be sequence of iterates generated by Algorithm ! and μ' be minimizer of the dual objective (8). Then μ( is feasible, and the dual objective value converges at least linearly to that of 4², that is, there exist 0 < 1) < 1 and qo such that for all q 2 qo-
((1 L(v²) ((((()))) L(x²)).
Related work
Table [
summarizes the existing methods in solving the primal or dual problem of II. Note that Column #ITERATIONS is the best-fit theoretical results of convergence analysis for each algorithm based on our literature search, which may be improved with further developments of convergence analysis. These algorithms can be broadly divided into the following categories.
Generic QP solvers. For instance, projected sub-gradient descent is readily applicable to the primal form (1) by computing the sub-gradient of each loss and the projection under affine constraints, although the projection operator itself is already difficult. As for the box-QP dual problem 77. CD Box-QP can also be efficiently solved using IPM, which leverages second-order cone programming techniques [3 [4] Moreover, ADMM was employed to iteratively minimize the augmented Lagrangian of II
converges linearly to the optimal solution [27 31 35 42
ERM solvers. There are multiple existing works developed tackle ERM problems. such as SAG 37 SVRG (25]. SAGA [8). proximal SDCA 44) and SDCA (39]. However, SAG, SVRG, SAGA and proximal SDCA can only handle smooth loss functions, with an optional non-smooth term that has an easy proximal operator. SDCA applies to general loss functions, but it requires the convex conjugate of loss functions, which is not necessarily simple to compute, also see discussion in Section 3.1 in 44 Moreover, it only guarantees sublinear convergence rate for non-smooth loss functions. In contrast, ReHLine supports all convex PLQ loss functions with optional general linear constraints, which could be non-smooth by construction, and they all enjoy linear computational complexity and provable linear convergence.
For non-smooth loss functions, another existing method to solve ERM is the smoothing technique 40 3 30), which approximates the non-smooth terms by smooth functions, and then uses smooth- based algorithms to solve the smooth problem. However, the choice of the smoothing function and smoothing parameter typically requires additional knowledge, and the linear convergence rate not always guaranteed.
7

---

Despite the successful use of the existing algorithms in solving (), these methods either suffer from slow convergence (P-GD), or high computational costs (CD, IPM, and ADMM). This is partly because these general methods do not fully exploit the inherent structure of the problem. In contrast, ReHLiie leverages the KKT condition ( to substantially improve numerical performance, thus achieving a linear computational complexity per iteration with respect to the sample size, as well as linear convergence rate.
In addition to generic algorithms, specialized algorithms have been developed for specific types of T For example, the dual coordinate descent of LIBLINEAR in solving SVMs [21] the semismooth Newton coordinate descent algorithm of HQREG in solving the ridge regularized Huber minimization
[4]. and SAGA. SAG, proximal SDCA, and SVRG of LIGHTNING solving smooth ERMs 44.
such as smoothed SVMs. Our algorithm is also compared to these specialized methods in Section
Software. Below are some commonly-used software for solving ① based the algorithms available.
CVX/CVXPY [7 [0] modeling language for convex optimization problems and supports wide range of solvers such as ECOS, SCS, OSQP, and MOSEK.
MOSEK
a commercial convex optimization solver based on IPM.
ECOS
an open-source solver of second-order cone programming based on IPM.
SCS an open-source convex quadratic cone solver based on ADMM.
DCCP a disciplined convex-concave programming solver built on top of CVXPY, which has been used solve classification problems with fairness constraints 47
LIBLINEAR 21: library for large linear classification.
HQREG 46 R package solving Lasso or elastic net penalized Huber loss regression and quantile regression.
LIGHTNING 4 a Python library solving various ERMs based on primal SDCA, SGD, AdaGrad,
SAG, SAGA, SVRG.
5 Experiments and benchmarks
In this section we demonstrate the performance of ReHLine compared with the existing state-of- the-art (SOTA) solvers various machine learning tasks. The source code of ReHLiine along with Python/R interface, is readily accessible on our GitHub repository (https: //gi thub com/ and our project page hhttps://rehlii: involve five prominent machine learning tasks that draw from diverse sets of data. Specifically, we focus four classification datasets and five regression datasets sourced from OpenML. (httpa: www.woopeml..r... org/). //), with various scales and dimensions. The experiment settings are summarized well-organized toolset and framework, the BENCHOPP library 29 to implement optimization benchmarks for all the SOTA solvers.
httptgttt....... github io/). Our experiments
in Table To achieve fair comparison,
Due to representation limits, only present the results using variant of the ReHLine algo- rithm that supports variable shrinkage (refer to Algorithm 2). Results utilizing ReHLine with or without shrinkage reported on our benchmark repository ReHLine- ben chmark). Note that there is no substantial difference between these variations, and they do not influence the conclusions drawn from experiments. Simultancously, in the benchmark repository, we provide a more detailed "time versus objective" plot for each task/dataset and each solver produced by BENCHOPT.
com sof tmin/
Table 4: The classification and regression datasets with different scales used in our experiments.
CLASSIFICATION DATASET (n X d)
REGRESSION DATASET (n x d)
Steel-plates-fault (SPF) 1941 34
Liver-disorders (LD) 345 6
Philippine 3832 309
Kin8nm 8191
Sviva-prior 14395 x 109
Topo-2-1 8885 x 267
Creditcard 284807 X 31
Ho 22784 x
Buzzin Set-zinn-Wiiteer Twitter (BTT): 583250 x 78
For the FairSVM task, ReHLine compared with the original solver 48 based on DCCP and other generic solvers in CVXPY. In this case, both feasibility and optimality examined, where the

---

optimality is measured by the objective function, and the feasibility is measured by the violation of constraints: max (n 1 x, 2 ElasticQR by considering the model defined in (A.2) with A X² 1, RidgeHuber of (A.4) with A1 = 0, A2 = 1, and SVM of (A.T) and sSVM of (A.5) with C = Table presents the running times in seconds of all solvers that converge based on a relative/absolute tolerance of 10-5.
p. 10 8). Moreover, we examine the performance of
Table 5: The averaged running times (+: standard deviation) of SOTA solvers machine learning tasks. "X" indicates cases where the solver produced an invalid solution or exceeded the allotted time limit. Speed-up refers to the speed-up in the averaged running time (on the largest dataset) achieved by ReHLiine where "oo" indicates that the solver fails to solve the problem.
TASK
DATASET
ECOS
MOSEK
SCS
DCCP
REHLINE
FairSVM
SPF (xle-4)
X
X
4.25 (+:0.5)
Philippine (ele=2)
1550 (1006)
87.4 (1:02)
130 (1:2)
1137 (19.2)
1.03 (#:0.2)
Sylva-prior (x le-2)
0.47 (4:0.1)
Creditcard (xle-1)
175 (10.2)
64.2 (100.1
161 (1:405)
0.64 (1:0.2)
Fail/Succeed
2/2
2/2 100x
2/2
3/1
0/4
Speed-up (on Crediccad)
273x
252x
00
TASK
DATASET
ECOS
MOSEK
SCS
REHLINE
ElasticQR
LD (xle-4)
106 (#:7)
34.9 (1:25.0)
2.60 (+0.30)
Kin8nm (xle-3)
92.0 (11.0)
63.1 (155.5)
4.12 (1.0.55)
House-&L (xle-3)
887 (+161)
277 (+34)
X
7.21 (#11.99)
Topo-2-1 (x1e:2)
4752 (4:2015)
X
3.04 (+:0.49)
BT (x1e-0)
7079 (+:2517)
2.49 (=0.56)
Fail/Succeed
3/2
2/3
3/2
0/5
Speed-up (on BT)
2843x
TASK
DATASET
ECOS
MOSEK
SCS
HQREG
REHLINE
RidgeHuber
1.1) (x1e-4)
X
X
X
4.90 (++0.00)
1.40 (++0220)
Kin8nm (x1c-3)
X
X
X
1.58 (#0.211)
2.04 (4:0.0))
House-8L (xle-3)
925 (1:2)
2.42 (1.0.44)
0.80 (1:0.21)
Topo-2-1 (x1e-2)
2620
(=104)
267 (#:1)
213 (#:2)
3.53 (200667
1.78 (+:.33)
BT (xle-1)
X
2384
(1-433)
X
12.5 (11.8)
5.28 (11.31)
Fail/Succeed
4/1
2/3
4/1
0/5
0/5
Speed-up up (on BT)
30
452x
00
2.37x
TASK
DATASET
ECOS
LIBLINEAR
REHLINE
SVM
SPF (x1c-4)
X
372 (#1)
237 (1:27)
12.7 (#00.1)
3.90 (+0.10)
Philippine (x le-2)
1653 (1411)
86.5 (1.0.2
153 (=146)
1.80 (1:0.02)
0.82 (1:0.02)
Sylva-prior (x
X
731 (+-2)
843 (+:1006)
16.0 (+006)
4.08 (+-0884)
Creditcard (x 1e-2)
2111 (4:84)
1731 (1:4510)
23.1 (4:2.5)
5.08 (11.455
Fail/Succeed
2/2
1/3
0/4
0/4
0/4
Speed-up Creditcard)
415x
340x
4.5x
TASK
DATASET
SAGA
SAG
SDCA
SVRG
REHLINE
sSVM
SPF (xle-4)
39.9 (114.6)
28.3 (15.00
15.0 (11..4)
41.4 (113.9)
4.80 (1:1.20
Philippine (x 1e-2)
24,3 (+:27.8)
5.53 (+-9.)
1.47 (++0.19)
15.8 (++6.8
0.89 (++0.10)
Sylva-prior (xle-2)
3.37 (+9.81)
3.00 (4-056)
(+0.23)
3.40 (+00844
0.86 (#0.14)
Creditcard (x
10.4 (+1.44
(1:2.0
14.0 (+1.9)
11.2 (+1.4)
6.36 (t:1922)
Fail/Succeed
0/4
0/4
0/4
0/4
Speed-up (on Creditcard)
1.6x
2.3x
2.2x
1.7x
The results in Table indicate that the proposed ReHLine algorithm/software achieves visible improvements over existing solvers. The major empirical findings are as follows. First, the amount of improvement achieved by ReHLine in terms of running time is substantial over generic solvers, including ECOS, SCS, and MOSEK based on CVX, with the largest improvement 100x-10000 speed-up on the largest-scale dataset. The improvement continues to expand even further as the scale of the problem grows. Second, ReHLine also yields considerable improvements for specialized algorithms, including LIBLINEAR in SVM, HQREG in Huber minimization, and LIGHTNING in sSVM. Third, given that the generic solvers fail or time out in multiple tasks and datasets, ReHLine has shown remarkable flexibility when it comes to solving various domain-specine problems.

---

6 Conclusion
In this paper, we present new algorithm ReHLine designed efficiently solve general regularized ERM with the convex PLQ loss function and optional linear constraints. Through the transformation of PLQ loss into a composite ReL.-RRHHU function, we have developed a novel CD algorithm that updates primal and dual variables simultaneously. This approach has proven to be highly effective in reducing the computational complexity while at the same time achieving a fast linear convergence rate. Based on our experiments, ReHLine demonstrates remarkable flexibility and computational efficiency, outperforming existing generic and specialized methods range of problems.
Limitation. Although our approach has exhibited substantial versatility in addressing diverse range of problems, it is important to acknowledge that it may not be universally applicable to all convex optimization problems. The main limit of our approach stems from the restriction imposed by the formulation (T. For instance, (T) requires a convex PLQ loss function and a strict Ie-reuularization We recognize this limitation and intend to address it in future research pursuits. Moreover, interested in exploiting additional loss function properties, such as symmetry.
7
Acknowledgements
Ben Dai's work supported in part by HK GRF-24302422 and GRF-14304823. Yixuan Qiu's work was supported in part by National Natural Science Foundation China (12101389), Shanghai Pujiang Program (21PJC056). MOE Project of Key Research Institute of Humanities and Social Sciences (22JJDI 10001), and Shanghai Research Center for Data Science and Decision Technology. The authors confirm contribution to the paper as follows. The main algorithm: Ben Dai and Yixuan Qiu; theory and proofs: Ben Dai; C++ implementation: Yixuan Qiu; Python interface: Ben Dai and Yixuan Qiu; R interface: Yixuan Qiu; manuscript preparation: Ben Dai and Yixuan Qiu, equally.
References
[1] Alekseeva, E., Brotcorne, L., Lepaul, S., and Montmeat, A. (2018). bilevel approach to optimize electricity prices. Yugoslav Journal of Operations Research, 29())9-30.
[2] ApS, M. (2019). The MOSEK Optimization Toolbox MATLAB Manual. Version 9.0.
[3] Beck, A. and Teboulle, M. (2012). Smoothing and first order methods: A unified framework. SIAM Journal on Optimization, 22(2):557-580
[4] Blondel, M. and Pedregosa, F. (2016). Lightning: large-scale linear classification, regression and ranking in Python.
[5] Boyd, S., Parikh, N., Chu, E.,, Peleato, B., Eckstein, J., al. (2011). Distributed optimization and statistical learning via the alternating direction method of multipliers. Foundations and Trends in Machine Learning, 3(1):1-122.
[6] Bubeck, S. al. (2015). Convex optimization: algorithms and complexity. Foundations and Trends@ in Machine Learning, 8(3-4):231-35
[7] Cortes, C. and Vapnik, V. (1995). Support-vector networks. Machine Learning, 20(3):273-297.
[8] Defazio, A., Bach, F., and Lacoste-Julien S. (2014). Saga: A fast incremental gradient method with support for non-strongly convex composite objectives. Advances in Neural Information Processing Systems, 27.
[9] Deng, W., Lai, M.-.,, Peng. Z., and Yin, W. (2017). Parallel multi-block ADMM with o(1/k) convergence. Journal of Scientific Computing, 71:712-736.
[10] Diamond, S. and Boyd, S. (2016). CVXPY: A Python-embedded modeling language for convex optimization. Journal of Machine Learning Research, 17(1):2909 2913.
[11] Domahidi, A., Chu, E.., and Boyd, (2013). Beos: socp solver for embedded systems. In 2013 European Control Conference (ECC), pages 3071-3076. IEEE.
10

---

[12] Fan, R.-E., Chang, K.-W., Hsieh, C.-J., Wang. X. and Lin, C.-J. (2008). LIBLINEAR: A library for large linear classification. Journal of Machine Learning Research, 9:1871-1874.
[13] Ferris, M. C. and Munson, T. S. (2002). Interior-point methods for massive support vector machines. SIAM Journal on Optimization, 13(3):783-804
[14] Fine, S. and Scheinberg, K. (2001). Efficient SVM training using low-rank kernel representa- tions. Journal of Machine Learning Research, 2(Dec):243-264.
[15] Fukushima, K. (1969). Visual feature extraction by multilayered network of analog threshold elements. IEEE Transactions on Systems Science and Cybernetics, 5(4):322-333.
[16] Gérard, M., Clautiaux, F,, and Sadykov, R. (2016). Column generation based approaches for a tour scheduling problem with a multi-skill heterogeneous workforce. European Journal of Operational Research, 20193000
[17] Grant, M. and Boyd, (2014). CVX: Matlab software for disciplined convex programming, version 2.1.
[18] Griva, I,, Shanno, D. F., and Vanderbei, R. J. (2004). Convergence analysis of primal-dul interior-poimt method for nonlinear programming. Optimization Online.
[19] Gu, Y., Fan, J., Kong, L., Ma. S., and Zou, H. (2018). ADMM for high-dimensional sparse penalized quantile regression. Technometrics, 60933311--33.
[20] He, B. and Yuan, X. (2012). On the o(1/n) convergence rate of the Douglas-Rachford alternating direction method. SIAM Journal Numerical Analysis, 50(2):700-709.
[21] Hsieh, C.-J,, Chang, K.-W., Lin, C.-J,, Keerthi, S. S., and Sundararajan, S. (2008). A dual coordinate descent method for large-scale linear SVM. International Conference on Machine Learning, pages 408-415
[22] Huber, P. J. (1964). Robust estimation of a location parameter. The Annals of Mathematical Statistics, 35(1):73-101.
[23] Huber, P. J. (1965). robust version of the probability ratio test. The Annals of Mathematical Statistics, pages 1753-1758.
[24] Jagannathan, R. and Ma, T. (2003). Risk reduction in large portfolios: Why imposing the wrong constraints helps. Journal of Finance, 18(4):1651-1663.
[25] Johnson, R. and Zhang, T. (2013). Accelerating stochastic gradient descent using predictive variance reduction. Advances in Neural Information Processing Systems, 26.
[26] Koenker, R. and Bassett Jr, G. (1978). Regression quantiles. Econometrica, 46(1):33-50,
[27] Luo, Z.-Q. and Tseng, P. (1992). On the convergence of the coordinate descent method for convex differentiable minimization. Journal of Optimization Theory and Applications, 72(1):7-35.
[28] Michel Valladão, D., Veiga, Á., and Street, (2018). A linear stochastic programming model for optimal leveraged portfolio selection. Computational Economics, 51(4): 10211:222
[29] Moreau, T., Massias, M., Gramfort, A., Ablin, P., Bannier, P.-A., Charlier, B., Dagréou, M., Dupre la Tour, T., Durif, G., Dantas, C. F., al. (2022). Benchopt: Reproducible, efficient and collaborative optimization benchmarks. In Advances in Neural Information Processing Systems, pages 25404-25421.
[30] Nesterov, Y. (2005). Smooth minimization of non-smooth functions. Mathematical Program- ming, 103:127-152.
[31] Nesterov, Y. (2012). Efficiency of coordinate descent methods on huge-scale optimization problems. SIAM Journal Optimization, 22(2):341-362.
[32] O'donoghue, B., Chu, E.,, Parikh, N., and Boyd, S. (2016). Conic optimization via operator splitting and homogeneous delf-dual embedding. Journal of Optimization Theory and Applications, 1042
11

---

[33] Rennie, J. D. and Srebro, N. (2005). Loss functions for preference levels: Regression with discrete ordered labels. Proceedings of the IJCAI Multidisciplinary Workshop on Advances in Preference Handling, volume 1. Citeseer.
[34] Rockafellar, R. T. and Wets, R. J.-B. (2009). Variational analysis. Springer Science & Business Media.
[35] Saha, A. and Tewari, A. (2013). On the nonasymptotic convergence of eyclic coordinate descent methods. SIAM Journal Optimization, 23(1):576-60)
[36] Salas, V. and Saurina, (2002). Credit risk in two institutional regimes: Spanish commercial and savings banks. Journal of Financial Services Research, 22(3):203-224,
[37] Schmidt, M., Le Roux, N.. and Bach, F. (2017). Minimizing finite sums with the stochastic average gradient. Mathematical Programming. 162:83-112.
[38] Schumaker, L. (2007). Spline Functions: Basic Theory. Cambridge university press.
[39] Shalev-Shwartz, S. and Zhang, T. (2013). Stochastic dual coordinate ascent methods for regularized loss minimization. Journal of Machine Learning Research, 14(1).
[40] Shalev-Shwartz, S. and Zhang. TT (2014). Accelerated proximal stochastic dual coordinate ascent for regularized loss minimization. In International Conference on Machine Learning, pages 64-72. PMLR.
[41] Shen, X., Diamond, S., Gu, Y., and Boyd, S. (2016). Disciplined convex-concave programming. In 2016 IEEE 55th Conference on Decision and Control (CDC), pages 1009-1014, IEFE.
[42] Tseng, P. and Yun, S. (2010). A coordinate gradient descent method for linearly constrained smooth optimization and support vector machines training. Computational Optimization and Applications, 47(2):179-206.
[43] Vapnik, V. (1991). Principles of risk minimization for learning theory. In Advances in Neural Information Processing Systems, pages 831-838.
[44] Vapnik, V. (1999). The Nature of Statistical Learning Theory. Springer science & business media.
[45] Wright, S. (2015). Coordinate descent algorithms. Mathematical Programming, 151(1):3-34.
[46] Yi, C. and Huang, (2017). Semismooth newton coordinate descent algorithm for elastic-net penalized Huber loss regression and quantile regression. Journal of Computational and Graphical Statistics, 26(3):547-557.
[47] Zafar, M. B., Valera, I., Gomez- M., and Gummadi, K. P. (2019). Fairness constraints: A flexible approach for fair classification. Journal of Machine Learning Research, 20(1):273) 2778.
[48] Zafar, M. B., Valera, I., Rogriguez, M. G., and Gummadi, K. P. (2017). Fairness constraints: Mechanisms for fair classification. In Artificial Intelligence and Statistics,, pages 962-970.
[49] Zou, H. and Hastie, T. (2005). Regularization and variable selection via the elastic net. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 6712):001-320.
12

---

A
Example Models
SVMs with fairness constraints (FairSVM).
The formulation for Va 48 is:
C
min BER
(1 x,)
+ 115112
i=1
1 Σ ==1
S.L. 2 x; S
2 x; 2 -p
(A.1)
1=1
where X; E Rd a feature vector, and Y/ E {-1,1} a binary label, 2, E IR a collection of centered sensitive features (E 2; 0), such as gender and/or race. The constraints in (A.1)
constrain the correlation between the sensitive features 2; and the decision function 3 x;, and the constants E T trade off predictive accuracy and fairness. Note that the FairSVM A.1 can be rewritten as ReHLine optimization of with
ZTX/n U - CCyT/ V C1 A +
b (
-ZXX/n)
where 1n (1, (1,.... T)T the length-n one vector, gaxd the feature matrix, and y (y1, ... (yx) is the response vector.
Elastic net regularized quantile regression (ElasticQR). The formulation for the elastic net penalized quantile regression 19 is:
min BERU+1
Σ
Pa(si - xI - Ba-1) + A 1|1 +
d 11312 2
(A.2)
where
Pi (u) = u(x 1(u < 0))
(A.3)
the check loss, 10, is a prespecified sample quantile, x; E Rd a feature vector, y/ E RR is response, and Åι, Xo 2 0 are weights of lasso and ridge penalties. Then, the ElasticQR (A.2) can be rewritten as a ReHLine optimization of with
K 1 RAg 12
d+1
n
YT
O d4+
1-x 1=x1T nAg n
O d+11 1 d+1 1
Y
X 1n X
) Id+1
U
OT
A1
V
nAg
OT
O_
A2 d+1
OT
0 d+1
where I is the identity matrix.
Elastic net regularized Huber minimization. minimization 23 is:
The formulation for the elastic net penalized Huber
Xa 11111
min 8
(M - xT + A111B111 +
(A.4)
71
1=1
where Hx(-) is the Huber loss with given parameter
//2
<< |2| S >
Hx(=)
x([2| x/2),
|2|
this case, A.4
can be rewritten as a ReHLine optimization with
11 0
O
Ri
ndg
1I
0
S
md²
n
T
nA
It
11
03 1
ndg Y
0
nd
Rig
nAa
11
AL 1
X
U
V+ 0,
X +
L
0g
A1 17
Au
13

---

Smoothed SVM (sSVM). The formulation for the smoothed SVMs 33 is:
min 1 ( x;) 11112
(A.5)
where E Rd is feature vector, and Vi E {-1,1] is binary label, and Vz(+) is the modified Huber loss the smoothed hinge loss:
0, (1 (1-2)//2,
= 1, < 1 2 <<0.
V(=)
(1/2- -12-2),
2),
In this case,
A.5 can be rewritten as ReHLine optimization with
S -V-NN/ T VC/n1/ It VC/n1.
B
Update Rules for (
B.1
Canonical CD updates
Update & given others.
For k 1,.., K,
enew argmin Sk
- a
a + Σ EwaffarEE Σ NuumsEaaf - E h
InnsonEnaaf + baEK
2
20
k'tk
15
a (Ex EwaN
Ei, Anaii;
En. Whi8hix;) + bk
max
(B.6)
|axl2
Update du given others. For l = 1,... and = =;,,.... -,....m.
A/s
new argmin 12
x +
Σ
Σ
E Uu a +
Σ
k=1
K Ex ExaR
E(IRDA Avp
Ene MC8SH08))) uii
4 xx:1
(B.7)
Update "hi given others.
For h = 1,... H and 1,...
Thi
new
argmin The 2aa
(5) + 1)7h- +
Σ
2
(1h,))
k
a.snaaxim + Σ
Sna x themii
A=1
P.8
SNix]
Ek Exaa E (1
NWPSNTX Eve Aeruvyx)) + the
Po.ra
x +1
(B.8)
The canonical CD updates in B.6 B.8) for one full iteration require computation complexity of O(K nd + nL + (K nL nH)).
14

---

B.2 ReHLine updates
Consider the Lagrangian function C in 7 Clearly,
VeL a
Ewalaa ak Σ
Auuuaxx Σ
Mishia + by.
Σ Ew Ewa ak
Σ
Σ
/hss + by h
l,
(44) + by-
Similarly, derive that
Vx, =
- (uxx] ((p) mu),
Vn Mhi (shixTB(pp) + tta)).
Then the CD updates [.6) can be simplified as:
VgL.ie ((Cold) max 0.
aT + bk
enew
P(o,+0c)
cold
Sk
|
|a 12
guew
Gold
Sk
enld
Jak.
(B.9)
The CD updates in (B.7 simplified as:
VA, C(2old) max 0. (0.in min
uxx! f U1i
A new = P/1111
1/
1
u X;
Aold , The CD updates in (B.8) are simplified as:
guew 3old
(A
(B.
Thi
Plo,tnil
Thi
old
V The ((yold)
xx +1
max 0, min
This thi old
8hix]
B
x 112
+ Lns
+
Thi
old
@new -glld
(The This 88
(B.11)
The computational complexity of the new ReHLine iterations is O(K nL nH) dd), shown in (B.9) through (B.1I). Notably, this complexity is linear with respect to both the sample size n and the dimensionality d, making ReHLine highly efficient algorithm.
C Screening rules for ReHLine
We extend the shrinking strategy in Algorithm 3 of 12 to the ReHLine solver. To begin with, the following theorem advocates utilizing "sparse" updates optimize the performance of ReHLine.
Theorem Let A"" ,33 be the convergent point of (g(1), A(4), ( /g()) based on (B.9.
(B.1I, then
(A) For & updates, let 9% aT by.
if > then 3y0 such that Vq 40. {(() 5h 0.
(B) For A updates, gi unx] 9 + vii
(9)
if gh 0, then 3go such that Vq 90. A Xi; 0;
g% 0, then Equ such that Vq 90. A X) 1.
(C) For I updates, denote Ph. =s.x. 1 + thi-
if 7hi = O and gh. < O, then Ego such that 4q 2 G0.
7he
(9) = 0;
(9) Thi-
1/s Thi and 9%. > This then 3go such that Vq qo 7im
Building upon the insights from Theorem 4 we have formulated Algorithm 2 that involves shrinking certain "boundary" variables during the implementation of ReHLine.
15

---

Algorithm 2: ReHLine solver minimizing
with shrinking
Input :X E Rnxd
U,
V
E _LxM; S, T
IR H xn A E 1xd
A
b I2 tol > 0
Initialize 0. E
A
0,
0:
2 Compute 3 Σ 1 Eąa Σi= x; EE/ Anuii + Σ Thisks
K
3 Compute and store = (r:) R", ri = ||x: and p = (px) R Pk
| ::
Initialize M, = oc, m, 00, s = (5. 7),, K= [K], C = [1] 2 [m]. H = H x
[n]:
while maximum iteration do Initialize M,,- = 00, m, 00, s (5. A,7); 7}: // CD updates for
for & E (EC do
gs: - a + ba;
(is_shrink, 9k) +- PG_ feessiaalllesseeeee Ik, Mc)
10
if' is shrink then
11
K K, [k]: continue:
12
Mg max
mg min Jk};
m -Pk gh): Ea Ex + B 3 a
// CD updates for A
14
for (1, E PERM (C) do
Uui = ( [3 vz):
(is_shrink,, (ts) PG ReLU(Ai Iii, My, m,);
17
if iss shrānk then
C C {1,i); continue;
19
Mr max Ju); my minfma, Iul:
max(-İ4, min(1 duo un Pri" gu)); du Xo + e 33 e" eu
// CD updates for r for i) E PERM(H) do
21
Ihi = Thi (sax) + 4n): (is_shrink, 9hi) t- PG ReHU Ihi, My, my); if fris_she shrink then
24
25
L H
ih continue:
26
Mr max(MM, Ji my f min(my, Init:
s" 1maa((-hi. min(Thi The (shara + 1) -h)))):
27
28
The + 7h + + B E" S"ShiKji
29
M, m, < and [M < S and m. S (E, then
30
if KC| [. + H7| = K + n(T + H) then
31
break:
32
else
33
Take all shrunken variables back,
34
i.e., K. [K], C = [L] x |n]. H = [1]] {[n]
35
M, f- 00, m, +- Xo, (5, 7}
if M, < 0 then M,
else M, M,, {, X, 7): if m, 2 0 then m. + co else m. + m. 1.13.7
37
O A, g
16

---

Algorithm 3: Projected gradients in ReHLine.
Function PG "d_feasiblee$ (Ck. Ok. .):
Ok 9k: fi__shrrnk <- False;
if ε 0 and 9k > .M then
is j___hrink < True
if & 0 and ge 2 0 then
Dh = Output: iishiin 9k
7 Function PG_ PG_ReLU
M, m):
It + gwv is Les__shrinn False;
(Au and gu M) or (Ax 1 and gi m) then
10
is shrinkk <- True
11
(As and Oii 0) or (Au: and Iü S 0) then
12
St = 0
Output:is.. Output:is_shrink 9t
Function PG_ReHU (Thi, Gho. This M, m)):
14
Ihe < 9hi is s__rriinnk- 4- False;
15
(Th = and Ohi M) or (The = Thi and Thi < m) then
16
is_shrink < True
17
(Ths = and ghi 2 0) or (The = Thi and ghi I 0) then
Ihi thrr ghs
Output: :is_ssisssisssiee
D Technical proofs
Lemma 1. Suppose fz) is nonnegative convex PLQ with knots and extreme points (tx) Then, there exists by E f1. K}, such that f(2) 2 f(tra) for any E R.
Proof. We prove it by contradiction. Suppose that for there exists such that f(z)) < (ta) for all {1. (1..... K}, then the minimum can only be obtained outside [t2, tx], that is
Ezo < <t,f(20)) f(te), or 320 tk. f(20) < f(t)) for any (1,-..,444 K).
For example, suppose that 20 (tx, oc). f(20) < f(4) for any & E (1,..., K}. Note that the quadratic function in (zx, oo) can be written as f(z) = 1f (tx)(=-tx)+ff1 (x
(t)) fttx). Based on the fact that (=0) < f(xx) and the convexity of f(=), we have ('x(x) < 0. Then, since (t is the set of all knots and extreme points, we have f'(=) < for
> tx.
and +om (2) 1om (f" (tx)(z tk) + (tx)) < O. This yields that (tx) = O. alternatively, f(z) is linear function in (tK, 00) as f(z) = (t tk) f(tx). Hence, when - 00, f(=) + -00, which contradicts with the fact that f(2) is nonnegative function. The same argument goes for the case when 20 (-o0, t1). This completes the proof.
Proof of Proposition I
Proof. Note that cReLU(z) cmax(z,0) max(cz, 0), and
0.
50
(ReHU,(x)
(2cz)2/2.
0< 0< Vcz s Ver = ReHUU,et(Vss)
-cT/2))
Vez Ver
17

---

for c > 0. Then,
H cReHU RReHU,,,s (sh + + + tn.
cL(pz + Σ
cReLU(u((2 +) + 1)) +
h=1
H
ReL.U(cu(pz + +) + co1) +
ReHU Von, (Vesa(pz + + vel,)
1=1
ReLU(u12 + v/) +
ReHU, 73 (6/,2 + tzy)
h=1
This completes the proof.
Proof of Theorem 1
Proof. The sufficiency (-) directly follows from the definitions of ReLU and ReHU. Next, the necessity (=)) of the theorem is established using induction. Without loss of generality, for nonnegative convex PLQ function f(2), we assume that (A1) f' (tx)f" (fx+1) 2 0, otherwise we can add the extreme point between (1x. tR+1] into the set of knots; (A2) ming=1,... f((x =
otherwise we can add ReLU (minx=1,.. K [ts)) to the final decomposition I shows that the minimum of fz) can only be obtained within the knots). Next, we conduct the mathematical induction.
K 1. Note that f(2) can be rewritten
{ 97 (z)
(x
2 25h,
f(z)
gt (2) ff(t))
t1). 2> f1
where 9i (z) and qjj (z) are left and right parts of f(z) separated
According to Lemma ] and Assumption A2, we have f(t) = 0 the minimum of f(2). Combined with Assumption AT, have (t) S for & t1, and (t) 2 0 for t f1. Then, f(z) can be decomposed as:
vi (2) ReHU 50 (-v/
4)) ReI.U(_" 4)), = = l1 0i (z) ReHU ( (4) 4)) +
f(:)
- 41),,
2 >> t
Hence, f(2) = ReHU (f t1)) + ReLURRL (t)) t1) + ReHU (VF
(2)) + ReLU(f (+2)(2 11).. The necessity of the theorem for = 1 then follows.
K == Ko + 1. Assume that convex PLQ ʃ 2 0 with K(K S Ko) knots is ReI.U-RHHU composite, now we turn to prove the statement still holds for f(z) with Ko + 1.
CASE 1. Suppose the minimum f(t)) is obtained at (12,.................. .t.)) (not in the boundary). In
this case, we can decompose f(z) as f(z) = fi(z) + fx) where fi(=) and f+(z) are defined as:
f()). 2 2st, t>t.
[+(=) {
0.
2 zst,
Si(=) O,
f(2),
2> t
Both f((2) and f+( are nonnegative convex PI.Qs with K S Kg knots, thus can be represented as ReHUS and ReLUs.
CASE 2. Suppose the minimum f(t") is obtained when t" t, or t For example, when !* t,, we conduct the same decomposition f(z) f((=== fx(x) f1. where (i()) is nonnegative convex PLQ with one knot, which ReHU-ReLU composite. For fr(z), we further decompose it as:
f+(2) ReLUU 41) ReHUU-
4)) R(2),
where 7 (T (1), and R(=) (with at most Ko knots) is defined as:
{0. {f(() f4(t)(( -t1 (1)tt)(tz f1)(z (t1 + t2)/2), z > tz. R(=)
2 S tz
18

---

Now, we show that R(2) is nonnegative convex. Note that f(2) is convex, then for 2 t2,
f(=) z [ -4) f+(t)))= - t2)
ReTUU(( 41)
where the first equality follows from the fact that f(2) is quadratic function in [t,, 12], and thus (1() = (f( - f1), and the second equality follows from f(t2) = (2(t1) (t2 11 (12
Thus, R(z) is nonnegative convex PLQ with Ko knots {/2, 1o+1). which can be decom- posed ReHUs and ReLUs. The same argument goes for the case when t tko+1.
This completes the proof, and the desirable results follow from the induction.
Proof of Proposition 2
Proof. The lemma follows from the fact that:
d
||3|| = 1181 = Σ
(ReLU(e]3)) ReLU(-e]3)).
j=1
where ey (0,... ,0,1,0 the d'tuple with all components equal to except the jth.
Proof of Theorem
Proof.
Note that the augmented Lagrangian of the primal problem (6) is:
H The +
K
Cp
Tii +
6
111112 2
Σ Am
Eu + bx)
Iml
1=1
hm
i=1
Aml
Ate(TT unx vii)
(Oni + ohe spir. B (xx)
h=1
H
Hhithi
i=1 t=1
k=1
elIle - 1101 Tr(TET)) || (TA
Tr -Ar(A 1U diag (X,3) V) - Tr (r(O + - -Saaag XB) - T)))
-Tr(AITT) TTT,
(D.12)
where E = (Ex E R A = (AB) E B4xa, r= ( E R A = -6xx E R and
V (Whi) &R Hx are dual variables. The KKT conditions for the augmented Lagrangian
OCp
n
H Thisha) 0.
B
ERa +
Σ
x;
03
1=1
1=1
h=1
3Lr
dLr
dLr Oops
Av
81
0,
Oni
Thi 0,
The 7hi 4/ O,
19

---

for all (i, h). Then, by substituting the KKT conditions into the primal Lagrangian, the (negative) dual objective function is derived as:
K K
LD =
EuEN a
+
+
Shi xxx
1a P.
R W
K
Σ
Σ>
1,8
k=1 A
LA W
Th
Σ
Eubu Σ Aucri Σ
Mithi
1
STAATE+ + -vec(A)TU 0 (3) (ec(A) vec((T) S(a) + +I) vec(T)
2 2
- AETAU (3) (e) EA 3) vec I) + vec(A) S (e))e 1
Tr(TTT)
where (s) and (g)) are the mode-3 unfolding of the tensors Ü E Lxxxd and S Hxxx and is an identity matrix.
Proof of Theorem 3
Proof. The treatment for the proof based on result of coordinate descent on boxed constrained convex optimization (27). To proceed, denote f (ET, vec(A)T, TT as the collection of all dual variables, then the dual objective can be rewritten as:
CD(μ) (Ep) c. μ.
where g(z) = = 2 and E and c are defined as:
E
A
0 (3)
S S(3)
E R BRd++)))((K
0
0
IH xn
(bT, vec( T vee
R +o+
Now, we tend verify assumptions of Theorem 2.1 (27. Specifically, (a) the dual problem is convex box-constrained QP, thus the set of optimal solutions is nonempty; (b) g(z) = 2 5| is strictly convex twice continuously differentiable function: (c) V²g(z) = Ia+ is always positive definite; (d) E has no zero column based on its definition. The desirable result then follows from Theorem 2.1 of (27].
Proof of Theorem
Proof. Note that the primal problem (6 is strictly convex, thus has the unique minimizer B Therefore, the sequence { provided by the proposed Algorithm converges to the unique minimizer B*. Next, we check the gradients each coordinate.
For É-coordinates, with & 2 0.
old = V. (p(
old = a. ay cold ld
+
Σ
cold sa
Σ 1,
2 oold una
Σ
Thi Shaa,x,, + bk
old
h,
n
a. cold a
old The
smi) + by = af + bb.
X
A
ups +
h=1
Suppose S the gradient of & based on ()) then lim Ir = =TT + bb =: gi-
Therefore, if of 0, then there exists qo for all g To we have g e(9+1) 0. Sk
(9) 0. which implies that
20

---

In the same manner, for A c-coordinates, with 1 2 \ui 2 0, old = VAn C1. (11 old) - u1:x B Dli-
Then, 1om Ih = u 3r - -11 =: =/:
if 9ú < 0, then 3go such that Vq 2 go. A/1 (g+1) = Xi: = 0;
if g/. > 0, then 3yo such that Vq 2 90. A/6
(g+1) X/s = 1.
For T-coordinates, with Thi Mhi 2 0, thus
old V Tn: LD
old Thi
ShiX. gold Ihi-
Thi
Then, 1o 9hi Thi ux 3 UN Th. which implies that
(4)
(4) 0;
if The and gh. 0, then 3go such that Vg q0. Thi
The The and ghs < 0, then 3go such that Vg 2 go. k Thi-
(q)
21


Reviews:
[{"review_id": "LUnAR0VOVo", "reviewer": "Reviewer_xSEV", "summary": "This paper introduced a new function class called composite ReLU-ReHU which is shown to be equivalent to the class of convex PLQ functions. Based on the ReLU-ReHU decomposition of convex PLQ functions, the authors then formulated a new box-constrained quadratic programming optimization problem, named ReHLine optimization which can be solved efficiently by the ReHLine algorithm. The proposed solver has a provable linear convergence rate and a linear per-iteration computation complexity. Benchmarking on various tasks and datasets, the ReHLine showed significant improvement over both generic and specific solvers, especially on large-scale datasets.", "questions": "1. Are the experimental results using Algorithm 1 or Algorithm 2 in the Appendix? Could the authors provide the running time of both algorithms for comparison?\n2. Could the authors provide an ablation study to demonstrate the objective gaps of different solvers when increasing the number of iterations?\n\nMinor:\n- There is a typo in Proposition 2: $T \\leftarrow \\begin{pmatrix} \\sqrt{\\frac{2}{\\lambda_2}} \\boldsymbol{T} & \\boldsymbol{0}_d^T \\end{pmatrix}$.", "limitations": "The authors have discussed the limitations of their work in Section 6. ", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- This paper is well-written and easy to follow.\n- To the best of my knowledge, the composite ReLU-ReHU functions and the ReHLine algorithm are novel.\n- The effectiveness of the proposed method is adequately backed by both theoretical results and empirical results.\n- The ReHLine algorithm successfully tackles any convex PLQ loss functions, thus it shows great potential for many machine learning tasks.\n", "weaknesses": "- While authors claimed that the simplification of Canonical CD updates is highly non-trivial, it seems to be a straight application of [18] to the ReHLine optimization problem.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "## Weaknesses\n\n> While authors claimed that the simplification of Canonical CD updates is highly non-trivial, it seems to be a straight application of [18] to the ReHLine optimization problem.\n\n**Reply**: Thanks for the comments. When the ReHLine optimization problem (1) and (3) is given, the paper seems to be a direct application of LibLinear [18]. However, the main contribution of our paper is actually the induction/construction of a class of problems that allows for the generalization of the LibLinear [18] algorithm, achieving both linear convergence and linear computational complexity.\n\nSpecifically, we would like to point out that whether the canonical CD can be simplified largely depends on whether the special relationship between the primal and dual variables as in equation (9) exists, and meanwhile we request the algorithm to have linear convergence. [18] is an important reference as it discovers this relationship in SVMs and proves its linear convergence. However, to the best of our knowledge, prior to ReHLine it is unknown whether these two properties hold for other more general cases. In this article, we have greatly expanded the class of models that possess these properties. In particular, we show that they are available not only for the hinge loss as in SVMs, but also for all convex PLQ functions. We think that this generalization is novel and non-trivial.\n\n## Questions\n\n> Are the experimental results using Algorithm 1 or Algorithm 2 in the Appendix? Could the authors provide the running time of both algorithms for comparison?\n\n**Reply**: Thank you for the comments. We are using Algorithm 2 (ReHLine solver with shrinking) for the experiments. Following your suggestion, we additionally conducted experiments to compare the ReHLine with and without shrinkage in the following tables.\n\nThe experimental results suggest that our algorithm's superiority (compared to other methods) is not influenced by shrinkage. In fact, in most cases, ReHLine without shrinkage may even outperform the one with shrinkage.\n\n| FairSVM | ReHLine(shrink=False) | ReHLine(shrink=True) |\n|:-----------|----------:|----------:|\n|SPF         | 3.951E-04 | 4.130E-04 |\n|philippine  | 1.186E-01 | 1.602E-02 |\n|sylva_prior | 9.234E-03 | 1.010E-02 |\n|creditcard  | 3.318E-01 | 2.157E-01 |\n\n\n| ElasticQR | ReHLine(shrink=False) | ReHLine(shrink=True) |\n|:---------------|----------:|----------:|\n|liver-disorders | 8.111E-05 | 2.165E-04 |\n|kin8nm          | 6.109E-04 | 3.491E-03 |\n|house_8L        | 1.618E-03 | 9.481E-03 |\n|topo_2_1        | 8.242E-03 | 4.549E-02 |\n|BT              | 1.454E-01 | 2.918E+00 |\n\n\n| SVM | ReHLine(shrink=False) | ReHLine(shrink=True) |\n|:-----------|----------:|----------:|\n|SPF         | 8.351E-04 | 3.968E-04 |\n|philippine  | 1.193E-01 | 1.083E-02 |\n|sylva_prior | 4.841E-03 | 5.709E-03 |\n|creditcard  | 2.973E-02 | 5.842E-02 |\n\n\n| ElasticHuber | ReHLine(shrink=False) | ReHLine(shrink=True) |\n|:---------------|----------:|----------:|\n|liver-disorders | 1.121E-04 | 1.210E-04 |\n|kin8nm          | 8.201E-04 | 2.176E-03 |\n|house_8L        | 7.009E-04 | 1.004E-03 |\n|topo_2_1        | 1.025E-1  | 2.084E-02 |\n|BT              | 1.737E-01 | 4.263E-01 |\n\n\n> Could the authors provide an ablation study to demonstrate the objective gaps of different solvers when increasing the number of iterations?\n\n**Reply**: Following your suggestion, we have included the figures depicting the optimization progress, specifically the objective function value over time or step, for all the benchmarks. Kindly refer to the attached PDF document for a comprehensive overview.\n\n> Minor: There is a typo in Proposition 2: $\\mathbf{T} \\leftarrow\\left(\\begin{array}{cc}\\sqrt{\\frac{2}{\\lambda_2}} \\mathbf{T} & \\mathbf{0}^\\intercal_d\\end{array}\\right)$\n\n**Reply**: Thanks for pointing out the typo. We will fix it in the revision."}, {"review_id": "2cOLIshRgd", "reviewer": "Reviewer_kbfD", "summary": "The authors consider tackling the optimization fo finite sum of piecewise linear quadratic functions (PLQ) arising from e.g. robust empirical risk minimization by means of a reformulation of PLQ functions as sums of ReLU and smoothed ReLU functions combined with a stochastic dual ascent algorithm. The authors start by presenting how any PLQ function can be formulated as sums of of ReLU and smoothed ReLU functions with a table summarizing some common ones. Then the proposed algorithm is presented whose convergence is ensured by previous work. Experimental results advocate for the superiority of the approach compared to old all-purpose solvers in various classification tasks. ", "questions": "- Can the authors do comparisons with state of the art algorithms for empirical risk minimization for each of the proposed problems rather than using all-purpose algorithms? A list of potential candidates has been given in the weakness section. \n- In particular, can the authors discuss and compare empirically their algorithm against SDCA in the absence of constraints?\n- How can the proposed algorithms handle nonlinear models such as kernel methods?", "limitations": "- No support for sparsity inducing norms\n", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "- The generality of the approach is well defended by the flexibility of the parameterization with ReLU and smoothed ReLU functions.\n- The proposed algorithm captures efficiently the underlying structure of the algorithm. In particular, compared to primal-based approaches (which could easily be implemented using the proximal operators of the losses), the approach can handle constraints. \n- The proposed algorithm is hyper-parameter free. \n- The experiments could be a strength, although they are simply misleading in their current state.", "weaknesses": "- The paper lacks some important related work. In particular, the proposed algorithm appears to simply be a form of stochastic coordinate dual ascent. Such a note would help gain perspective on the contributions of the authors.\n- A \n- Numerous algorithms have been developed to tackle empirical risk minimization problems. For example, s-SVM or SVM^2 losses are smooth and could be tackled by fast incremental solvers such as SAGA, SVRG or SDCA or even their accelerated versions such as [1]. see also [2]. Non-smooth losses can also generally easily be smoothed, see e.g. [3], to be amenable theoretically and in practice to fast resolution by fast incremental algorithms. True, most of these algorithms require stepsizes, which make the proposed approach friendlier at first. However, SDCA for example does not require any and could a priori be formulated like the proposed algorithm (using that the proximal operator of the conjugate of the losses is available here). A thorough discussion and comparison to all of the aforementioned algorithms is necessary to held this paper to the claims the authors make. Additional \n- Non-negativity, box and monotonicity constraints can be handled by simple projections and primal algorithms could also be considered in those cases.\n- The current approach cannot handle non-strongly convex penalty such as sparsity inducing penalties. \n- The experimental results and claims are misleading: no one uses all-purpose solvers such as MOSEK. In particular no one uses interior point methods for empirical risk minimization. Claiming 1000 times improvement is overselling the method which in turns is detrimental to the contributions of the paper. Usual comparisons consist in claiming improvements over state of the art method, which appear to hold but at much smaller scale. The proposed algorithm may provide quantitative benefits, but a comprehensive comparison would help the paper.\n\nMinor details:\n- E appears not have been defined in Theorem 2. Similarly recalling what a mode-3 unfolding of a tensor is, would improve the readability of the paper. \n\n[1] A. Defazio. A Simple Practical Accelerated Method for Finite Sums. NIPS 2016\n[2] S Shalev-Shwartz, T Zhang. Accelerated Mini-Batch Stochastic Dual Coordinate Ascent. NIPS 2013\n[3] A Beck. Smoothing and First Order Methods: A Unified Framework. SIOPT 2012", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "## Weaknesses\n\n> (related work)\n\n**Reply**: Thanks for the comments. We emphasize that our contribution is not only applying CD (and its extensions such as SDCA) to a specific optimization problem, but also to understand the class of problems such that the CD variant can be implemented with linear computational complexity and provable linear convergence.\n\nIn our revised manuscript, we have added more related works and discussed their connections with ReHLine. See our global response.\n\n> (s-SVM example)\n\n**Reply**: Following your suggestion, we have now added a new benchmark for smoothed SVM (s-SVM) to compare the running times of various algorithms, including SAGA, SAG, SDCA, SVRG, and the proposed ReHLine. In this experiment, the specialized implementations of the algorithms for s-SVM are implemented using the Python library `lightning`. It is worth mentioning that computationally demanding parts in `lightning` are implemented via `Cython` (C-Extensions for Python), so we believe that the comparison is roughly fair.\n\n|s-SVM|SAGA|SAG|SDCA|SVRG|ReHLine|\n|-|-|-|-|-|-|\n|SPF|3.154E-03|2.462E-03|1.316E-03|3.095E-03|5.298E-04|\n|philippine|1.292E-01|7.322E-02|2.210E-02|1.325E-01|1.086E-02|\n|sylva_prior|5.084E-02|4.055E-02|2.371E-02|3.140E-02|1.008E-02|\n|creditcard|1.146E-01|1.625E-01|1.438E-01|1.164E-01|6.451E-02|\n|**speed-up**|1.8~11.9x|2.5~6.7x|2.0~2.5x|1.8~12.2x|--|\n\nAs indicated in the table, even in smooth problems like s-SVM, ReHLine has shown reasonable improvement (or at least demonstrated comparable performance) of these specialized implementations.\n\n> (other algorithms for ERM)\n\n**Reply**: Thanks for providing these references. As we have shown in the response to the first question, in the planned revision we have added discussions on such related works. Some key points are again summarized below:\n\n1. SAG, SVRG, and SAGA require smooth loss functions, with an optional non-smooth term that has an easy proximal operator.\n2. Smoothing methods as in [BT12] do not have a linear convergence guarantee.\n3. ReHLine can handle linear constraints, whereas SDCA does not consider such cases.\n4. SDCA requests the convex conjugate of the loss function, but computing the convex conjugate of a general convex PLQ function is not straightforward.\n5. SDCA is only guaranteed a sub-linear convergence rate for **non-smooth** losses.\n\n> (special constraints)\n\n**Reply**: Of course, these methods can make use of the special structures of the constraints, but in ReHLine all linear constraints can be handled in a unified way, without sacrificing the computational efficiency and convergence rate.\n\n> (sparsity inducing penalties)\n\n**Reply**: Thank you for bringing this up. While strong convexity is necessary for our theoretical analysis, sparsity can still be achieved by incorporating elastic net penalties (L1+L2), as described in Proposition 2. Additionally, for L1 penalty, the technique used in Equation (5) of Section 4.2 in [SZ14] can also be applied to ReHLine.\n\n> (experimental results)\n\n**Reply**: Thank you for the comments. Following your suggestion, we have added a new benchmark to compare our proposed method with other specialized implementations in s-SVM. Please refer to our reply to the previous points for more details.\n\nHowever, we would like to emphasize that there is still a certain usefulness for all-purpose solvers. For example, even in simple problems like FairSVM [39, 40], the authors opted to use an all-purpose solver but less efficient software such as DCCP for implementation. \n\nSince our algorithm possesses some attributes similar to all-purpose solvers, including custom losses and constraints, to clarify our experimental conclusions, we would like to present them in two parts: (i) When compared to all-purpose solvers, ReHLine has achieved significant improvements; (ii) In the case of specialized algorithms/implementations, ReHLine can compete with and even outperfom specialized implementations in specific scenarios.\n\n> (minor details)\n\n**Reply**: Thanks for pointing out. $\\mathbf{E}$ stands for a matrix of all ones, and the mode-3 unfolding of a tensor $\\mathbf{U}=(u_{ijk})\\in\\mathbb{R}^{m\\times n\\times l}$ is the matrix $V=(v_{ks})\\in\\mathbb{R}^{l\\times (mn)}$, where the $k$-th row of $V$ is the vectorization of the $k$-th slice of $\\mathbf{U}$, i.e., $V_{k\\cdot}=\\mathbf{vec}(U_{\\cdot\\cdot k})^T$.\n\n## Questions\n\n> (first question already covered in the Weakness section)\n\n> (comparison with SDCA)\n\n**Reply**: The differences between ReHLine and SDCA are summarized in the global response.\n\n> (kernel methods)\n\n**Reply**: Thank you for the insightful comments. Our method can be extended to the following kernel learning problem:\n\n$$\\min_{\\beta}\\sum_{i=1}^n L_i(\\beta^TK_i)+\\frac{1}{2}\\|\\beta\\|_K^2,\\quad\\text{s.t. }A\\beta+b\\geq 0,$$\n\nwhere $K$ is a p.d. kernel matrix, $\\beta$ becomes a length-$n$ vector, and the prediction function is now $f(x_j)=\\beta^T K_j$.\n\nUsing Cholesky factorization $K=Q^TQ$, and denoting $\\alpha=Q\\beta$, it can be rewritten as:\n\n$$\\min_{\\alpha}\\sum_{i=1}^n L_i(K^T_i Q^{-1}\\alpha)+\\frac{1}{2}\\|\\alpha\\|_{2}^2,\\quad\\text{s.t. }AQ^{-1}\\alpha+b\\geq 0,$$\n\nwhich follows the form of (1) in our paper, with $x_i \\leftarrow K^T_i Q^{-1}$, $A\\leftarrow AQ^{-1}$, thus can be solved by the proposed ReHLine method.\n\nHowever, one major contribution of this article is to reduce the per-iteration cost from $\\mathcal{O}(n^2)$ to $\\mathcal{O}(nd)$. In the kernel learning $d = n$, yielding the computational complexity re-increase to $\\mathcal{O}(n^2)$. In this light, when it comes to solving problems related to kernel learning or high-dimensional problems where d is greater than n, we do not recommend using the ReHLine algorithm. We will add more discussion in the revision.\n\n[BT12] Beck and Teboulle (2012). Smoothing and First Order Methods: A Unified Framework.\n\n[SZ14] Shalev-Shwartz and Zhang (2014). Accelerated Proximal Stochastic Dual Coordinate Ascent for Regularized Loss Minimization."}, {"review_id": "LpdcyIPG8u", "reviewer": "Reviewer_7RWz", "summary": "This paper introduces an algorithm, ReHLine, which has a linear convergence rate on minimizing convex piecewise linear-quadratic loss functions. Experiments on several tasks show great performance gains over existing algorithms.", "questions": "See above.", "limitations": "See above.", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "General: The paper is clean well-written and easy to follow. Definitions are Proposition are well explained with examples.\n\n\nPros:\n\n* Introduces ReLU-ReHU decomposition and show commonly-used loss functions have their ReLU-ReHU representation.\n* Finds the dual formulation and shows that proposed algorithm has linear convergence rate by Coordinates Descent.\n* Experimental results greatly outperform existing algorithms.\n\nCons:\n\n* The main contribution of this paper is proposal of ReLU-ReHU decomposition and its resulting formulation. Its convergence result is shown by classical result. I feel the contribution is limited and marginally above the borderline.\n\n\nIn all, this paper is well written and the framework ReLU-ReHU is promising. Although the contribution seems limited but I am happy to vote for accept.\n\n----------\n\nMain body checked but did not throughly go through appendix.", "weaknesses": "See above.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "## Weaknesses\n\n> The main contribution of this paper is proposal of ReLU-ReHU decomposition and its resulting formulation. Its convergence result is shown by classical result. I feel the contribution is limited and marginally above the borderline.\n\n**Reply**: Thanks for the comments. We would like to point out that our contribution has different aspects, and the convergence result is only one of them. As can be seen in Table 1, there exist other algorithms that achieve the same convergence rate as ReHLine. However, the total computational cost is affected by both the convergence rate and the per-iteration cost. Making one of the components efficient is possible, but achieving **both** is nontrivial.\n\nTo make the linear convergence and linear computational complexity **simultaneously hold**, we need a special relation between the primal and dual variables, as is characterized by equation (9), and the objective function needs to have some specific structures to support the convergence rate. One major finding in this article is that the relation (9) is met by a very general loss function class, the convex PLQ functions, and the resulting algorithm has provable linear convergence. To the best of our knowledge, prior to ReHLine, these properties were only studied for SVMs in the LibLinear solver, and now we have greatly expanded the class of models that enjoy the linear convergence and linear computation properties."}, {"review_id": "qMsng2WHnA", "reviewer": "Reviewer_xA17", "summary": "This paper studies the Empirical Risk Minimization(ERM), which is an important framework of machine learning, and focus on a general regularized ERM based on a convex PLQ loss with linear constraints. According to the authors, the existing algorithms are faced with the problem of slow convergence or high computational cost. The authors leverage the linear property of Karush-Kuhn-Tucker conditions and propose ReHLine, whose total complexity is lower than existing algorithms.\n\n\n## Comments\nIt is recommended to indicate the unit of running time on Table 5.", "questions": "Please see my question in the \"Weakness\" section.", "limitations": "Not applicable.", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. The authors cleverly decompose convex PLQ into a series of ReLU and the so-called rectified Huber units to deal with the problem that the standard definition of convex PLQ is not convenient to optimization algorithms.When solving the decomposed optimization problem, the authors prove that their ReHLine algorithm can achieve linear convergence rate.\n2. After comparing with other state-of-the-art algorithms on multiple datasets, the experimental results show that the ReHLine algorithm can solve a variety of domain-specific problems, presenting excellent flexibility, and on large-scale datasets, the algorithm shows excellent acceleration ratio.", "weaknesses": "Actually, rather than those solvers tested in the experiment part, there are a number of general optimization solvers, including CPLEX, GUROBI, SCIP, etc. Those optimization solvers have exhibited great success in solving linear/nonlinear continuous/discrete optimization problems. Have you compared your proposed method against those general optimization solvers?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "## Weaknesses\n\n> Actually, rather than those solvers tested in the experiment part, there are a number of general optimization solvers, including CPLEX, GUROBI, SCIP, etc. Those optimization solvers have exhibited great success in solving linear/nonlinear continuous/discrete optimization problems. Have you compared your proposed method against those general optimization solvers?\n\n**Reply**: Thank you for drawing our attention to the availability of more general solvers. Solving Constraint Integer Programs (SCIP) is specifically designed for integer programming, so it does not align well with the nature of our problem. CPLEX and GUROBI are commercial solvers which require licenses to run. At present, we only obtain free licenses for conducting small-scale experiments, and they have failed for running the datasets described in our paper. To illustrate their effectiveness, we run the solvers in simulated datasets with RidgeQR compared with all other solvers (`fail` indicates that we are unable to obtain results under the current licenses). As indicated in the following table, the proposed algorithm continues to exhibit remarkable superiority over CPLEX and GUROBI.\n\n| RidgeQR | CPLEX | ReHLine | GUROBI | MOSEK | SCS | ReHLine | \n|:------|----------:|----------:|----------:|----------:|----------:|----------:|\n|n=100  | 5.073E-02 | 7.370E-05 | 1.337E-02 | 8.348E-03 | 2.082E-03 | 1.322E-04 | \n|n=300  | 1.437E-01 | 8.979E-05 | fail      | 1.024E-02 | 1.207E-02 | 4.093E-04 | \n|n=500  | fail      | 1.091E-04 | fail      | 1.260E-02 | 1.952E-02 | 4.160E-04 | \n|n=700  | fail      | 1.365E-04 | fail      | 1.522E-02 | 2.865E-02 | 5.973E-04 | \n|n=1000 | fail      | 1.870E-04 | fail      | 1.828E-02 | 2.663E-02 | 1.459E-03 |"}, {"review_id": "5ppoEEnHqN", "reviewer": "Reviewer_2E2D", "summary": "This paper proposes a new optimization algorithm for convex piecewise linear-quadratic objectives with $L_2$ regularization and linear constraints, which achieves the best known iteration complexity simultaneously with the best known per-iteration computational cost, resulting in a smaller total computational complexity than any previous algorithm. Experiments show that the proposed algorithm achieves ~1000x speedup compared to previous solvers and improves over problem-specific solvers despite its generality.", "questions": "1. Line 123 states that the proposed algorithm is inspired by LibLinear, but there is no further discussion of the relationship between ReHLine and LibLinear. What components of ReHLine are inspired by LibLinear, and what components of ReHLine are novel compared to LibLinear?\n2. The introduction mentions that this work focuses on the case that the dataset size is much larger than the input dimension and the number of constraints, but this large-scale condition is not discussed further. Why is this large-scale condition significant to your approach, and how does ReHLine compare to baseline algorithms when $n$ is small?\n3. Please clarify the questions I posed about the details of experimental evaluation in the Weaknesses section.", "limitations": "The authors provided a reasonable discussion of the limitations of their approach. A discussion of potential negative societal impact is, in my opinion, not necessary for this work.", "rating": 7, "confidence": 3, "soundness": 4, "presentation": 3, "contribution": 4, "strengths": "1. The considered optimization problem is broad, encompassing training for many linear models (SVM, LAD, SVR, QR) and constraint sets, including fairness.\n2. The technique appears to be novel.\n3. The comparison of total complexity (Table 2) demonstrates a significant theoretical advantage over previous algorithms.\n4. The empirical comparison (Table 5) shows that the proposed algorithm enjoys a significant practical advantage over previous algorithms.", "weaknesses": "1. The experimental evaluation is missing a few details that need clarifying. For example, line 212 states that the BenchOpt framework is used \"to implement optimization benchmarks for all the SOTA solvers\". I find this statement a little unclear. Did the authors implement all of the baseline solvers from scratch? If not, what do they mean by \"implement optimization benchmarks\"? Also, is it common to implement these solvers in Python? My understanding is that the production solvers are most commonly implemented in C/C++, and using a slower language like Python might amplify the difference in speed between algorithms.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "## Weaknesses\n\n> The experimental evaluation is missing a few details that need clarifying. For example, line 212 states that the BenchOpt framework is used \"to implement optimization benchmarks for all the SOTA solvers\". I find this statement a little unclear. Did the authors implement all of the baseline solvers from scratch? If not, what do they mean by \"implement optimization benchmarks\"?\n\n**Reply**: We apologize for the confusion. To clarify, we would like to mention that we do not develop baseline  solvers from scratch in our study. Instead, we utilize well-established  software for benchmarking purposes. We appreciate your suggestion and will revise the manuscript accordingly in order to clarify any potential  misleading description.\n\n> Also, is it common to implement these solvers in Python? My understanding is that the production solvers are most commonly implemented in C/C++, and using a slower language like Python might amplify the difference in speed between algorithms.\n\n**Reply**: Thank you for the thorough comments. All the methods that we compared, including our own, are fundamentally built upon C/C++ backends. Python just serve as a facilitator by providing an API interface that enhances their usability. See more description of backends in CVXPY and CVXCORE. Python acts as a bridge between the users and the underlying C/C++ backends, enabling a seamless and more user-friendly experience. Therefore, our Python experiments do not significantly impact the execution speed of the software itself. As a result, the experimental results are reliable and reproducible. We will add more descriptions in our manuscript.\n\n## Questions\n\n> Line 123 states that the proposed algorithm is inspired by LibLinear, but there is no further discussion of the relationship between ReHLine and LibLinear. What components of ReHLine are inspired by LibLinear, and what components of ReHLine are novel compared to LibLinear?\n\n**Reply**: The development of ReHLine was largely inspired by the linear computational complexity and linear convergence property of LibLinear, which contribute to the great success of LibLinear in practice. However, in LibLinear, these properties are developed specifically ONLY for SVMs.\n\nOne major merit of ReHLine is that these good properties actually apply to a much broader loss function class (with or without constraints), characterized by the convex PLQ functions. Therefore, with ReHLine, the efficiency of LibLinear can be extended to many other ERM problems, and LibLinear can be viewed as a special case of ReHLine.\n\n> The introduction mentions that this work focuses on the case that the dataset size is much larger than the input dimension and the number of constraints, but this large-scale condition is not discussed further. Why is this large-scale condition significant to your approach, and how does ReHLine compare to baseline algorithms when $n$ is small?\n\n**Reply**: As is shown in Table 1, one major contribution of this article is to reduce the per-iteration cost from $\\mathcal{O}(n^2)$ to $\\mathcal{O}(n)$, so that ReHLine demonstrates great advantages when $n$ is large and forms the bottleneck of computation. It can be seen that equation (8) is a quadratic function of $(K+nL+nH)$ variables. If $n$ is small, then the scale of the quadratic programming problem (7) is also small, which means that a general-purpose solver may already be sufficient."}, {"review_id": "S84DxGggLw", "reviewer": "Reviewer_gb3q", "summary": "This paper proposes a dual-coordinate descent solver for a class of ERM (Empricial Risk Minimization) problem with general linear inequality constraint, which is of linear convergence rate and efficient coordinate-update cost of O(n) (n is the number of samples). The main contribution is it extends existing method of dual-coordinate descent to handle problem with additional linear inequality constraints.", "questions": "To my knowledge, Liblinear is a special case of solver proposed in this paper, speicificaly designed for the SVM optimization problem.\nThen why could solver proposed in this paper perform more efficiently than Liblinear on SVM problem? Is it implementation issue? or their objective functions are different? The author should have made this clear to avoid confusion. ", "limitations": "The paper does not need to limit itself to the case of Relu-Rehu function.\nInstead, any smooth loss (in the primal) yields strongly-convex dual problem. Therefore can apply same optimization method convergence rate proposed in the paper. (ref: https://www.jmlr.org/papers/volume14/shalev-shwartz13a/shalev-shwartz13a.pdf)", "rating": 5, "confidence": 5, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "The paper has solid theoretical motivation, derivation and analysis, and has presented the idea clearly with solid experiments support on pratical optimization problems.", "weaknesses": "The paper completely missed discussion on existing works of dual-coordinate ascent on general ERM problem. There is a thread of works such as:\n\nShalev-Shwartz, Shai, and Tong Zhang. \"Stochastic dual coordinate ascent methods for regularized loss minimization.\" Journal of Machine Learning Research 14.1 (2013).\n\non dual-coordinate ascent method for general ERM. The novelty of problem setup proposed in this reviewed paper is it adds an additional linear inequality constraints. The \"piecewise-quadratic loss\" is simply a special case of \"smooth loss\" of the above mentioned paper and SAME convergence rate applies to all this type of problem, so the author chould have made this work more general (by extending the formulation of the above work by simply adding linear inequliaty).", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "## Weaknesses\n\n> The paper completely missed discussion on existing works of dual-coordinate ascent on general ERM problem. There is a thread of works such as:\n>\n>> Shalev-Shwartz, Shai, and Tong Zhang. \"Stochastic dual coordinate ascent methods for regularized loss minimization.\" Journal of Machine Learning Research 14.1 (2013).\n>\n> on dual-coordinate ascent method for general ERM. The novelty of problem setup proposed in this reviewed paper is it adds an additional linear inequality constraints. The \"piecewise-quadratic loss\" is simply a special case of \"smooth loss\" of the above mentioned paper and SAME convergence rate applies to all this type of problem, so the author chould have made this work more general (by extending the formulation of the above work by simply adding linear inequliaty).\n\n**Reply**: Thank you for bringing the reference for our attention. We will add the reference into our revision; see our global response for the discussions. In terms of novelty, we have the following important differences compared with the algorithm in the reference.\n\n1. **Convergence rates are different**. In general, the \"convex piecewise-quadratic loss\" or the proposed ReHLine loss is NOT \"smooth loss\" but only \"Lipschitz loss\", thus the results in [SZ13] yielding that a sub-linear convergence (which is suboptimal compared to the linear convergence of Theorem 3 in our manuscript). Although refined analysis for \"almost smooth function\" is provided in [SZ13], the refined results, even for linear SVM, depend on the assumptions of underlying data distribution (see discussion in Section 5 and the definition of $N(u)$ in Theorem 16 in their paper).\n2. **The technical details of the two methods are different**. [SZ13] and ours do share some similarities at a structual level. However, they differ in the details of the iterative steps. Specifically, [SZ13] directly solves the convex conjugate of the loss function, while our method further decomposes this loss function into simpler components for multi-step updates. Intuitively, our approach further leverages the decomposability and linearity of the loss function to perform multiple iterative updates on the decomposed components.\n3. **Improve the practical feasibility of convex conjugate**. In fact, in each iteration, solving and implementing the convex conjugate of a general function is not an easy task, especially when the function is a composition of multiple different functions, also see the discussion in [SZ14]: \"In general, this optimization problem is still not necessarily simple to solve because $\\phi^*$ may also be complex.\" In comparison, our proposed ReLU-ReHU decomposition offers better implementation capabilities and property (e.g., Proposition 1).\n4. As you mentioned, they cannot handle constrained optimization problems.\n\n## Questions\n\n> To my knowledge, Liblinear is a special case of solver proposed in this paper, speicificaly designed for the SVM optimization problem. Then why could solver proposed in this paper perform more efficiently than Liblinear on SVM problem? Is it implementation issue? Or their objective functions are different? The author should have made this clear to avoid confusion.\n\n**Reply**: Thank you for the comments. It is indeed correct that the algorithm used in LibLinear to solve SVM is a special case of our method. We make this claim in order to highlight the efficiency of our implementation/software, as we regard software as a significant contribution of our work. We will revise the manuscript to clarify the confusion that this advantage is caused by implementation. \n\n## Limitations:\n\n> The paper does not need to limit itself to the case of Relu-Rehu function. Instead, any smooth loss (in the primal) yields strongly-convex dual problem. Therefore can apply same optimization method convergence rate proposed in the paper. (ref: https://www.jmlr.org/papers/volume14/shalev-shwartz13a/shalev-shwartz13a.pdf)\n\n**Reply**: Thank you for the comments. As we mentioned in the previous point, our optimization objective is not a smooth function, and therefore does not readily yield linear convergence results based on the theoretical results in [SZ13]. In addition to theoretical considerations, we also place great emphasis on the algorithm's computational efficiency. Specifically, in the proposed framework, each coordinate can be solved with a simple analytic solution. This feature has demonstrated excellent efficacy and practical feasibility in our experiments.\n\n[SZ13] Shalev-Shwartz and Zhang (2013). Stochastic Dual Coordinate Ascent Methods for Regularized Loss Minimization.\n\n[SZ14] Shalev-Shwartz and Zhang (2014). Accelerated Proximal Stochastic Dual Coordinate Ascent for Regularized Loss Minimization."}]

Author rebuttal:
# To All Reviewers

Thank you all reviewers for the encouraging and insightful comments. We appreciate the time and effort the reviewers have dedicated to providing valuable feedback on our manuscript. In this round, we have made every effort to address all the comments of the reviewers. **The point-to-point responses are provided in the reply**.

Additionally, we would like to take this opportunity to highlight our contributions in software development. We have dedicated a significant amount of time to optimize our C/C++ implementation, Python API, quick routines of various losses and constraints, and provide comprehensive documentation. Our motivation behind these efforts stems from the fact that there are currently not many practical software solutions available. For instance, in simple problems like FairSVM [39,40], the authors opted to use generic but less efficient software such as DCCP for implementation. Given that our software not only operates efficiently but also offers a high level of flexibility and practicality in various applications, we sincerely hope that the reviewers will take into consideration our contribution to software when evaluating this article.

As suggested by the reviewers, one point we want to highlight globally is the added discussions on the existing literature. In our revised manuscript, we have included more related works and discussed their connections with ReHLine. Below is an excerpt of the planned revision:

> There are multiple existing works developed to tackle ERM problems, such as SAG, SVRG, SAGA, and SDCA. However, SAG, SVRG, and SAGA can only handle smooth loss functions, with an optional non-smooth term that has an easy proximal operator. SDCA applies to general loss functions, but it requests the convex conjugate of loss functions, which is not necessarily simple to compute. Moreover, it only guarantees a sub-linear convergence rate for non-smooth loss functions. In contrast, ReHLine supports all convex PLQ loss functions with optional general linear constraints, which are non-smooth by construction, and they all enjoy a linear computational complexity and provable linear convergence.

>
> For non-smooth loss functions, another existing method to solve ERM is the smoothing technique (Beck and Teboulle, 2012), which approximates the non-smooth terms by smooth functions, and then uses gradient-based algorithms to solve the smooth problem. However, the choice of the smoothing function and smoothing parameter typically requires additional knowledge, and the linear convergence rate is not always guaranteed.

Furthermore, we would like to supplement the following important differences between ReHLine and SDCA.

1. **Convergence rates are different for non-smooth losses**. In general, the PLQ or the proposed ReHLine loss is non-smooth but Lipschitz function, thus the theoretical results of SDCA [SZ13] yield a *sub-linear* convergence (which is suboptimal compared to the linear convergence of Theorem 3 in our manuscript). Although refined analysis for "almost smooth function" is provided in [SZ13], the refined results, even for linear SVM, depend on the assumptions of underlying data distribution (see the discussion in Section 5 and the definition of $N(u)$ in Theorem 16 of [SZ13]).
2. **The technical details of the two methods are different**. [SZ13] and ours do share some similarities at a structual level. However, they differ in the details of the iterative steps. Specifically, [SZ13] directly solves the convex conjugate of the loss function, while our method further decomposes this loss function into simpler components for multi-step updates. Intuitively, our approach further leverages the decomposability and linearity of the loss function to perform multiple iterative updates on the decomposed components.
3. **Improve the practical feasibility of convex conjugate**. In fact, solving and implementing the convex conjugate of a general function is not an easy task, especially when the function is a composition of multiple different functions. In comparison, our proposed ReLU-ReHU decomposition offers better implementation capabilities and properties (e.g., Proposition 1).
4. SDCA cannot handle **constrained** optimization problems.

[SZ13] Shalev-Shwartz and Zhang (2013). Stochastic Dual Coordinate Ascent Methods for Regularized Loss Minimization.

Meta-review: At the end of the discussion phase, all reviewers recommended an acceptance. I would recommend that the authors carefully consider the various comments of the reviewers when preparing the final version, such as carefully discussing the prior work (such as SDCA) and toning down the exaggerated performance claim (there are various standard low-iteration-cost methods that could be applied instead of MOSEK so the 1000x faster claim is not believable).
Decision: Accept (poster)