An active-set method for sparse approximations. Part II: General piecewise-linear terms

TL;DR

This paper introduces an efficient active-set method tailored for convex quadratic programming problems with general piecewise-linear terms, significantly improving memory efficiency and robustness for applications like risk minimization and sparse approximation.

Contribution

It develops a novel active-set algorithm that exploits piecewise-linear structure to reduce memory use and enhance computational efficiency in convex quadratic programming.

Findings

Successfully applied to risk-averse portfolio selection

Effective in quantile regression and SVM classification

Outperforms general-purpose solvers on real datasets

Abstract

In this paper we present an efficient active-set method for the solution of convex quadratic programming problems with general piecewise-linear terms in the objective, with applications to sparse approximations and risk-minimization. The method exploits the structure of the piecewise-linear terms appearing in the objective in order to significantly reduce its memory requirements, and thus improve its efficiency. We showcase the robustness of the proposed solver on a variety of problems arising in risk-averse portfolio selection, quantile regression, and binary classification via linear support vector machines. We provide computational evidence to demonstrate, on real-world datasets, the ability of the solver of efficiently handling a variety of problems, by comparing it against an efficient general-purpose interior point solver as well as a state-of-the-art alternating direction method of multipliers. This work complements the accompanying paper [``An active-set method for sparse approximations. Part I: Separable $\ell_1$ terms", S. Pougkakiotis, J. Gondzio, D. S. Kalogerias], in which we discuss the case of separable $\ell_1$ terms, analyze the convergence, and propose general-purpose preconditioning strategies for the solution of its associated linear systems.

Tables

Table 1. Table 1 : Portfolio optimization datasets.

Name	# of assets	# of data points	Timeline
DowJones	28	1363	Feb. 1990–Apr. 2016
NASDAQ100	82	596	Nov. 2004–Apr. 2016
FTSE100	83	717	Jul. 2002–Apr. 2016
FF49Industries	49	2325	Jul. 1969–Jul. 2015
SP500	442	595	Nov. 2004–Apr. 2016
NASDAQComp	1203	685	Feb. 2003–Apr. 2016

Table 2. Table 2 : CVaR portfolio selection: varying confidence level ( tol = 10 − 5 superscript 10 5 10^{-5} , τ = 10 − 2 𝜏 superscript 10 2 \tau=10^{-2} ).

Dataset	$𝜶$	Iterations			Time (s)
missingmissingmissing  3-5
missingmissingmissing  6-8
	PMM(SSN)[Fact.]	IP–PMM	OSQP	AS	IP–PMM	OSQP
DowJones	$0.05$	35(169)[142]	18	32,575	0.65	0.71	10.94
	$0.10$	36(166)[145]	21	42,625	0.75	0.83	13.61
	$0.15$	36(144)[125]	13	43,100	0.67	0.53	14.45
NASDAQ100	$0.05$	31(156)[138]	19	30,575	0.66	0.60	7.92
	$0.10$	33(168)[150]	16	34,700	0.65	0.50	8.81
	$0.15$	33(163)[142]	15	40,500	0.62	0.57	10.46
FTSE100	$0.05$	32(179)[164]	17	21,825	0.77	0.90	10.96
	$0.10$	34(171)[147]	19	35,825	0.80	0.99	17.70
	$0.15$	35(179)[156]	11	36,550	0.93	0.59	18.24
FF49Industries	$0.05$	41(215)[186]	27	$50,{000}^{‡}$1	1.77	4.65	${59.16}^{‡}$
	$0.10$	42(224)[183]	22	$50,{000}^{‡}$	2.37	3.62	${57.92}^{‡}$
	$0.15$	39(205)[157]	11	$50,{000}^{‡}$	2.66	2.18	${51.25}^{‡}$
SP500	$0.05$	30(203)[193]	20	21,075	5.93	12.82	65.70
	$0.10$	34(199)[187]	11	32,375	5.88	6.98	99.39
	$0.15$	34(173)[160]	17	39,225	5.78	10.17	119.95
NASDAQComp	$0.05$	28(198)[192]	10	16,100	22.07	35.52	143.82
	$0.10$	31(171)[169]	14	31,275	23.31	48.99	291.97
	$0.15$	34(173)[167]	19	44,775	21.67	67.14	403.75

Table 3. Table 3 : CVaR portfolio selection: varying confidence level ( tol = 10 − 5 superscript 10 5 10^{-5} , τ = 10 − 1 𝜏 superscript 10 1 \tau=10^{-1} ).

Dataset	$𝜶$	Iterations			Time (s)
missingmissingmissing  3-5
missingmissingmissing  6-8
	PMM(SSN)[Fact.]	IP–PMM	OSQP	AS	IP–PMM	OSQP
DowJones	$0.05$	35(168)[144]	21	24,975	0.57	0.85	8.31
	$0.10$	36(154)[133]	20	21,925	0.64	0.79	7.32
	$0.15$	36(143)[128]	17	29,625	0.68	0.69	9.95
NASDAQ100	$0.05$	34(157)[124]	19	17,225	0.46	0.57	4.50
	$0.10$	34(155)[133]	16	18,100	0.59	0.49	4.80
	$0.15$	35(157)[132]	16	20,775	0.56	0.48	5.49
FTSE100	$0.05$	34(166)[137]	21	18,900	0.65	1.17	10.10
	$0.10$	34(174)[157]	16	24,750	0.92	0.86	13.15
	$0.15$	36(166)[147]	17	28,850	0.95	0.92	15.53
FF49Industries	$0.05$	43(215)[171]	19	48,000	1.63	3.00	59.33
	$0.10$	41(184)[146]	22	40,475	1.57	3.53	48.48
	$0.15$	41(160)[136]	19	35,525	1.45	3.06	40.91
SP500	$0.05$	32(163)[144]	20	21,425	4.68	12.26	66.29
	$0.10$	33(165)[143]	18	22,100	5.86	12.43	73.56
	$0.15$	35(169)[148]	19	32,225	5.61	11.51	98.25
NASDAQComp	$0.05$	32(190)[182]	26	17,200	15.55	90.78	157.83
	$0.10$	34(180)[172]	15	13,925	22.75	51.81	127.26
	$0.15$	34(170)[167]	18	14,675	24.69	62.40	137.42

Table 4. Table 4 : MAsD portfolio selection: varying regularization ( tol = 10 − 5 superscript 10 5 10^{-5} ).

Dataset	$𝝉$	Iterations			Time (s)
missingmissingmissing  3-5
missingmissingmissing  6-8
	PMM(SSN)[Fact.]	IP–PMM	OSQP	AS	IP–PMM	OSQP
DowJones	$0.01$	51(142)[139]	10	32,350	1.05	0.56	10.47
	$0.05$	52(142)[139]	13	14,050	0.72	0.54	4.60
NASDAQ100	$0.01$	47(167)[160]	11	18,800	0.67	0.26	5.02
	$0.05$	48(176)[161]	10	40,550	0.67	0.27	10.72
FTSE100	$0.01$	46(153)[150]	10	18,325	0.79	0.27	9.56
	$0.05$	46(153)[150]	14	47,125	0.74	0.35	23.29
FF49Industries	$0.01$	53(141)[136]	11	$50,{000}^{‡}$1	1.49	1.72	${57.59}^{‡}$
	$0.05$	52(137)[132]	9	12,050	1.35	1.51	14.24
SP500	$0.01$	42(165)[163]	17	28,375	5.33	2.91	87.10
	$0.05$	41(157)[153]	17	41,275	4.39	3.07	128.09
NASDAQComp	$0.01$	44(143)[133]	17	34,175	9.75	9.91	310.12
	$0.05$	44(143)[132]	18	41,000	10.22	9.72	367.32

Table 5. Table 5 : Quantile regression datasets.

Name	# of data points	# of features
space_ga	3107	6
abalone	4177	8
cpusmall	8192	12
cadata	20,640	8
E2006-tfidf	16,087	150,360

Table 6. Table 6 : Quantile regression: varying confidence level ( tol = 10 − 4 superscript 10 4 10^{-4} , λ = 10 − 2 𝜆 superscript 10 2 \lambda=10^{-2} , τ = 0.5 𝜏 0.5 \tau=0.5 ).

Problem	$𝜶$	Iterations			Time (s)
missingmissingmissing  3-5
missingmissingmissing  6-8
	PMM(SSN)[Fact.]	IP–PMM	OSQ	AS	IP–PMM	OSQP
space_ga	$0.50$	15(43)[30]	19	2350	0.38	0.35	0.51
	$0.65$	15(50)[35]	20	16,300	0.37	0.39	3.36
	$0.80$	15(62)[48]	19	16,975	0.46	0.35	3.46
	$0.95$	15(78)[77]	20	7400	0.48	0.36	1.51
abalone	$0.50$	14(68)[64]	35	6350	1.02	1.29	2.15
	$0.65$	14(68)[63]	35	6300	0.89	1.36	2.17
	$0.80$	14(79)[74]	27	12,425	0.89	1.01	4.13
	$0.95$	14(87)[79]	18	2875	0.82	0.69	0.97
cpusmall	$0.50$	15(64)[64]	29	$⩞$1	2.01	2.57	$⩞$
	$0.65$	16(74)[74]	30	$⩞$	2.33	2.81	$⩞$
	$0.80$	16(86)[85]	26	$⩞$	2.64	2.37	$⩞$
	$0.95$	15(110)[110]	22	$⩞$	2.85	1.96	$⩞$
cadata	$0.50$	3(66)[65]	49	7125	2.96	14.63	16.10
	$0.65$	4(85)[83]	42	9050	3.65	12.43	20.27
	$0.80$	3(77)[76]	45	7275	3.45	13.13	16.45
	$0.95$	3(225)[224]	80	7450	13.14	23.32	17.08
E2006-tfidf	$0.50$	14(27)[20]	$†$2	$†$	84.82	$†$	$†$
	$0.65$	15(34)[26]	$†$	$†$	95.08	$†$	$†$
	$0.80$	16(46)[35]	$†$	$†$	112.24	$†$	$†$
	$0.95$	17(79)[71]	$†$	$†$	165.78	$†$	$†$

Table 7. Table 7 : Quantile regression: varying regularization ( tol = 10 − 4 superscript 10 4 10^{-4} , α = 0.8 𝛼 0.8 \alpha=0.8 ).

Problem	$(𝝉,𝝀)$	Iterations			Time (s)
missingmissingmissing  3-5
missingmissingmissing  6-8
	PMM(SSN)[Fact.]	IP–PMM	OSQP	AS	IP–PMM	OSQP
space_ga	$(0.2,5\cdot {10}^{-2})$	15(58)[49]	16	37,125	0.51	0.30	8.00
	$(0.4,1\cdot {10}^{-2})$	15(61)[48]	26	12,050	0.46	0.49	2.48
	$(0.6,5\cdot {10}^{-3})$	15(67)[46]	25	13,025	0.44	0.46	2.68
	$(0.8,1\cdot {10}^{-3})$	15(82)[62]	27	20,575	0.61	0.49	4.21
abalone	$(0.2,5\cdot {10}^{-2})$	14(70)[65]	25	16,225	0.88	1.32	5.49
	$(0.4,1\cdot {10}^{-2})$	14(71)[67]	34	15,575	0.86	1.48	5.16
	$(0.6,5\cdot {10}^{-3})$	15(80)[76]	39	18,675	0.92	1.50	6.26
	$(0.8,1\cdot {10}^{-3})$	15(116)[109]	30	9975	1.44	1.16	3.34
cpusmall	$(0.2,5\cdot {10}^{-2})$	15(83)[81]	20	$⩞$1	1.98	2.14	$⩞$
	$(0.4,1\cdot {10}^{-2})$	15(81)[78]	26	$⩞$	2.10	2.39	$⩞$
	$(0.6,5\cdot {10}^{-3})$	16(101)[94]	21	$⩞$	2.79	1.94	$⩞$
	$(0.8,1\cdot {10}^{-3})$	15(106)[103]	20	$⩞$	3.16	1.85	$⩞$
cadata	$(0.2,5\cdot {10}^{-2})$	9(147)[143]	62	5875	15.42	18.05	13.18
	$(0.4,1\cdot {10}^{-2})$	4(63)[62]	50	7350	2.76	14.68	17.41
	$(0.6,5\cdot {10}^{-3})$	5(121)[118]	45	9775	5.49	13.32	21.70
	$(0.8,1\cdot {10}^{-3})$	7(276)[276]	33	875	14.01	9.91	2.08
E2006-tfidf	$(0.2,5\cdot {10}^{-2})$	16(47)[38]	$†$2	$†$	115.87	$†$	$†$
	$(0.4,1\cdot {10}^{-2})$	16(43)[36]	$†$	$†$	118.89	$†$	$†$
	$(0.6,5\cdot {10}^{-3})$	16(43)[34]	$†$	$†$	115.90	$†$	$†$
	$(0.8,1\cdot {10}^{-3})$	17(48)[40]	$†$	$†$	128.38	$†$	$†$

Table 8. Table 8 : Binary classification datasets.

Name	# of training points	# of features
rcv1	20,242	47,236
real-sim	72,309	20,958
news20	19,996	1,355,191

Table 9. Table 9 : Binary classification via elastic-net SVMs: varying regularization ( tol = 10 − 5 superscript 10 5 10^{-5} , λ = 10 − 2 𝜆 superscript 10 2 \lambda=10^{-2} ).

Problem	${𝝉}_{\mathrm{𝟏}}$	${𝝉}_{\mathrm{𝟐}}$	Iterations	Time (s)
missingmissingmissing  4-4
missingmissingmissing  5-5
	PMM(SSN)[Fact.]	AS
rcv1	0.2	0.2	48(139)[105]	43.36
	0.8	0.2	47(100)[58]	23.28
	0.2	0.8	47(165)[133]	54.81
	5	5	39(85)[45]	20.57
real-sim	0.2	0.2	50(130)[100]	252.25
	0.8	0.2	47(85)[48]	132.20
	0.2	0.8	48(101)[68]	201.56
	5	5	47(85)[48]	128.75
news20	0.2	0.2	50(142)[102]	213.12
	0.8	0.2	41(85)[43]	133.94
	0.2	0.8	70(212)[161]	74.42
	5	5	41(85)[43]	130.33

Equations

$$\begin{array}{rl}\underset{x\in {\mathbb{R}}^n}{\min }& \left\{c^{\top }x+\frac{1}{2}{x}^{\top }Qx+\sum _{i=1}^l{\left({\left(Cx+d\right)}_i\right)}_{+}+\Vert Dx{\Vert }_1+{\delta }_{\mathcal{K}}(x)\right\},\\ \text{s.t.}& Ax=b,\end{array}$$

$$\begin{array}{rl}\underset{(x,w)\in {\mathbb{R}}^n\times {\mathbb{R}}^l}{\min }& \left\{c^{\top }x+\frac{1}{2}{x}^{\top }Qx+\sum _{i=1}^l{\left(w_i\right)}_{+}+\Vert Dx{\Vert }_1+{\delta }_{\mathcal{K}}(x)\right\},\\ \text{s.t.}& Cx+d-w=0_l,\\ & Ax=b.\end{array}$$

$$\Vert Cx+d{\Vert }_1\equiv -\text{mathds}{1}_l^{\top }(Cx+d)+\sum _{j=1}^l{\left(2(Cx+d{)}_j\right)}_{+},$$

$$\sum _{i=1}^l\max \left\{{\left(C_{1}x+d_1\right)}_i,{\left(C_{2}x+d_2\right)}_i\right\},$$

$$\text{mathds}{1}_l^{\top }\left(C_{2}x+d_2\right)+\sum _{i=1}^l{\left({\left(C_{1}x+d_1-C_{2}x-d_2\right)}_i\right)}_{+}\equiv \sum _{i=1}^l\max \left\{{\left(C_{1}x+d_1\right)}_i,{\left(C_{2}x+d_2\right)}_i\right\}.$$

$$\begin{array}{rl}{\mathcal{L}}_{\beta }(x,w;y,z):=& c^{\top }x+\frac{1}{2}{x}^{\top }Qx+\sum _{i=1}^l{\left(w_i\right)}_{+}+\Vert Dx{\Vert }_1-y^{\top }\left(\left[\begin{array}{c}Cx+d-w\\ Ax-b\end{array}\right]\right)\\ & +\frac{\beta }{2}{\Vert \left[\begin{array}{c}Cx+d-w\\ Ax-b\end{array}\right]\Vert }^2-\frac{1}{2\beta }\Vert z{\Vert }^2+\frac{1}{2\beta }\Vert z+\beta x-\beta {\Pi }_{\mathcal{K}}({\beta }^{-1}z+x){\Vert }^2.\end{array}$$

$$\begin{split}\scalebox{1.0}{\mbox{$\displaystyle\phi(x,w)\equiv\ \phi_{\beta_{k},\rho_{k}}(x,w;x_{k},y_{k},z_{k})\coloneqq$}}&\ \scalebox{1.0}{\mbox{$\displaystyle\mathcal{L}_{\beta_{k}}(x,w;y_{k},z_{k})+\frac{1}{2\rho_{k}}\left\|x-x_{k}\right\|^{2}-g\left(x,w\right)$}}.\end{split}$$

$$\underset{x,w}{\min }\left\{{\mathcal{L}}_{{\beta }_k}(x,w;y_k,z_k)+\frac{1}{2{\rho }_k}{\Vert x-x_k\Vert }^2\right\}\equiv \underset{x,w}{\min }\left\{\varphi (x,w)+g(x,w)\right\},$$

$$\scalebox{1.0}{\mbox{$\displaystyle\left(\nabla\phi(x^{*},w^{*})\right)^{\top}\left((x,w)-(x^{*},w^{*})\right)+g(x,w)-g(x^{*},w^{*})\geq 0,\quad\forall\ (x,w)\in\mathbb{R}^{n}\times\mathbb{R}^{l}$}}.$$

$$\begin{array}{rl}{\nabla }_x\varphi (x,t,w)=& c+Qx-\left[\begin{array}{cc}{C}^{\top }& A^{\top }\end{array}\right]y_k+{\beta }_k\left[\begin{array}{cc}{C}^{\top }& A^{\top }\end{array}\right]\left(\left[\begin{array}{c}Cx+d-w\\ Ax-b\end{array}\right]\right)\\ & +(z_k+{\beta }_{k}x)-{\beta }_k{\Pi }_{\mathcal{K}}({\beta }_k^{-1}{z}_k+x)+{\rho }_k^{-1}(x-x_k),\\ {\nabla }_w\varphi (x,t,w)=& (y_k{)}_{1:l}-{\beta }_k\left(Cx+d-w\right).\end{array}$$

$$(0_{n+l},0_{l+m})\in {\mathcal{M}}_{{\beta }_k,{\rho }_k}(x,w,y;x_k,y_k,z_k)\equiv {\mathcal{M}}_k(x,w,y),$$

$$\begin{split}\mathcal{M}_{k}(x,w,y)\coloneqq\left\{\scalebox{1.0}{\mbox{$\displaystyle(u^{\prime},v^{\prime})\colon u^{\prime}\in r_{k}(x,y)+\begin{bmatrix}\partial g_{1}(x)\\ \partial g_{2}(w)\end{bmatrix},\ v^{\prime}=\begin{bmatrix}Cx+d-w\\ Ax-b\end{bmatrix}+\beta_{k}^{-1}(y-y_{k})$}}\right\},\end{split}$$

$$r_{k}(x,y)\coloneqq\begin{bmatrix}\scalebox{1.0}{\mbox{$\displaystyle c+Qx-\begin{bmatrix}C^{\top}&A^{\top}\end{bmatrix}y+(z_{k}+\beta_{k}x)-\beta_{k}\Pi_{\mathcal{K}}(\beta_{k}^{-1}z_{k}+x)$}}+\rho_{k}^{-1}(x-x_{k})\\ y_{1:l}\end{bmatrix}.$$

$${\epsilon }_k:=\text{dist}\left(0_{n+2l+m},{\mathcal{M}}_k\left(x_{k+1},w_{k+1},y_{k+1}\right)\right)\le {ϵ}_k,$$

$${\epsilon }_k=\Vert \left[\begin{array}{c}\hat{r}+{\Pi }_{\mathcal{U}}\left(-\hat{r}\right)\\ \left[\begin{array}{c}C\\ A\end{array}\right]x_{k+1}+\left[\begin{array}{c}d\\ -b\end{array}\right]-\left[\begin{array}{c}{I}_l\\ 0_{m,l}\end{array}\right]w_{k+1}+{\beta }_k^{-1}(y_{k+1}-y_k)\end{array}\right]\Vert .$$

$$z_{k+1}$$

$${\beta }_{k+1}$$

$${\hat{\mathcal{M}}}_k\left(x,w,y\right):={\zeta }_k\left[\begin{array}{c}{\zeta }_k^{-1}\left(\left[\begin{array}{c}x\\ w\end{array}\right]-{\text{prox}}_{{\zeta }_{k}g}\left(x-{\zeta }_k{r}_k(x,y)\right)\right)\\ \left[\begin{array}{c}C\\ A\end{array}\right]x+\left[\begin{array}{c}d\\ -b\end{array}\right]-\left[\begin{array}{c}{I}_l\\ 0_{m,l}\end{array}\right]w+{\beta }_k^{-1}(y-y_k)\end{array}\right]=\left[\begin{array}{c}{0}_{n+l}\\ 0_{l+m}\end{array}\right].$$

$$\underbrace{\scalebox{1.0}{\mbox{$\displaystyle\begin{bmatrix}H_{k_{j}}&0_{n,l}&-\zeta_{k}B_{g_{1},k_{j}}\begin{bmatrix}C^{\top}&A^{\top}\end{bmatrix}\\ 0_{l,n}&\left(I_{l}-B_{g_{2},k_{j}}\right)&\zeta_{k}\begin{bmatrix}B_{g_{2},k_{j}}&0_{l,m}\end{bmatrix}\\ \zeta_{k}\begin{bmatrix}C\\ A\end{bmatrix}&\zeta_{k}\begin{bmatrix}-I_{l}\\ 0_{m,l}\end{bmatrix}&\zeta_{k}\beta_{k}^{-1}I_{l+m}\end{bmatrix}$}}}_{{M_{k_{j}}}}\begin{bmatrix}d_{x}\\ d_{w}\\ d_{y}\end{bmatrix}=-\widehat{\mathcal{M}}_{k}\left(x_{k_{j}},w_{k_{j}},y_{k_{j}}\right),$$

$$H_{k_j}:=I_n-B_{g_1,k_j}+{\zeta }_k{\beta }_k{B}_{g_1,k_j}\left(\left(1+{\rho }_k^{-1}{\beta }_k^{-1}\right)I_n-B_{\delta ,k_j}+{\beta }_k^{-1}Q\right),$$

$$\begin{array}{rl}& B_{\delta ,k_j}\in {\partial }_x^C{\Pi }_{\mathcal{K}}\left({\beta }_k^{-1}{z}_k+x_{k_j}\right),\\ & B_{g_1,k_j}\in {\partial }_x^C\left({\text{prox}}_{{\zeta }_k{g}_1}\left(x_{k_j}-{\zeta }_k(r_k(x_{k_j},y_{k_j}){)}_{1:n}\right)\right),\\ & B_{g_2,k_j}\in {\partial }_w^C\left({\text{prox}}_{{\zeta }_k{g}_2}\left(w_{k_j}-{\zeta }_k(y_{k_j}{)}_{1:l}\right)\right).\end{array}$$

$${\partial }_{u_i}^C\left({\Pi }_{[a_i,b_i]}(u_i)\right)=\{\begin{array}{ll}\{1\},& \text{if}{u}_i\in (a_i,b_i),\\ \{0\},& \text{if}{u}_i\notin [a_i,b_i],\\ [0,1],& \text{if}{u}_i\in \{a_i,b_i\}.\end{array}$$

$${\left({\text{prox}}_{{\zeta }_k{g}_1}\left(u\right)\right)}_i={\text{soft}}_{\left[-{\zeta }_k{D}_{(i,i)},{\zeta }_k{D}_{(i,i)}\right]}\left(u_i\right)\equiv \max \left\{|u_i|-{\zeta }_k{D}_{(i,i)},0\right\}\text{sign}(u_i).$$

$${\left({\text{prox}}_{{\zeta }_k{g}_2}\left(u\right)\right)}_i={\text{soft}}_{\left[0,{\zeta }_k\right]}\left(u_i\right)\equiv \max \{u_i-{\zeta }_k,0\}+\min \{u_i,0\}.$$

$${\left({\partial }_u^C\left({\text{prox}}_{{\zeta }_k{g}_1}\left(u\right)\right)\right)}_i=\{\begin{array}{ll}\{1\},& \text{if}|u_i|>{\zeta }_k{D}_{(i,i)}\text{or}{D}_{(i,i)}=0,\\ \{0\},& \text{if}|u_i|<{\zeta }_k{D}_{(i,i)},\\ [0,1],& \text{if}{u}_i={\zeta }_k{D}_{(i,i)},\end{array}\text{for all }i\in \{1,\dots ,n\},$$

$${\left({\partial }_u^C\left({\text{prox}}_{{\zeta }_k{g}_2}\left(u\right)\right)\right)}_i=\{\begin{array}{ll}\{1\},& \text{if}{u}_i>{\zeta }_k\text{or}{u}_i<0,\\ \{0\},& \text{if}0<u_i<{\zeta }_k,\\ [0,1],& \text{if}{u}_i={\zeta }_k\text{or}{u}_i=0,\end{array}\text{for all }i\in \{1,\dots ,l\}.$$

$${\Theta }_k(x,w,y):={\Vert {\hat{M}}_k\left(x,w,y\right)\Vert }^2.$$

$$M_{k_j}\left[\begin{array}{c}{d}_x\\ d_w\\ d_y\end{array}\right]=-{\hat{M}}_k\left(x_{k_j},w_{k_j},y_{k_j}\right).$$

$${\Theta }_k\left(\bar{x},\bar{w},\bar{y}\right)\le \left(1-2\mu {\delta }^{m_j}\right){\Theta }_k\left(x_{k_j},w_{k_j},y_{k_j}\right)$$

$${\left(B_{\delta ,k_j}\right)}_{(i,i)}:=\{\begin{array}{ll}1,& \text{if}{\left({\beta }_k^{-1}{z}_k+x_{k_j}\right)}_i\in \left(a_{l_i},a_{u_i}\right),\\ 0,& \text{otherwise},\end{array},$$

Full text

1 Introduction

In this paper we are interested in the solution of the following optimization problem

[TABLE]

where $c\in\mathbb{R}^{n}$, $Q\in\mathbb{R}^{n\times n}$ is a positive semi-definite matrix, $C\in\mathbb{R}^{l\times n}$, $d\in\mathbb{R}^{l}$, $A\in\mathbb{R}^{m\times n}$ is a linear constraint matrix with $b\in\mathbb{R}^{m}$ a given right-hand side, $D\in\mathbb{R}^{n\times n}$ is a diagonal positive semi-definite (“weight”) matrix, and $\mathcal{K}$ is a closed convex set $\mathcal{K}\coloneqq\{x\in\mathbb{R}^{n}\colon x\in\left[a_{l},a_{u}\right]\}$ with $a_{l},\ a_{u}\in\mathbb{R}^{n}$, such that $(a_{l})_{i}\in\mathbb{R}\cup\{-\infty\},\ (a_{u})_{i}\in\mathbb{R}\cup\{+\infty\}$ for all $i=1,\ldots,n$. Additionally, $(\cdot)_{+}\equiv\max\{\cdot,0\}$, while $\delta_{\mathcal{K}}(\cdot)$ is an indicator function for the set $\mathcal{K}$, that is, $\delta_{\mathcal{K}}(x)=0$ if $x\in\mathcal{K}$ and $\delta_{\mathcal{K}}(x)=\infty$, otherwise.

We introduce an auxiliary variable $w\in\mathbb{R}^{l}$, and reformulate (1.1) in the following form:

[TABLE]

Remark 1**.**

Let us notice that the model in (P) can readily accommodate terms of the form of $\|Cx+d\|_{1}$ where $C\in\mathbb{R}^{l\times n}$, and $d\in\mathbb{R}^{l}$. Indeed, letting $c=-C^{\top}\mathds{1}_{l}$ and adding the constant term $-\mathds{1}_{l}^{\top}d$ in the objective of (P), we notice that

[TABLE]

where $\mathds{1}_{l}\coloneqq(1,\ldots,1)^{\top}\in\mathbb{R}^{l}$. Similarly, any piecewise-linear term of the form

[TABLE]

where $C_{1},\ C_{2}\in\mathbb{R}^{l\times n}$ and $d_{1},\ d_{2}\in\mathbb{R}^{l}$, can also be readily modeled. Indeed, setting $c=C_{2}^{\top}\mathds{1}$ and adding the term $d_{2}^{\top}\mathds{1}_{l}$ in the objective yields

[TABLE]

Finally, it is important to note that model (P) allows for multiple piecewise-linear terms of the form $\max\{Cx+d,0_{l}\}$, $\|Cx+d\|_{1}$ or $\max\left\{C_{1}x+d_{1},C_{2}x+d_{2}\right\}$, since we can always adjust $l$ to account for more than one terms. Hence, one can observe that (P) is quite general and can be used to model a plethora of very important problems that arise in practice.

In light of the discussion in Remark 1, it is easily observed that problem (P) can model a plethora of very important problems arising in several application domains spanning, among others, operational research, machine learning, data science, and engineering. More specifically, various lasso and fussed lasso instances (with applications to sparse approximations for classification and regression [19, 54], portfolio allocation [1], or medical diagnosis [23], among many others) can be readily modeled by (P). Additionally, various risk-minimization problems with linear random cost functions can be modeled by (P) (e.g. see [33, 44, 48]). Furthermore, even risk-minimization problems with nonlinear random cost functions, which are typically solved via Gauss-Newton schemes (e.g. see [10]), often require the solution of sub-problems like (P). Finally, continuous relaxations of integer programming problems with applications to operational research (e.g. [32]) often take the form of (P). Given the multitude of problems requiring easy access to (usually accurate) solutions of (P), the derivation of efficient, robust, and scalable solution methods is of paramount importance.

Problem (P) can be solved by various first- or second-order methods. In particular, using a standard reformulation, by introducing several auxiliary variables, (P) can be written as a convex quadratic programming (QP) one and efficiently solved by, among others, an interior point method (IPM; e.g. [19, 38]), an alternating direction method of multipliers (ADMM; e.g. [20]), or a proximal point method (e.g. [28, 43]). However, the reformulation of (P) into a convex QP is not expected to lead to scalable solution methods, since the dimension of the problem is significantly increased, and hence an already large-scale instance might be very difficult to solve in this way.

Alternatively, ADMM (or general splitting) schemes can be developed for the solution of (P) without the need of additional auxiliary variables (see Section 3). However, no first-order scheme would be able to consistently yield sufficiently accurate solutions (i.e. of 4-, 5- or 6-digit accuracy). If such a solution is sought, we have to employ a semismooth Newton method (SSN; e.g. [12, 24, 37]), or a combination of a proximal point method with an SSN scheme utilized for the solution of its associated sub-problems. SSN methods have been successfully applied to a plethora of problems, however, their success is heavily reliant on the properties of the problem at hand (e.g. the rank of the linear constraints, or the conditioning of the Hessian). On the other hand, the combination of the proximal point method with the SSN can circumvent these issues, since the associated nonsmooth sub-problems can be guaranteed to be well-defined and well-conditioned.

Various such solvers have been developed and analyzed in the literature. For example, the authors in [53] developed a dual augmented Lagragian scheme combined with an SSN method for the solution of semidefinite programming problems, and obtained very promising results. This scheme was then utilized for the solution of linear programming problems in [35], and for lasso-regularized problems in [34]. A similar primal-dual approach for $\ell_{1}$-regularized convex quadratic programming problems was developed and analyzed in our accompanying paper [39] and was shown to be especially efficient for the solution of elastic-net linear regression and $L^{1}$-regularized partial differential equation constrained optimization problems. In fact, the proposed active set method developed in this work is a direct extension of the method given in [39], altered in a specific way so that it can efficiently handle most piecewise-linear terms that appear in practice, via restricting its memory requirements.

Indeed, we showcase that each of the nonsmooth terms in the objective of (P) can be utilized for reducing the memory requirements of the proposed method. In particular, when computing the Clarke subdifferential ([13]) of an appropriate augmented Lagrangian penalty function, we can show that the Clarke derivatives of such piecewise-linear terms can act as projectors. As a result, we obtain an active-set scheme that reduces the sub-problems’ dimensions significantly, leading to better performance and reduced memory requirements. In particular, we observe that a thresholding operation (originating from the presence of $\ell_{1}$ terms in the objective) determines which of the variables $x$ are inactive, allowing us to throw away entire columns from matrices $A$ and $C$ when solving the associated sub-problems. Furthermore, the $\max\{\cdot,0\}$ terms in the objective determine which of the rows of $C$ are non-important, allowing us to further eliminate such rows.

We showcase the robustness and the efficiency of the resulting active-set scheme on various optimization problems arising in risk-averse portfolio selection, quantile regression, and binary classification via linear support vector machines. In each of these cases the proposed scheme is compared against the robust polynomially convergent regularized interior point method developed and analyzed in [38], as well as the well-known ADMM-based OSQP solver ([49]). We demonstrate the reduced memory requirements of the active set scheme (and hence its improved scalability), as compared to interior point and ADMM alternatives applied to QP reformulations, while showcasing its efficiency and robustness.

Structure of the paper

In Section 2 we briefly present the proposed inner-outer active-set method. In particular, in Section 2.1 we derive a proximal method of multipliers (outer scheme) for the solution of (P), assuming that we can find an $\epsilon$-optimal solution of the associated sub-problems. Then, in Section 2.2, we briefly present a semismooth Newton method (inner scheme) for finding approximate solutions to sub-problems arising from the proximal method of multipliers. Focusing on the structure of problem (P), by selecting appropriate Clarke derivatives, we show that the proposed inner-outer method is in fact an active-set scheme, the associated linear systems of which are well-conditioned and stable. In Section 2.3, we discuss an extension of the method for dealing with problems having arbitrary piecewise-linear terms in the objective (with applications to, among others, robust optimization) that are not currently covered within the paper.

Subsequently, in Section 3, we derive a proximal alternating direction method of multipliers for the approximate solution of (P) in order to obtain good initial estimates for the primal and dual variables of the problem. This can then be used to warm-start the proposed second-order algorithm. A good starting point for the algorithm could mean that only a small portion of the problem data matrices is used at each inner-outer iteration of the scheme, while the outer method is expected to achieve its local linear (and potentially superlinear) convergence rate in a small number of inner-outer iterations.

In Section 4, we showcase the efficiency and robustness of the approach on several important real-life applications arising in risk-averse portfolio optimization, statistical regression, and binary classification via linear support vector machines. Finally, we discuss our conclusions in Section 5.

Notation

Given a vector $x$ in $\mathbb{R}^{n}$, $\|x\|$ denotes its Euclidean norm. Given a closed set $\mathcal{K}\subset\mathbb{R}^{n}$, $\Pi_{\mathcal{K}}(x)$ denotes the Euclidean projection onto $\mathcal{K}$, that is $\Pi_{\mathcal{K}}(x)\coloneqq\arg\min\{\|x-z\|\colon z\in\mathcal{K}\}$. Given a closed set $\mathcal{K}$, we write $\textnormal{dist}(z,D)\coloneqq\inf_{z^{\prime}\in\mathcal{K}}\|z-z^{\prime}\|$. Given a convex function $p\colon\mathbb{R}^{n}\mapsto\mathbb{R}$, we define the proximity operator as $\textbf{prox}_{p}(u)\coloneqq\arg\min_{x}\left\{p(x)+\frac{1}{2}\|u-x\|^{2}\right\}$. Given an index set $\mathcal{D}$, $|\mathcal{D}|$ denotes its cardinality. Given a rectangular matrix $A\in\mathbb{R}^{m\times n}$ and an index set $\mathcal{B}\subseteq\{1,\ldots,n\}$, we denote the columns of $A$, the indices of which belong to $\mathcal{B}$, as $A_{\mathcal{B}}$. Given a matrix $A\in\mathbb{R}^{m\times n}$, and two index sets $\mathcal{B}_{1}\subseteq\{1,\ldots,m\}$ and $\mathcal{B}_{2}\subseteq\{1,\ldots,n\}$, we denote the subset of rows and columns of $A$, the indices of which belong to $\mathcal{B}_{1},\ \mathcal{B}_{2}$ respectively, as $A_{\left(\mathcal{B}_{1},\mathcal{B}_{2}\right)}$. Finally, given an arbitrary vector $d$ with $n$ components as well as some indexes $1\leq i_{1}\leq i_{2}\leq n$, we denote by $d_{i_{1}:i_{2}}$ the vector containing the $i_{1}$-st up to the $i_{2}$-nd component of this vector. To avoid confusion, indexes $i$ are always used to denote entries of vectors or matrices, while $k$ and $j$ are reserved to denote iteration counters (outer and inner, respectively).

2 An active-set method

In this section we derive an active-set method for the solution of (P). The algorithm is an inner-outer scheme, which results by combining an outer proximal method of multipliers, and an inner semismooth Newton method for the solution of the PMM sub-problems. Following the discussion in [39, Section 2], we briefly derive the outer PMM scheme. Then, we briefly present the inner semismooth Newton method and discuss the solution of its associated linear systems, which is the main bottleneck of the algorithm.

2.1 Outer scheme: A primal-dual proximal method of multipliers

In this section, following and extending the developments in [39], we derive a primal-dual proximal method of multipliers for the approximate solution of (P). A convergence analysis is not provided, and the reader is referred to [39, Section 2], for an outline of such an analysis (since it applies directly to the case under consideration).

Given a penalty parameter $\beta>0$, and some dual multiplier estimates $y,\ z$, we follow [39, Section 2] to obtain the augmented Lagrangian corresponding to (P), which reads

[TABLE]

Indeed, this can be verified utilizing certain standard properties of Fenchel duality and of the proximity operator (see [39, Section 2]).

During iteration $k\geq 0$ of the proximal method of multipliers, we have the estimates $(x_{k},y_{k},z_{k})$ as well as the penalty parameters $\beta_{k},\ \rho_{k}$, such that $\rho_{k}\coloneqq\frac{\beta_{k}}{\tau_{k}}$, where $\{\tau_{k}\}_{k=0}^{\infty}$ is a non-increasing positive sequence. For simplicity of exposition, let $g\left(x,w\right)\coloneqq g_{1}(x)+g_{2}(w)$, where $g_{1}(x)\coloneqq\|Dx\|_{1},$ and $g_{2}(w)\coloneqq\sum_{i=1}^{l}\left(w_{i}\right)_{+}.$ We consider the following regularized continuously differentiable function:

[TABLE]

Notice that we introduce primal proximal regularizer only for variable $x$, and not for $w$. This is a very important algorithmic choice that departs from the developments in [39]. We want to develop a memory-efficient active-set method, and for that reason we avoid introducing a proximal term for the auxiliary variables $w$. This choice, which does not hinder the stability of the proposed approach, leads to better memory efficiency of the resulting active-set scheme, by promoting a sparsification of the associated linear systems solved at each inner-outer iteration (this point will become clear in Section 2.2.1).

The minimization of the proximal augmented Lagrangian function can be written as

[TABLE]

and thus we need to find $(x^{*},w^{*})\in\mathbb{R}^{n}\times\mathbb{R}^{l}$ such that

[TABLE]

To that end, we observe that

[TABLE]

By introducing the dual variable $y\coloneqq y_{k}-\beta_{k}\left(\begin{bmatrix}Cx+d-w\\ Ax-b\end{bmatrix}\right)\in\mathbb{R}^{l+m}$, the optimality conditions of $\underset{x,w}{\min}\ \left\{\phi(x,w)+g(x,w)\right\}$ can be written as

[TABLE]

where

[TABLE]

[TABLE]

We now describe the proposed proximal method of multipliers in Algorithm PMM.

Let us notice that given a triple $(\tilde{x},\tilde{w},\tilde{y})$, we can easily evaluate the distance in (2.3), due to the piecewise linear structure of the associated function $g(\cdot)$. A trivial extension of the analysis in [39, Section 3] yields that Algorithm PMM is globally convergent under the standard assumption of primal-dual feasibility. Furthermore, given some additional conditions on the error sequence $\{\epsilon_{k}\}$ one can show that in fact Algorithm PMM achieves a local linear convergence rate (which becomes superlinear if $\beta_{k}\rightarrow\infty$ at a suitable rate). Finally, the algorithm exhibits a global linear convergence rate, assuming that the starting point is chosen properly. For more details, we refer the reader to [39, Theorems 2.2, 2.3].

2.2 Inner scheme: A semismooth Newton method

Next, we employ a semismooth Newton (SSN) method to solve problem (2.3) appearing in Algorithm PMM, and verify that the resulting inner-outer scheme admits an active-set interpretation. Given the estimates $(x_{k},y_{k},z_{k})$ as well as the penalties $\beta_{k},\ \rho_{k}$, we apply SSN to approximately solve (2.2). Given any bounded positive penalty $\zeta_{k}>0$, the optimality conditions in (2.2) can be written as

[TABLE]

We set $x_{k_{0}}=x_{k}$, $y_{k_{0}}=y_{k}$, and at every iteration $k_{j}$ of SSN, we solve

[TABLE]

where we have introduced the notation

[TABLE]

assuming that we are given some arbitrary matrices

[TABLE]

The symbol $\partial_{x}^{C}(\cdot)$ denotes the Clarke subdifferential of a function (see [13]). In our case, any element of the Clarke subdifferential is a Newton derivative (see [14, Chapter 13]), since $r_{k}(\cdot,\cdot)$ and $g(\cdot,\cdot)\equiv g_{1}(\cdot)+g_{2}(\cdot)$ are piecewise continuously differentiable and regular functions. For any $u\in\mathbb{R}^{n}$ and any $i\in\{1,\ldots,n\}$, it holds that

[TABLE]

Since $g_{1}(x)=\|Dx\|_{1}$ and $D\succeq 0_{n}$ is diagonal, we have that for all $u\in\mathbb{R}^{n}$ and any $i\in\{1,\ldots,n\}$,

[TABLE]

Finally, since $g_{2}(w)=\max\{w,0\}$, we have that for all $u\in\mathbb{R}^{l}$ and any $i=1,\ldots,l$,

[TABLE]

Then, we can show (e.g. see [14, Example 14.9]) that

[TABLE]

and

[TABLE]

We complete the derivation of the SSN, by defining a primal-dual merit function for globalizing the semismooth Newton method via a backtracking line-search scheme. Following [39], we employ the following merit function

[TABLE]

This function is often used for globalizing SSN schemes applied to nonsmooth equations of the form of (2.6). In Algorithm SSN, we outline a locally superlinearly convergent semismooth Newton method for the approximate solution of (2.3). The associated linear systems can be solved approximately (e.g. by Krylov subspace methods as in [39]), although in this work a suitable factorization scheme is utilized.

Local superlinear convergence of Algorithm SSN follows directly from [37, Thereom 3], since the equation in (2.6) is BD-regular (i.e. the Bouligand subdifferential at the optimum contains nonsingular matrices) by construction, as it models “regularized” sub-problems arising from the proximal method of multipliers. On the other hand, if we assume that the directional derivative of (2.9) is continuous at the optimum, then we can mirror the analysis in [25] to show that Algorithm SSN is globally convergent. There is a wide literature on similar semismooth Newton schemes, and we refer the reader to [12, 25, 37, 40] and the references therein for more details.

2.2.1 The SSN linear systems

The major bottleneck of the previously presented inner-outer scheme is the solution of the associated linear systems given in (2.7). One could alter Algorithm SSN so that it does not require an exact solution. In turn this would allow for the utilization of preconditioned Krylov subspace solvers for the efficient solution of such systems (e.g. as in [39, Section 3]). In particular, any preconditioner derived in [21] can be utilized for the proposed solver. However, for simplicity of exposition, we employ a standard factorization approach. As will become evident in Section 4, the active-set nature of the proposed algorithm enables the use of factorization even for large-scale problems of interest, since most of the data are inactive at each inner-outer iteration $k_{j}$. Indeed, in the presence of multiple piecewise-linear terms in the objective, and assuming that the method is appropriately warm-started, one can ensure that most of the problem data are not incorporated when forming the active Hessian at each inner-outer iteration.

In what follows we derive the associated linear systems. Let $k\geq 0$ be an arbitrary iteration of Algorithm PMM, and $j\geq 0$ an arbitrary iteration of Algorithm SSN. Firstly, let us notice that any element $B_{\delta}\in\partial_{x}^{C}\left(\Pi_{\mathcal{K}}\left(\cdot\right)\right)$ yields a Newton derivative (see [14, Theorem 14.8]). The same applies for any element $B_{g_{1}}\in\partial_{x}^{C}\left(\textbf{prox}_{\zeta_{k}g_{1}}\left(\cdot\right)\right)$, and $B_{g_{2}}\in\partial_{w}^{C}\left(\textbf{prox}_{\zeta_{k}g_{2}}\left(\cdot\right)\right)$. Thus, using (2.8) we can choose $B_{\delta,k_{j}},\ B_{g_{1},k_{j}},\ B_{g_{2},k_{j}}$ from the Bouligand subdifferential to improve computational efficiency, by reducing the active variables and constraint rows. To that end, we set $B_{\delta,k_{j}}$ as a diagonal matrix with

[TABLE]

for all $i\in\{1,\ldots,n\}$, $B_{g_{1},k_{j}}$ as the following diagonal matrix

[TABLE]

for all $i\in\{1,\ldots,n\}$, where $\widehat{u}_{k_{j}}\coloneqq x_{k_{j}}-\zeta_{k}(r_{k}(x_{k_{j}},y_{k_{j}}))_{1:n}$, and $B_{g_{2},k_{j}}$ as

[TABLE]

for all $i\in\{1,\ldots,l\}$.

Given the aforementioned choices for the projection matrices, we can now explicitly eliminate certain variables from the system in (2.7), in order to obtain a saddle-point system. To that end, from the second block-equation in (2.7), we have

[TABLE]

Let $\mathcal{B}_{g_{2},k_{j}}\coloneqq\left\{i\in\{1,\ldots,l\}\colon\left(B_{g_{2},k_{j}}\right)_{(i,i)}=1\right\}$. Then, we obtain

[TABLE]

where the right-hand side can be evaluated easily. Letting $\mathcal{N}_{g_{2},k_{j}}\coloneqq\{1,\ldots,l\}\setminus\mathcal{B}_{g_{2},k_{j}}$, we have

[TABLE]

On the other hand, from the third block-equation of (2.7), we observe that

[TABLE]

which can be computed after solving the reduced system. We define the following index sets

[TABLE]

From the first block equation of (2.7), we obtain

[TABLE]

After a straightforward (if tedious) calculation, we can pivot $(d_{y})_{\mathcal{B}_{g_{2},k_{j}}}$, $d_{w}$, and $(d_{x})_{\mathcal{N}_{g_{1},k_{j}}}$ in system (2.7). This results in the following saddle-point system:

[TABLE]

where $\widehat{C}\coloneqq C_{\left({\mathcal{N}_{g_{2},k_{j}}},\mathcal{B}_{g_{1},k_{j}}\right)},\ \widehat{A}\coloneqq A_{\mathcal{B}_{g_{1},k_{j}}}$, and

[TABLE]

Notice that the coefficient matrix of (2.15) is symmetric and quasi-definite (see [50]). As such, it is strongly factorizable, that is, each symmetric permutation of it admits an $L\Delta L^{\top}$ decomposition, with $\Delta$ diagonal. Alternatively, one could eliminate further variables in order to obtain a positive definite linear system. This could be beneficial in certain applications, but is omitted here for the sake of generality. We note that $\widehat{C}$ contains only a subset of the columns and rows of $C$. Similarly, $\widehat{A}$ contains only a subset of the columns of $A$. As we will verify in practice, relatively close to the solution of (P) the active-set matrices $\widehat{C}$, $\widehat{A}$, can be significantly smaller than $C$ and $A$, respectively, allowing the proposed approach to solve large-scale instances, without requiring excessive memory. Finally, we note that if we had included a proximal term for the auxiliary variables $w$ (i.e. if we were to directly apply the algorithm given in [39]), then $\widehat{C}$ would only contain a subset of the columns of $C$ but all of its rows (which, in general, are expected to be much more problematic, since it is often the case that $l\gg n$).

2.3 Extension to arbitrary piecewise-linear terms

Before closing this section we would like to notice that there are certain piecewise-linear terms that are not captured by the model in (P). In particular, given a set of $K$ pairs $(C_{r},d_{r})\in\mathbb{R}^{l\times n}\times\mathbb{R}^{l}$ (with $K\geq 3$), where $r\in\{1,\ldots,K\}$, we could consider problems of the following form

[TABLE]

In order to do so, we would have to introduce $K$ auxiliary vectors $w_{r}\in\mathbb{R}^{l}$, $r\in\{1,\ldots,K\}$, and reformulate the problem in the following form

[TABLE]

Subsequently, in the derivation of the semismooth Newton method, we would have to evaluate the proximity operator of $\tilde{g}(u)\coloneqq\max_{i\in\{1,\ldots,K\}}(u_{i})$, $u\in\mathbb{R}^{K}$, which admits a closed-form solution:

[TABLE]

Then, the Clarke derivative of the previous could be computed by utilizing [14, Theorem 14.7]. This was not considered in this work for brevity of presentation, and is left as a future research direction. Indeed, this extension paves the way for the generalization of the proposed active-set scheme to a plethora of robust optimization problems that appear in practice (e.g. see [4]), as well as delivering an alternative to standard cutting-plane methods (e.g. see [29]), decomposition methods (e.g. see [3, 5, 18, 46]), or specialized interior point (e.g. see [22]) approaches appearing in the literature.

3 Warm-starting

Following the developments in [35, 39], we would like to find a starting point for Algorithm PMM that is relatively close to a primal-dual solution. Indeed, this is crucial since, on the one hand, then we can expect to observe early linear convergence of Algorithm PMM, while on the other hand, we can expect to be close to identifying the correct active-sets which in turn implies that the memory requirements of Algorithm SSN are significantly reduced. To that end, we employ a proximal alternating direction method of multipliers (pADMM; e.g. see [20]). We reformulate (P) by introducing an artificial variable $u\in\mathbb{R}^{n+l}$, as

[TABLE]

Given a penalty $\sigma>0$, we associate the following augmented Lagrangian to (P’)

[TABLE]

where $y\in\mathbb{R}^{m+n+2l}$ denotes the dual multipliers associated to the linear equality constraints of (P’). Let an arbitrary positive definite matrix $R\in\mathbb{R}^{(n+l)\times(n+l)}$ be given, and assume the notation $\|v\|^{2}_{R}=v^{\top}Rv$, for any $v\in\mathbb{R}^{n+l}$. Also, as in Section 2.1, we denote $g(x,w)=\|Dx\|_{1}+\sum_{i=1}^{l}\left(w_{i}\right)_{+}$. Algorithm pADMM describes a proximal ADMM for the approximate solution of (P’).

Let us notice that under certain standard assumptions on (P’), Algorithm pADMM converges globally, potentially at a linear rate (see [20]). We can employ the regularization matrix $R$ as a means of ensuring memory efficiency of Algorithm pADMM. For example, we can recover the prox-linear ADMM [20, Section 1.1], where given some sufficiently large constant $\hat{\sigma}>0$, one defines

[TABLE]

Given this choice of $R$, the second step of Algorithm pADMM consists of only matrix-vector multiplications with $A,\ C$ and $Q$, and thus no (non-diagonal) matrix inversion is required. If memory is not an issue, one could use a positive-definite diagonal regularizer $R$, yielding a standard regularized ADMM. In our implementation, the user can specify which of the two strategies should be employed. In either case, the first step of Algorithm pADMM is trivial to solve, since we know that the proximity operator of $g(\cdot)$ can be computed analytically as in Section 2.2.

Finally, once an approximate solution $(\tilde{x},\tilde{w},\tilde{u},\tilde{y})$ is computed, we set the starting point of Algorithm PMM as $(x_{0},w_{0},y_{0},z_{0})=\left(\tilde{x},\tilde{w},\left(\tilde{y}\right)_{(1:m+l)},\tilde{z}\right)$, where

[TABLE]

Indeed, an optimal primal-dual solution of (P’) satisfies

[TABLE]

thus the characterization of $\tilde{z}$ in Algorithm PMM follows from the Appendix, where one can also find the termination criteria of Algorithm pADMM (see (A.2)).

4 Applications and numerical results

In this section we present various applications that can be modeled by problem (P), focusing on portfolio optimization, quantile regression, and binary classification. In particular, we first discuss (and numerically demonstrate) the effectiveness of the approach for the solution of single-period mean-risk portfolio optimization problems, where risk is measured via the conditional value at risk or the mean absolute semi-deviation. Subsequently, we apply the proposed scheme for the solution of quantile regression problems, demonstrating its scalability. Finally, we apply the active-set scheme for the solution of binary classification problems via linear support vector machines on certain large-scale datasets.

Implementation details

Before proceeding to the numerical results, we mention certain implementation details of the proposed algorithm. The implementation follows closely the developments in [39], and can be found on GitHub111https://github.com/spougkakiotis/Active_set_method_for_CQP_piecewise_LP. The code is written in MATLAB, and the experiments are run on a PC with a 2.2GHz Intel core i7-8750H processor (hexa-core), 16GB of RAM, using the Windows 10 operating system.

We run Algorithm pADMM (warm-start) for at most 100 iterations (or until a 3-digit accurate solution is found). The user can choose whether the warm-starting scheme should be matrix-free or not. In the presented experiments, the matrix-free scheme was only used for the largest quantile regression and classification problems, for which standard ADMM crashed due to excessive memory requirements. Then, starting with $\beta_{0}=10$, $\rho_{0}=50$, we run Algorithms PMM-SSN. Following [39], when solving the PMM sub-problems using Algorithm SSN we use a predictor-corrector-like heuristic in which the first iteration is accepted without line-search and then line-search is activated for subsequent iterations. Algorithm PMM is allowed to run for at most 200 iterations, while Algorithm SSN is stopped after 20 inner iterations. An iterate is accepted as optimal if the conditions given in (A.1) are satisfied for the tolerance specified by the user. Most other implementation details follow directly from the developments in Sections 2.1–2.2. We refer the reader to the implementation on GitHub for additional details.

In the presented experiments we compare the proposed approach against the interior point solver given in [38], the implementation of which can be found on GitHub222https://github.com/spougkakiotis/IP_PMM, as well as the ADMM-based OSQP solver given in [49], the implementation of which can also be found on GitHub333https://github.com/osqp/osqp-matlab. We note that most problems considered herein were especially challenging for OSQP (mostly due to requesting highly accurate solutions), and for that reason we allow it to run for 50,000 iterations (unlike its default iteration threshold, which is 4000).

4.1 Portfolio optimization

We first consider the mean-risk portfolio selection problem (originally proposed in [36]), where we minimize some convex risk measure, while keeping the expected return of the portfolio above some desirable level. A variety of models for this problem have been intensely analyzed and solved in the literature (e.g. see [1, 31, 44]). The departure from the variance as a measure of risk often allows for great flexibility in the decision making of investors, enabling them to follow potential regulations as well as to better control the risk associated with an “optimal” portfolio. Optimality conditions and existence of solutions for several general deviation measures have been characterized in [45]. The comparison of portfolios obtained by minimizing different risk measures has also been considered (e.g. see [41, 42, 48]).

The method presented in this paper is general and not related to a specific model choice. Indeed, we would like to showcase the efficiency of our approach for obtaining accurate solutions to portfolio optimization problems with various risk measures of practical interest. Thus, we focus on the solution of a standard portfolio selection model that has been used in the literature. All numerical results are obtained on real-world datasets. As a result the problems are of medium-scale. Nonetheless, even for such medium-scale problems, we will be able to demonstrate the efficiency of the proposed approach, when compared to the efficient interior point method employed in the literature [19] for similar problems. We also note that both second-order solvers (i.e. IPM and active-set) significantly outperform OSQP on these instances, however it was included in the comparison for completeness. Some large-scale instances will be tackled in the following subsections, where the method will be applied to quantile regression and binary classification instances.

Let $x\in\mathbb{R}^{n}$ represent a portfolio of $n$ financial instruments, such that

[TABLE]

This requirement indicates that short positions for each stock are restricted by some percentage ($a_{l_{i}}\%$) of the available wealth (assuming $a_{l_{i}}<0$), and no more than $a_{u_{i}}\%$ of the total wealth can be invested to instrument $i$. Let $\bm{\xi}\in\mathbb{R}^{n}$ denote a random vector, the $i$-th entry of which represents the random return of the $i$-th instrument. Then, the random loss (i.e. the negative of the random return) of a given portfolio $x$ is given by $f(x,\bm{\xi})\coloneqq-x^{\top}\bm{\xi}$. In this paper we assume that $\bm{\xi}$ follows some continuous distribution $p(\bm{\xi})$, as well as that there is a one-to-one correspondence between percentage return and monetary value (as in [44, Section 3]). Additionally, given some expected benchmark return $r$ (e.g. the market index), we only consider portfolios that yield an expected return above a certain threshold, i.e. $\mathbb{E}[-f(x,\bm{\xi})]\geq r$.

Finally, given the previously stated constraints, we would like to minimize some convex risk measure $\varrho(\cdot)$ of interest. However, in order to make sure that the problem is well-posed, while the transaction costs are not excessive, we include an $\ell_{1}$ term in the objective. This is a well-known modeling choice that yields sparse portfolios, and thus regulates the transaction costs in the single-period portfolio setting (e.g. see [1, 16, 31]). Additionally, with an appropriate tuning of the $\ell_{1}$ regularization parameter $\tau>0$, one could also control the amount of short positions (see [16]). It should be mentioned here that in the multi-period setting such an $\ell_{1}$ term does not guarantee a reduction in the transaction costs (but mostly in the holding costs), and an additional total variation term should also be added in the objective (see [15, 19]). This is omitted here, since we focus on the single-period case, but the model in (P) could easily incorporate such an additional term (see Remark 1). By putting everything together, the model reads as

[TABLE]

There are several methods for solving such stochastic problems; let us mention two important variants. There is the parametric approach (e.g. as in [1, 44]), where one assumes that the returns follow some known distribution which is subsequently sampled to yield finite-dimensional optimization problems, and the sampling approach (e.g. as in [31]), where one obtains a finite number of samples (without assuming a specific distribution). Such samples are often obtained by historical observations, and this approach is also followed in this paper. It is well-known that historical data cannot fully predict the future (see [31]), however it is a widely-used practice. The reader is referred to [27] for an extensive discussion onprobabilistic models for portfolio selection problems.

Additional soft or hard constraints can be included when solving a portfolio selection problem. Such constraints can either be incorporated directly via the use of a model (e.g. see [7, Section 2]) as hard constraints or by including appropriate penalty terms in the objective (soft constraints). It is important to note that the model given in (P) is quite general and as a result has great expressive power, allowing one to incorporate various important risk measures (and their combinations), as well as modeling constraints of interest.

Real-world datasets:

In what follows, we solve two different instances of problem (4.1). In particular, we consider two potential risk measures; the conditional value at risk (e.g. see [44]), as well as the mean absolute semi-deviation (e.g. see [47, Section 6.2.2.], noting that this is in fact equivalent to the mean absolute deviation originally proposed in [30]). In each of the above cases problem (4.1) has the form of (P) and thus Algorithm PMM can directly be applied.

We showcase the effectiveness of the proposed approach on 6 real datasets taken from [9]. Each dataset contains time series for weekly asset returns and market indexes for different major stock markets, namely, DowJones, NASDAQ100, FTSE100, SP500, NASDAQComp, and FF49Industries. In the first 5 markets the authors in [9] provide the market indexes, while for the last dataset the uniform allocation strategy is considered as a benchmark. Additional information on the datasets is collected in Table 1. We note that stocks with less than 10 years of observations have been disregarded.

4.1.1 Conditional value at risk

First, we consider portfolio optimization problems that seek a solution minimizing the conditional value at risk; a measure which is known to be coherent (see [2] for a definition of coherent risk measures). In particular using the notation introduced earlier, we consider the following optimization problem

[TABLE]

where $f(x,\bm{\xi})$ is the random cost function, $A\in\mathbb{R}^{m\times n}$ models the linear constraint matrix of problem (4.1) (where an auxiliary variable has been introduced to transform the inequality constraint involving the expectation into an equality), and $\mathcal{K}\coloneqq[a_{l},a_{u}]$. In the above $1-\alpha\in(0,1)$ is the confidence level. It is well-known ([44]) that given a continuous random variable $X$, the conditional value at risk can be computed as

[TABLE]

We can write problem (4.2) in the following equivalent form:

[TABLE]

where the expectation has been substituted by summation since we assume the availability of a dataset $\{\xi_{1},\ldots,\xi_{l}\}$. Introducing an auxiliary variable $w\in\mathbb{R}^{l}$, the previous can be re-written as

[TABLE]

We solve problem (4.2) using the proposed active-set method (AS), the interior point-proximal method of multipliers (IP-PMM) given in [38], and OSQP ([49]). We allow any short position, that is $a_{l_{i}}=-1$, and restrict investing more that $60\%$ of the available wealth on a single stock (i.e. $a_{u_{i}}=0.6$). We set $\texttt{tol}=10^{-5}$ and $\tau=10^{-2}$, and run the three methods for each of the datasets described in Table 1 for varying confidence level. We report the confidence parameter $\alpha$, the number of PMM, SSN, IP-PMM and OSQP iterations, the CPU time needed by each of the three schemes, as well as the number of factorizations used within the active set (PMM-SSN) scheme. Indeed, it often happens that the active-set is not altered from one iteration to the next, and the factorization does not need to be re-computed. The results are collected in Table 4.1.1. In all the numerical results that follow, the lowest running time exhibited by a solver, assuming it successfully converged, is presented in bold.