An active-set algorithm for norm constrained quadratic problems

TL;DR

This paper introduces an active-set algorithm for solving nonconvex quadratic problems with norm constraints, leveraging efficient trust-region subproblem solutions and eigenpair computations, demonstrating superior scalability and performance.

Contribution

It extends trust-region subproblem solutions to nonconvex quadratic problems with norm constraints, including efficient computation of nonglobal minimizers.

Findings

Algorithm outperforms state-of-the-art methods in experiments.

Effective in high-dimensional problems like Sparse PCA with 70 million nonzeros.

Scalable approach applicable to large-scale optimization tasks.

Abstract

We present an algorithm for the minimization of a nonconvex quadratic function subject to linear inequality constraints and a two-sided bound on the 2-norm of its solution. The algorithm minimizes the objective using an active-set method by solving a series of Trust-Region Subproblems (TRS). Underpinning the efficiency of this approach is that the global solution of the TRS has been widely studied in the literature, resulting in remarkably efficient algorithms and software. We extend these results by proving that nonglobal minimizers of the TRS, or a certificate of their absence, can also be calculated efficiently by computing the two rightmost eigenpairs of an eigenproblem. We demonstrate the usefulness and scalability of the algorithm in a series of experiments that often outperform state-of-the-art approaches; these include calculation of high-quality search directions arising in Sequential Quadratic Programming on problems of the CUTEst collection, and Sparse Principal Component Analysis on a large text corpus problem (70 million nonzeros) that can help organize documents in a user interpretable way.

Tables

Table 1. Table 1: Results of Richtarik2012 on applying sparse PCA on the NYtimes dataset. Each column lists words that correspond to nonzero entries in the respective principal vectors. The resulting principal vectors have clear associated meanings that correspond to sports, economics, elections, education, and the US.

Vector #1	Vector #2	Vector #3	Vector #4	Vector #5
game	companies	campaign	children	attack
play	company	president	program	government
player	million	zzz_al_gore	school	official
season	percent	zzz_bush	student	zzz_u_s
team	stock	zzz_george_bush	teacher	zzz_united_states

Table 2. Table 2: Results of applying sparse PCA on the NYtimes dataset with TPower , allowing weights with mixed signs. Each column lists words that correspond to nonzero entries in the generated principal vector. Note that increasing the cardinality of the sparse principal vector has severely affected the interpretability of the principal vector as compared to the results of Table 1 that have clear associated meaning.

Words	Weights	Words (continued)	Weights (continued)
teamed	$-\mathrm{0.451\hspace{0.17em}97}$	hitch	$-\mathrm{0.083\hspace{0.17em}649}$
gameday	$-0.3698$	playland	$-\mathrm{0.081\hspace{0.17em}834}$
seasonable	$-\mathrm{0.346\hspace{0.17em}98}$	leaguer	$-\mathrm{0.076\hspace{0.17em}07}$
playful	$-\mathrm{0.281\hspace{0.17em}05}$	nightcap	$-\mathrm{0.072\hspace{0.17em}675}$
playa	$-\mathrm{0.243\hspace{0.17em}74}$	percentage	$\mathrm{0.272\hspace{0.17em}63}$
gamesmanship	$-\mathrm{0.209\hspace{0.17em}57}$	companywide	$\mathrm{0.234\hspace{0.17em}43}$
coachable	$-\mathrm{0.139\hspace{0.17em}91}$	companion	$\mathrm{0.145\hspace{0.17em}81}$
pointe	$-\mathrm{0.133\hspace{0.17em}15}$	stockade	$\mathrm{0.127\hspace{0.17em}27}$
runaround	$-\mathrm{0.112\hspace{0.17em}66}$	marketability	$\mathrm{0.123\hspace{0.17em}68}$
win95	$-0.1104$	millionaire	$\mathrm{0.110\hspace{0.17em}05}$
yardage	$-\mathrm{0.105\hspace{0.17em}03}$	governmental	$\mathrm{0.092\hspace{0.17em}267}$
guzzle	$-\mathrm{0.093\hspace{0.17em}052}$	billionaire	$\mathrm{0.090\hspace{0.17em}433}$
player	$-\mathrm{0.089\hspace{0.17em}249}$	businesses	$\mathrm{0.076\hspace{0.17em}009}$
wonder	$-\mathrm{0.086\hspace{0.17em}247}$	zzz_bush_administration	$\mathrm{0.072\hspace{0.17em}731}$
ballad	$-\mathrm{0.085\hspace{0.17em}514}$	fundacion	$\mathrm{0.070\hspace{0.17em}655}$

Table 3. Table 3: Results of applying non-negative sparse PCA on the NYtimes dataset with the suggested active-set algorithm. Each column lists words that correspond to nonzero entries in the respective principal vectors. The words are sorted in decreasing magnitude of the respective weights. Given a set of computed principal vectors { x i } subscript 𝑥 𝑖 \{x_{i}\} , the next principal vector is computed on the deflated data matrix D − ∑ ( D ​ x i ) ​ x i T 𝐷 𝐷 subscript 𝑥 𝑖 superscript subscript 𝑥 𝑖 𝑇 D-\sum(Dx_{i})x_{i}^{T} . The algorithm is initialized with the 30 30 30 most positive, or most negative, entries (depending on which one leads to higher explained variance) of the first singular vector of the respective matrix. Furthermore we perform a standard post-processing polishing step at each resulting principal vector, as described in ( Journee2010, , §4.2) . Note that, unlike Table 2 our algorithm results in principal vectors that have clear, user interpretable associated meanings that corresponds to sports, economics, politics, education, and US foreign policy. These words could be used for document classification Zhang2011 .

Vector #1	Vector #2	Vector #3	Vector #4	Vector #5
team	percent	zzz_al_gore	official	school
game	company	zzz_bush	government	student
season	million	zzz_george_bush	attack	children
player	companies	campaign	zzz_u_s	program
play	stock	president	zzz_united_states	teacher
games	market	election	zzz_bush	parent
point	billion	political	military	high
coach	business	zzz_white_house	zzz_american	percent
run	fund	republican	palestinian	public
win	analyst	voter	zzz_taliban	education
yard	money	vote	zzz_afghanistan	kid
guy	zzz_enron	presidential	terrorist	district
won	firm	democratic	war	test
played	investor	tax	leader	family
hit	sales	zzz_republican	zzz_israel	college
ball	industry	administration	group	home
playing	plan	zzz_clinton	bin	child
league	cost	ballot	laden	class
left	investment	votes	country	group
home	quarter	plan	administration	job
fan	deal	democrat	security	boy
shot	financial	zzz_washington	forces	money
field	customer	zzz_congress	zzz_israeli	help
playoff	economy	poll	zzz_pakistan	mother
night	chief	support	american	system
goal	price	candidate	zzz_washington	friend
final	executive	zzz_florida	nation	girl
start	growth	candidates	terrorism	private
quarterback	earning	governor	weapon	voucher
football	share	vice	foreign	grade

Table 4. Table 4: Comparison of the suggested active-set algorithm against Ipopt in calculating SQP search directions on 63 nonconvex problems of the CUTEst library. These include all the problems with linear inequality constraints, at most n = 2000 𝑛 2000 n=2000 variables and an indefinite Hessian at the initial point. The initial point is computed with GUROBI as the nearest feasible point (in the two norm) to the initial point suggested by CUTEst . The column m 𝑚 m denotes the number of inequality constraints, while columns t , δ ​ f ∗ , error 𝑡 𝛿 superscript 𝑓 error t,\delta f^{*},\text{error} denote the solution time (in seconds), the difference in objective value f ​ ( δ ​ x ∗ ) − f ​ ( 0 ) 𝑓 𝛿 superscript 𝑥 𝑓 0 f(\delta x^{*})-f(0) achieved and the residual error defined by ( 30 ) of the suggested active-set solver. Similarly, δ ​ f ipopt ∗ , error ipopt 𝛿 subscript superscript 𝑓 ipopt subscript error ipopt \delta f^{*}_{\text{ipopt}},\text{error}_{\text{ipopt}} denote respectively identical metrics for Ipopt . t ipopt sparse superscript subscript 𝑡 ipopt sparse t_{\text{ipopt}}^{\text{sparse}} ( t ipopt dense ) superscript subscript 𝑡 ipopt dense (t_{\text{ipopt}}^{\text{dense}}) denotes Ipopt ’s timings when the matrices are passed in sparse (dense) format. Hardware used: Intel E5-2640v3, 64GB memory.

Problem	Parameter	$n$	$m$	$t$	$t_{\text{ipopt}}^{\text{dense}}$	$t_{\text{ipopt}}^{\text{sparse}}$	$\delta f^{\ast }$	$\delta f_{\text{ipopt}}^{\ast }$	error	${\text{error}}_{\text{ipopt}}$
AVION2		$49$	$113$	$\mathrm{4.444\hspace{0.17em}538}\text{×}{10}^{-02}$	$\mathrm{4.670\hspace{0.17em}914}\text{×}{10}^{-02}$	$\mathrm{1.723\hspace{0.17em}294}\text{×}{10}^{-02}$	$-\mathrm{3.645\hspace{0.17em}564}\text{×}{10}^{02}$	$-\mathrm{3.645\hspace{0.17em}789}\text{×}{10}^{02}$	$\mathrm{1.479\hspace{0.17em}577}\text{×}{10}^{-10}$	$\mathrm{5.426\hspace{0.17em}607}\text{×}{10}^{-04}$
BLOCKQP1	N=100	$210$	$521$	$\mathrm{7.184\hspace{0.17em}425}\text{×}{10}^{-03}$	$\mathrm{1.308\hspace{0.17em}548}\text{×}{10}^{00}$	$\mathrm{5.072\hspace{0.17em}931}\text{×}{10}^{-02}$	$-\mathrm{6.745\hspace{0.17em}688}\text{×}{10}^{00}$	$-\mathrm{6.745\hspace{0.17em}688}\text{×}{10}^{00}$	$\mathrm{9.992\hspace{0.17em}007}\text{×}{10}^{-16}$	$\mathrm{9.472\hspace{0.17em}623}\text{×}{10}^{-09}$
BLOCKQP2	N=100	$210$	$521$	$\mathrm{1.333\hspace{0.17em}553}\text{×}{10}^{-02}$	$\mathrm{1.234\hspace{0.17em}218}\text{×}{10}^{00}$	$\mathrm{5.287\hspace{0.17em}435}\text{×}{10}^{-02}$	$-\mathrm{6.745\hspace{0.17em}688}\text{×}{10}^{00}$	$-\mathrm{6.745\hspace{0.17em}688}\text{×}{10}^{00}$	$\mathrm{6.661\hspace{0.17em}338}\text{×}{10}^{-16}$	$\mathrm{9.374\hspace{0.17em}904}\text{×}{10}^{-09}$
BLOCKQP3	N=100	$210$	$521$	$\mathrm{1.551\hspace{0.17em}006}\text{×}{10}^{-02}$	$\mathrm{7.109\hspace{0.17em}685}\text{×}{10}^{-01}$	$\mathrm{3.531\hspace{0.17em}169}\text{×}{10}^{-02}$	$-\mathrm{4.018\hspace{0.17em}526}\text{×}{10}^{00}$	$-\mathrm{4.018\hspace{0.17em}526}\text{×}{10}^{00}$	$\mathrm{2.164\hspace{0.17em}935}\text{×}{10}^{-15}$	$\mathrm{9.224\hspace{0.17em}197}\text{×}{10}^{-09}$
BLOCKQP4	N=100	$210$	$521$	$\mathrm{9.190\hspace{0.17em}950}\text{×}{10}^{-03}$	$\mathrm{7.812\hspace{0.17em}505}\text{×}{10}^{-01}$	$\mathrm{4.003\hspace{0.17em}907}\text{×}{10}^{-02}$	$-\mathrm{4.018\hspace{0.17em}526}\text{×}{10}^{00}$	$-\mathrm{4.018\hspace{0.17em}526}\text{×}{10}^{00}$	$\mathrm{1.015\hspace{0.17em}854}\text{×}{10}^{-14}$	$\mathrm{8.997\hspace{0.17em}483}\text{×}{10}^{-09}$
BLOCKQP5	N=100	$210$	$521$	$\mathrm{9.125\hspace{0.17em}178}\text{×}{10}^{-03}$	$\mathrm{7.250\hspace{0.17em}970}\text{×}{10}^{-01}$	$\mathrm{3.328\hspace{0.17em}339}\text{×}{10}^{-02}$	$-\mathrm{4.018\hspace{0.17em}526}\text{×}{10}^{00}$	$-\mathrm{4.018\hspace{0.17em}526}\text{×}{10}^{00}$	$\mathrm{1.160\hspace{0.17em}183}\text{×}{10}^{-14}$	$\mathrm{9.325\hspace{0.17em}761}\text{×}{10}^{-09}$
BLOWEYA	N=100	$202$	$304$	$\mathrm{2.057\hspace{0.17em}082}\text{×}{10}^{-01}$	$\mathrm{6.153\hspace{0.17em}741}\text{×}{10}^{-01}$	$\mathrm{4.804\hspace{0.17em}298}\text{×}{10}^{-02}$	$-\mathrm{2.902\hspace{0.17em}540}\text{×}{10}^{-03}$	$-\mathrm{2.903\hspace{0.17em}560}\text{×}{10}^{-03}$	$\mathrm{1.103\hspace{0.17em}075}\text{×}{10}^{-12}$	$\mathrm{4.457\hspace{0.17em}334}\text{×}{10}^{-05}$
BLOWEYB	N=100	$202$	$304$	$\mathrm{3.518\hspace{0.17em}280}\text{×}{10}^{-01}$	$\mathrm{5.780\hspace{0.17em}202}\text{×}{10}^{-01}$	$\mathrm{4.473\hspace{0.17em}963}\text{×}{10}^{-02}$	$-\mathrm{2.574\hspace{0.17em}329}\text{×}{10}^{-03}$	$-\mathrm{2.575\hspace{0.17em}152}\text{×}{10}^{-03}$	$\mathrm{1.452\hspace{0.17em}606}\text{×}{10}^{-12}$	$\mathrm{4.834\hspace{0.17em}189}\text{×}{10}^{-05}$
BLOWEYC	N=100	$202$	$304$	$\mathrm{3.704\hspace{0.17em}172}\text{×}{10}^{-02}$	$\mathrm{5.043\hspace{0.17em}310}\text{×}{10}^{-01}$	$\mathrm{4.545\hspace{0.17em}440}\text{×}{10}^{-02}$	$-\mathrm{2.714\hspace{0.17em}783}\text{×}{10}^{-03}$	$-\mathrm{2.715\hspace{0.17em}942}\text{×}{10}^{-03}$	$\mathrm{3.039\hspace{0.17em}952}\text{×}{10}^{-15}$	$\mathrm{3.156\hspace{0.17em}699}\text{×}{10}^{-05}$
EQC		$9$	$21$	$\mathrm{2.374\hspace{0.17em}110}\text{×}{10}^{-04}$	$\mathrm{3.362\hspace{0.17em}073}\text{×}{10}^{-02}$	$\mathrm{3.607\hspace{0.17em}733}\text{×}{10}^{-02}$	$-\mathrm{1.653\hspace{0.17em}931}\text{×}{10}^{00}$	$-\mathrm{4.903\hspace{0.17em}961}\text{×}{10}^{03}$	$\mathrm{0.000\hspace{0.17em}000}\text{×}{10}^{00}$	$\mathrm{5.729\hspace{0.17em}870}\text{×}{10}^{00}$
EXPFITA		$5$	$22$	$\mathrm{9.989\hspace{0.17em}510}\text{×}{10}^{-04}$	$\mathrm{5.855\hspace{0.17em}086}\text{×}{10}^{-03}$	$\mathrm{7.309\hspace{0.17em}660}\text{×}{10}^{-03}$	$-\mathrm{1.529\hspace{0.17em}569}\text{×}{10}^{01}$	$-\mathrm{1.529\hspace{0.17em}569}\text{×}{10}^{01}$	$\mathrm{7.813\hspace{0.17em}195}\text{×}{10}^{-12}$	$\mathrm{1.159\hspace{0.17em}203}\text{×}{10}^{-07}$
EXPFITB		$5$	$102$	$\mathrm{8.978\hspace{0.17em}320}\text{×}{10}^{-04}$	$\mathrm{1.018\hspace{0.17em}311}\text{×}{10}^{-02}$	$\mathrm{1.349\hspace{0.17em}212}\text{×}{10}^{-02}$	$-\mathrm{7.615\hspace{0.17em}627}\text{×}{10}^{01}$	$-\mathrm{7.615\hspace{0.17em}627}\text{×}{10}^{01}$	$\mathrm{1.547\hspace{0.17em}296}\text{×}{10}^{-11}$	$\mathrm{2.935\hspace{0.17em}301}\text{×}{10}^{-06}$
EXPFITC		$5$	$502$	$\mathrm{2.641\hspace{0.17em}202}\text{×}{10}^{-03}$	$\mathrm{4.953\hspace{0.17em}891}\text{×}{10}^{-02}$	$\mathrm{6.193\hspace{0.17em}263}\text{×}{10}^{-02}$	$-\mathrm{3.810\hspace{0.17em}092}\text{×}{10}^{02}$	$-\mathrm{3.810\hspace{0.17em}092}\text{×}{10}^{02}$	$\mathrm{8.583\hspace{0.17em}134}\text{×}{10}^{-13}$	$\mathrm{1.278\hspace{0.17em}743}\text{×}{10}^{-04}$
FERRISDC	n=100	$1100$	$2110$	$\mathrm{1.632\hspace{0.17em}799}\text{×}{10}^{01}$	$\mathrm{4.871\hspace{0.17em}856}\text{×}{10}^{02}$	$\mathrm{2.109\hspace{0.17em}886}\text{×}{10}^{00}$	$-\mathrm{2.314\hspace{0.17em}228}\text{×}{10}^{-01}$	$-\mathrm{2.314\hspace{0.17em}213}\text{×}{10}^{-01}$	$\mathrm{9.975\hspace{0.17em}523}\text{×}{10}^{-09}$	$\mathrm{4.939\hspace{0.17em}870}\text{×}{10}^{-05}$
GOULDQP1		$32$	$81$	$\mathrm{4.174\hspace{0.17em}306}\text{×}{10}^{-03}$	$\mathrm{2.652\hspace{0.17em}730}\text{×}{10}^{-02}$	$\mathrm{1.418\hspace{0.17em}419}\text{×}{10}^{-02}$	$-\mathrm{8.790\hspace{0.17em}888}\text{×}{10}^{02}$	$-\mathrm{8.790\hspace{0.17em}888}\text{×}{10}^{02}$	$\mathrm{9.783\hspace{0.17em}285}\text{×}{10}^{-13}$	$\mathrm{7.851\hspace{0.17em}076}\text{×}{10}^{-05}$
HIMMELBJ		$45$	$61$	$\mathrm{4.420\hspace{0.17em}100}\text{×}{10}^{00}$	$\mathrm{3.329\hspace{0.17em}924}\text{×}{10}^{-02}$	$\mathrm{2.115\hspace{0.17em}145}\text{×}{10}^{-02}$	$-\mathrm{7.390\hspace{0.17em}332}\text{×}{10}^{01}$	$-\mathrm{7.391\hspace{0.17em}036}\text{×}{10}^{01}$	$\mathrm{3.586\hspace{0.17em}903}\text{×}{10}^{00}$	$\mathrm{1.144\hspace{0.17em}493}\text{×}{10}^{-05}$
HS105		$8$	$17$	$\mathrm{1.351\hspace{0.17em}519}\text{×}{10}^{-03}$	$\mathrm{4.024\hspace{0.17em}936}\text{×}{10}^{-03}$	$\mathrm{4.477\hspace{0.17em}112}\text{×}{10}^{-03}$	$-\mathrm{8.081\hspace{0.17em}024}\text{×}{10}^{01}$	$-\mathrm{8.081\hspace{0.17em}024}\text{×}{10}^{01}$	$\mathrm{1.452\hspace{0.17em}172}\text{×}{10}^{-12}$	$\mathrm{2.669\hspace{0.17em}379}\text{×}{10}^{-07}$
HS24		$2$	$5$	$\mathrm{3.546\hspace{0.17em}430}\text{×}{10}^{-04}$	$\mathrm{6.277\hspace{0.17em}818}\text{×}{10}^{-03}$	$\mathrm{6.520\hspace{0.17em}859}\text{×}{10}^{-03}$	$-\mathrm{1.311\hspace{0.17em}008}\text{×}{10}^{-01}$	$-\mathrm{1.311\hspace{0.17em}008}\text{×}{10}^{-01}$	$\mathrm{2.220\hspace{0.17em}446}\text{×}{10}^{-16}$	$\mathrm{1.623\hspace{0.17em}253}\text{×}{10}^{-08}$
HS36		$3$	$7$	$\mathrm{2.946\hspace{0.17em}660}\text{×}{10}^{-04}$	$\mathrm{5.016\hspace{0.17em}933}\text{×}{10}^{-03}$	$\mathrm{5.453\hspace{0.17em}507}\text{×}{10}^{-03}$	$-\mathrm{1.832\hspace{0.17em}051}\text{×}{10}^{02}$	$-\mathrm{1.832\hspace{0.17em}051}\text{×}{10}^{02}$	$\mathrm{7.105\hspace{0.17em}427}\text{×}{10}^{-14}$	$\mathrm{9.974\hspace{0.17em}060}\text{×}{10}^{-09}$
HS37		$3$	$8$	$\mathrm{2.950\hspace{0.17em}220}\text{×}{10}^{-04}$	$\mathrm{5.691\hspace{0.17em}497}\text{×}{10}^{-03}$	$\mathrm{4.989\hspace{0.17em}736}\text{×}{10}^{-03}$	$-\mathrm{1.832\hspace{0.17em}051}\text{×}{10}^{02}$	$-\mathrm{1.832\hspace{0.17em}051}\text{×}{10}^{02}$	$\mathrm{9.947\hspace{0.17em}598}\text{×}{10}^{-14}$	$\mathrm{9.974\hspace{0.17em}059}\text{×}{10}^{-09}$
HS41		$4$	$9$	$\mathrm{1.630\hspace{0.17em}200}\text{×}{10}^{-04}$	$\mathrm{4.818\hspace{0.17em}620}\text{×}{10}^{-03}$	$\mathrm{5.279\hspace{0.17em}881}\text{×}{10}^{-03}$	$-\mathrm{1.562\hspace{0.17em}500}\text{×}{10}^{-02}$	$-\mathrm{1.562\hspace{0.17em}500}\text{×}{10}^{-02}$	$\mathrm{5.551\hspace{0.17em}115}\text{×}{10}^{-17}$	$\mathrm{1.329\hspace{0.17em}137}\text{×}{10}^{-08}$
HS44		$4$	$10$	$\mathrm{4.904\hspace{0.17em}870}\text{×}{10}^{-04}$	$\mathrm{8.343\hspace{0.17em}742}\text{×}{10}^{-03}$	$\mathrm{6.473\hspace{0.17em}175}\text{×}{10}^{-03}$	$-\mathrm{1.304\hspace{0.17em}760}\text{×}{10}^{00}$	$-\mathrm{1.304\hspace{0.17em}760}\text{×}{10}^{00}$	$\mathrm{2.442\hspace{0.17em}491}\text{×}{10}^{-15}$	$\mathrm{3.577\hspace{0.17em}725}\text{×}{10}^{-08}$
HS44NEW		$4$	$10$	$\mathrm{6.557\hspace{0.17em}700}\text{×}{10}^{-04}$	$\mathrm{4.320\hspace{0.17em}227}\text{×}{10}^{-03}$	$\mathrm{4.416\hspace{0.17em}368}\text{×}{10}^{-03}$	$-\mathrm{1.847\hspace{0.17em}553}\text{×}{10}^{00}$	$-\mathrm{1.847\hspace{0.17em}553}\text{×}{10}^{00}$	$\mathrm{2.220\hspace{0.17em}446}\text{×}{10}^{-15}$	$\mathrm{8.037\hspace{0.17em}391}\text{×}{10}^{-09}$
HS55		$6$	$14$	$\mathrm{2.517\hspace{0.17em}480}\text{×}{10}^{-04}$	$\mathrm{5.743\hspace{0.17em}211}\text{×}{10}^{-03}$	$\mathrm{5.443\hspace{0.17em}675}\text{×}{10}^{-03}$	$-\mathrm{1.000\hspace{0.17em}000}\text{×}{10}^{00}$	$-\mathrm{2.236\hspace{0.17em}068}\text{×}{10}^{00}$	$\mathrm{4.000\hspace{0.17em}000}\text{×}{10}^{00}$	$\mathrm{5.616\hspace{0.17em}085}\text{×}{10}^{-05}$

Table 5. Table 5: Continuation of Table 4

Problem	Parameter	$n$	$m$	$t$	$t_{\text{ipopt}}^{\text{dense}}$	$t_{\text{ipopt}}^{\text{sparse}}$	$\delta f^{\ast }$	$\delta f_{\text{ipopt}}^{\ast }$	error	${\text{error}}_{\text{ipopt}}$
LEUVEN7		$360$	$1605$	$\mathrm{5.521\hspace{0.17em}594}\text{×}{10}^{-01}$	$\mathrm{6.015\hspace{0.17em}317}\text{×}{10}^{00}$	$\mathrm{3.374\hspace{0.17em}051}\text{×}{10}^{-01}$	$-\mathrm{4.858\hspace{0.17em}191}\text{×}{10}^{02}$	$-\mathrm{4.858\hspace{0.17em}192}\text{×}{10}^{02}$	$\mathrm{2.160\hspace{0.17em}050}\text{×}{10}^{-12}$	$\mathrm{5.712\hspace{0.17em}026}\text{×}{10}^{-05}$
LINCONT		$1257$	$2435$	$\mathrm{6.700\hspace{0.17em}989}\text{×}{10}^{01}$	$\mathrm{4.709\hspace{0.17em}799}\text{×}{10}^{02}$	$\mathrm{6.443\hspace{0.17em}721}\text{×}{10}^{-01}$	$-\mathrm{1.113\hspace{0.17em}381}\text{×}{10}^{-18}$	$-\mathrm{1.558\hspace{0.17em}091}\text{×}{10}^{-05}$	$\mathrm{2.899\hspace{0.17em}090}\text{×}{10}^{-10}$	$\mathrm{3.035\hspace{0.17em}288}\text{×}{10}^{-07}$
MPC10		$1530$	$3571$	$\mathrm{2.982\hspace{0.17em}143}\text{×}{10}^{01}$	$\mathrm{5.949\hspace{0.17em}118}\text{×}{10}^{02}$	$\mathrm{9.300\hspace{0.17em}218}\text{×}{10}^{-01}$	$-\mathrm{3.155\hspace{0.17em}428}\text{×}{10}^{03}$	$-\mathrm{3.155\hspace{0.17em}438}\text{×}{10}^{03}$	$\mathrm{1.636\hspace{0.17em}501}\text{×}{10}^{-11}$	$\mathrm{9.900\hspace{0.17em}918}\text{×}{10}^{-05}$
MPC11		$1530$	$3571$	$\mathrm{3.138\hspace{0.17em}670}\text{×}{10}^{01}$	$\mathrm{6.484\hspace{0.17em}969}\text{×}{10}^{02}$	$\mathrm{9.605\hspace{0.17em}708}\text{×}{10}^{-01}$	$-\mathrm{3.155\hspace{0.17em}339}\text{×}{10}^{03}$	$-\mathrm{3.155\hspace{0.17em}349}\text{×}{10}^{03}$	$\mathrm{2.258\hspace{0.17em}621}\text{×}{10}^{-11}$	$\mathrm{1.353\hspace{0.17em}352}\text{×}{10}^{-04}$
MPC12		$1530$	$3571$	$\mathrm{3.015\hspace{0.17em}642}\text{×}{10}^{01}$	$\mathrm{6.625\hspace{0.17em}879}\text{×}{10}^{02}$	$\mathrm{9.935\hspace{0.17em}781}\text{×}{10}^{-01}$	$-\mathrm{3.155\hspace{0.17em}068}\text{×}{10}^{03}$	$-\mathrm{3.155\hspace{0.17em}078}\text{×}{10}^{03}$	$\mathrm{1.451\hspace{0.17em}272}\text{×}{10}^{-11}$	$\mathrm{9.204\hspace{0.17em}114}\text{×}{10}^{-05}$
MPC13		$1530$	$3571$	$\mathrm{3.139\hspace{0.17em}020}\text{×}{10}^{01}$	$\mathrm{6.049\hspace{0.17em}481}\text{×}{10}^{02}$	$\mathrm{8.909\hspace{0.17em}481}\text{×}{10}^{-01}$	$-\mathrm{3.154\hspace{0.17em}808}\text{×}{10}^{03}$	$-\mathrm{3.154\hspace{0.17em}819}\text{×}{10}^{03}$	$\mathrm{1.516\hspace{0.17em}229}\text{×}{10}^{-11}$	$\mathrm{1.890\hspace{0.17em}977}\text{×}{10}^{-04}$
MPC14		$1530$	$3571$	$\mathrm{3.036\hspace{0.17em}349}\text{×}{10}^{01}$	$\mathrm{6.181\hspace{0.17em}861}\text{×}{10}^{02}$	$\mathrm{9.324\hspace{0.17em}590}\text{×}{10}^{-01}$	$-\mathrm{3.154\hspace{0.17em}810}\text{×}{10}^{03}$	$-\mathrm{3.154\hspace{0.17em}821}\text{×}{10}^{03}$	$\mathrm{2.793\hspace{0.17em}515}\text{×}{10}^{-11}$	$\mathrm{1.047\hspace{0.17em}976}\text{×}{10}^{-04}$
MPC15		$1530$	$3571$	$\mathrm{3.230\hspace{0.17em}341}\text{×}{10}^{01}$	$\mathrm{5.987\hspace{0.17em}057}\text{×}{10}^{02}$	$\mathrm{9.142\hspace{0.17em}717}\text{×}{10}^{-01}$	$-\mathrm{3.153\hspace{0.17em}846}\text{×}{10}^{03}$	$-\mathrm{3.153\hspace{0.17em}859}\text{×}{10}^{03}$	$\mathrm{2.210\hspace{0.17em}518}\text{×}{10}^{-11}$	$\mathrm{1.069\hspace{0.17em}349}\text{×}{10}^{-04}$
MPC16		$1530$	$3571$	$\mathrm{3.241\hspace{0.17em}830}\text{×}{10}^{01}$	$\mathrm{6.266\hspace{0.17em}897}\text{×}{10}^{02}$	$\mathrm{9.186\hspace{0.17em}500}\text{×}{10}^{-01}$	$-\mathrm{3.153\hspace{0.17em}841}\text{×}{10}^{03}$	$-\mathrm{3.153\hspace{0.17em}853}\text{×}{10}^{03}$	$\mathrm{1.616\hspace{0.17em}848}\text{×}{10}^{-11}$	$\mathrm{1.070\hspace{0.17em}100}\text{×}{10}^{-04}$
MPC2		$1530$	$3571$	$\mathrm{2.824\hspace{0.17em}377}\text{×}{10}^{01}$	$\mathrm{5.600\hspace{0.17em}366}\text{×}{10}^{02}$	$\mathrm{8.767\hspace{0.17em}303}\text{×}{10}^{-01}$	$-\mathrm{3.157\hspace{0.17em}144}\text{×}{10}^{03}$	$-\mathrm{3.157\hspace{0.17em}151}\text{×}{10}^{03}$	$\mathrm{1.161\hspace{0.17em}907}\text{×}{10}^{-11}$	$\mathrm{1.056\hspace{0.17em}836}\text{×}{10}^{-04}$
MPC3		$1530$	$3571$	$\mathrm{2.783\hspace{0.17em}808}\text{×}{10}^{01}$	$\mathrm{6.103\hspace{0.17em}022}\text{×}{10}^{02}$	$\mathrm{9.163\hspace{0.17em}411}\text{×}{10}^{-01}$	$-\mathrm{3.157\hspace{0.17em}321}\text{×}{10}^{03}$	$-\mathrm{3.157\hspace{0.17em}329}\text{×}{10}^{03}$	$\mathrm{1.224\hspace{0.17em}183}\text{×}{10}^{-11}$	$\mathrm{1.047\hspace{0.17em}022}\text{×}{10}^{-04}$
MPC4		$1530$	$3571$	$\mathrm{3.018\hspace{0.17em}797}\text{×}{10}^{01}$	$\mathrm{6.347\hspace{0.17em}866}\text{×}{10}^{02}$	$\mathrm{9.388\hspace{0.17em}815}\text{×}{10}^{-01}$	$-\mathrm{3.156\hspace{0.17em}964}\text{×}{10}^{03}$	$-\mathrm{3.156\hspace{0.17em}972}\text{×}{10}^{03}$	$\mathrm{1.548\hspace{0.17em}417}\text{×}{10}^{-11}$	$\mathrm{1.033\hspace{0.17em}535}\text{×}{10}^{-04}$
MPC5		$1530$	$3571$	$\mathrm{2.943\hspace{0.17em}601}\text{×}{10}^{01}$	$\mathrm{6.344\hspace{0.17em}659}\text{×}{10}^{02}$	$\mathrm{9.331\hspace{0.17em}457}\text{×}{10}^{-01}$	$-\mathrm{3.157\hspace{0.17em}075}\text{×}{10}^{03}$	$-\mathrm{3.157\hspace{0.17em}083}\text{×}{10}^{03}$	$\mathrm{1.574\hspace{0.17em}789}\text{×}{10}^{-11}$	$\mathrm{1.044\hspace{0.17em}170}\text{×}{10}^{-04}$
MPC6		$1530$	$3571$	$\mathrm{3.035\hspace{0.17em}928}\text{×}{10}^{01}$	$\mathrm{5.720\hspace{0.17em}201}\text{×}{10}^{02}$	$\mathrm{9.081\hspace{0.17em}594}\text{×}{10}^{-01}$	$-\mathrm{3.151\hspace{0.17em}828}\text{×}{10}^{03}$	$-\mathrm{3.151\hspace{0.17em}837}\text{×}{10}^{03}$	$\mathrm{1.558\hspace{0.17em}678}\text{×}{10}^{-11}$	$\mathrm{1.905\hspace{0.17em}436}\text{×}{10}^{-04}$
MPC7		$1530$	$3571$	$\mathrm{2.897\hspace{0.17em}083}\text{×}{10}^{01}$	$\mathrm{6.455\hspace{0.17em}158}\text{×}{10}^{02}$	$\mathrm{9.593\hspace{0.17em}935}\text{×}{10}^{-01}$	$-\mathrm{3.155\hspace{0.17em}357}\text{×}{10}^{03}$	$-\mathrm{3.155\hspace{0.17em}366}\text{×}{10}^{03}$	$\mathrm{1.857\hspace{0.17em}497}\text{×}{10}^{-11}$	$\mathrm{1.770\hspace{0.17em}294}\text{×}{10}^{-04}$
MPC8		$1530$	$3571$	$\mathrm{2.761\hspace{0.17em}454}\text{×}{10}^{01}$	$\mathrm{5.973\hspace{0.17em}750}\text{×}{10}^{02}$	$\mathrm{9.020\hspace{0.17em}797}\text{×}{10}^{-01}$	$-\mathrm{3.155\hspace{0.17em}408}\text{×}{10}^{03}$	$-\mathrm{3.155\hspace{0.17em}418}\text{×}{10}^{03}$	$\mathrm{1.359\hspace{0.17em}729}\text{×}{10}^{-11}$	$\mathrm{1.503\hspace{0.17em}306}\text{×}{10}^{-04}$
MPC9		$1530$	$3571$	$\mathrm{2.949\hspace{0.17em}440}\text{×}{10}^{01}$	$\mathrm{6.456\hspace{0.17em}766}\text{×}{10}^{02}$	$\mathrm{9.470\hspace{0.17em}979}\text{×}{10}^{-01}$	$-\mathrm{3.155\hspace{0.17em}392}\text{×}{10}^{03}$	$-\mathrm{3.155\hspace{0.17em}402}\text{×}{10}^{03}$	$\mathrm{1.521\hspace{0.17em}752}\text{×}{10}^{-11}$	$\mathrm{1.044\hspace{0.17em}509}\text{×}{10}^{-04}$
NASH		$72$	$138$	$\mathrm{2.023\hspace{0.17em}284}\text{×}{10}^{-02}$	$\mathrm{9.083\hspace{0.17em}129}\text{×}{10}^{-02}$	$\mathrm{2.347\hspace{0.17em}952}\text{×}{10}^{-02}$	$-\mathrm{1.335\hspace{0.17em}068}\text{×}{10}^{-16}$	$-\mathrm{2.562\hspace{0.17em}575}\text{×}{10}^{-05}$	$\mathrm{1.101\hspace{0.17em}341}\text{×}{10}^{-13}$	$\mathrm{2.652\hspace{0.17em}532}\text{×}{10}^{-07}$
NCVXQP1	N=1000	$1000$	$2500$	$\mathrm{3.468\hspace{0.17em}481}\text{×}{10}^{-01}$	$\mathrm{5.037\hspace{0.17em}756}\text{×}{10}^{01}$	$\mathrm{2.530\hspace{0.17em}511}\text{×}{10}^{-01}$	$-\mathrm{7.824\hspace{0.17em}962}\text{×}{10}^{04}$	$-\mathrm{7.824\hspace{0.17em}962}\text{×}{10}^{04}$	$\mathrm{9.094\hspace{0.17em}947}\text{×}{10}^{-12}$	$\mathrm{3.340\hspace{0.17em}764}\text{×}{10}^{-07}$
NCVXQP2	N=1000	$1000$	$2500$	$\mathrm{3.615\hspace{0.17em}064}\text{×}{10}^{-01}$	$\mathrm{5.018\hspace{0.17em}531}\text{×}{10}^{01}$	$\mathrm{2.535\hspace{0.17em}988}\text{×}{10}^{-01}$	$-\mathrm{5.736\hspace{0.17em}680}\text{×}{10}^{04}$	$-\mathrm{5.736\hspace{0.17em}680}\text{×}{10}^{04}$	$\mathrm{1.273\hspace{0.17em}293}\text{×}{10}^{-11}$	$\mathrm{2.509\hspace{0.17em}192}\text{×}{10}^{-07}$
NCVXQP3	N=1000	$1000$	$2500$	$\mathrm{3.681\hspace{0.17em}914}\text{×}{10}^{-01}$	$\mathrm{4.789\hspace{0.17em}657}\text{×}{10}^{01}$	$\mathrm{2.396\hspace{0.17em}625}\text{×}{10}^{-01}$	$-\mathrm{4.998\hspace{0.17em}258}\text{×}{10}^{04}$	$-\mathrm{4.998\hspace{0.17em}258}\text{×}{10}^{04}$	$\mathrm{1.273\hspace{0.17em}293}\text{×}{10}^{-11}$	$\mathrm{2.732\hspace{0.17em}257}\text{×}{10}^{-07}$
NCVXQP4	N=1000	$1000$	$2250$	$\mathrm{3.495\hspace{0.17em}667}\text{×}{10}^{-01}$	$\mathrm{3.453\hspace{0.17em}693}\text{×}{10}^{01}$	$\mathrm{2.044\hspace{0.17em}327}\text{×}{10}^{-01}$	$-\mathrm{7.824\hspace{0.17em}962}\text{×}{10}^{04}$	$-\mathrm{7.824\hspace{0.17em}962}\text{×}{10}^{04}$	$\mathrm{5.184\hspace{0.17em}120}\text{×}{10}^{-11}$	$\mathrm{3.340\hspace{0.17em}765}\text{×}{10}^{-07}$
NCVXQP5	N=1000	$1000$	$2250$	$\mathrm{3.553\hspace{0.17em}980}\text{×}{10}^{-01}$	$\mathrm{3.450\hspace{0.17em}249}\text{×}{10}^{01}$	$\mathrm{2.057\hspace{0.17em}896}\text{×}{10}^{-01}$	$-\mathrm{5.736\hspace{0.17em}680}\text{×}{10}^{04}$	$-\mathrm{5.736\hspace{0.17em}680}\text{×}{10}^{04}$	$\mathrm{1.714\hspace{0.17em}398}\text{×}{10}^{-09}$	$\mathrm{2.509\hspace{0.17em}195}\text{×}{10}^{-07}$

Table 6. Table 6: Continuation of Table 4

Problem	Parameter	$n$	$m$	$t$	$t_{\text{ipopt}}^{\text{dense}}$	$t_{\text{ipopt}}^{\text{sparse}}$	$\delta f^{\ast }$	$\delta f_{\text{ipopt}}^{\ast }$	error	${\text{error}}_{\text{ipopt}}$
NCVXQP6	N=1000	$1000$	$2250$	$\mathrm{3.642\hspace{0.17em}895}\text{×}{10}^{-01}$	$\mathrm{3.295\hspace{0.17em}662}\text{×}{10}^{01}$	$\mathrm{1.933\hspace{0.17em}023}\text{×}{10}^{-01}$	$-\mathrm{4.998\hspace{0.17em}258}\text{×}{10}^{04}$	$-\mathrm{4.998\hspace{0.17em}258}\text{×}{10}^{04}$	$\mathrm{5.820\hspace{0.17em}766}\text{×}{10}^{-11}$	$\mathrm{2.732\hspace{0.17em}277}\text{×}{10}^{-07}$
NCVXQP7	N=1000	$1000$	$2750$	$\mathrm{3.584\hspace{0.17em}502}\text{×}{10}^{-01}$	$\mathrm{7.042\hspace{0.17em}072}\text{×}{10}^{01}$	$\mathrm{2.941\hspace{0.17em}084}\text{×}{10}^{-01}$	$-\mathrm{7.824\hspace{0.17em}962}\text{×}{10}^{04}$	$-\mathrm{7.824\hspace{0.17em}962}\text{×}{10}^{04}$	$\mathrm{5.129\hspace{0.17em}550}\text{×}{10}^{-10}$	$\mathrm{3.340\hspace{0.17em}764}\text{×}{10}^{-07}$
NCVXQP8	N=1000	$1000$	$2750$	$\mathrm{3.669\hspace{0.17em}888}\text{×}{10}^{-01}$	$\mathrm{7.371\hspace{0.17em}150}\text{×}{10}^{01}$	$\mathrm{2.933\hspace{0.17em}839}\text{×}{10}^{-01}$	$-\mathrm{5.736\hspace{0.17em}680}\text{×}{10}^{04}$	$-\mathrm{5.736\hspace{0.17em}680}\text{×}{10}^{04}$	$\mathrm{8.305\hspace{0.17em}506}\text{×}{10}^{-09}$	$\mathrm{2.509\hspace{0.17em}189}\text{×}{10}^{-07}$
NCVXQP9	N=1000	$1000$	$2750$	$\mathrm{3.835\hspace{0.17em}907}\text{×}{10}^{-01}$	$\mathrm{6.705\hspace{0.17em}865}\text{×}{10}^{01}$	$\mathrm{2.787\hspace{0.17em}119}\text{×}{10}^{-01}$	$-\mathrm{4.998\hspace{0.17em}258}\text{×}{10}^{04}$	$-\mathrm{4.998\hspace{0.17em}258}\text{×}{10}^{04}$	$\mathrm{2.119\hspace{0.17em}577}\text{×}{10}^{-08}$	$\mathrm{2.732\hspace{0.17em}247}\text{×}{10}^{-07}$
PENTAGON		$6$	$15$	$\mathrm{1.087\hspace{0.17em}885}\text{×}{10}^{-03}$	$\mathrm{9.631\hspace{0.17em}484}\text{×}{10}^{-03}$	$\mathrm{9.260\hspace{0.17em}053}\text{×}{10}^{-03}$	$-\mathrm{2.699\hspace{0.17em}815}\text{×}{10}^{-02}$	$-\mathrm{2.266\hspace{0.17em}403}\text{×}{10}^{-02}$	$\mathrm{8.267\hspace{0.17em}692}\text{×}{10}^{-15}$	$\mathrm{2.092\hspace{0.17em}105}\text{×}{10}^{-06}$
QC		$9$	$22$	$\mathrm{2.120\hspace{0.17em}020}\text{×}{10}^{-04}$	$\mathrm{1.348\hspace{0.17em}843}\text{×}{10}^{-02}$	$\mathrm{1.216\hspace{0.17em}467}\text{×}{10}^{-02}$	$-\mathrm{1.114\hspace{0.17em}234}\text{×}{10}^{02}$	$-\mathrm{1.114\hspace{0.17em}236}\text{×}{10}^{02}$	$\mathrm{0.000\hspace{0.17em}000}\text{×}{10}^{00}$	$\mathrm{1.875\hspace{0.17em}249}\text{×}{10}^{-05}$
QCNEW		$9$	$21$	$\mathrm{2.276\hspace{0.17em}460}\text{×}{10}^{-04}$	$\mathrm{1.657\hspace{0.17em}338}\text{×}{10}^{-02}$	$\mathrm{1.513\hspace{0.17em}761}\text{×}{10}^{-02}$	$-\mathrm{1.653\hspace{0.17em}931}\text{×}{10}^{00}$	$-\mathrm{4.112\hspace{0.17em}503}\text{×}{10}^{01}$	$\mathrm{0.000\hspace{0.17em}000}\text{×}{10}^{00}$	$\mathrm{1.990\hspace{0.17em}100}\text{×}{10}^{-04}$
QPNBLEND		$83$	$157$	$\mathrm{1.501\hspace{0.17em}935}\text{×}{10}^{-02}$	$\mathrm{6.937\hspace{0.17em}614}\text{×}{10}^{-02}$	$\mathrm{2.064\hspace{0.17em}464}\text{×}{10}^{-02}$	$-\mathrm{9.282\hspace{0.17em}172}\text{×}{10}^{-03}$	$-\mathrm{9.283\hspace{0.17em}588}\text{×}{10}^{-03}$	$\mathrm{2.504\hspace{0.17em}663}\text{×}{10}^{-13}$	$\mathrm{4.679\hspace{0.17em}297}\text{×}{10}^{-05}$
QPNBOEI1		$384$	$891$	$\mathrm{6.116\hspace{0.17em}935}\text{×}{10}^{-01}$	$\mathrm{1.115\hspace{0.17em}413}\text{×}{10}^{01}$	$\mathrm{1.653\hspace{0.17em}756}\text{×}{10}^{-01}$	$-\mathrm{4.366\hspace{0.17em}064}\text{×}{10}^{01}$	$-\mathrm{4.366\hspace{0.17em}092}\text{×}{10}^{01}$	$\mathrm{7.786\hspace{0.17em}588}\text{×}{10}^{-12}$	$\mathrm{8.287\hspace{0.17em}324}\text{×}{10}^{-05}$
QPNBOEI2		$143$	$363$	$\mathrm{4.073\hspace{0.17em}949}\text{×}{10}^{-02}$	$\mathrm{1.011\hspace{0.17em}505}\text{×}{10}^{00}$	$\mathrm{6.225\hspace{0.17em}512}\text{×}{10}^{-02}$	$-\mathrm{1.003\hspace{0.17em}864}\text{×}{10}^{03}$	$-\mathrm{1.004\hspace{0.17em}486}\text{×}{10}^{03}$	$\mathrm{1.224\hspace{0.17em}766}\text{×}{10}^{-09}$	$\mathrm{8.420\hspace{0.17em}519}\text{×}{10}^{-05}$
QPNSTAIR		$467$	$905$	$\mathrm{7.413\hspace{0.17em}311}\text{×}{10}^{-01}$	$\mathrm{1.500\hspace{0.17em}511}\text{×}{10}^{01}$	$\mathrm{2.324\hspace{0.17em}631}\text{×}{10}^{-01}$	$-\mathrm{1.150\hspace{0.17em}055}\text{×}{10}^{02}$	$-\mathrm{1.150\hspace{0.17em}130}\text{×}{10}^{02}$	$\mathrm{4.352\hspace{0.17em}719}\text{×}{10}^{-11}$	$\mathrm{2.310\hspace{0.17em}892}\text{×}{10}^{-04}$
SOSQP1	N=1000	$2000$	$5001$	$\mathrm{1.025\hspace{0.17em}841}\text{×}{10}^{00}$	$\mathrm{3.387\hspace{0.17em}020}\text{×}{10}^{02}$	$\mathrm{2.187\hspace{0.17em}173}\text{×}{10}^{-01}$	$-\mathrm{5.000\hspace{0.17em}000}\text{×}{10}^{-01}$	$-\mathrm{5.000\hspace{0.17em}000}\text{×}{10}^{-01}$	$\mathrm{3.108\hspace{0.17em}624}\text{×}{10}^{-14}$	$\mathrm{4.572\hspace{0.17em}343}\text{×}{10}^{-09}$
SOSQP2	N=1000	$2000$	$5001$	$\mathrm{1.189\hspace{0.17em}768}\text{×}{10}^{00}$	$\mathrm{2.171\hspace{0.17em}972}\text{×}{10}^{03}$	$\mathrm{1.061\hspace{0.17em}519}\text{×}{10}^{00}$	$-\mathrm{5.000\hspace{0.17em}000}\text{×}{10}^{-01}$	$-\mathrm{5.000\hspace{0.17em}000}\text{×}{10}^{-01}$	$\mathrm{4.940\hspace{0.17em}492}\text{×}{10}^{-15}$	$\mathrm{4.988\hspace{0.17em}599}\text{×}{10}^{-09}$
STATIC3		$434$	$240$	$\mathrm{9.991\hspace{0.17em}469}\text{×}{10}^{-01}$	$\mathrm{7.633\hspace{0.17em}038}\text{×}{10}^{-01}$	$\mathrm{5.708\hspace{0.17em}679}\text{×}{10}^{-02}$	$-\mathrm{1.473\hspace{0.17em}959}\text{×}{10}^{04}$	$-\mathrm{1.473\hspace{0.17em}959}\text{×}{10}^{04}$	$\mathrm{4.599\hspace{0.17em}613}\text{×}{10}^{-06}$	$\mathrm{1.455\hspace{0.17em}835}\text{×}{10}^{-04}$
STNQP1	P=10	$1025$	$2560$	$\mathrm{3.909\hspace{0.17em}290}\text{×}{10}^{-01}$	$\mathrm{5.022\hspace{0.17em}181}\text{×}{10}^{01}$	$\mathrm{1.634\hspace{0.17em}204}\text{×}{10}^{-01}$	$-\mathrm{1.629\hspace{0.17em}919}\text{×}{10}^{03}$	$-\mathrm{1.629\hspace{0.17em}919}\text{×}{10}^{03}$	$\mathrm{2.557\hspace{0.17em}954}\text{×}{10}^{-13}$	$\mathrm{1.049\hspace{0.17em}186}\text{×}{10}^{-08}$
STNQP2	P=10	$1025$	$2560$	$\mathrm{3.919\hspace{0.17em}327}\text{×}{10}^{-01}$	$\mathrm{4.356\hspace{0.17em}643}\text{×}{10}^{01}$	$\mathrm{1.601\hspace{0.17em}439}\text{×}{10}^{-01}$	$-\mathrm{1.629\hspace{0.17em}919}\text{×}{10}^{03}$	$-\mathrm{1.629\hspace{0.17em}919}\text{×}{10}^{03}$	$\mathrm{4.547\hspace{0.17em}474}\text{×}{10}^{-13}$	$\mathrm{1.047\hspace{0.17em}931}\text{×}{10}^{-08}$

Equations

$$\begin{array}[]{ll}\text{minimize}&f(x):=\frac{1}{2}x^{T}Px+q^{T}x\\ \text{subject to}&r_{\min}\leq\mathopen{}\mathclose{{}\left\lVert x}\right\rVert_{2}\leq r_{\max}\\ &Ax\leq b\end{array}$$

$$\begin{array}[]{ll}\text{minimize}&\frac{1}{2}x^{T}Px+q^{T}x\\ \text{subject to}&\mathopen{}\mathclose{{}\left\lVert x}\right\rVert_{2}=r\\ &Ax=b.\end{array}$$

$$M:=\left[\begin{array}{cc}-P& q{q}^T/r^2\\ I& -P\end{array}\right].$$

$$\begin{array}[]{ll}\text{minimize}&\frac{1}{2}x^{T}Px+q^{T}x\\ \text{subject to}&\mathopen{}\mathclose{{}\left\lVert x}\right\rVert_{2}=r,\\ \end{array}$$

$$Px+q+\mu x=0\Rightarrow x=-(P+\mu I{)}^{-1}q,$$

$$\mathopen{}\mathclose{{}\left\lVert x}\right\rVert_{2}^{2}=r^{2}\Leftrightarrow\mathopen{}\mathclose{{}\left\lVert(P+\mu I)^{-1}q}\right\rVert_{2}^{2}-r^{2}=0.$$

$$\displaystyle\mathopen{}\mathclose{{}\left\lVert(P+\mu I)^{-1}q}\right\rVert_{2}^{2}$$

$$\displaystyle=\mathopen{}\mathclose{{}\left\lVert(\Lambda+\mu I)^{-1}W^{T}q}\right\rVert_{2}^{2}=\sum_{i=1}^{n}{\frac{(w_{i}^{T}q)^{2}}{(\lambda_{i}+\mu)^{2}}},$$

$$s(\mu ):=\sum _{i=1}^n\frac{(w_i^{T}q{)}^2}{({\lambda }_i+\mu {)}^2}-r^2=0.$$

$$\underset{:=M}{\underbrace{\left[\begin{array}{cc}-P& q{q}^T/r^2\\ I& -P\end{array}\right]}}\left[\begin{array}{c}{z}_1\\ z_2\end{array}\right]=\mu \left[\begin{array}{c}{z}_1\\ z_2\end{array}\right].$$

$$x^{*}=-\text{sign}(q^{T}z_{2}^{*})r\frac{z_{1}^{*}}{\mathopen{}\mathclose{{}\left\lVert z_{1}^{*}}\right\rVert_{2}},$$

$$-P{z}_1^{\ast }+q{q}^T{z}_2^{\ast }/r^2={\mu }^{\ast }{z}_1^{\ast }\Rightarrow (z_2^{\ast }{)}^T(P+{\mu }^{\ast }I)z_1^{\ast }=(z_2^{\ast }{)}^{T}q{q}^T{z}_2^{\ast }/r^2,$$

$$\mathopen{}\mathclose{{}\left\lVert z_{1}^{*}}\right\rVert^{2}_{2}=(q^{T}z_{2}^{*}/r)^{2}\Rightarrow\mathopen{}\mathclose{{}\left\lVert-\frac{r^{2}}{q^{T}z_{2}^{*}}z_{1}^{*}}\right\rVert_{2}^{2}=r^{2}\Rightarrow\mathopen{}\mathclose{{}\left\lVert y^{*}}\right\rVert_{2}=r.$$

$$-P{z}_1^{\ast }+q{q}^T{z}_2^{\ast }/r^2={\mu }^{\ast }{z}_1^{\ast }\Rightarrow P\frac{-z_1^{\ast }{r}^2}{q^T{z}_2^{\ast }}+q+{\mu }^{\ast }\frac{-z_1^{\ast }{r}^2}{q^T{z}_2^{\ast }}=0\Rightarrow P{y}^{\ast }+q+{\mu }^{\ast }{y}^{\ast }=0.$$

$$y^{*}=y^{*}\frac{r}{\mathopen{}\mathclose{{}\left\lVert y^{*}}\right\rVert_{2}}=\frac{-r^{2}z_{1}^{*}}{q^{T}z_{2}^{*}}{\frac{|q^{T}z_{2}^{*}|}{\mathopen{}\mathclose{{}\left\lVert-r^{2}z_{1}^{*}}\right\rVert}}_{2}r=-\text{sign}(q^{T}z_{2}^{*})r\frac{z_{1}^{*}}{\mathopen{}\mathclose{{}\left\lVert z_{1}^{*}}\right\rVert_{2}}=x^{*}.$$

$$(P-{\lambda }_{1}I)x+q=0$$

$$\displaystyle\mathopen{}\mathclose{{}\left\lVert x}\right\rVert_{2}=r.$$

$$\begin{array}[]{ll}\text{minimize}&\frac{1}{2}x^{T}Px+q^{T}x\\ \text{subject to}&\mathopen{}\mathclose{{}\left\lVert x}\right\rVert_{2}=r\\ &Ax=0\end{array}$$

$$\begin{array}[]{ll}\text{minimize}&\frac{1}{2}y^{T}\tilde{P}y+\tilde{q}^{T}y\\ \text{subject to}&\mathopen{}\mathclose{{}\left\lVert y}\right\rVert_{2}=r\end{array}$$

$$\tilde{P}:=Z^{T}PZ,\tilde{q}:=Z^{T}q.$$

$$\tilde{M}:=\left[\begin{array}{cc}-\tilde{P}& \tilde{q}{\tilde{q}}^T/r^2\\ I& -\tilde{P}\end{array}\right].$$

$$x^{\ast }=Z{y}^{\ast }.$$

$$\tilde{M}=\left[\begin{array}{cc}-\tilde{P}& \tilde{q}{\tilde{q}}^T/r^2\\ I& -\tilde{P}\end{array}\right]=\left[\begin{array}{cc}{Z}^T& 0\\ 0& Z^T\end{array}\right]\left[\begin{array}{cc}-P& q{q}^T/r^2\\ I& -P\end{array}\right]\left[\begin{array}{cc}Z& 0\\ 0& Z\end{array}\right].$$

$$\left[\begin{array}{cc}Z{Z}^T& 0\\ 0& Z{Z}^T\end{array}\right]\underset{=M}{\underbrace{\left[\begin{array}{cc}-P& q{q}^T/r^2\\ I& -P\end{array}\right]}}\left[\begin{array}{cc}Z{Z}^T& 0\\ 0& Z{Z}^T\end{array}\right]$$

$$s^{\prime }(\alpha )=-2\sum _{j=1}^n\frac{(w_j^{T}q{)}^2}{({\lambda }_j+\alpha {)}^3}$$

$$\text{Im}(s(a+bi))=-\sum _{j=1}^n\frac{2({\lambda }_j+a)b}{(({\lambda }_j+a{)}^2+b^2{)}^2}(w_j^{T}q{)}^2,$$

$$\frac{1}{2}{x}^T(P-\alpha I)x+q^{T}x$$

$$x^{\ast }$$

$$\displaystyle=-\text{sign}(q^{T}z_{2}^{*})r\frac{z_{1}^{*}}{\mathopen{}\mathclose{{}\left\lVert z_{1}^{*}}\right\rVert_{2}}$$

$$(\tilde{P}+{\mu }^{\ast }I)y+\tilde{q}=0$$

$$(P+{\mu }^{\ast }I)x+q=0$$

Full text

††journal: Mathematical Programming

1 Introduction

Optimization with spherical \mathopen{}\mathclose{{}\left\lVert x}\right\rVert_{2}=1 or norm constraints r_{\min}\leq\mathopen{}\mathclose{{}\left\lVert x}\right\rVert_{2}\leq r_{\max} is used in a number of scientific fields. Many optimization problems, such as eigenvalue problems Saad2011 , dimensionality reduction Jolliffe2011 and compressed sensing Boufounos2008 , are naturally posed on or in the unit sphere, while norm constraints are often useful for regularization, e.g. in general nonlinear Conn2000 or robust optimization Jeyakumar2014 . In this paper we consider the solution of norm constrained problems in the form

[TABLE]

where $x\in\mathbb{R}^{n}$ is the decision variable, $P\in\mathbb{S}^{n}$, i.e. an $n\times n$ symmetric matrix, $A=[a_{1}^{T}\dots a_{m}^{T}]\in\mathbb{R}^{m\times n}$, $q\in\mathbb{R}^{n}$, $b=[b_{1}\dots b_{m}]^{T}\in\mathbb{R}^{m}$ and $r_{\min},r_{\max}$ non-negative scalars.

It is generally difficult to solve problems of the form (P) due to the nonconvexity of the lower bounding norm constraint and the potential indefiniteness of the matrix $P$, which renders many, but certainly not all Jeyakumar2014 , of these problems intractable. Even finding a feasible point for (P) or testing if a first-order critical point of (P) is a local minimizer can be intractable (see Proposition 4 and Gould2005 ). As a result, we restrict our attention to the search for first order critical points and we assume that a feasible point is given or can be computed.

A specific tractable variation of (P) is the famous Trust-Region Subproblem111Note that the TRS is typically defined in the ball \mathopen{}\mathclose{{}\left\lVert x}\right\rVert\leq r rather than the sphere \mathopen{}\mathclose{{}\left\lVert x}\right\rVert=r. We consider only the boundary solution as this fits better the needs of this paper.:

[TABLE]

Perhaps surprisingly, the TRS can be solved to global optimality despite its nonconvexity, and global solution of the TRS has been widely studied in the literature. Specific approaches to its solution include exact Semidefinite Reformulations, Gradient Steps, Truncated Conjugate Gradient Steps, Truncated Lanczos Steps and Newton root-finding, with associated software packages for its solution. The reader can find details in Conn2000 and (Nocedal2006, , Chapter 4).

We suggest an active-set algorithm for the solution of (P) that, as we proceed to describe, takes advantage of the existing efficient global solution methods for the TRS. Using such an approach, (P) is optimized by solving a series of equality constrained subproblems. Each subproblem is a modification of (P) where a subset of its inequalities (called the working-set) is replaced by equalities while the rest are ignored. If a subproblem does not include a norm constraint then it is called an Equality-constrained Quadratic Problem (EQP); otherwise it is a Trust Region Subproblem (TRS). As described earlier, the global solution for the TRS, and also for the EQP Gould2001 , has been widely studied in the literature. Thus we can efficiently address the global solution of each subproblem. However, a complication exists for the TRS because, unlike EQP, it can exhibit local-nonglobal minimizers i.e. local minimizers that are not globally optimal, which complicate the analysis. Indeed, if the optimal working-set includes a norm constraint and the global solutions of the respective TRS subproblem are infeasible for (P), then the solution of (P) must be obtained from a local-nonglobal optimizer of the TRS subproblem.

Thus, an algorithm for computing local-nonglobal minimizers of the TRS is required for our active-set algorithm. Unfortunately, algorithms for detecting the presence of/computing local-nonglobal solutions of (TRS) are significantly less mature than their global counterparts. A notable exception is the work of Martinez Martinez1994 , which proved that there exists at most one local-nonglobal minimizer and gave conditions for its existence. Moreover, Martinez1994 presented a root-finding algorithm for the computation of the Lagrange Multiplier associated with the local-nonglobal minimizer that entailed the computation of the two smallest eigenvalues of the Hessian and the solution of a series of indefinite linear systems.

More recent work is based on a result from Adachi2017 , which shows that each KKT point of the TRS, in the absence of linear equality constraints $Ax=b$ (i.e. with a norm constraint on $x$ only), can be extracted from an eigenpair of

[TABLE]

This result was used in Salahi2017 to calculate the local-nonglobal minimizer under the assumption that it exists. We extend these results by proving that both checking if the local-nonglobal minimizer exists and computing it, in the case where it exists, can be performed at the same time by simply calculating the two rightmost eigenpairs of $M$. Crucially, and similarly to the results of Adachi2017 for the global minimizer, this allows the use of efficient factorization-free methods such as the Arnoldi method for the detection and computation of the local-nonglobal minimizer. Furthermore, we show that this approach can efficiently handle equality constraints $Ax=b$ in the Trust Region Subproblem, without the need for variable reduction, by means of a projected Arnoldi method.

The paper is organised as follows. In Section 2 we present an introduction to the TRS and present novel results for the detection/computation of its local-nonglobal minimizer and for the incorporation of linear equality constraints. In Section 3.1 we introduce an active-set algorithm for solving (P) starting from the special case where $r_{\min}=r_{\max}$, and then generalizing to any $r_{\min},r_{\max}$. Finally, in Section 4, we conclude with a series of numerical results.

Notation used:

$\mathbb{S}^{n}$ denotes the set of symmetric $n\times n$ real matrices. Given a vector (or matrix) $x$, $x^{T}$ denotes its transpose and $x^{H}$ its conjugate transpose.

2 The Trust-Region Subproblem

This section concerns the computation of global as well as local-nonglobal minimizers for (TRS). Although this section only considers solutions on the boundary \mathopen{}\mathclose{{}\left\lVert x}\right\rVert=r, the results can be trivially extended to the interior \mathopen{}\mathclose{{}\left\lVert x}\right\rVert\leq r. Indeed, if an interior solution exists then the TRS is essentially convex ($P$ is positive semidefinite in the nullspace of $A$ (Gould1999, , §1)), thus no local-nonglobal solution exists and the interior solution is obtained by solving an equality constrained quadratic problem. The presence of a (necessarily global) interior solution can be detected by checking the sign of the Lagrange Multiplier of the global boundary solution(s) (Martinez1994, , Lemma 2.2).

For clarity of exposition, we will first assume that there are no linear equality constraints in (TRS), resulting in the following problem:

[TABLE]

where $x\in\mathbb{R}^{n}$ is the decision variable, and $P\in\mathbb{S}^{n}$, $q\in\mathbb{R}^{n},r\in\mathbb{R}_{+}$ are the problem data. For the rest of the section we will assume that $n>1$, excluding the trivial one dimensional case where the feasible set consists of two points, and that $r>0$. We will then show in §2.1 how to extend the results of this section to allow for the inclusion of linear equality constraints. Finally, in §2.A we prove the main result of this section, which is presented in a separate subsection to facilitate the presentation.

In the sequel we will make frequent use of the eigendecomposition of $P:=W\Lambda W^{T}$ defined by an orthonormal matrix $W:=[w_{1}\dots w_{n}]$ and $\Lambda:=\mathrm{diag}(\lambda_{1},\dots,\lambda_{n})$ where $\lambda_{1}\leq\dots\leq\lambda_{n}$.

Every KKT point of (T) satisfies

[TABLE]

where $\mu$ is a Lagrange multiplier associated with the norm constraint222Technically, $\mu$ is the Lagrange multiplier of the equivalent constraint \frac{1}{2}\mathopen{}\mathclose{{}\left\lVert x}\right\rVert_{2}^{2}=\frac{1}{2}r^{2}. We define (T) with \mathopen{}\mathclose{{}\left\lVert x}\right\rVert_{2}=r for simplicity and to match the notation of Adachi2017 .. Equation (1) is well defined when $\mu\neq-\lambda_{i},\;i=1,\dots,n$. For the purposes of clarity of the subsequent introduction, which follows (Nocedal2006, , §4), we will assume that this holds. However, the conclusions of this section are independent of this assumption.

Using the feasibility of $x$, we have

[TABLE]

Then, noting that

[TABLE]

we can express the rightmost condition in (2) as

[TABLE]

Determining the KKT points of (T) is therefore equivalent to finding the roots of $s$, which is often referred to as the secular equation Nocedal2006 . We depict $s$ for a particular choice of $P,q,r$ in Figure 1. As might be expected from the tight connection between polynomial root-finding problems and eigenproblems, solving $s(\mu)=0$ is equivalent to the following eigenproblem (Adachi2017, , Equation (22)):

[TABLE]

This elegant result was first noted more than 50 years ago (Forsythe1965, , Equation (2.21)), only to be later disregarded as inefficient by some of the same authors Gander1989 , and then rediscovered by Adachi2017 who highlighted its great applicability owing to the remarkable efficacy of modern eigensolvers Lehoucq1998 .

The relation between the spectrum of $M$ and the Lagrange multipliers of (T) is formally stated below:

Lemma 1

The spectrum of $M$ includes the Lagrange multiplier of every KKT point of (T).

Proof

This is a consequence of Proposition 2 of subsection 2.A and (Forsythe1965, , Theorem 4.1). ∎

Note that the eigenvalues of $M$ are complex in general since $M$ is asymmetric. Furthermore, Lemma 1 and all the subsequent results allow for potential “degenerate” KKT points with $\mu=-\lambda_{i}$.

We will focus on the KKT points that are minimizers of (T) rather than maximizers or saddle points. The characterization of global minimizers is widely known in the literature. The following theorem shows that the Lagrange multiplier $\mu^{g}$ of the global minimizer corresponds to the rightmost eigenvalue of $M$ in $\mathbb{C}$, which is always real:

Theorem 2.1

A KKT point of (T) with Lagrange multiplier $\mu^{g}$ is a global minimizer if and only if $\mu^{g}$ is the rightmost eigenvalue of $M$. Furthermore, we necessarily have that $\mu^{g}\in[-\lambda_{1},\infty)$.

Proof

See (Adachi2017, , Theorem 3.2). ∎

In general, (T) is guaranteed to possess a unique global optimizer, except perhaps in a special case that has deservedly been given a special name:

Definition​​ (Hard case)

Problem (T) belongs to the hard case when $\mu^{g}=-\lambda_{1}$.

We will see that in the hard case (T) is still guaranteed to have a global minimizer, but it is not necessarily unique. Furthermore, the computation of the minimizer(s) can be more challenging than in the “standard” case.

Moreover, due to the nonconvexity of (T) there can exist a local minimizer that is not global. We derive an analogous result for this minimizer that is called “local-nonglobal” Martinez1994 :

Theorem 2.2

Problem (T) has a second-order sufficient local-nonglobal minimizer if and only if it does not belong to the hard case and the second-rightmost eigenvalue of $M$ is real and simple. Furthermore, if such a minimizer exists, then its Lagrange multiplier $\mu^{\ell}$ is equal to this second-rightmost eigenvalue.

Proof

See §2.A. ∎

The Lagrange multipliers for the global and local-nonglobal minimizers of (T) can therefore be identified by calculating the two rightmost eigenvalues333Except perhaps when the local-nonglobal minimizer is not second-order sufficient. of $M$, which can be done efficiently e.g. with the Arnoldi method.

Furthermore, for each of the respective multipliers a corresponding minimizer can be extracted immediately by the respective eigenvectors $[z_{1}^{*};\;z_{2}^{*}]$ of $M$ using as

[TABLE]

unless $z_{1}^{*}=0$. Indeed, this follows from the next Proposition:

Proposition 1

Every eigenpair $(\mu^{*},[z_{1}^{*};\;z_{2}^{*}])$ of $M$ gives a KKT point $(\mu^{*},x^{*})$ for (T) where $x^{*}$ is defined by (5), unless $z_{1}^{*}=0$.

Proof

This proof is based on (Adachi2017, , §3.3) that shows the result only for the rightmost eigenpair of $M$.

We will first show the result for $y^{*}=\frac{-r^{2}}{q^{T}z_{2}^{*}}z_{1}^{*}$ and then show that $y^{*}=x^{*}$. Note that using $x^{*}$ is preferred over $y^{*}$ because it can avoid unnecessary numerical errors according to (Adachi2017, , §3.3). We also emphasize that all of the proof assumes that $z_{1}^{*}$ is nonzero.

For $y^{*}$, we have to show that it is well defined, feasible and stationary. Regarding the first two points, note that (4) (first row) gives

[TABLE]

or since $(P+\mu^{*}I)z_{2}^{*}=z_{1}^{*}$ due to (4) (second row), we have \mathopen{}\mathclose{{}\left\lVert z_{1}^{*}}\right\rVert^{2}={(q^{T}z_{2}^{*}/r)}^{2}, which implies that $q^{T}z_{2}^{*}\neq 0$ because $z_{1}^{*}\neq 0$ by assumption. Thus:

[TABLE]

To show stationarity, note that the first row of (4) also gives:

[TABLE]

Finally, we show that $y^{*}=x^{*}$:

[TABLE]

∎

Dealing with the hard case:

Extracting a solution $x^{*}$ via (5) is not possible when $z_{1}^{*}=0$, since (5) is not well-defined in that case. In this case, $(\mu^{*},z_{2}^{*})$ is an eigenpair of $-P$ due to (4). Note that, according to (Martinez1994, , Lemma 3.3), this can never happen for the local-nonglobal minimizer; thus we can always extract the local-nonglobal minimizer with (5). However, it is possible for global minimizers in the hard case, i.e. the case where $\mu^{*}=-\lambda_{1}$. In the hard case, the necessary and sufficient conditions for global optimality are, due to the KKT conditions of (T) and Theorem (2.1), the following:

[TABLE]

Note that $P-\lambda_{1}I$ is singular, and the above system of equations represents the intersection of an affine subspace with a sphere. This intersection must be non-empty, since (T) necessarily has a global minimizer as it arises from the minimization of a smooth function over a compact subset of $\mathbb{R}^{n}$. One can then solve (6)–(7), by computing the minimum length solution $y_{\text{min}}$ of $(P-\lambda_{1}I)x+q=0$ and return $x^{*}=y_{\text{min}}+\alpha\upsilon$ where $\alpha$ is a scalar such that \mathopen{}\mathclose{{}\left\lVert x^{*}}\right\rVert_{2}=\mathopen{}\mathclose{{}\left\lVert y_{\text{min}}+\alpha\upsilon}\right\rVert_{2}=r and $\upsilon$ is any null-vector of $P-\lambda_{1}I$. Interestingly, when $z_{1}^{*}=0$, $z_{2}^{*}$ is a null vector of $P-\lambda_{1}I$ due to (4), thus the only additional computation for extracting a solution when $z_{1}^{*}=0$ is finding the minimum-length solution of a symmetric linear system. This can be achieved e.g. with MINRES-QLP Choi2011 , or with the Conjugate Gradient method (Adachi2017, , Theorem 4.3).

2.1 Equality-constrained Trust-Region Subproblems

We will now extend the preceding results of this section to also account for the presence of linear equality constraints, $Ax=b$, in the TRS. We will show that, given an operator that projects an $n-$dimensional vector to the nullspace of $A$, we can calculate global as well as local-nonlocal minima of (TRS) by applying a projected Arnoldi method to $M$.

For simplicity, and without loss of generality444Problems with linear constraints of the form $Ax=b$ can be transformed to (8) via a change of variables $\tilde{x}=x+x_{0}$ where $x_{0}\perp\mathcal{N}(A)$ such that $Ax_{0}=b$, resulting in a quadratic cost and constraints: A\tilde{x}=0,\;\mathopen{}\mathclose{{}\left\lVert\tilde{x}}\right\rVert_{2}=\sqrt{r^{2}-\mathopen{}\mathclose{{}\left\lVert x_{0}}\right\rVert_{2}^{2}}., we will assume that $b=0$, thus solving the following problem:

[TABLE]

where $x\in\mathbb{R}^{n}$ is the decision variable, $P\in\mathbb{S}^{n}$, i.e. an $n\times n$ symmetric matrix, $A$ is a $m\times n$ matrix of full row rank and $q$ is an $n-$dimensional vector. Similarly to the preceding results for (T), we assume that $n-m>1$ and $r>0$.

In principle, (8) can be solved with the “eigenproblem approach” described in this section if we reduce (8) into a smaller TRS subproblem via a nullspace elimination procedure (Nocedal2006, , §15.3), obtaining

[TABLE]

where we define $Z$ to be a $n\times(n-m)$ orthonormal matrix with $\mathcal{R}(Z)=\mathcal{N}(A)$ and

[TABLE]

According to the preceding results of this section, global and local-nonglobal minimizer(s) $y^{*}$ of (9) can be computed by calculating the two rightmost eigenvectors of

[TABLE]

The respective solution(s) $x^{*}$ of (P) can then be recovered as

[TABLE]

Some disadvantages of this approach are that a nullspace basis is required and that the structure or sparsity of the problem might be destroyed.

We will present an alternative method that avoids these issues, mirroring standard approaches for solving equality constrained quadratic programs Gould2001 . To this end, note that

[TABLE]

Since Z is an orthonormal matrix, every eigenvalue-eigenvector pair $\{\mu,[\tilde{z}_{1};\;\tilde{z}_{2}]\}$ of $\tilde{M}$ corresponds to an eigenvalue-eigenvector pair $\{\mu,[z_{1};\;z_{2}]\}$ of the matrix

[TABLE]

with $Z[\tilde{z}_{1};\;\tilde{z}_{2}]=[z_{1};\;z_{2}]$. Note that (11) also has an extra eigenvalue at zero with multiplicity $2m$.

Recalling that $Z$ is an orthonormal matrix spanning $\mathcal{N}(A)$, we note that multiplication with $ZZ^{T}$ is equivalent to projection onto the nullspace of $A$. In practice, eigenvalues of (11) can therefore be calculated using the Arnoldi method on $M$ where we project the starting vector and every requested matrix-vector product with $M$ onto the nullspace of $\begin{bmatrix}A&0\\ 0&A\end{bmatrix}$.

Note that the two rightmost eigenvalues of $\tilde{M}$ correspond to the two rightmost eigenvalues of (11) if at least two eigenvalues of $\tilde{M}$ have positive real parts. We will show that we can always transform (8) so that this holds. Indeed if we shift $P$ to $P-\alpha I$ with an appropriate $\alpha\in\mathbb{R}$ so that $\tilde{P}$ become negative definite555Such an $\alpha$ can be obtained by applying e.g., the Lanczos method, or the Gershgorin circle Theorem on $P$., then $\tilde{M}$ has at least two eigenvalues with positive real part:

Lemma 2

If $\tilde{P}\in\mathbb{S}^{n-m}$ is negative definite then $\tilde{M}$ has at least two eigenvalues with positive real part, unless we trivially have $n-m=1$.

Proof

To avoid redefining the notation used in the introduction of §2 we will show the equivalent result for $(P,M)$, i.e. that $M$ has at least two eigenvalues with positive real part when $P$ is negative definite and $n>1$.

To do so, it suffices to show that $s$ has at most one root with non-positive real part. This is because $M$ has $2n\geq 4$ eigenvalues and the spectrum of $M$ matches the roots of $s$ on the left-hand complex plane (shown in Proposition 2 that we will prove on §2.A, given that $P$ is negative definite).

To show this, note that $s$ has at most one (simple) real root with nonpositive real part since

[TABLE]

is positive for any $\alpha\leq 0$ since $\lambda_{1},\dots,\lambda_{n}<0$ by assumption. To conclude the proof, it suffices to show that all the roots of $s$ with nonpositive real part are real. Indeed, note that for any $a,b\in\mathbb{R}$ we have

[TABLE]

where $i$ is the imaginary unit. A detailed derivation of the above equation is provided in Lemma 4. For any $a\leq 0$ we have $-a\geq 0>\lambda_{1},\dots,\lambda_{n}\Rightarrow 2(\lambda_{j}+a)<0$, $j=1,\dots,n$ by assumption. Thus $\text{Im}(s(a+bi))\neq 0\Rightarrow s(a+bi)\neq 0$ for any $a\leq 0$ unless $b=0$. ∎

Such a shift would not affect the (local/global) optimizers of (P) as

[TABLE]

where $\alpha\in\mathbb{R}$. Thus, we can always transform (8) so that the two rightmost eigenvalues of $\tilde{M}$ are equal to the two rightmost eigenvalues of (11).

Having computed the two rightmost eigenvalues/vectors with a projected Arnoldi method applied to $M$, global/local solutions of (P) can be obtained from the respective eigenvectors. Indeed, recall that according to (5) and (10), an appropriate (rightmost/second-rightmost) eigenvector $[\tilde{z}_{1}^{*};\;\tilde{z}_{2}^{*}]$ of (11) gives a solution $x^{*}$ to (P) as follows

[TABLE]

where $[z_{1}^{*};\;z_{2}^{*}]$ is an eigenvector of $M$. Thus, $x^{*}$ can be extracted solely by the eigenvector of $M$ unless $z_{1}^{*}=0$.

Dealing with the hard case:

In the case where $z_{1}^{*}=0$, which can appear in the hard case, one would have to calculate the minimum length solution of

[TABLE]

and then construct a global solution $x^{*}$ as $x^{*}=Z(y+\alpha\tilde{z}_{2}^{*})=Zy+\alpha z_{2}^{*}$, where $\alpha$ is a scalar such that \mathopen{}\mathclose{{}\left\lVert x^{*}}\right\rVert_{2}=r. However, the minimum length solution of (15) can be obtained by the minimum length solution of

[TABLE]

which can be calculated via projected MINRES-QLP or a projected CG algorithm (Adachi2017, , Theorem 4.3) Gould2001 .

2.2 An algorithm for computing global and local solutions of the TRS

Using the ideas presented in this section, we now present Algorithm 1, a practical algorithm for computing global and second-order sufficient local-nonglobal solutions of equality constrained trust-region subproblems. In the remainder of this section, we establish the correctness of Algorithm 1.

Algorithm 1 uses the results of this section (in particular Theorems 2.1 and 2.2, Equation (13) and the paragraph “dealing with the hard case” of §2.1) to compute the minimizers of (8). The only novel point is the way we detect the presence of the local-nonglobal minimisers, which we proceed to show to be valid.

If the TRS is in the hard case, then the algorithm should only return global minimiser(s). Indeed, this is the case as

If $z_{1}^{g}=0$, then the algorithm terminates at line 16, according to the remarks of the paragraph “Dealing with the hard case” of §2.1. 2. 2.

If $z_{1}^{g}\neq 0$, then there must exist another eigenvector associated with $\mu^{g}$ besides $[z_{1}^{g};\;z_{2}^{g}]$. This is because $\mu^{g}=-\lambda_{1}$ thus according to Proposition 2 of §2.A there exists a vector $u\in\mathcal{N}(P-\lambda_{1}I)$ that is orthogonal to $q$ which makes $(\mu^{g},[0;\;u])$ an eigenpair of $M$ due to (4).

Thus $\mu^{g}=\mu^{l}$, which implies that $\mu^{l}$ is not simple. As a result the algorithm will only return a global minimizer using (5) at line 11.

If the TRS is not in the hard case, then we must have $z_{1}^{g}\neq 0$ (see remarks of paragraph “Dealing with the hard case” of §2). Thus the algorithm will return either in lines 9 or 11, detecting correctly the presence of the local-nonglobal minimiser according to Theorem 2.2.

Finally, we emphasize that in the hard case the TRS can have an infinite number of solutions, but only one or two of them are returned by Algorithm 1.

2.A Proof of Theorem 2.2

In this subsection we provide a proof of our main theoretical result, Theorem 2.2, which shows that (T) has a second-order sufficient local-nonglobal minimizer if and only if it does not belong to the hard case and the second-rightmost eigenvalue of the matrix

[TABLE]

with counted algebraic multiplicities, is real and simple. Furthermore, we show that if such a minimizer exists, its Lagrange multiplier $\mu^{\ell}$ is equal to this second-rightmost eigenvalue.

The proof is based on the following Lemma, which follows from Martinez1994 :

Lemma 3

A KKT point $x^{\ell}$ of (T) with Lagrange multiplier $\mu^{\ell}$ is a second-order sufficient local-nonglobal minimizer iff $q\not\perp w_{1}$, $\mu^{\ell}\in(-\lambda_{1},-\lambda_{2})$ and $s^{\prime}(\mu^{\ell})>0$.

Proof

See (Martinez1994, , Theorem 3.1). ∎

Before we use Lemma 3 we will need the following two ancillary results:

Lemma 4

For any complex number $a+bi\in\mathbb{C}$ where $a,b$ are real numbers such that $a$ is in the domain of $s$ and $b\neq 0$, we have:

[TABLE]

Proof

Recall the definition of $s$:

[TABLE]

and define, for convenience, $\kappa_{j}=w_{j}^{T}q$, $j=1,\dots,n$. Thus

[TABLE]

or, equivalently, $\mathrm{Re}(s(a+bi))<s(a)$. ∎

Lemma 5

If $s(x)<0$ for some real $x$ in the domain of $s:\mathbb{C}\mapsto\mathbb{C}$ then $s$ has an equal number of poles and zeros with real parts in $[x,\infty)$.

Proof

We will prove this result with the argument principle. We propose encircling the half-space

[TABLE]

with the contour of Figure 2 and show that it maps to a contour that lies strictly in the left half plane. The suggested closed contour $C$ intersects the real axis at $x$ and at $\infty$. It consists of two segments: $C_{1}=\{\omega+\beta i\mid\beta\in\mathbb{R}\}$, i.e. a segment that is parallel to the imaginary axis, and $C_{2}$, a half-circle with infinite radius.

By the argument principle Ahlfors1966 we have $\frac{1}{2\pi i}\int_{C}\frac{s^{\prime}(z)}{s(z)}dz=Z-P$, where $Z,P$ are the number of roots and poles of $s$ in the domain enclosed by $C$. It follows that the number of poles in $S$ is equal to the number of zeros in $S$ if $\int_{C}\frac{s^{\prime}(z)}{s(z)}dz=0$, which holds if there are no poles or zeros of $s$ on $C$ and $s(C)$ does not enclose the origin. By examining (3) we conclude that $s$ has only real, finite poles. Thus, the contour $C$ does not intersect any pole of $s$. It remains to show that there exist no zeros of $s$ on $C$ and $s(C)$ does not enclose the origin. We show both by proving that every point in $s(C)=s(C_{1})\cup s(C_{2})$ has negative real part. For $s(C_{2})$, it is easy to see that, since $s(C_{2})=\{-r^{2}\}$ as $\lim_{|w|\to\infty}s(w)=-r^{2}$. Regarding $s(C_{1})$ note that $s(x)<0$, and thus Lemma 4 gives $\mathrm{Re}(s(x+\beta i))\leq s(x)<0\;\forall\beta\in\mathbb{R}$, i.e., any point in $s(C_{1})$ has negative real part. ∎

We will now connect the second-rightmost root of the secular equation $s$ with the results of Lemma 3. Note that, due to Lemma 4, if $s$ has a real root $\mu$ then it cannot also have an unreal root with real part $\mu$ and vice versa. Thus, the concept of “the second-rightmost root of $s$” is well defined as long as we assume (or are prepared to check) that this root is real.

Lemma 6

1. (i)

Suppose that $q\not\perp w_{1}$ and $\mu\in(-\lambda_{2},-\lambda_{1})$ is a simple root of $s$ with $s^{\prime}(\mu)>0$. Then, $\mu$ is the second-rightmost root of $s$ in $\mathbb{C}$. 2. (ii)

Suppose that $q\not\perp w_{1},w_{2}$, $\lambda_{1}\neq\lambda_{2}$, $\mu\in\mathbb{R}$ is the second-rightmost root of $s$ in $\mathbb{C}$ and $\mu$ is simple. Then, $\mu\geq-\lambda_{2}$.

Proof

Note that in both cases $q\not\perp w_{1}$, thus (T) does not belong to the hard case Adachi2017 . Hence $s$ has a real root in $(-\lambda_{1},\infty)$ (Theorem 2.1) which we denote by $\mu^{g}$. Moreover, $\mu^{g}$ is simple because $s^{\prime}(\mu^{g})=-2\sum_{j=1}^{n}(w_{j}^{T}q)^{2}/(\lambda_{j}+\mu^{g})^{3}$ is negative since $\mu^{g}>-\lambda_{1}>\dots>-\lambda_{n}$.

In the subsequent analysis, we will frequently refer to the poles of $s$ which are

[TABLE]

and have double multiplicity, and the following parametric set

[TABLE]

defined for any $x\in\mathbb{R}$.

We will begin with claim (i). Since $\mu<-\lambda_{1}<\mu^{g}$ we conclude that $\mu$ is not the rightmost root of $s$. In order to show that it is the second rightmost, it suffices to show that there are at most two roots of $s$ in $S(\mu)$.

Since $s(\mu)=0$ and $s^{\prime}(\mu)>0$, we have $s(\mu-\epsilon)<0$ for a sufficiently small $\epsilon>0$. According to Lemma 5, $s$ has equal number of zeros and poles in $S(\mu-\epsilon)$. By assumption $q\not\perp w_{1}$ and $\mu>-\lambda_{2}$, thus $-\lambda_{1}$ is the only (double) pole of $s$ in $S(\mu-\epsilon)$. Hence, the number of zeros in $S(\mu-\epsilon)\supset S(\mu)$ must be two. This concludes the proof of claim (i).

We will now prove claim (ii). Assume the contrary, i.e. that $\mu<-\lambda_{2}$. Consider first the case where $s^{\prime}(\mu)>0$. Since $s(\mu)=0$ we can choose a sufficiently small $\epsilon>0$ such that $s(\mu-\epsilon)<0$. Thus $s$ has equal number of zeros and poles in $S(\mu-\epsilon)$ (Lemma 5). By assumption $\mu<-\lambda_{2}$ and $q\not\perp w_{1},w_{2}$, and thus the double poles $-\lambda_{1}$ and $-\lambda_{2}$ are in $S(\mu-\epsilon)$. Thus $s$ has (at least) four zeros in $S(\mu-\epsilon)$, and since $\epsilon$ can be arbitrarily small, the same holds for $S(\mu)$.

We conclude that there should be extra roots for $s$ in $S(\mu)$, besides the simple roots $\mu^{g}$ and $\mu$. Moreover, these extra roots cannot have real part equal to $\mu$ as according to Lemma 4 $\mathrm{Re}(s(\mu+bi))<s(\mu)=0$ for any real $b\neq 0$. Thus their real parts must be in $(\mu,\infty)$. This means that $\mu$ is not the second-rightmost root of $s$, which is a contradiction.

Finally, if $s^{\prime}(\mu)<0$, then, using the same arguments as before, we conclude that there exist four roots in $S(\mu+\epsilon)$ for any sufficiently small $\epsilon>0$. Thus, there exist extra roots (besides the simple root $\mu^{g}$) with real part in $(\mu,\infty)$, which, again, contradicts the assumption that $\mu$ is the second-rightmost root of $s$. ∎

Lemma 6, in combination with Lemma 3, provides an “almost” iff relationship between the second-rightmost root of $s$ and the local-nonglobal minimizer of (T). However, the quantity of primary interest is the second-rightmost eigenvalue of $M$ instead of the second-rightmost root of $s$. The following result characterizes the spectrum of $M$ and its tight connection with the roots of $s$:

Proposition 2

The spectrum of $M$ consists of:

• •

the roots of $s$ with matched algebraic multiplicity; and

• •

*all of the eigenvalues $-\lambda_{i}$ of $-P$ that are non-simple for $-P$ or have $w_{i}$ (the *$i-$th eigenvector of $P$) orthogonal to $q$.

Furthermore, every $-\lambda_{i}$ that is in the spectrum of $M$ is non-simple for $M$.

Proof

We will prove the result by deriving a closed form expression for the characteristic polynomial of $M$. Since $-I$ and $P+\mu I$ commute, we have (Silvester2000, , Theorem 3):

[TABLE]

for all $\mu\neq-\lambda_{1},\dots-\lambda_{n}$. Thus, recalling (2)-(3) we have

[TABLE]

for all $\mu\neq-\lambda_{1},\dots-\lambda_{n}$. Thus, (20) and the characteristic polynomial of $M$ coincide in all $\mathbb{C}$ except perhaps on $2n$ points. It follows from continuity arguments that (20) is in fact the characteristic polynomial of $M$.

By examining (19) (2nd leftmost) and (20), we conclude that the eigenvalues of $M$ include all the roots of $s$ with matched algebraic multiplicity. The secular equation $s$ has $2n$ roots whenever the eigenvalues of $P$ are distinct and $w_{i}^{T}q\neq 0$ for all $i=1,\dots,n$. Otherwise, i.e. for every $i$ with $\lambda_{i}$ non-simple or $w_{i}^{T}q=0$, $M\in\mathbb{R}^{2n\times 2n}$ has extra non-simple eigenvalues at these $-\lambda_{i}$, which constitute the only differences between the spectrum of $M$ and the roots of $s$. ∎

The following result follows directly from the above Proposition and will be useful for the rest of this section.

Corollary 1

No simple eigenvalue of $M$ is in the spectrum of $-P$.

We are now ready to prove the main result of this section, Theorem 2. The proof will occupy the rest of this subsection. Note, again, that, due to Lemma 4 and Proposition 2, if $M$ has a real eigenvalue $\mu$ then it cannot also have an unreal eigenvalue with real part $\mu$ and vice versa. Thus, the concept of “the second-rightmost eigenvalue of $M$” is well defined as long as we assume (or are prepared to check) that this eigenvalue is real.

“only-if” In this case we know that $\mu^{\ell}$ is the Lagrange multiplier of a second-order sufficient local-nonglobal minimizer and we have to show that

1. (i)

(T) is not in the hard case. 2. (ii)

$\mu^{\ell}$ is a simple eigenvalue of $M$. 3. (iii)

$\mu^{\ell}$ is the second-rightmost eigenvalue of $M$.

Points (i), and (ii) follow directly from Lemma 3 and Proposition 2 (note that $q\not\perp w_{1}$ implies that (T) is not in the hard case Adachi2017 ). In order to prove (iii) we will first show that the spectrum (counted with algebraic multiplicities) of $M$ with real part larger than $-\lambda_{2}$ coincides with the roots of $s$ in the same region. This follows from Proposition 2, because $\lambda_{1}$ is simple (due to $\mu^{\ell}\in(-\lambda_{2},-\lambda_{1})$) and $q^{T}w_{1}\neq 0$. Point (iii) then follows from Lemma 6(i) which shows that $\mu^{\ell}$ is the second-rightmost root of $s$, and thus the second-rightmost eigenvalue of $M$.

“if” In this case, we know that (T) is not in the hard case and that $\mu^{\ell}$ is the second-rightmost eigenvalue of $M$ which is real and simple, and we want to show that $\mu^{\ell}$ is the Lagrange multiplier of the second-order sufficient local-nonglobal minimizer of (T).

Note that, by assumption, $\mu^{\ell}$ is in the spectrum of $M$ but, due to Corollary 1, not in that of $-P$. Thus $\mu^{\ell}$ is a root of $s$ due to Proposition 2. As a result, it suffices to show that $\mu^{\ell}$ is in $(-\lambda_{2},-\lambda_{1})$ and $s^{\prime}(\mu^{\ell})>0$ (Lemma 3).

We will first show that $\mu^{\ell}\in(-\lambda_{2},-\lambda_{1})$. Since we are not in the hard case, $M$ has a real eigenvalue in $(-\lambda_{1},\infty)$ (Theorem 2.1) which is simple and unique in $(-\lambda_{1},\infty)$ due to Proposition 2 and the fact that $s^{\prime}(\mu)=-2\sum_{j=1}^{n}(w_{j}^{T}q)^{2}/(\lambda_{j}+\mu)^{3}$ is negative in that region. Thus, $-\lambda_{1}$ is not in the spectrum of $M$ as otherwise it would be its second-rightmost eigenvalue (which by assumption is real) and it would imply that $\mu^{\ell}$ is in the spectrum of $-P$, which we have already excluded. Thus, $\mu^{\ell}\in(-\infty,-\lambda_{1})$, and $\mu^{\ell}\neq-\lambda_{2}$.

It remains to show that $\mu^{\ell}\geq-\lambda_{2}$. If $q^{T}w_{2}=0$ then, $(-\lambda_{2},[0;w_{2}])$ is an eigenpair of $M$, thus we must have $\mu^{\ell}>-\lambda_{2}$ so as to avoid $\mu^{\ell}$ (the second-rightmost eigenvalue of $M$) being in the spectrum of $-P$. Consider now the case where $q^{T}w_{2}\neq 0$. Note that $\lambda_{1}\neq\lambda_{2}$ and $q\not\perp w_{1}$ as otherwise $-\lambda_{1}$ is in the spectrum of $M$ (Proposition 2), which, according to the paragraph above is a contradiction. Thus, $\lambda_{1}\neq\lambda_{2}$ and $q\not\perp w_{1},w_{2}$ resulting in $\mu^{\ell}\geq-\lambda_{2}$ due to Lemma 6(ii).

Finally we show that $s^{\prime}(\mu^{\ell})>0$, thereby concluding the proof. Note that $s^{\prime}(\mu^{\ell})\neq 0$ since by assumption $\mu^{\ell}$ is a simple root of $s$. Assume the contrary, i.e. $s^{\prime}(\mu^{\ell})<0$. Since $s$ is convex in $(-\lambda_{2},-\lambda_{1})$ (Martinez1994, , (3.12)) and $\lim_{\mu\to-\lambda_{1}}s(\mu)=\infty$ we conclude that there must exist another root of $s$ in $(\mu^{\ell},-\lambda_{1})$. This contradicts the assumption that $\mu^{\ell}$ is the second-rightmost root of $s$.

3 An active-set algorithm for (P)

We are now in a position to present the main contribution of this paper, an active-set algorithm for solving (P). We will first present an algorithm for the special that $r_{\min}=r_{\max}:=r$ and then describe how to generalize for any $r_{\min},r_{\max}$.

3.1 Solving (P) when $r_{\min}=r_{\max}:=r$

In this section we introduce a primal active-set approach for the optimization of (P) when $r_{\min}=r_{\max}:=r$. It will be useful in the subsequent analysis to recall the Karush-Kuhn-Tucker (KKT) conditions of (P), which are

[TABLE]

As is typical from a primal active-set approach, our algorithm starts from a given feasible point of (P) and generates iterates that remain feasible for (P) and have non-increasing objective values. At every iteration we treat a subset of the inequality constraints $Ax\leq b$ as equalities. We refer to this subset as the working set $\mathcal{W}_{k}$ and define

[TABLE]

where $[\cdot]$ denotes vertical (row-wise) concatenation. To simplify the analysis, we will assume that one of the simplest and most common constraint qualification holds for every iterate of the algorithm:

Assumption 1

*Linear Independence Constraint Qualification [LICQ]

The LICQ holds for every iterate of the suggested Algorithm, i.e. $\begin{bmatrix}\bar{A}_{k}\\ \bar{x}_{k}^{T}\end{bmatrix}$ is full row rank.*

We will first give a brief, schematic description of our algorithm. At every iteration $k$ of our active-set based minimization procedure, we use Algorithm 1 to compute minimizers of the subproblem

[TABLE]

We then attempt to move towards those minimizers in that hope that we either (i) “hit” a new constraint or (ii) set $x_{k+1}$ to a (potentially local) minimizer of (SP). Unfortunately due to the nonconvexity of the problem, we might achieve neither (i) nor (ii) by simply moving towards the minimizer of (SP), as the objective of (P) is not guaranteed to be non-increasing along these moves. In that case, which we will show can happen only when (SP) has a unique minimizer, we perform a projected gradient descent followed (if necessary) by a move towards the global minimizer of (SP) that is guaranteed to achieve either (i) or (ii)666Technically, in this case $x_{k+1}$ will be a second order necessary stationary point of (SP).. Finally, at the end of each iteration, we check for termination and update the working set.

Our overall active-set procedure is shown in Algorithm 2. We now proceed in the following subsections that describe in detail each of the Algorithm’s steps. For the rest of this section, we will say that a point is “feasible” when it is in the feasible region of (P). Also, given a point $x$ we will call the projection of $\nabla f(x)$ to the nullspace of constraint normals of (SP), $\begin{bmatrix}\bar{A}_{k}\\ x^{T}\end{bmatrix}$, as “projected gradient of $f(\cdot)$ at $x$”.

3.1.1 Move towards the minimizers of (SP)

Notice that Algorithm 1 always returns a global minimizer $x^{g}$ of (SP) and possibly a second minimizer $x^{*}$.

Let us consider first the case where two distinct minimizers are returned which is treated in lines 5-6 of Algorithm 2. In this case it is always possible to either set $x_{k+1}$ to a minimizer of (SP) or “hit” a new constraint as follows. Consider the circle defined by the intersection of the sphere \mathopen{}\mathclose{{}\left\lVert x}\right\rVert_{2}=r and the plane defined by $x_{k}$ and the two returned minimizers. When (SP) is constrained in this circle, then it can be reduced to a two-dimensional TRS. As such, we will show that it has at most 2 maximizers and 2 minimizers, except perhaps when $x_{k}$ is a global minimizer of (SP) in which case we can trivially set $x_{k+1}=x_{k}$.

Proposition 3

A two dimensional TRS has at most 2 minimizers and 2 maximizers except when its objective is constant across its feasible region.

Proof

It suffices to show the result only for the minimizers, as a negation of the TRS’s objective would show the same result for its maximizers. Suppose that the TRS does not belong to the hard case. Then it has a unique global minimizer and at most one local-nonglobal minimizer Martinez1994 . If it belongs to the hard case, then only global minimizers can exist. These are intersections of the affine subspace (6) and the sphere (7). Note that this intersection is either a distinct point, two points or a circle. In the latter case, every feasible point of (T) is a global minimizer. This completes the proof. ∎

We can identify the two minimizers of this two-dimensional TRS as $x^{g}$ and $x^{*}$. It follows that at least in one of the circular arcs connecting $x_{k}$ with $x^{g}$ and $x^{*}$ the objective value is always less than $f(x_{k})$. Thus, by moving into that circular arc one will either end up with $x^{g}$ or $x^{*}$ or “hit” a new constraint.

On the other hand, if a single global minimizer $x^{g}$ was returned, which corresponds to lines 7-8 of Algorithm 2, then we consider the circle defined by the intersection of the sphere \mathopen{}\mathclose{{}\left\lVert x}\right\rVert_{2}=r and the plane defined by $x_{k}$, $x^{g}$ and the projected gradient of $f(\cdot)$ at $x_{k}$. Obviously, if $x^{g}$ is feasible for (P) then we choose $x_{k+1}=x^{g}$, but otherwise we are not guaranteed to “hit” a new constraint. This is because the associated two-dimensional TRS might possess a second minimizer, which is not necessarily a stationary point of (SP), but might be the best feasible solution on this circle.

3.1.2 Perform projected gradient descent (if necessary)

This subsection concerns the case where the procedure described in the previous subsection could neither set $x_{k+1}$ to a minimizer of SP nor “hit” a new constraint. In this case (lines 10-14 of Algorithm 2), we proceed by performing projected gradient descent (PGD) on (SP) starting from the current iterate of the active-set algorithm. This is guaranteed to converge to a KKT point $x^{s}$ with $f(x^{s})\leq f(x_{k})$ (Beck2018, , Theorem 4.5). Suppose that $x^{s}$ is not feasible for (P). Since the feasible region is a closed set, PGD exits the feasible region in finitely many steps. By stopping the projected gradient descent just before it exits the feasible region, we can find a point $x_{k+1}$ with a newly “hit” constraint.

Consider now the case where $x^{s}$ is feasible for (P). If the minimum eigenvalue $\lambda_{\min}$ of the projected Hessian of $f(\cdot)$ at $x^{s}$ is nonnegative (lines 12-14 of Algorithm 2), then we set $x_{k+1}=x^{s}$ and we proceed to the steps outlined in the next subsection. Otherwise $\lambda_{\min}<0$, and there exists a feasible sequence $\{z_{k}\}$ (w.r.t. to (SP)) converging to $x^{s}$ with an appropriate limiting direction $d$

[TABLE]

such that

[TABLE]

In practice, we can choose $d$ equal to some projected eigenvector associated with $\lambda_{\min}<0$. Consider now the circle defined by the intersection of the sphere \mathopen{}\mathclose{{}\left\lVert x}\right\rVert_{2}=r and the plane defined by $x^{s},x^{g}$ and the limiting direction $d$. As before, this can be reduced to a two-dimensional TRS that possesses at most 2 minimizers and 2 maximizers (Proposition 3). We can identify two of them: $x^{g}$ which is a global minimum and $x^{s}$ which, according to (26), is a local maximum. It follows that in at least one of the two circular arcs connecting $x^{s}$ with $x^{g}$ the objective value is always less than $f(x^{s})$. Thus, by moving into that circular arc, and since $x^{s}$ is feasible for (P) and $x^{g}$ infeasible, we can identify a suitable feasible point $x_{k+1}$ such that $f(x^{g})\leq f(x_{k+1})\leq f(x^{s})$ for which a new constraint, not in the current working set, becomes active.

3.1.3 Update the working set and check for termination

If a new constraint was identified in the steps above, then we update the working set and proceed to the next iteration (lines 15-16 of Algorithm 2). Otherwise, $x_{k+1}$ was set to a second-order necessary stationary point of (SP) and we proceed as follows (lines 17-18 of Algorithm 2). Stationarity of $x_{k+1}$ gives

[TABLE]

for some Lagrange multipliers $\mu$ and $\kappa_{i}\;\forall i\in\mathcal{W}_{k}$. Thus, $(x_{k+1},\mu,\kappa)$ satisfy the first KKT condition (21) if we define $\kappa_{i}=0\;\forall\not\in\mathcal{W}_{k}$. It follows from feasibility of $x_{k+1}$ and the definition of $\kappa$ that the third and fourth KKT conditions (23)-(24) also hold. If we also have that $\kappa_{i}\geq 0\;\forall i\in\mathcal{W}_{k}$, then $(x_{k+1},\mu,\kappa)$ is a KKT pair for (P) and we terminate the algorithm (lines 23-24 of Algorithm 2). If, on the other hand, $\kappa_{i}<0$ for some $i\in\mathcal{W}_{k}$ (lines 25-27 of Algorithm 2), then (22) is not satisfied and $x_{k+1}$ is not a local minimizer for (P). In fact, the objective $f(\cdot)$ may be decreased by dropping one of the constraints corresponding to a negative Lagrange multiplier (Nocedal2006, , Section 12.3).

3.1.4 Finite termination

We will show that Algorithm 2 terminates in finitely many outer iterations, assuming that $x_{k+1}\neq x_{k}$ for every $k=1,2.\dots$. Indeed at every iteration the algorithm either

• •

Moves to a stationary point of the current TRS subproblem; or

• •

Activates a new constraint.

Since there can be at most $m$ constraints in the working set it follows that $x_{k}$ visits a stationary point of the $k-$th TRS subproblem periodically (at least every m iterations). Furthermore, note that every TRS subproblem, resulting from a particular working set, has at most $2n$ sets of stationary points of the same objective value (Forsythe1965, , Theorem 4.1 and subsequent comments). Since there are a finite number of different working set configurations, it follows that there exist finitely many such sets. On the other hand, every constraint that is deleted from the working set has an associated negative Lagrange multiplier, thus every iterate $x_{k}$ is not a stationary point of the next iteration. Furthermore, note that if $x_{k}$ is not a stationary point for the $k+1$ subproblem then our algorithm generates $x_{k+1}$ with $f(x_{k+1})<f(x_{k})$ unless $x_{k}=x_{k+1}$ which is excluded by assumption. This means that once the algorithm visits one of these sets of stationary points with equal objective value, it can never visit it again. Hence, the algorithm terminates after finitely many iterations.

We believe that it is not possible to bound the number of iterations by a polynomial in $m,n$. Indeed, even finding an initial feasible point is NP-complete, as we show in §3.2.1.

3.1.5 Further Remarks

Similarly to active-set algorithms for quadratic programs, we can always update the working set such that $\bar{A}_{k}$ is full row rank. However, LICQ might still fail to be satisfied when the gradient of the spherical constraint lies on $\mathcal{R}(\bar{A}^{T}_{k})$. In these cases we terminate the algorithm without guarantees about local optimality.

3.2 Solving (P) for any $r_{\min},r_{\max}$

We are now ready to present an active-set algorithm that solves (P) for any $r_{\min}$ and $r_{\max}$. At every iteration of the suggested active-set algorithm, the spherical inequality constraint will either be in the working set or not. If it is, then we iterate as described in Algorithm 2; otherwise we proceed with a generic (nonconvex) quadratic programming active-set algorithm (Nocedal2006, , §16.8), Gill1991 , Gould2002 . We switch between the two algorithms when an iterate of the QP algorithm hits the spherical boundary or when the Lagrange Multiplier $\mu$ of the norm constraint in Algorithm 2 is negative and less than any other Lagrange Multiplier $\kappa_{i}$.

In the special the case where $r_{\min}=0$, the norm constraint reduces to \mathopen{}\mathclose{{}\left\lVert x}\right\rVert_{2}\leq r_{\max}. In this case, when a switch from Algorithm 2 to the QP Algorithm happens, the projected Hessian of the Lagrangian of $f(\cdot)$ at $x_{k}$ is positive semidefinite, i.e.

[TABLE]

where $Q:=[Q_{1}\;q_{2}]$ is defined by the “thin” QR decomposition of $[\bar{A}_{k}^{T}\;x_{k}]$ and $q_{2}$ is an appropriate vector. Recall that for a switch to happen, we need $\mu<0$; hence $Q^{T}PQ$ is positive definite. Thus, $Q_{1}^{T}PQ_{1}$, i.e. the projected Hessian of the next iteration that is to be handled by the generic QP algorithm, has at most one nonpositive eigenvalue (Saad2011, , Corollary 4.1). This means that even the popular, but less flexible, class of “inertia controlling” QP algorithms can be used as part of the suggested active-set algorithm for solving (P) when $r_{\min}=0$.

3.2.1 Finding an initial point

Algorithm 2 is a primal active-set algorithm and, as such, it requires an initial feasible point. Unfortunately, finding such a point is, in general, NP-complete as we formally establish at the end of this subsection. Nevertheless, we can find a feasible point for (P) with standard tools, albeit in exponential worst-case complexity, as we proceed to show. Indeed, define the following problems:

[TABLE]

[TABLE]

i.e. a convex minimization and a convex maximization problem. Denote with $x^{*}_{\min}$ and $x^{*}_{\max}$ the solutions to the convex minimization and maximization problems respectively. Then, (P) is feasible iff \mathopen{}\mathclose{{}\left\lVert x^{*}_{\min}}\right\rVert_{2}\leq r_{\max} and \mathopen{}\mathclose{{}\left\lVert x^{*}_{\max}}\right\rVert_{2}\geq r_{\min}, in which case a feasible point for (P) can be found by interpolating between $x^{*}_{\min}$ and $x^{*}_{\max}$.

The convex minimization problem above can be solved in polynomial time with e.g. interior point or first order Stellato2017 methods. On the other hand, the convex maximization problem can have exponential worst-case complexity, but it can be solved to local or even global optimality with standard commercial solvers e.g. with IBM CPLEX.

Finally, we formally show that determining if (P) is feasible is NP-complete:

Proposition 4

Determining if there is a solution to

[TABLE]

is NP-complete.

Proof

Determining if (F) has a solution is NP since we can easily check whether a candidate point $x$ is feasible for (F). Furthermore, it can be decomposed into the following two independent problems: (i) “is there a point in the polytope $Ax\leq b$ such that \mathopen{}\mathclose{{}\left\lVert x}\right\rVert_{2}\leq r_{\max}?” and (ii) “is there a point in the polytope $Ax\leq b$ such that \mathopen{}\mathclose{{}\left\lVert x}\right\rVert_{2}\geq r_{\min}?” Problem (i) can be answered in polynomial time with e.g. interior point methods. Thus determining if (F) has a solution is reducible to problem (ii) which is known to be NP-complete (Murty1987, , Problem 3). ∎

4 Applications and Experiments

In this section we present numerical results of the suggested active-set algorithms on a range of numerical optimization problems. A Julia implementation of the algorithm, along with code for the generation of the results of this section is available at:

https://github.com/oxfordcontrol/QPnorm.jl

As described in the previous sections, the algorithm is based on a TRS solver, and a general (nonconvex) Quadratic Programming solver. We implemented these dependencies as separate packages that we also release publicly as listed below.

• •

TRS.jl: A Julia package for the computation of global and local-nonglobal minimizers of

[TABLE]

essentially implementing in Julia the theoretical results of Section 2 and Adachi2017 . Available at:

https://github.com/oxfordcontrol/TRS.jl

• •

GeneralQP.jl: A Julia implementation of Gill1978 , i.e. an “inertia controlling” active-set solver for nonconvex, dense quadratic problems, with efficient and numerically stable routines for updating QR decompositions of the working set and LDLt factorizations of the projected Hessian. Available at:

https://github.com/oxfordcontrol/GeneralQP.jl

This solver is useful as a part of the suggested algorithm for solving (P) according to the remarks of §3.2.

For simplicity our implementations of the active-set algorithms are based on dense linear algebra that uses a QR factorization to compute/update a nullspace basis for every working set. Thus the presented results are limited to dense problems, except for §4.3 where the special structure of the constraints of (32) results in trivial, sparse QR factorizations of the nullspace bases.

Before presenting the results we will discuss some practical considerations regarding the suggested active-set algorithm. The eigenproblems (4) associated with the trust-region subproblems of each working set are solved with ARPACK Lehoucq1998 . In the occasional cases where ARPACK fails, a direct eigensolver is used. Finally, before solving the subproblem (SP) of every iteration, we perform a few projected gradient steps with the hope to quickly activate a new constraint.

We now proceed with the results, starting with the simplest case of dense random problems.

4.1 Random Dense Constant Norm QPs

We compared the performance of our algorithm against the state-of-the-art non-linear solver Ipopt Wachter2006 with its default parameters. We use the Julia interface of Ipopt, Ipopt.jl, which exhibits negligible interface overhead times due to the excellent interfacing capabilities of Julia with C++.

We consider a set of these randomly generated problems with varying number of variables $n$. A feasible point for each problem instance is calculated for our Algorithm as described in Section 3.2.1. The time required to compute the initial feasible points is included in the subsequent results.

Figure 3 (left) shows the execution times. We observe that our algorithm is significantly faster than Ipopt by a factor of up to 50. Both of the solvers achieve practically identical objectives (relative difference less than $10^{-9}$) in all of the problem instances, except in two cases where there is a considerable difference due to the fact that the solvers ended up converging in different local minimizers. Finally, Figure 3 (right) shows the infeasibility of the returned solution, where we observe that our solver returns solutions of significantly smaller (i.e. better) infeasibility.

4.2 Computing search directions for Sequential Quadratic Programming

Sequential Quadratic Programming (SQP) is a powerful and popular algorithm that aims to solve the problem

[TABLE]

where $g:\mathbb{R}^{n}\mapsto\mathbb{R}$ and $h:\mathbb{R}^{n}\mapsto\mathbb{R}^{m}$ are general (nonlinear) functions. At every iteration, SQP minimizes the nonlinear objective over some search directions. These search directions can, for example, be obtained by the optimization of a quadratic approximation of the original objective subject to some linear approximation of the original constraints. It is common to introduce a norm constraint \mathopen{}\mathclose{{}\left\lVert\delta x}\right\rVert_{2}\leq r that ensures that the solution will be inside a confidence or trust region of this approximation thus resulting in a problem of the form (P). This is particularly useful at points where the Hessian of the nonlinear objective is singular or indefinite as the search directions computed by the respective QPs might be unrealistically large or unbounded.

Typically, SQP algorithms solve problems simpler than (P) by e.g. not including the inequality constraints (P) in explicit form for the calculation of the search directions (Nocedal2006, , §18.2). We will show, however, that our solver is capable of computing search directions from (P) directly, without the need for these simplifications.

We demonstrate this on all the problems of the CUTEst collection Gould2015 that have linear constraints and number of variables less than $2000$. For each of these problems we consider the quadratic approximation

[TABLE]

where $x_{0}$ is the closest feasible point (in the 2-norm) to the initial point suggested by CUTEst, which we compute with GUROBI gurobi , and $\delta x:=x-x_{0}$. We discard problems where GUROBI fails to calculate an initial feasible point. Furthermore, we do not consider problems for which (29) is convex, since these problems can be solved to global optimality in polynomial time with standard solvers such as MOSEK or COSMO.jl Garstka2019 . Finally, some problems in the CUTEst collection are parametric; if the default parameters lead to a problem with more than 2000 variables then we choose a parameter that, if possible, leads to a number of variables close to, but not exceeding, 2000. Tables 4-6 list all the 63 problems considered, along with any non-default parameters.

We compare the performance of our algorithm against Ipopt with its default options. For some problems, Ipopt terminates without indication of failure but returns a low quality solution. We consider a problem as “solved” using the same criteria as Ferreau2014 . That is, we require that the overall error of the KKT conditions is less than $10^{-4}$, defined as

[TABLE]

and $\kappa,\mu$ are the Lagrange multipliers corresponding to the linear inequality and norm constraints of (29). Figure 4 shows the performance profile for the problems considered. The performance profile suggests that our algorithm significantly outperforms Ipopt on this set of problems especially in reliability, in the sense that we consistently obtain high quality solutions. Note that unlike the dense implementation of our solver, Ipopt can exploit the sparsity of the problems efficiently. Further speedups can be brought to our algorithm with a sparse implementation which might allow its use in large scale sparse problems. Moreover, as an active-set solver, our algorithm can efficiently exploit warm starting, which might be highly desirable in SQP where repeated solution of problems of the form (29) is required as part of the SQP procedure for minimizing (28). This is in contrast to Ipopt whose Interior Point nature makes warm starting very difficult to exploit.

4.3 Sparse Principal Component Analysis

Principal Component Analysis (PCA) is a standard, widely used dimensionality reduction algorithm. It has applications in numerous fields, including statistics, machine learning, bioinformatics, genetics, meteorology and others. Given a $k\times n$ data matrix $D$ that consists of $k$ points of $n$ variables, PCA suggests a few linear combinations of these variables, which we call principal vectors, that explain as much variance of the data as possible. Standard PCA amounts to the following problem

[TABLE]

where $\Sigma$ is the covariance matrix of the data. We will assume for the rest of this section that the data matrix $D$ is “centered”, i.e. that it has column-wise zero mean. We can then consider $\Sigma$ to be the empirical covariance matrix $\Sigma:=\frac{D^{T}D}{k-1}\in\mathbb{S}_{+}^{n}$. The above problem is essentially an eigenvalue problem on $\Sigma$ or a singular value problem on $D$ and, as such, can be solved with standard linear algebra tools.

In general, each principal vector is a linear combination of all the variables, i.e. typically all components of the principal vectors are non-zero. This can pose an issue of interpretability of the reduced dimensions, as it is often desirable to express the principal vector as a combination of a few variables, especially when the variables are associated with a user interpretable meaning. To alleviate this problem, sparsity has to be enforced in the original PCA problem, resulting in a new optimization problem that aims to identify a small set of variables the linear combination of which will hopefully still be able to explain a significant proportion of the variance of the data.

Unfortunately, enforcing sparsity in PCA results in a combinatorial optimization problem. Various remarkably efficient and scalable heuristics have been suggested to avoid this intractability Journee2010 , Yuan2013 . We will focus on one of the most popular heuristics that is based on a lasso type constraint originally introduced by Jolliffe2003 , resulting in the following optimization problem:

[TABLE]

where $x\in\mathbb{R}^{n}$ is the decision variable and $\gamma\geq 1$ is a parameter that controls the sparsity of the solution.

We will now show how Problem (31) can be solved with our algorithm. Note that (31) is not in the standard norm-constrained form (P) that we address. However, we show in the Appendix that it is equivalent to

[TABLE]

where $[w_{1}^{T}\;w_{2}^{T}]:=w^{T}\in\mathbb{R}^{2n}$, in the sense that (31) and (32) have the same optimal value and every minimizer of (32) $\bar{w}$ defines a minimizer $\bar{x}=\bar{w}_{1}-\bar{w}_{2}$ for (31). Also, the optimal $\bar{w}$ has the same cardinality as the respective optimal $\bar{x}=\bar{w}_{1}-\bar{w}_{2}$.

The reader might think that solving (32) with an active-set method requires a solve time that is polynomial in the number of variables, thus limiting its scalability as compared to simpler gradient based methods that are scalable but might have inferior accuracy and, potentially, weaker convergence guarantees. This is not true however, since our algorithm takes advantage of the sparsity of the iterates. Indeed, excluding the lasso constraint, the working set $\mathcal{W}_{k}$ of Algorithm 2 applied in Problem (32) can be interpreted as the set of variables $w_{i}$ that are fixed to zero. Thus every subproblem can be trivially reduced to $k$ dimensions, where $k$ is the number of variables $w_{i}$ that are not included in the current working set. This property is particularly attractive when the user is interested in a very sparse solution and has been recognized in the literature on $\ell_{1}$-penalized active-set algorithms (Bach2012, , §6).

We will demonstrate the scalability and flexibility of our algorithm by applying it in a large “bag of words” dataset from articles of the New York Times newspaper obtained by the University of California, Irvine (UCI) Machine Learning Repository777Available at http://archive.ics.uci.edu/ml/datasets/Bag+of+Words. In this dataset, every row of the associated data matrix $D$ corresponds to a word used in any of the documents and every column to a document. The entries of the matrix are the number of occurrences of a word in a document. The dataset contains $m\approx 300,000$ articles that consist of $n\approx 102,660$ words (or keywords) with $\approx 70$ million non-zeros in the dataset matrix $D$. We seek to calculate a few sparse principal vectors, each of which will hopefully correspond to a user-interpretable category, like politics, economics, the arts etc.

The usefulness of Sparse PCA on this dataset has been already demonstrated in the literature. Both Richtarik2012 and Zhang2011 have generated a set of five sparse principal vectors, each of which is a linear combination of five words. The results of Richtarik2012 are listed in Table 1. The five resulting principal vectors have distinct, user interpretable associated meanings that correspond to sports, economics, politics, education, and US foreign policy. The resulting words can be used to create user interpretable categories of the documents Zhang2011 . However, this interpretability is in general lost when we increase the cardinality of the sparse principal vectors, i.e. when we require more representative words for each category. For example, the result of applying TPower Yuan2013 (a state-of-the art algorithm for Sparse PCA) in the NYtimes dataset to generate a sparse principal vector that consists of 30 words, listed in Table 2, does not have a distinct associated meaning. For reference, TPower runs in $4.91$ seconds on a standard Laptop with an Intel i7-5557U CPU and 8GB of memory and has a resulting variance of $58.79$. Another popular algorithm, $\ell_{1}-$GPower, runs in $355.86$ seconds and results in a variance of $56.66$. The results of our active-set algorithm are better than $\ell_{1}-$GPower888Note that $\ell_{1}$ penalized algorithms (such as our algorithm and $\ell_{1}-$GPower) do not allow for the specification of the desired sparsity directly. Instead, their sparsity is controlled by a penalty parameter. We find a penalty parameter that results in 30 nonzeros via means of a binary search, while exploiting warm starting of the solvers. and comparable (slightly inferior) to TPower: our algorithm runs in $5.22$ seconds and has a resulting variance equal to $58.23$. However, our algorithm is more general, in the sense that it can solve any problem of the form (P) which, as we proceed to show, allows for variants of (31) that can generate principal vectors with more meaningful interpretation.

Indeed, improving interpretability of the results of Sparse PCA can be achieved by the incorporation of additional constraints in the underlying optimization problem. Unfortunately, most of the algorithms that are specific for sparse PCA do not allow for the presence of extra constraints. An exception is perhaps GPower of Journee2010 , which is a generic framework for minimizing a concave function $f$ over a compact subset $\mathcal{F}$. The iterates of GPower are generated as

[TABLE]

where $f^{\prime}(\cdot)$ is a subgradient of $f(\cdot)$. Journee2010 showed that for the case of Sparse PCA of (31) the above optimization problem reduces to two matrix multiplications with $D$ and $D^{T}$. However, in the presence of general linear inequality constraints, solving (33) can be at least as hard as solving a linear program (LP), which is obviously considerably more computationally demanding than two matrix multiplications.

In contrast, due to the generality of our approach, imposing any linear constraint is straightforward for our algorithm. To demonstrate this, we impose a non-negativity constraint on the elements of the principal vector to avoid cases where the principal vectors consist of elements that have opposing meanings which happens, for example, in the results of TPower listed in Table 2. Table 3 lists five sparse principal vectors obtained with this approach which take 24.37 seconds to compute on a standard Laptop with an Intel i7-5557U CPU and 8GB of memory. Note that, unlike the results of Table 2, the principal vectors of Table 3 have a clear associated meaning that can be labeled as sports, economics, politics, US foreign policy and education. These words could for example be used for organizing the documents in a user interpretable way Zhang2011 .

5 Conclusion

This paper introduced an active-set algorithm for the solution of quadratic functions subject to linear constraints and a single norm constraint. The suggested algorithm is based on repeated solutions of the TRS, for which we derived novel theoretical results regarding its local-nonglobal minimizer. The usefulness of the suggested algorithm was demonstrated in a range of real world experiments.

Acknowledgements:

We are grateful to the detailed feedbacks from an anonymous reviewer that significantly improved the quality of this document. Furthermore, we thank Nick Trefethen for providing valuable comments. This work was supported by the EPSRC AIMS CDT grant EP/L015987/1 and Schlumberger. We would also like to acknowledge the support of the NII International Internship Program, which funded the first author’s visit to NII, during which this collaboration was launched. We thank Ron Estrin for valuable discussions. We used the University of Oxford Advanced Research Computing (ARC) facility in carrying out this work. http://dx.doi.org/10.5281/zenodo.22558.

Appendix

In this Section we show that the suggested active-set algorithm of §3.2 can perform $\ell_{1}$ sparse PCA, that is, to solve the problem

[TABLE]

Theorem .1

Problems (31) and (32) have the same optimal value and every (local) minimizer $\bar{w}$ of (32) defines a (local) minimizer $\bar{x}=\bar{w}_{1}-\bar{w}_{2}$ for (31).

Proof

First note that both Problems (31) and (32) have a minimizer as their objective is smooth and their feasible set compact. For any feasible $x$ of (31), the choice

[TABLE]

where $x_{+}$ and $x_{-}$ are defined element-wise as follows

[TABLE]

is feasible for (32) since \mathopen{}\mathclose{{}\left\lVert w_{1}}\right\rVert_{2}^{2}+\mathopen{}\mathclose{{}\left\lVert w_{2}}\right\rVert_{2}^{2}=\mathopen{}\mathclose{{}\left\lVert x}\right\rVert_{2}^{2}=1 and \mathbf{1}^{T}w_{1}+\mathbf{1}^{T}w_{2}=\mathopen{}\mathclose{{}\left\lVert x}\right\rVert_{1}\leq\gamma, while $f(x)=g([w_{1};\;w_{2}])$ since

[TABLE]

Thus the optimal value of (32) is less than that of (31).

Moreover, for any minimizer $(\bar{w}_{1},\bar{w}_{2})$ of (32), choosing $\bar{x}=\bar{w}_{1}-\bar{w}_{2}$ gives $f(\bar{x})=g([\bar{w}_{1};\;\bar{w}_{2}])$ with $\bar{x}$ feasible for (31) as we proceed to show. Indeed, note that $\bar{w}_{1}^{T}\bar{w}_{2}\geq 0$ and \mathopen{}\mathclose{{}\left\lVert\bar{w}_{1}}\right\rVert_{2}^{2}+\mathopen{}\mathclose{{}\left\lVert\bar{w}_{2}}\right\rVert_{2}^{2}=1 thus

[TABLE]

and

[TABLE]

Assume that $\bar{x}$ is not a minimizer for (31). Then, there exists an $\tilde{x}$ with \mathopen{}\mathclose{{}\left\lVert\tilde{x}-\bar{x}}\right\rVert_{\infty} arbitrarily small such that $f(\bar{x})>f(\tilde{x})$. Thus $(\tilde{w}_{1},\tilde{w}_{2}):=(\tilde{x}_{+},-\tilde{x}_{-})$ is feasible for (32) with $g([\tilde{w}_{1};\;\tilde{w}_{2}])=f(\tilde{x})<f(\bar{w})=g([\bar{w}_{1};\;\bar{w}_{2}]))$ and

[TABLE]

i.e. $(\tilde{w}_{1},\tilde{w}_{2})$ is arbitrarily close to $(\bar{w}_{1},\bar{w}_{2})$. This contradicts the assumption that $(\bar{w}_{1},\bar{w}_{2})$ is a local minimizer for (32).

Finally, since the optimal value of (32) is less than the one of (31) and the global optimum $(\bar{w}_{1},\bar{w}_{2})$ of (32) defines a feasible $\bar{x}$ for (31) with $g([\bar{w}_{1};\;\bar{w}_{2}])=f(\bar{x})$ we conclude that the optimal values of (31) and (32) coincide. ∎

Furthermore, the following complementarity condition holds which implies that the sparsity of a solution $\bar{w}$ of (32) is equal to the sparsity of the solution $\bar{x}=\bar{w}_{1}-\bar{w}_{2}$ for Problem (31).

Lemma 7

Every (local) minimizer $\bar{w}$ of (32) has $\bar{w}_{1}^{T}\bar{w}_{2}=0$.

Proof

Assume the contrary, i.e. that $\bar{w}_{1}^{T}\bar{w}_{2}\neq 0$. Due to feasibility of $\bar{w}=[\bar{w}_{1};\;\bar{w}_{2}]$ we have $\bar{w}\geq 0$. Thus, there exists an $i\in\{1,\dots,n\}$ such that $\bar{w}_{(i)},\bar{w}_{(i+n)}>0$. Define $\tilde{w}$ element-wise as

[TABLE]

for $\epsilon>0$ sufficiently small so that $\tilde{w}>0$. We will show that for a sufficiently small, positive $\alpha$ the point $(1+\alpha)\tilde{w}$ is feasible and of smaller objective value than the local minimizer $\bar{w}$ while \mathopen{}\mathclose{{}\left\lVert\bar{w}-(1+\alpha)\tilde{w}}\right\rVert can become arbitrarily small for an appropriate choice of $\alpha$ and $\epsilon$. This will be a contradiction of (local) optimality of $\bar{w}$.

Indeed, note that $\mathbf{1}^{T}\tilde{w}<\mathbf{1}^{T}\bar{w}\leq\gamma$ and \mathopen{}\mathclose{{}\left\lVert\tilde{w}}\right\rVert_{2}<\mathopen{}\mathclose{{}\left\lVert\bar{w}}\right\rVert_{2}\leq 1, thus $(1+\alpha)\tilde{w}$ is feasible for any $\alpha>0$ sufficiently small. Moreover, $g(\tilde{w})=g(\bar{w})$ giving $g((1+\alpha)\tilde{w})=(1+\alpha)^{2}g(\bar{w})<g(\bar{w})$, since $g(\bar{w})<0$ as otherwise $\bar{w}$ would be a global maximizer of $\eqref{eqn:scotlass_reformulated}$. ∎

Finally, we present detailed comparison results of our algorithm against Ipopt on the problems described in Section 4.2.

[FIGURE:]

[FIGURE:]

[FIGURE:]