論文 / 著書情報
Article / Book Information
題目(和文)
Title(English)

Study on Cutting-plane Algorithms for Mixed-integer Semidefinite
Optimization

著者(和文)

小林健

Author(English)

Ken Kobayashi

出典(和文)

学位:博士（工学）,学位授与機関:東京工業大学,
報告番号:甲第11762号,
授与年月日:2022年3月26日,
学位の種別:課程博士,
審査員:中田 和秀,水野 眞治,松井 知己,塩浦 昭義,梅室 博行

Citation(English)

Degree:Doctor (Engineering),
Conferring organization: Tokyo Institute of Technology,
Report number:甲第11762号,
Conferred date:2022/3/26,
Degree Type:Course doctor,
Examiner:,,,,

学位種別(和文)

博士論文

Type(English)

Doctoral Thesis

Powered by T2R2 (Tokyo Institute Research Repository)

Study on Cutting-plane Algorithms
for Mixed-integer Semidefinite
Optimization

Ken Kobayashi
Department of Industrial Engineering and Economics
School of Engineering
Tokyo Institute of Technology

A thesis submitted in fulfillment of the requirements
for the degree of Doctor of Engineering

Supervisor: Professor Kazuhide Nakata

February 2022

ABSTRACTThis thesis focuses on mixed-integer semidefinite optimization problems, which minimize or maximize a linear objective function subject to constraints wherein a given matrix
formed from the decision variables is positive semidefinite, and some of the variables areinteger-valued. Since this problem includes nonlinearity and discreteness, various practical optimization problems can be formulated as a mixed-integer optimization problem.
However, studies on algorithms for mixed-integer semidefinite optimization problems are
relatively few, and solving large-sized mixed-integer semidefinite optimization problems
has been difficult in practice. With this background in mind, this thesis reviews some
cutting-plane algorithms for efficiently solving mixed-integer semidefinite optimization
problems. Cutting-plane algorithms are optimization techniques to handle complex constraints or an objective function. Because of the extendability and ease of implementation, the cutting-plane algorithms have the potential to efficiently solve mixed-integer
semidefinite optimization problems by exploiting the structure of the individual problems.
First, we consider the standard form of the mixed-integer semidefinite optimization
problems and propose a general-purpose cutting-plane algorithm for solving it. We also
devise a branch-and-cut algorithm to improve computational efficiency, which integrates
the cutting-plane algorithm and a branch-and-bound algorithm.
Second, we focus on the best subset selection problem for eliminating multicollinearity
from linear regression models. We give a mixed-integer semidefinite optimization formulation of the problem, in which a subset of explanatory variables is selected under anupper bound on the condition number of the correlation matrix of the selected variables.
Then, we propose a specialized cutting-plane algorithm that exploits the structure of thecondition number constraint.
Finally, we focus on portfolio selection problems formulated as a mixed-integer semidefinite optimization problem. Specifically, we consider a distributionally robust portfolio
optimization problem with limiting the number of invested assets. By exploiting the
problem structure, we develop a scalable cutting-plane algorithm for solving the modelwith the technique of positive semidefinite matrix completion. In addition, we extend thecutting-plane algorithm for the cardinality-constrained distributionally robust portfolio
optimization so as to deal with another portfolio optimization problem that minimizes
the conditional value-at-risk under a cardinality constraint. For this problem, we propose a cutting-plane method that incorporates into its subroutine another cutting-planealgorithm that efficiently minimizes the conditional value-at-risk.

Acknowledgements
First of all, I would like to express my deepest gratitude to Professor Kazuhide
Nakata. He introduced me to the study of mathematical optimization and operations research and always gave me insightful comments. The stimulating discussions with him have been invaluable experiences to me. Without his support,
I could not complete my research.
I would also like to thank Professor Yuichi Takano. His passionate and energetic attitude towards research always motivated me. Through the collaboration
with him, I learned what a researcher should be.
I am grateful to my co-authors, Mr. Ryuta Tamura, Professor Ryuhei Miyashiro,and Professor Tomomi Matsui for their valuable comments and advices.
I would also like to thank my superiors, Dr. Hirokazu Anai, Dr. Hiromich
Kobashi, Dr. Yuhei Umeda, and Dr. Hiroya Inakoshi for giving me a chance topursue a doctoral degree. They carefully coordinated my work so that I could
carry out my research comfortably.
Last but not least, I wish to thank my family for their continuous encouragement and support.
Ken Kobayashi, February 2022

Contents

1 Introduction
1.1 Background and motivation . . . . . . . . . . . . . . . . . . . . . .
1.2 Objective and results . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3 Mixed-integer semidefinite optimization problems . . . . . . . . . .
1.4 Cutting-plane algorithms . . . . . . . . . . . . . . . . . . . . . . .
1.5 Thesis overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.6 Notation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1
1
3
5
6
8
9

2 Branch-and-cut Algorithm for Mixed-integer Semidefinite Optimization
10
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.2 Problem formulation . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.3 Cutting-plane algorithm . . . . . . . . . . . . . . . . . . . . . . . . 12
2.4 Branch-and-cut algorithm . . . . . . . . . . . . . . . . . . . . . . . 14
2.5 Numerical experiments . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.5.1 Experimental design . . . . . . . . . . . . . . . . . . . . . . 15
2.5.2 Random instances . . . . . . . . . . . . . . . . . . . . . . . 16
2.5.3 Computing restricted isometry constants . . . . . . . . . . . 17
2.5.4 Robust truss topology design . . . . . . . . . . . . . . . . . 19
2.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
3 Cutting-plane Algorithm for Best Subset Selection for Eliminating Multicollinearity
25
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.2 Linear regression and multicollinearity . . . . . . . . . . . . . . . . 27
3.2.1 Linear regression . . . . . . . . . . . . . . . . . . . . . . . . 27
3.2.2 Subset selection for eliminating multicollinearity . . . . . . 27
3.3 Mixed-integer semidefinite optimization Approach . . . . . . . . . 28
3.3.1 Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.3.2 Normal-equation-based constraints . . . . . . . . . . . . . . 30
3.4 Cutting-plane algorithm . . . . . . . . . . . . . . . . . . . . . . . . 31
3.5 Numerical experiments . . . . . . . . . . . . . . . . . . . . . . . . . 33
3.5.1 Computational performance of MISDO approaches . . . . . 33
3.5.2 Computational performance of cutting-plane algorithms . . 36
3.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
4 Cutting-plane Algorithm for Cardinality-constrained Distributionally Robust Portfolio Optimization
39
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
4.2 Problem formulation . . . . . . . . . . . . . . . . . . . . . . . . . . 41
4.2.1 Cardinality constraint . . . . . . . . . . . . . . . . . . . . . 41
4.2.2 Moment-based uncertainty set . . . . . . . . . . . . . . . . 42

iv

4.3

4.4

4.5

4.6

4.2.3 Piecewise-linear utility and loss functions . . . . . . . . . .
4.2.4 Portfolio optimization model . . . . . . . . . . . . . . . . .
Cutting-plane algorithm . . . . . . . . . . . . . . . . . . . . . . . .
4.3.1 Bilevel optimization reformulation . . . . . . . . . . . . . .
4.3.2 Algorithm description . . . . . . . . . . . . . . . . . . . . .
Problem reduction of the dual lower-level problem . . . . . . . . .
4.4.1 Reduced problem formulation . . . . . . . . . . . . . . . . .
4.4.2 Equivalence of the original and reduced problems . . . . . .
Numerical experiments . . . . . . . . . . . . . . . . . . . . . . . . .
4.5.1 Computational performance . . . . . . . . . . . . . . . . . .
4.5.2 Out-of-sample investment performance . . . . . . . . . . . .
Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

42
42
43
43
45
46
46
47
50
50
53
55

5 Bilevel Cutting-plane Algorithm for Cardinality-constrained MeanCVaR Portfolio Optimization
59
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
5.2 Problem formulation . . . . . . . . . . . . . . . . . . . . . . . . . . 60
5.2.1 Conditional value-at-risk . . . . . . . . . . . . . . . . . . . . 60
5.2.2 Portfolio optimization model . . . . . . . . . . . . . . . . . 61
5.3 Cutting-plane algorithms . . . . . . . . . . . . . . . . . . . . . . . 62
5.3.1 Bilevel optimization reformulation . . . . . . . . . . . . . . 62
5.3.2 Upper-level cutting-plane algorithm . . . . . . . . . . . . . 64
5.3.3 Lower-level cutting-plane algorithm . . . . . . . . . . . . . 65
5.3.4 Efficient cut generation for the upper-level problem . . . . . 66
5.3.5 Bilevel cutting-plane algorithm . . . . . . . . . . . . . . . . 68
5.4 Numerical experiments . . . . . . . . . . . . . . . . . . . . . . . . . 70
5.4.1 Problem instances . . . . . . . . . . . . . . . . . . . . . . . 70
5.4.2 Methods for comparison . . . . . . . . . . . . . . . . . . . . 70
5.4.3 Evaluation metrics . . . . . . . . . . . . . . . . . . . . . . . 71
5.4.4 Results for various numbers of scenarios . . . . . . . . . . . 72
5.4.5 Sensitivity to hyperparameter values . . . . . . . . . . . . . 73
5.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
6 Conclusion and Prospects
80
6.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
6.2 Future directions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
A Proofs
95
A.1 Proof of Theorem 4.2 . . . . . . . . . . . . . . . . . . . . . . . . . . 95
A.2 Proof of Lemma 4.6 . . . . . . . . . . . . . . . . . . . . . . . . . . 97
A.3 Proof of Theorem 5.1 . . . . . . . . . . . . . . . . . . . . . . . . . . 99
A.4 Proof of Theorem 5.4 . . . . . . . . . . . . . . . . . . . . . . . . . . 100
B Detailed Experimental Results

102vChapter 1
Introduction
Mixed-integer semidefinite optimization problems involve minimizing or maximizing a linear objective function subject to the constraints in which a given
matrix formed from the decision variables is positive semidefinite and some ofthe variables are integer-valued. Due to the nonlinearity of the positive semidefinite constraint and discreteness of the integer constraints, this problem includes
various practical optimization problems. This thesis concerns cutting-plane algorithms for solving mixed-integer semidefinite optimization problems. In thischapter, we give an introduction to this thesis.
We describe the background and motivation of this thesis in Section 1.1. In
Section 1.2, we present the objective of this thesis and summarize our results. In
Section 1.3, we define the mixed-integer semidefinite optimization problem andreview its applications. In Section 1.4, we introduce the general framework of
the cutting-plane algorithm and give a brief review of its recent development,
including our results. In Section 1.5, we give an overview of this thesis, and in
Section 1.6, we define the notation used in this thesis.

1.1

Background and motivation

In recent years, we have faced the problems of a shrinking workforce and exhaustion of fossil fuels. Therefore, there is a growing need to efficiently use limited
management resources to realize a sustainable society. In addition, with the development of high-speed information and communication technologies, we often
have to make decisions quickly in complex and uncertain situations.
A key technology to assist such complex decisions making is mathematicaloptimization. Mathematical optimization refers to methods of finding a solutionthat minimizes or maximizes a given objective function from a set of feasible
solutions. This is one of the most important fields in applied mathematics andhas been used in various fields of science and engineering, such as operations
research, statistics, finance engineering, control theory, and chemistry.
When we formulate a certain problem as an optimization problem, we often
encounter nonlinear objective functions, constraints, and discrete variables. For
example, in statistics, we estimate the parameters of a model by minimizing aloss function. The loss function is often nonlinear, and such loss minimization
problem is formulated as a nonlinear optimization problem. Also, if the datacontains a lot of irrelevant features, it is necessary to exclude these features fromthe model when we estimate the model parameters. We can achieve such exclusion
by introducing discrete variables that indicate whether we use the features forthe estimation or not.
This thesis focuses on solving mixed-integer semidefinite optimization (MISDO)

1

problems. MISDO problems involve minimizing a linear objective function subject to constraints in which a given matrix formed from the decision variables ispositive semidefinite (PSD) and some of the variables are integer-valued. Various nonlinear constraints can be expressed as PSD constraints [155, 161], andsome important combinatorial optimization problems can be reformulated as (orapproximated by) semidefinite optimization (SDO) problems [134, 147]. Also,
integrality constraints in the MISDO problem can naturally describe the discrete
nature of practical decision-making (e.g., a yes-or-no decision) [163]. Because of
these advantages, MISDO problems have found uses in diverse fields, such as architecture [40, 171], control systems [89, 111, 153, 158], graph theory [10, 11, 135],
signal processing [61, 130], surgery planning [172], machine learning [152], andstatistical data analysis [8, 127, 151]. Therefore, solving MISDO problems efficiently is of great value in the sense of providing one of the tools for solving
various kinds of problems in these fields.
In recent years, there has been remarkable progress in developing commercial
and non-commercial solvers for mixed-integer optimization (MIO) problems andsemidefinite optimization (SDO) problems. For MIO problems, we can use powerful commercial solvers such as Gurobi and CPLEX. These solvers are based ona branch-and-bound algorithm, which explores an enumeration tree while solving a continuous relaxation problem at each node. With these powerful solvers,
large-scale mixed-integer linear optimization (MILO) problems and mixed-integer
convex quadratic optimization (MIQO) problems can be solved efficiently. For
SDO problems, on the other hand, interior-point methods are powerful algorithms to solve in practice. In fact, most of the available SDO solvers, such as
MOSEK [120], SDPA [169], SDPT3 [156], and SeDuMi [148] are based on theinterior-point method, and improving their scalability has been an important andactive area of research.
While algorithms and solvers for MIO and SDO problems have been wellstudied, there are relatively few treatments for MISDO problems. The current
MIO solvers (i.e., Gurobi and CPLEX) cannot handle semidefinite constraints
by themselves, and methods of solving MISDO problems are currently in an immature stage. A natural way of solving an MISDO problem is to combine abranch-and-bound algorithm for MIO problems and an interior-point algorithm
for SDO problems. Recently, Gally et al. [62] developed a general-purpose MISDO
solver named SCIP-SDP; it is regarded as the most promising solver at present.
SCIP-SDP combines the branch-and-bound framework of SCIP [1] with SDO
solvers that use interior-point methods. Also, most of the specialized algorithms
for specific MISDO applications are based on the branch-and-bound algorithm[8, 10, 11, 40, 130, 135, 139, 171]. However, there is a serious drawback tothis approach, wherein if we just combine a branch-and-bound algorithm and aninterior-point method, we cannot implement a warm-starting strategy for efficiently solving a series of SDO problems because the interior-point methods do
not receive a given initial solution. Indeed, Gleixner et al. [63] implemented several warm-starting techniques [39, 68, 81, 143, 146] in SCIP-SDP, but their effects
are marginal or sometimes negative in some problem instances. Therefore, thesize of the MISDO problems that we can solve still has been limited.
Against this background, this thesis focuses on cutting-plane algorithms asalternatives to the branch-and-cut algorithm and reviews the cutting-plane algorithms for solving MISDO problems. Cutting-plane algorithms are optimization
techniques to handle complex constraints or an objective function. A cuttingplane algorithm first relaxes the constraints or the objective function and searches
for an optimal solution by solving the resultant relaxed problems while adding
2

linear constraints called cutting planes or valid inequalities. The cutting-planealgorithm was devised by Gomory [67] to solve MILO problems, and it was later
extended by Kelly [86] to handle convex optimization problems. The cuttingplane algorithm has the advantages of extensibility and ease of implementation.
Because of its tractability, the cutting-plane algorithm has been actively studied
in nonsmooth convex optimization (see, for example, [64, 73, 100, 110]). In addition, if the relaxation problem is easy to solve, we can improve the performance ofthe cutting-plane algorithm by using this property. Thus, if an MISDO problem
can be relaxed into a tractable one, a cutting-plane algorithm will be a promising
way of solving it efficiently.

1.2

Objective and results

From this perspective, we aim to design efficient cutting-plane algorithms to solve
MISDO problems. Specifically, we firstly propose a general-purpose cutting-planealgorithm for solving the standard form of MISDO problems. This method allows us to handle the general MISDO problem with existing state-of-the-art MIO
solvers. Secondly, we devise specialized cutting-plane algorithms for some important MISDO problems. In statistics and financial engineering, we sometimes
encounter large-sized MISDO problem instances due to the increasing size ofthe data sets in these fields, and such problems are difficult to solve even withgeneral-purpose algorithms. To handle large-sized problem instances, we design specialized cutting-plane algorithms that exploit the individual structures ofthese problems. In particular, this thesis focuses on a variable selection and portfolio selection problem formulated as MISDO problems and proposes specialized
cutting-plane algorithms for solving them efficiently.
The results of this thesis are summarized as follows:
Branch-and-cut algorithm for mixed-integer semidefinite optimization(Chapter 2)First, we study general-purpose cutting-plane algorithms for solving the standard form of MISDO problems. We develop a cutting-plane algorithm for solving
MISDO problems by extending the existing algorithm for SDO problems proposedby Konno et al. [100]. In the proposed algorithm, we relax the PSD constraint
and solve the relaxed MILO problem repeatedly while adding a cutting plane tothe constraints at each iteration so as to satisfy the PSD constraint. We prove theconvergence properties of the algorithm. In addition, to speed up the computation, we devise a branch-and-cut algorithm, where cutting planes are dynamically
added during a branch-and-bound procedure. Our experimental results demonstrate that, for many problem instances, our branch-and-cut algorithm delivered
superior performance compared with other general-purpose MISDO solvers interms of computational efficiency and stability.
Cutting-plane algorithm for best subset selection for eliminating multicollinearity (Chapter 3)
Second, we consider a subset selection of explanatory variables for eliminating
multicollinearity from linear regression models. Specifically, we select the best
subset of the variables subject to an upper bound on the condition number of thecorrelation matrix of selected variables. We first formulate the best subset selection under the condition number constraint as an MISDO problem. Since the3current general-purpose MISDO solver only handled very small-sized instances,which we will show in the corresponding chapter, we also developed a specialized
cutting-plane algorithm. To approximate the condition number constraint, thiscutting-plane algorithm iteratively appends cutting planes to the relaxed MIQO
problem. Here, we also propose to generate strong cutting planes by effectively
using heuristic search to improve the computational efficiency of our algorithm.
Computational results demonstrate that subset selection with our MISDO formulation succeeds when the number of candidate explanatory variables is small.
Additionally, our specialized cutting-plane algorithm frequently provides solutions of better quality than those obtained by local search algorithms for subset
selection.
Cutting-plane algorithm for cardinality-constrained distributionally robust portfolio optimization (Chapter 4)
Third, we focus on portfolio selection problems. In particular, we study a distributionally robust portfolio optimization model with a cardinality constraint, which
limits the number of invested assets. This problem is formulated as a MISDO
problem, and thus, solving it exactly is computationally challenging when thenumber of investable assets is large. To overcome this issue, we propose a specialized cutting-plane algorithm to solve the cardinality-constrained distributionally robust optimization problem. We first reformulate the problem as a bilevel
optimization problem and design a cutting-plane algorithm for solving the upperlevel problem. To generate cutting planes efficiently, we apply the technique of
positive semidefinite matrix completion to the lower-level problem and show that
we can calculate subgradients efficiently regardless of the number of investableassets. The numerical experiments demonstrate that our cutting-plane algorithm
was very effective compared with the existing MISDO solver, especially when thenumber of investable assets is large. In addition, the out-of-sample investment
performances given by the cardinality-constrained distributionally robust model
were better than those of the conventional mean-variance model.
Bilevel cutting-plane algorithm for cardinality-constrained mean-CVaR
optimization (Chapter 5)
Finally, we focus on cardinality-constrained conditional value-at-risk (CVaR)
minimization problems and extend the cutting-plane algorithm for the distributionallyrobust portfolio optimization problem. While the cardinality-constrained CVaR
minimization problem is formulated as an MILO problem, which is a special caseof MISDO problem and can be handled by state-of-the-art MIO solvers, its problem size depends not only on the number of investable assets but also on thenumber of asset return scenarios. Thus, the computational efficiency decreases
when the number of scenarios is large. To overcome this challenge, we propose
a specialized cutting-plane algorithm named the bilevel cutting-plane algorithm
for exactly solving the cardinality-constrained mean-CVaR portfolio optimization problem. In our proposal, we extend the cutting-plane algorithm discussed
in Chapter 4 so that it incorporates another cutting-plane algorithm that efficiently minimizes CVaR. Numerical experiments demonstrate that, comparedwith other MIO approaches, our algorithm can provide optimal solutions to large
problem instances faster.

4

1.3

Mixed-integer semidefinite optimization problems

In this section, we define MISDO problems as an extension of SDO problems andgive a brief overview of their applications.
We begin by introducing linear optimization (LO) problems. An LO problem
is the most basic optimization problem. It aims to minimize or maximize a linear
function subject to finite linear constraints. The standard form of LO problems
is formulated as:(LO-P)

minimize c> x

maximize b> y

(LO-D)y∈RM

x∈RN

subject to A> y ≤ c,subject to Ax = b,(1.1)

x ≥ 0,
where (LO-P) and (LO-D) denote the primal and dual LO problem, respectively;
A ∈ RM ×N , b ∈ RM , and c ∈ RN are given constants; (x, y) is a decision variable
and x ≥ 0 represents that all the entries of x are nonnegative.
SDO is an extension of LO problems to an optimization problem over the
space of N × N real symmetric matrices. Here, we denote the set of N × N real
symmetric matrices by S N , and define the standard inner product between two
matrices X = (Xnm ), Y = (Ynm ) ∈ S N as
X • Y = Tr(XY ) =

N X
N
X

Xnm Ynm .

n=1 m=1

Let A1 , A2 , . . . , AM ∈ S N , b ∈ RM , and C ∈ S N be given constants. Then, thestandard primal and dual forms of SDO problems are given by(SDO-P)(SDO-D)

minimize C • X
X∈S Nmaximize b> yy∈RMsubject to A> (y)  C,

subject to A(X) = b,

(1.2)

X  O,
where (SDO-P) and (SDO-D) denote the primal and dual SDO problem, respectively; (X, y) is a decision variable, and X  O denotes that X is positive
semidefinite i.e., v > Xv ≥ 0 for all v ∈ RN ; A : S N → RM is a linear mapping
defined by A(X) := (A1 • X,P
A2 • X, . . . , AM • X)> and A> : RM → S N is its>
adjoint defined by A (y) := M
m=1 ym Am .
It is easy to see that the SDO problem (1.2) includes the LO problem (1.1)
as a special case. For a vector x ∈ RN , let Diag(x) ∈ S N be the diagonal matrix
whose diagonal elements are x. Then, the nonnegativity of x is equivalent to thepositive semidefiniteness of Diag(x) as follows:x ≥ O ⇐⇒ Diag(x)  O.
In addition, the inner product of two vectors x, y ∈ RN is expressed as that ofDiag(x) and Diag(y):x> y = Diag(x) • Diag(y).
Thus, we can reformulate the LO problem (1.1) as an SDP problem (1.2).
The MISDO problem is defined as an SDO problem with integer constraints
and its dual form is given by:(MISDO-D)

maximize b> yy∈RMsubject to A> (y)  C,y ∈ {0, 1}Mb × RMc .
5

(1.3)

where y := (ym ) ∈ RM ; Mb and Mc are, respectively, the numbers of binary
and continuous variables satisfying M = Mb + Mc . Alternatively, we may also
consider MISDO problems in the primal form: .
(MISDO-P)

minimize C • X
X∈S N

subject to A(X) = b,
X  O,
Xnm ∈ {0, 1}, (∀(n, m) ∈ B),
where X := (Xnm ) ∈ S N , and the index set B ⊆ {1, 2, . . . , N } × {1, 2, . . . , N }imposes binary-valued entries in X.
Since the MISDO problem involves nonlinearity and discreteness because of
the semidefinite and binary constraints, it includes various practical optimization
problems. One of the main fields where MISDO problems appear is combinatorial optimization. For example, graph partitioning problems are formulated as an
MISDO problem [10, 11]. Also, Rendl and Rinaldi [135] proposed an algorithm
for exactly solving max-cut problems with an MISDO formulation. In structuraldesign, Yonekura and Kannno [171] considered a robust truss topology optimization and formulated this problem as an MISDO problem. In control theory,
Joshi and Boyd [89] focused on a sensor selection problem that minimizes theerror in estimated parameters and showed that this problem can be expressed as
an MISDO problem. In data mining, MISDO problems are used to obtain an
optimal clustering in the sense of a given metric [8, 127]. In statistics, Tamura
et al. [151] examined the best subset selection for eliminating multicollinearity
from a linear regression model and proposed a MISDO formulation with a condition number constraint. Moreover, MISDO problems have recently been used inrobust optimization. In particular, Zhang et al. [172] focused on a robust optimization for allocating surgery blocks and formulated the problem as an MISDO
problem.

1.4

Cutting-plane algorithms

In this section, we explain the general framework of cutting-plane algorithms
for solving convex optimization problems and give a brief overview of recent
developments on cutting-plane algorithms, including our results.
The cutting-plane algorithm was first introduced by Gomory [67] for solving
MILO problems, and it was extended by Kelly [86] to general convex optimization
problems. Here, we explain the general framework proposed by Kelly [86]. Let
us consider the following convex optimization problem:
minimize
x∈RN

f (x)

subject tox ∈ F,(1.5)

where F ⊆ RN is a convex set, f : RN → R is a convex function, and we assume
that the above problem has an optimal solution.
The main ingredient of cutting-plane algorithms is the hyperplane separation
theorem [36, 82, 136]. For a convex set F ⊆ RN and a point x̂ ∈/ F, the hyperplane
separation theorem ensures the existence of a 6= 0 and b such that a> x̂ > b anda> x ≤ b for all x ∈ F. We call this inequality separating F and x̂ a cutting
plane. Also, we call the inequality satisfying a> x ≤ b for all x ∈ F valid or avalid inequality for F.
The cutting-plane algorithm solves Problem (1.5) by refining convex polyhedrons that contain F while cutting off infeasible solutions with the hyperplane6

separation theorem. Here, let us define F1 to be the initial feasible region that
contains F. At the tth iteration (t ≥ 1), the cutting-plane algorithm solves thefollowing optimization problem:
minimize
x∈RN

f (x)

subject tox ∈ Ft ,(1.6)

where Ft is a relaxed feasible region at the tth iteration satisfying Ft ⊆ F1 . Let
xt be an optimal solution to Problem (1.6). Then, we check the feasibility of
xt . If xt ∈ F, xt is optimal for the original problem (1.5), and we terminate thealgorithm and output xt as an optimal solution. If xt ∈/ F, we find a cutting
plane a> x ≤ b separating xt from Ft and update the relaxed feasible region byadding the cutting plane to the constraints as follows:
Ft+1 ← Ft ∩ {x ∈ RN | a> x ≤ b},(1.7)

After updating the feasible region, we set t ← t + 1 and solve Problem (1.6)
again. We repeat this procedure until xt ∈ F (possibly up to some tolerance).
We summarize the general scheme of the cutting-plane algorithm for Problem(1.5) by Algorithm 1.1.
Algorithm 1.1 General scheme of cutting-plane algorithms for solving Problem (1.5)
Step 0 (Initialization) Set t ← 1 and the initial feasible region F1 ⊇ F.
Step 1 (Relaxed Problem) Solve Problem (1.6). Let xt be an optimal solution.
Step 2 (Cut Generation) If xt ∈ F, terminate the algorithm. Otherwise,find a cutting plane a> x ≤ b separating xt from F.
Step 3 (Update Feasible Region) Update the feasible region as in Eq.
(1.7).
Step 4 Set t ← t + 1 and return to Step 1.
Cutting-plane algorithms have the following advantages:• They do not require the objective and constraint functions to be differentiable. As long as we can determine whether a solution belongs to F andgenerate cutting planes for F, we can apply the cutting-plane algorithm tothe problem to be solved.
• If the relaxed problem has a nice solvable structure, we can use it to solve
the relaxed problem and accelerate the entire computation. For example,
if state-of-the-art solvers can handle the relaxed problem, we can use themin the cutting-plane algorithm.
• Cutting-plane algorithms have good extendability. If we can generate acutting plane for F, Algorithm 1.1 can also be applied to the following
mixed-integer optimization problem:
minimize
x∈RN

f (x)

subject tox ∈ F ∩ ({0, 1}p × RN −p ),(1.8)

where p denotes the number of binary variables, and {0, 1}p is the set of pdimensional vectors whose all entries are 0 or 1.
7

Because of these advantages, cutting-plane algorithms have come to be applied to various optimization problems. Regrettably, we can not describe allof the studies on the cutting-plan algorithms, but we will mention the resultsthat are strongly relevant to this thesis. Konno et al. [100] proposed a cuttingplane algorithm for solving SDO problems. A similar algorithm was proposedby Krishnan et al. [102]. Recently, an MISDO solver named CUTSDP, which is
implemented in YALMIP [106], was made available; it employs the outer approximation algorithm, which is an extension of the cutting-plane algorithm. Coeyet al. [45] developed a general-purpose solver named Pajarito for mixed-integer
conic optimization problems. Kobayashi and Takano [94] extended the cuttingplane algorithm for SDO problems and proposed a branch-and-cut algorithm forsolving MISDO problems.
For specific MIO problems, Bertsimas and King [32] considered subset selection in linear regression, which is a classical problem in statistics, and employed
a cutting-plane strategy. In financial engineering, Bertsimas and Cory-Wrightproposed a cutting-plane algorithm for cardinality-constrained mean-variance
portfolio optimization [24]. This cutting-plane algorithm was extended to coversparse principal component analysis [22] and covariance selection problems [33].
Kobayashi et al. [96] considered a cardinality-constrained distributionally robust
portfolio optimization, which is formulated as an MISDO problem, and devised
an efficient cutting-plane algorithm for solving it with the technique of positive
semidefinite matrix completion [60, 121]. Künzi-Bay and Mayer [103] proposeda cutting-plane algorithm for minimizing CVaR, which is a risk measure in finance, and variants of this algorithm were developed [4, 80, 149]. Kobayashi
et al. [95] proposed an algorithm named the bilevel cutting-plane algorithm forsolving cardinality-constrained CVaR minimization problems.
Finally, we mention three points that should be kept in mind when we design
an efficient cutting-plane algorithm. First, the relaxed problem (1.6), which issolved at each iteration, should be tractable. Second, cutting planes should be
generated with a low computational cost. The last point is the strength of cutting>
planes. For two cutting planes a>
1 x ≤ b1 and a2 x ≤ b2 for F, we call the former
is stronger than the latter if the following condition holds:N>{x ∈ RN | a>1 x ≤ b1 } ⊆ {x ∈ R | a2 x ≤ b2 }.

In this case, the stronger cutting plane cuts off more infeasible solutions; so it
is preferred when we add a cutting plane to the constraints. Thus, when we
generate a cutting plane, it should be strong.

1.5

Thesis overview

The remainder of this thesis is structured as follows. In Chapter 2, we study
the standard form of MISDO problems and devise general-purpose cutting-planealgorithms for solving the problem. In Chapter 3, we consider the best subset
selection problem for eliminating multicollinearity from linear regression models. We formulate this problem as an MISDO problem and propose a specialized cutting-plane algorithm that exploits the structure of the condition number
constraint. In Chapters 4 and 5, we focus on cardinality-constrained portfolio
optimization problems. We first focus on distributionally robust portfolio optimization in Chapter 4. We devise a scalable cutting-plane algorithm where we
generate cutting planes efficiently by using the technique of positive semidefinite
matrix completion. In Chapter 5, we deal with the cardinality-constrained CVaR
8minimization problem and extend the cutting-plane algorithm discussed in Chapter 4 so that it can handle a large number of asset return scenarios. Chapter 6 isdevoted to concluding remarks. Note that the details of the existing research oneach problem are given in the corresponding chapter.

1.6

Notation

We define the set of real numbers, integers, and positive integers as R, Z, andN, respectively. For N ∈ N, [N ] represents the set {1, 2, . . . , N }. We denote by
x := (xn ) ∈ RN an N -dimensional column vector of real numbers. The zero andall-one vectors of appropriate size are written as 0 and 1, respectively. If all theentries of x ∈ RN are nonnegative, we write x ≥ 0. For x, y ∈ RN , we use x ≥ y
to mean x − y ≥ 0. The transpose of a matrix x ∈P
RN is denoted by x> . The
N>
standard inner product on R is given by x y = n∈[N ] xn yn for x, y ∈ RN .
P
1/p
p
The `p -norm of the vector x is denoted by kxkp :=|x|.
nn∈[N ]We denote by X := (Xnm ) ∈ RN ×M an N × M real matrix. The zero andidentity matrices of appropriate sizes are written as O and I, respectively. The
trace Tr(·) is the sum of the diagonal elements of a square matrix. The diag(·)
operator extracts the main diagonal from the matrix as a vector, and the Diag(·)
operator maps the vector to the diagonal matrix. The transpose of a matrix
X ∈ RN ×M is denoted by X > .
The set of all symmetric N × N real matrices
by S N . The standard
P is denotedPN
inner product on S is given by X • Y = m∈[N ] n∈[N ] Xmn Ymn for X, Y ∈S N . The minimum and maximum eigenvalues of matrix X ∈ S N are denoted by
λmin (X) and λmax (X), respectively.
A matrix X ∈ S N is positive semidefinite (PSD) if
v > Xv ≥ 0
for all v ∈ RN . We write X  O if the matrix X is PSD. For X, Y ∈ S N , weuse X  Y to mean X − Y  O. A matrix X ∈ S N is positive definite if
v > Xv > 0
for all v 6= 0. We write X  0 if the matrix X is positive definite. For asymmetric matrix X ∈ S N , its positive semidefiniteness and nonnegativity of its
eigenvalues are equivalent as follows:
X  O ⇐⇒ λmin (X) ≥ 0.

9Chapter 2
Branch-and-cut Algorithm for Mixed-integer
Semidefinite Optimization
In this chapter, we study the general-purpose cutting-plane algorithms for solving MISDO problems. In particular, we focus on the standard form of MISDO
problems and develop a cutting-plane algorithm for solving it by extending theexisting algorithm for solving SDO problems proposed by Konno et al. [100]. We
prove the convergence properties of the algorithm. In addition, to speed up thecomputation, we devise a branch-and-cut algorithm, in which cutting-planes are
dynamically added during a branch-and-bound procedure. The content of thischapter is included in Kobayashi et al. [94].
We give an introduction to this study and summarize our contribution in
Section 2.1. In Section 2.2, we show the definition of the MISDO problems again.
In Section 2.3, we present our cutting-plane algorithm and analyze its convergenceproperties. In Section 2.4, we describe our branch-and-cut algorithm for solving
MISDO problems. In Section 2.5, we report the computational results of tests of
our algorithms, and in Section 2.6, we conclude with a brief summary.

2.1

Introduction

In this chapter, we focus on the standard form of MISDO problems. As we
mentioned in Chapter 1, the standard MISDO problem includes various kinds of
optimization problems appearing in diverse fields, such as combinatorial optimization [134, 147], architecture [40, 171], control systems [89, 111, 153, 158], graphtheory [10, 11, 135], signal processing [61, 130], surgery planning [172], and statistical data analysis [8, 127, 151]. Thus, developing an efficient general-purpose
algorithm for solving the standard MISDO formulation is of great significance interms of providing a framework for handling these various optimization problems
in a unified manner.
A common way to solve MISDO problems is to use a branch-and-bound(B&B) algorithm. Gally et al. [62] developed a general-purpose MISDO solver,
SCIP-SDP, by combining the B&B framework of SCIP [1] with SDO solvers thatuse interior-point methods. Also, the existing specialized algorithms for specific
MISDO applications are mainly based on the B&B algorithm [8, 10, 11, 40, 130,
135, 139, 171]. These B&B algorithms involve solving a continuous relaxation(i.e., SDO) problem at each node of the enumeration tree. In this case, however,
we cannot implement a warm-starting strategy for efficiently solving a series of
SDO problems because the interior-point methods do not receive a given initialsolution. Moreover, the progress of the B&B algorithm is sometimes disrupted
and returns an incorrect solution due to numerical instability, especially when
Slater’s condition does not hold for some of the SDO problems.
10

As an alternative to the B&B framework, we focus on the cutting-plane algorithm [86, 162]. Specifically, we shall extend the cutting-plane algorithms [97, 102]that were originally developed for solving SDO problems. These algorithms canreadily be extended to solving MISDO problems by removing the PSD constraint
and repeatedly solving the resultant MILO problem in which a cutting plane is
imposed at each iteration to exclude infeasible solutions for the removed PSD
constraint. This approach has the advantage of being able to use state-of-the-art
optimization software (e.g., Gurobi or CPLEX) when solving MILO problems.
However, to solve an MISDO problem with a sufficient degree of accuracy, wemust deal with an exponentially large number of MILO problems [38]. A similar cutting-plane algorithm for solving MISDO problems has previously been
described [108] and implemented in YALMIP [106]. To the best of our knowledge, however, the theoretical properties of the algorithm have not been fully
described, and no detailed computational results have been reported.
The purpose of this chapter is to describe an effective computational framework for solving MISDO problems. First, we formulate a cutting-plane algorithm for solving MISDO problems and prove its convergence properties. We
then devise a high-performance branch-and-cut (B&C) algorithm on the basisof this cutting-plane algorithm. Specifically, we start by solving an MILO problem in which the PSD constraint has been relaxed, and then we include cutting
planes for the constraint dynamically during the B&B procedure. We implement
this B&C algorithm through the use of a callback function in the Gurobi Optimizer. We compare the computational performance of our algorithms with thatof SCIP-SDP [62] and CUTSDP for three types of MISDO problems: random instances, computing restricted isometry constants [61], and robust truss topology
design [171].
The main advantages of our B&C algorithm can be summarized as follows:• Our B&C algorithm solves an MILO problem subject to a series of cutting
planes. In contrast to the existing B&B algorithms that handle a series of
SDO problems, a warm-starting strategy using the dual simplex method
makes the B&C computation much faster.
• Our B&C algorithm can be implemented using sophisticated MILO software, such as Gurobi or CPLEX, and it deals with no SDO problems that
may cause numerical instability. Hence, computation using our algorithm
is very stable.
• Our B&C algorithm adds cutting planes dynamically during the MILO
computation. As a result, our algorithm needs to execute the B&B procedure only once, in contrast with the cutting-plane algorithm, which repeats
the B&B procedure.

2.2

Problem formulation

In this chapter, we focus on solving the dual standard form of MISDO problems.
Recall that the dual standard dual form of MISDO problems is expressed as11

follows:
maximizeyb> y(2.1a)

subject to A> (y) :=X

Am ym  C,

(2.1b)m∈[M ]y ∈ {0, 1}

Mb

× R Mc ,(2.1c)

where b := (bm ) ∈ RM and Am ∈ S N (m ∈ [M ]) are given constants, and y isa vector of binary and continuous decision variables. Note that M = Mb + Mc
holds, where Mb and Mc are, respectively, the numbers of binary and continuous
decision variables. Note that we may also consider MISDO problems in the primal
form in Eq. (1.4). In the following sections, we shall deal with the problem (2.1),
but our algorithms can readily be applied to the primal standard form.
We assume that the feasible region of the problem (2.1) is bounded even if
the PSD constraint (2.1b) is removed. Since
C − A> (y)  O ⇒ diag(C − A> (y)) ≥ 0,
we define the relaxed feasible region
Y := {y ∈ {0, 1}Mb × RMc | diag(C − A> (y)) ≥ 0}.
Assumption 2.1. The relaxed feasible region Y is bounded.
Let y ? be an optimal solution to the problem (2.1). Note that Assumption 2.1
is fulfilled if the following constraint is incorporated into the definition of thefeasible region:−U ≤ ym ≤ U, Mb + 1 ≤ ∀i ≤ M,
where U is a sufficiently large number such that this constraint is satisfied by y ? .

2.3

Cutting-plane algorithm

Our cutting-plane algorithm for solving MISDO problems is an extension of the
cutting-plane algorithms [97, 102] for solving SDO problems. It can also be
regarded as an improved version of the extended cutting-plane algorithm [162] inthe sense that our algorithm makes it possible to handle PSD constraints.
Our algorithm starts with defining the initial feasible region as F1 := Y
where the PSD constraint (2.1b) has been removed. At tth iteration (t ≥ 1), ouralgorithm solves the following relaxed MILO problem:
maximizeyb> ysubject toy ∈ Ft ,(2.2)

where Ft is a relaxed feasible region at the tth iteration such that Ft ⊆ F1 . If
Ft = ∅, the problem (2.1) is infeasible because its feasible region is contained in
Y. Otherwise, since the feasible region Ft is bounded from Assumption 2.1 andFt ⊆ F1 , there exists an optimal solution yt to the problem (2.2).
We next check whether the PSD constraint (2.1b) is satisfied by the relaxedsolution yt . If λmin (C−A> (yt )) ≥ 0, then C−A> (yt ) is PSD and yt is an optimal
solution to the original problem (2.1), so we terminate the algorithm. Otherwise,
let dt ∈ RN be the normalized eigenvector of C − A> (yt ) corresponding to theminimum eigenvalue. It follows that>>d>t (C − A (yt ))dt = λmin (C − A (yt )) < 0.

12

Accordingly, we can exclude the solution yt from the feasible region by incorporating the constraint
>
d>
(2.3)
t (C − A (y))dt ≥ 0.
Note that this constraint is a valid inequality for the PSD constraint (2.1b) becauseDE
C − A> (y)  O ⇐⇒ d> (C − A> (y))d ≥ 0, ∀d ∈ Rn .
After constraint (2.3) is included, we set t ← t + 1 and solve the MILO
problem (2.2) again. We repeat this process (i.e., solving the MILO problem andadding a cutting plane) until C − A> (yt ) becomes nearly PSD:
λmin (C − A> (yt )) ≥ −ε,(2.4)

where ε is a sufficiently small non-negative number representing a tolerance forfeasibility. In this case, constraints (2.1c) and (2.4) are satisfied by yt , and thevalue of the objective function is not less than the maximum (i.e., b> yt ≥ b> y ? )because yt solves the relaxed problem (2.2). Here, we call such a solution optimal
within the feasibility tolerance ε. Our cutting-plane algorithm is summarized
as Algorithm 2.1; if the algorithm terminates, it proves the infeasibility of the
problem (2.1) or yields an optimal solution with the feasibility tolerance ε.
Algorithm 2.1 Cutting-plane algorithm for solving MISDO Problems
Step 0 (Initialization) Let ε ≥ 0 be a tolerance for feasibility. Define theinitial feasible region as F1 := Y. Set t ← 1.
Step 1 (Relaxed Problem) Solve Problem (2.2). Let yt be an optimal solution.
Step 2 (Termination Condition)(a) If Ft = ∅, then the problem (2.1) is infeasible.
(b) If λmin (C − A> (yt )) ≥ −ε, then yt is an optimal solution within
the feasibility tolerance ε.
Step 3 (Cut Generation) Update the feasible region,
Ft+1 ← Ft ∩ y ∈ {0, 1}Mb × RMc | Eq.(2.3) .
Step 4 Set t ← t + 1 and return to Step 1.
Following the approach used in the case of solving SDO problems [97, 98], weprove convergence of the algorithm.
Theorem 2.2. For any ε > 0, Algorithm 2.1 terminates in a finite number of
iterations.
Proof. Suppose that the algorithm does not terminate. The infinite sequence{(dt , yt )} generated by the algorithm is bounded because kdt k = 1 and yt ∈ Y
for all t. Therefore, we can choose a subsequence {(dt , yt ) | t ∈ T } that converges¯ ȳ).
to an accumulation point (d,Since the termination condition is never satisfied, we have that for all t ∈ T ,>>d>t (C − A (yt ))dt = λmin (C − A (yt )) < −ε.

13

It follows thatlimt→∞ (t∈T )(dt )> (C − A> (yt ))dt = d¯> (C − A> (ȳ))d¯ ≤ −ε.

(2.5)

For each t < `, y` satisfies the tth cutting plane>d>t (C − A (y` ))dt ≥ 0.

It follows that
lim

lim

t→∞ (t∈T ) `→∞ (`∈T )

>
>
¯>
¯
d>
t (C − A (y` ))dt = d (C − A (ȳ))d ≥ 0.

(2.6)

Eqs. (2.5) and (2.6) imply that ε = 0, which contradicts the hypothesis.
Theorem 2.3. Suppose that ε = 0. Even if Algorithm 2.1 does not terminate,
every accumulation point of the sequence of solutions {yt } is an optimal solution
to the problem (2.1a)–(2.1c).
Proof. Suppose that {(dt , yt ) | t ∈ T } is a sequence that converges to an accu¯ ȳ). Note that λmin (C −A> (yt )) converges to λmin (C −A> (ȳ))
mulation point (d,as t → ∞ (t ∈ T ) (see, e.g., Theorem 2.4.9.2 [84]). Due to Eq. (2.6), we haveλmin (C − A> (ȳ)) ==limt→∞ (t∈T )limλmin (C − A> (yt ))(dt )> (C − A> (yt ))dt

t→∞ (t∈T )= d¯> (C − A> (ȳ))d¯ ≥ 0,
which means that the matrix C − A> (ȳ) is PSD. Since b> yt ≥ b> y ? for all t, wehave that b> ȳ ≥ b> y ? . This implies that the accumulation point ȳ is an optimal
solution to the problem (2.1).
When all of the decision variables are binary, we can prove finite convergenceof the algorithm even in the case ε = 0.
Theorem 2.4. If Mc = 0, Algorithm 2.1 terminates in a finite number of iterations.
Proof. Recall that Step 3 excludes the solution yt from the feasible region Ft in
the tth iteration. Since the number of possible solutions y ∈ {0, 1}M is 2M , thealgorithm terminates with Ft = ∅ after at most 2M iterations.

2.4

Branch-and-cut algorithm

Note that Step 1 of the cutting-plane algorithm solves an MILO problem at
every iteration, and therefore we must execute the B&B algorithm many times
from scratch. To resolve this issue, we develop a B&C algorithm that generates
cutting planes dynamically during the B&B procedure. Such dynamic constraint
generation is known as lazy constraint callback in the optimization literature.
This functionality is offered by modern optimization software (e.g., CPLEX or
Gurobi) as a state-of-the-art computational feature. Lazy constraints have been
employed to speed up the computation for various applications, including harvest
scheduling [159], energy supply systems [170], distribution network design [124],
robust optimization [28], and statistical model selection [31]. To the best of our
14knowledge, however, we are the first to use this approach to dealing with PSD
constraints.
Our method starts with the B&B algorithm for solving the relaxed MILO
problem (2.2). Once a feasible solution ŷ ∈ Y is found, the callback procedure
is initiated; if the PSD constraint (2.1b) is violated by ŷ, a cutting plane isappended to exclude ŷ from the feasible region. After that, we proceed with the
B&B algorithm from the middle of the MILO computation. We continue thisprocess until the optimality of the obtained solution is proved through the B&B
procedure. This B&C algorithm is summarized as Algorithm 2.2. Although ouralgorithms are designed for solving the problem (2.1) of the dual form, they canreadily be modified to solve the primal standard form.
Algorithm 2.2 B&C algorithm for solving MISDO problem (2.1)
Step 0 (Initialization) Let ε ≥ 0 be a tolerance for feasibility. Set the feasible
region as F := Y.
Step 1 (B&B Procedure) Start (or continue) the B&B algorithm for solving
the MILO problem,
maximizeyb> ysubject toy ∈ F.

Step 2 (Callback Procedure) Once a feasible solution ŷ ∈ F is found, go to
the following steps:
Step 2.1 (Cut Generation) If λmin (C − A> (ŷ)) < −ε, update the feasible region,on
F ← F ∩ y ∈ {0, 1}Mb × RMc | dˆ> (C − A> (y))dˆ ≥ 0 ,where dˆ is the normalized eigenvector of C − A> (ŷ) corresponding
to the minimum eigenvalue.
Step 2.2 Return to Step 1.

2.5

Numerical experiments

We performed several computational experiments to evaluate the effectiveness of
our algorithms for solving MISDO problems.
2.5.1

Experimental design

We compare the computational performance of the following methods:
SCIP:

MISDO solver SCIP-SDP1 [62],

CUTSDP:MISDO solver CUTSDP implemented in YALMIP2 [106],

CPA:

cutting-plane algorithm (Algorithm 2.1),

B&C:

B&C algorithm (Algorithm 2.2).

12

http://www.opt.tu-darmstadt.de/scipsdp/
https://yalmip.github.io/solver/cutsdp/

15

We also solved some instances by Pajarito [45] to evaluate its performance. Our
preliminary experiments showed that Pajarito was very slow and failed to complete the computations even for small-sized instances. Therefore, we omitted
Pajarito from our target solvers for the comparison in the following numericalexperiments.
We used SCIP-SDP 3.1.0 on the NEOS Server3 [46], where the B&B framework of SCIP4 5.0.0 [63] and the SDO solver DSDP5 5.8 [21] were combined in
the default configuration. When we ran CUTSDP, we used the MATLAB R2015
b and Gurobi Optimizer6 8.0 to solve relaxed linear optimization problems and
MILO problems. CPA and B&C were implemented in the Python language;
MILO problems were solved using the Gurobi Optimizer 8.0, and its callback
function was employed for Step 2 of Algorithm 2.2. The tolerance for feasibility
was set as ε := 10−5 · (1 + kbk1 ) on the basis of the DIMACS error [117], whereb is the coefficient vector of the objective function. These CUTSDP, CPA, andB&C computations were performed on a Windows 7 PC with an Intel Core i74790 CPU (3.60 GHz) and 16 GB of memory. The computation of each method
was terminated if it did not finish by itself within 7200 s. In these cases, theresults obtained at 7200 s were taken as the final outcome.
The column labels used in the tables of experimental results are defined asfollows.
Time:computation time in seconds,#Abort: number of times that the computation was aborted by a numericalissue,#Limit:number of times that the computation reached the time limit of 7200 s,#Cuts:

number of cutting planes generated by CUTSDP, CPA, and B&C,

#Nodes: number of nodes in the enumeration tree formed by SCIP and B&C.
Following Gally et al. [62], we use the shifted geometric mean to aggregate resultsof random instances. The shifted geometric mean of values x1 , x2 , . . . , xN is
defined as

1/N
Y

(xn + s)− s,n∈[N ]

where the shift s was set to 10, 1000, and 100 for “Time,” “#Cuts,” and“#Nodes,” respectively, as in Gally et al. [62].
2.5.2

Random instances

In this subsection, we report experimental results for solving random instancesof the problem (2.1a)–(2.1c).
Instance generation
Following Yamashita et al. [169], we generated random instances of the problem (2.1). First, we constructed ỹ as a feasible solution: ỹm was randomly chosen
3

https://neos-server.org/neos/
http://scip.zib.de/5
http://www.mcs.anl.gov/hs/software/DSDP/
6
http://www.gurobi.com/
416

from {0, 1} for m ≤ Mb and was drawn from a uniform distribution on the interval [0, 1] for m ≥ Mb + 1. The elements of each of the matrices Am (m ∈ [M ])were drawnPfrom a uniform distribution on the interval [−1, 1]. We next set
C := αI + m∈[M ] Am ỹm , where α is a positive-valued parameter ensuring thatthe matrix C − A> (ỹ) is PSD. We then set bm := Am • I for m ∈ [M ].
Experimental results
Table 2.1 gives the computational results for the case of the random instances.
Five instances of the problem (2.1a)–(2.1c) were generated for each tuple (N, Mb , Mc , α),and the shifted geometric means of “Time,” “#Cuts,” and “#Nodes” were computed.
We can see from Table 2.1 that our algorithm B&C was always the fastest
among the four methods. The number of generated nodes in an enumeration
tree was much larger for B&C than for SCIP; however, the series of continuous
relaxation problems at such nodes can be solved efficiently in the case of B&C
computation. CPA was often slower than SCIP, but it is noteworthy that the
SCIP computations were aborted in some instances by numerical issues. In contrast to SCIP, our algorithms do not handle SDO problems, so they can solve
MISDO problems without encountering numerical instability. A typical example is the case of (N, Mb , Mc , α) = (30, 30, 60, 1.0): the computations of CPA
and B&C finished on average in 2069.8 s and 18.8 s, respectively, but the SCIP
computations were aborted in all five instances.
Table 2.1 showed that B&C was always faster than CUTSDP. The main difference between B&C and CUTSDP is that CUTSDP does not implement the lazy
constraints callback function. Thus, CUTSDP does not employ a warm-starting
strategy unlike B&C, and this is considered to be the reason why CUTSDP
was not efficient compared with our B&C algorithm. In addition, comparing CPA and CUTSDP, CPA was faster than CUTSDP except for the case of(N, Mb , Nc , α) = (15, 30, 30, 1.0). While CPA solves the relaxed MILOs from the
beginning, CUTSDP first solves a sequence of relaxed continuous linear optimization problems iteratively with adding cuts at each iteration. Then, CUTSDP
keeps the added constraints and moves into the next phase, imposing integer
constraints and solving MILOs iteratively as our CPA does. The results that
CUTSDP was less efficient than CPA suggest that CUTSDP would add too many
constraints during the first phase and make solving the relaxed MILOs hard in
the second phase.
2.5.3

Computing restricted isometry constants

Here, we compare the performance of the algorithms for the problem of computing
restricted isometry constants [61].
Problem formulation
Let us consider a system of linear equations Ax = b, where the coefficient matrix A ∈ Rµ×N and the right-hand side vector b ∈ Rµ are given. Compressedsensing [13] focuses on finding the sparsest solution x̂0 to an undetermined linear
system with µ ≤ N :x̂0 := arg min{kxk0 | Ax = b},where kxk0 denotes the number of nonzero elements of vector x ∈ RN . This
problem is known to be NP-hard [13].
17Table 2.1: Results for random instances of the problem (2.1a)–(2.1c)
N

Mb

Mc

α

MethodTime#Abort#Limit#Cuts#Nodes

30

30

30

1.0

15

30

30

1.0

60

30

30

1.0

30

15

30

1.0

30

60

30

1.0

30

30

15

1.0

30

60

1.0

61.5
378.4
131.0
4.2
35.9
253.3
315.8
4.4
142.2
1086.7
186.3
4.7
34.6
158.1
78.0
2.1
128.7
978.8
439.8
5.3
30.6
83.5
14.4
0.7
—
2637.02069.8
18.8
22.8
282.5
145.1
2.0
75.3
144.0
97.8
55.1
164.8
3046.61123.811.8
558.0
7169.65069.8
31.0

0
0
0
0
1
0
0
0
1
0
0
0
2
0
0
0
0
0
0
0
0
0
0
0
5
0
0
0
3
0
0
0
2
0
0
0
0
0
0
0
0
0
0
0

0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0

30

SCIP
CUTSDP
CPA
B&C
SCIP
CUTSDP
CPA
B&C
SCIP
CUTSDP
CPA
B&C
SCIP
CUTSDP
CPA
B&C
SCIP
CUTSDP
CPA
B&C
SCIP
CUTSDP
CPA
B&C
SCIP
CUTSDP
CPA
B&C
SCIP
CUTSDP
CPA
B&C
SCIP
CUTSDP
CPA
B&C
SCIP
CUTSDP
CPA
B&C
SCIP
CUTSDP
CPA
B&C—
1541.2
543.6
556.2
—
1279.7
524.8
770.8
—
2235.7
633.1
628.8
—
1478.7
540.3
546.6
—
1624.6
553.8
584.2
—
690.7
185.7
199.1
—
3858.81559.31617.3—
1273.0
445.7
431.1
—
995.1
401.8
3947.6—3107.71066.51112.2—4468.71644.01663.525.4——
1071.763.4——
5270.123.0——
1010.8
15.7
—
—
629.2
51.6
—
—
5183.2
24.1
—
—
446.6
—
—
—
2977.7
1.0
—
—
474.2
47.2
—
—
13111.2
38.8——
2403.850.0——4125.5

30

30

30

0.1

30

30

30

10.0

45

45

45

1

60

60

60

1

18

—
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
4
0
0

The lower and upper restricted isometry constants (RICs) [55] of order κ aredefined asακ := min{kAxk22 | kxk22 = 1, kxk0 ≤ κ},(2.7)βκ := max{kAxk22 | kxk22 = 1, kxk0 ≤ κ}.

(2.8)

When certain conditions on RICs are fulfilled, the sparsest solution x̂0 coincides
with the following `1 -norm minimizer [55]:x̂1 := arg min{kxk1 | Ax = b}.
Gally and Pfetsch [61] proved that exact values of ακ and βκ can be computedby solving the following MISDO problem:
minimize/maximize

Tr(A> AX)(2.9a)

Tr(X) = 1,Xzn ≤ κ,(2.9b)

X,zsubject to(2.9c)n∈[N ]

1
1− zn ≤ Xnm ≤ zn ,22
X  O,(∀(n, m) ∈ [N ] × [N ]),Nz ∈ {0, 1} ,(2.9d)(2.9e)(2.9f)

where X := (Xnm ) ∈ S N and z := (zn ) ∈ {0, 1}N are decision variables.
Experimental results
Tables 2.2 and 2.3 give the results for computing lower and upper RICs, respectively. The size of the problem (2.9) is represented by N ∈ {10, 15, 20, 25}, andthe order of RICs is κ ∈ {3, 5}. Note that the number of linear equations in thesystem Ax = b was fixed to µ = 10 because it is independent of the size of the
problem (2.9). As in the case of Gaussian-distributed matrices [61], each elementof matrix A was drawn from the standard normal distribution. Five instances ofmatrix A were generated for each pair (N, κ), and the shifted geometric meansof “Time,” “#Cuts,” and “#Nodes” were computed. In Tables 2.2 and 2.3,
each column “Mb ” and “Mc ” are the number of integer variables and that ofcontinuous ones, respectively.
We can see from Tables 2.2 and 2.3 that, as was the case for random instances,
B&C was the fastest among the four methods. CPA was also faster than SCIP and
CUTSDP in the case of computing RICs. In addition, we notice that when the
problem size (i.e., N ) was large, the SCIP computation was frequently aborted
due to a numerical issue. Indeed, when N ≥ 20, SCIP failed to complete thecomputation for 19 out of 20 instances of computing lower RICs (Table 2.2) andfor 11 out of 20 instances of computing upper RICs (Table 2.3). Also, CUTSDP
did not solve all instances of computing lower RICs (Table 2.2) and upper RICs
(Table 2.3) within the time limit of 7200 s for N ≥ 25.
2.5.4

Robust truss topology design

In this subsection, we report the experimental results for our tests of the threealgorithms on the problem of robust truss topology design [171].

19

Table 2.2: Results for computing lower RICs (µ = 10)
N

κ

Mb

Mc

MethodTime#Abort#Limit#Cuts#Nodes

10

3

45

15

10

5

45

15

15

3

105

15

15

5

105

15

20

3

190

20

5

190

20

6.1
21.9
1.7
1.0
16.6
50.7
14.0
2.3
140.3
2468.1
11.5
3.8
693.4
7200.0237.9
14.1
1037.67200.058.7
8.9—
7200.01938.7112.5—
7200.0188.4
16.9—
7200.07200.0373.0—
7200.0616.4845.33—
7200.07200.01343.6—7200.02257.4102.1—
7200.07200.03802.6

0
0
0
0
2
0
0
0
0
0
0
0
0
0
0
0
4
0
0
0
5
0
0
0
5
0
0
0
5
0
0
0
5
0
0
0
5
0
0
0
5
0
0
0
5
0
0
0

0
0
0
0
0
0
0
0
0
1
0
0
0
5
0
0
0
5
0
0

20

SCIP
CUTSDP
CPA
B&C
SCIP
CUTSDP
CPA
B&C
SCIP
CUTSDP
CPA
B&C
SCIP
CUTSDP
CPA
B&C
SCIP
CUTSDP
CPA
B&C
SCIP
CUTSDP
CPA
B&C
SCIP
CUTSDP
CPA
B&C
SCIP
CUTSDP
CPA
B&C
SCIP
CUTSDP
CPA
B&C
SCIP
CUTSDP
CPA
B&C
SCIP
CUTSDP
CPA
B&C
SCIP
CUTSDP
CPA
B&C

—
730.5
55.4
479.2
—
945.9
183.1
1087.0—4369.6
107.6
943.5
—
4296.1342.4
1718.1—5478.4191.7
1155.7—5153.0
669.5
533.3
—
5149.2294.4
1407.6—5332.31178.91008.0—5116.3449.1
1123.2—5619.52104.32272.5—6635.7
652.3
609.7
—
6696.32776.73336.7

27.6
—
—
751.6
91.4
—
—
2827.2
187.6——
1864.8891.4——
6494.2373.0——
3031.0———21239.5———
5435.7———60996.7———
9828.7———136486.1———
16269.0———225854.9

25

25

30

30

35

35

3

5

3

5

3

5

300

300

435

435

595

595

25

25

30

30

35

35

20

—
5
0
0
—
5
0
0
—
5
5
0
—
5
0
0
—
5
5
0
—
5
0
0
—
5
5
0

Table 2.3: Results for computing upper RICs (µ = 10)
N

κ

Mb

Mc

MethodTime#Abort#Limit#Cuts#Nodes

10

3

10

45

10

5

10

45

15

3

15

105

15

5

15

105

20

3

20

190

20

5

20

190

25

3

25

300

5

25

300

5.4
13.5
1.5
0.6
9.7
53.9
4.8
2.0
41.3
231.3
3.8
2.7
54.5
2657.9
22.8
6.2
322.0
1479.7
8.1
3.8
605.9
3651.5129.2
14.6
1224.07200.023.9
7.9—
7200.01285.268.9—
7200.040.8
14.7—
7200.02380.4177.9—
7200.092.4
26.5—
7200.06426.1
318.0

0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
2
0
0
0
4
0
0
0
5
0
0
0
5
0
0
0
5
0
0
0
5
0
0
0
5
0
0
0

0
0
0
0
0
0
0
0
0
0
0
0
0
3
0
0
0
1
0
0
0
3
0
0
0
5
0
0

25

SCIP
CUTSDP
CPA
B&C
SCIP
CUTSDP
CPA
B&C
SCIP
CUTSDP
CPA
B&C
SCIP
CUTSDP
CPA
B&C
SCIP
CUTSDP
CPA
B&C
SCIP
CUTSDP
CPA
B&C
SCIP
CUTSDP
CPA
B&C
SCIP
CUTSDP
CPA
B&C
SCIP
CUTSDP
CPA
B&C
SCIP
CUTSDP
CPA
B&C
SCIP
CUTSDP
CPA
B&C
SCIP
CUTSDP
CPA
B&C

—
416.7
30.2
308.4
—
826.5
77.8
880.4
—
1091.734.2
683.0—
2028.393.7
1084.2—1177.4
34.8
511.9
—
922.2
121.9
997.3
—
1072.646.4
681.2—
1106.3
165.4
922.4
—
1185.444.6
654.5—
1424.9168.91199.2—1519.346.8
301.3—
1887.3
158.4
938.1

26.3
—
—
429.1
74.3
—
—
2233.651.7——
1465.076.6——
2608.2
93.1——
1662.8130.8——6100.2137.0——
3211.8———13854.4———4552.8———24939.4———
5424.4———56816.4

30

30

35

35

3

5

3

5

30

30

35

35

435

435

595

595

21

—
5
0
0
—
5
0
0
—
5
0
0
—
5
0
0
—
5
4
0

Problem formulation
Let us consider an undirected graph called the ground structure in d-dimensional
space. We suppose that some nodes are fixed, and d-dimensional forces act on theq unfixed nodes; so the external load vector is defined as f ∈ Rν with ν := d · q.
The index set in the external load vector is denoted by J := {1, 2, . . . , ν}. The
index set of edges in the ground structure is denoted by I := {1, 2, . . . , µ}, which
represents candidate members (or bars) of a truss. Truss topology design involves
choosing a subset of members and determining their cross-sectional areas so thatthe truss structure will be lightweight but sufficiently rigid.
We begin by introducing the decision variables. We denote by a := (ai )i∈I ∈Rµ the vector of member cross-sectional areas. Additionally, the vector p :=(pj )j∈J ∈ {0, 1}ν corresponds to the existence of each node, and the vector
t := (ti )i∈I ∈ {0, 1}µ represents the selection of members in the truss structure.
ν
P Given matrices Ki ∈ S (i ∈ I), the stiffness matrix is defined as K(a) :=i∈I ai Ki . The following matrix is used to express uncertainty about the external loads applied at each node in the truss structure:
Q := (f˜, rv1 , rv2 , . . . , rvν−1 ) ∈ Rν×ν ,where f˜ ∈ Rν is the nominal external load, r is the level of uncertainty, andvj ∈ Rν (j = 1, 2, . . . , ν − 1) are orthonormal basis vectors of the orthogonal
complement of f˜. When an ellipsoidal uncertainty set of external loads [19] isemployed, the worst-case compliance constraint is posed asτ̄ I(Diag(p)Q)>
 O,(2.10)
Diag(p)Q
K(a)where τ̄ is the upper-bound parameter for the compliance; see Yonekura andKanno [171] for the details.
We define Jf ⊆ J to represent the nodes that should remain in the truss
structure. Other nodes can be removed, but we must maintain the nodes thatare connected to at least one member. Such constraints can be expressed as
0 ≤ pj ≤ 1,(∀j ∈ J \ Jf ),(2.11)ti ≤ pj ,(∀i ∈ Ik̂(j) , ∀j ∈ J \ Jf ),(2.12)

pj = 1,(∀j ∈ Jf ),(2.13)

where Ik̂(j) ⊆ I is the set of members that are connected to the node associatedwith the jth element of f˜.
The following constraints are imposed on the member cross-sectional areas:amin ti ≤ ai ≤ amax ti ,ti ∈ {0, 1},(∀i ∈ I),(∀i ∈ I),

(2.14)(2.15)

where amin and amax are, respectively, lower and upper bounds on cross-sectional
areas.
The robust truss topology design minimizing the structural volume is formulated as the following MISDO problem:
Xminimize`i ai(2.16a)a, p, tsubject toi∈I

Eqs. (2.10)–(2.15),(2.16b)where `i is the length of the ith member. Although p is treated as a vector ofcontinuous variables, the binary constraint p ∈ {0, 1}ν can be satisfied.
22

Experimental results
We considered a two-dimensional lattice as the ground structure; nodes wereplaced on lattice points, and every pair of distinct nodes was connected by aunique edge. Note that the longer edge was removed if two edges overlapped
along the same straight line, and the remaining edges were treated as candidate
members.
As explained below, the settings of parameters were the same as those usedby Yonekura and Kanno [171]. Specifically, the lengths of horizontal and verticalmembers were, respectively, 100 cm and 50 cm, the lower and upper bounds oncross-sectional areas were, respectively, amin := 10 cm2 and amax := 1000 cm2 ,and the elastic modulus was 200 GPa. The leftmost nodes of trusses were fixed,
and the bottom-right node was loaded with a vertical force of 1 kN as a nominal
external load f˜. The level of uncertainty was r := 0.1, and the upper-boundparameter for compliance was τ̄ := 10 kN. Note that all lengths were measured
in centimeters, and that the value of each member length, `i , was divided by 1000
to reduce numerical error.
Table 2.4 gives the computational results for the robust truss topology design
problem. In Tables 2.4, each column “N”, “Mb ”, and “Mc ” are the matrix sizein Eq. (2.10), the number of integer variables, and that of continuous ones,
respectively. For instance, Lat(3,3) is a problem instance corresponding to 3 × 3
lattice points. Also, 22mem and 51mem stand for the 22- and 51-member truss
instances defined in Yonekura and Kanno [171]. The column labeled “Obj” shows
the value of the objective function (2.16a). If the computation reached the time
limit of 7200 s, the obtained lower and upper bounds on the objective function
are shown.
We can see from Table 2.4 that there is no clear inferior-to-superior relationship between B&C and SCIP for small-sized instances. For example, B&C was
three times faster than SCIP for the Lat(3,4) instance, whereas B&C was fourtimes slower than SCIP for the 22mem instances. For large-sized instances, thecomputation time of B&C was relatively long. Indeed, B&C took much longer
to solve the Lat(4,4) and 51mem instances than SCIP did. However, it is noteworthy that B&C successfully found solutions of good quality to such large-sized
instances. For the Lat(4,4) instance, although the B&C computation was terminated due to the time limit, its computed upper-bound value (i.e., the objective
value of the obtained feasible solution) was 71.01, which is very close to the
optimal value (i.e., 69.34) obtained by SCIP. For the Lat(5,3) instance, B&C
computed an upper-bound value of 294.49, which is very close to the optimalvalue (i.e., 292.68) provided by CPA and CUTSDP, whereas SCIP failed to solve
this instance due to a numerical issue. On the other hand, CPA showed poorer
computational performance than both SCIP and B&C for robust truss topology
design.
Comparing B&C and CUTSDP, Table 2.4 showed that B&C was faster than
CUSTDP except for the cases of the two instances, 22mem and Lat(5,3). For the
22mem instance, B&C was slower than CUTSDP, but the difference in the computational time between these methods was small. Moreover, as we mentioned
before, B&C obtained a good upper-bound value for the Lat(5,3) instance. In
terms of a comparison between CPA and CUTSDP, Table 2.4 showed that CPA
was faster then CUTSDP except for the cases of the three instances, Lat(3,3),
Lat(4,3) and 22mem. For each of the three instances, while CPA took longer
to solve than CUTSDP did, the computational time of CPA and CUTSDP werecomparable.
23

Table 2.4: Results for robust truss topology design
Instance

ν

µ

N

Mb

Mc

Method

Time

Obj

#Cuts

#NodesLat(3,3)

12

26

24

26

38

Lat(3,4)

16

46

32

46

65

Lat(4,3)

18

47

36

47

96

Lat(3,5)

20

70

40

70

90

Lat(4,4)

24

83

48

83

107

Lat(5,3)

24

72

48

72

96

22mem

12

22

24

22

34

51mem

18

51

36

51

69

SCIP
CUTSDP
CPA
B&C
SCIP
CUTSDP
CPA
B&C
SCIP
CUTSDP
CPA
B&C
SCIP
CUTSDP
CPA
B&C
SCIP
CUTSDP
CPA
B&C
SCIP
CUTSDP
CPA
B&C
SCIP
CUTSDP
CPA
B&C
SCIP
CUTSDP
CPA
B&C

18.6
35.2
55.1
24.4
45.2
1590.3
943.1
12.9
326.1
629.8
741.3
324.8
593.0
>7200
>7200
424.8
850.8
>7200
>7200
>7200
Abort
6852.76324.3
>7200
7.7
24.1
39.2
36.3
336.5
>7200
>7200
4308.7

32.40
32.40
32.40
32.40
21.24
21.24
21.24
21.24
114.07114.07114.07114.0717.32[15.92, —][16.14, —]
17.32
69.34
[64.82, —][64.61, —][69.31, 71.01]—
292.68292.68[292.67, 294.49]
33.27
33.27
33.27
33.27
75.05
[74.95, —][74.86, —]
75.05

—
497
347
2612
—
764
479
1111
—
1384
982
8855
—
>962
>433
3902
—
>1690
>996
>10059—2963
2127>22398
—
467
321
7507
—
>1999
>1375
17680

71
—
—
17115
101
—
—
14594
755
—
—
85147
354
—
—
149298525——>423279———>336055
23
—
—
15343
236
—
—
267973

2.6

Conclusion

In this chapter, we have described general-purpose algorithms to solve such problems. Firstly, we formulated a cutting-plane algorithm to solve MISDO problems
and proved its convergence properties. We then developed a B&C algorithm,
where cutting planes are added dynamically to the relaxed MILO problem during a B&B procedure. The experimental results show that our B&C algorithm
was effective in solving three different MISDO problems: random instances, computing restricted isometry constants, and robust truss topology design. Ouralgorithm is based on the B&B procedure for solving MILO problems, so it canuse a warm-starting strategy to solve a series of continuous relaxation problems
efficiently. Moreover, our algorithms do not solve any MISDO problems, whichmakes their computation very stable and practical.

24Chapter 3
Cutting-plane Algorithm for Best Subset
Selection for Eliminating Multicollinearity
In this chapter, we study the best subset selection problem for eliminating multicollinearity from linear regression models. Specifically, we focus on the problem
of selecting the best subset of explanatory variables subject to an upper bound on
the condition number of the correlation matrix of the selected variables. First, weformulate this problem as an MISDO problem. Second, we develop a specialized
cutting-plane algorithm that, to approximate the condition number constraint,
iteratively appends cutting planes to the mixed-integer quadratic optimization
problem. Here, to improve the computational efficiency of our algorithm, we
propose to generate strong cutting planes by effectively using a heuristic search.
The content of this chapter is included in Tamura et al. [151].
We give an introduction to this study and summarize our contribution in
Section 3.1. In Section 3.2, we give a brief review of linear regression and subset
selection for eliminating multicollinearity. In Section 3.3, we formulate the best
subset selection under the condition number constraint as an MISDO problem.
In Section 3.4, we describe our specialized cutting-plane algorithm for solving thebest subset selection problem. In Section 3.5, we show the computational results,
and in Section 3.6, we conclude with a brief summary.

3.1

Introduction

Multicollinearity, which exists when two or more explanatory variables in a regression model are highly correlated, is a frequently encountered problem in multiple
regression analysis [54, 76, 112]. Such an interrelationship among explanatory
variables obscures their relationship with the explained variable, leading to computational instability in model estimation. Moreover, the reliability of the regression analysis is decreased in the presence of multicollinearity by the low qualityof the resultant estimates.
Several approaches can be used to avoid the deleterious effects of multicollinearity [42]. One approach is orthogonal transformation through procedures
such as principal component regression [88, 115] and partial least squares regression [164, 165]. In this approach, a set of correlated variables is transformed
into a set of linearly uncorrelated variables (i.e., principal components) for use ina regression model. Orthogonal transformation can enhance the computational
stability of model estimation but often leads to worse predictive performance andresults that are strongly influenced by the presence of outliers [58, 78].
Another approach is penalized regression, such as ridge regression [83], lasso [154],and elastic net [174]. This approach introduces a penalty function to shrink regression coefficient estimates toward zero. Penalized regression helps prevent25

regression models from overfitting noisy datasets and, accordingly, is effective forachieving high predictive performance. However, the penalty functions produce
biased estimates, which are undesirable from the standpoint of model interpretation [30, 32].
In this chapter, we focus on subset selection, which is a simple but effective approach for eliminating multicollinearity. Conventionally in this approach,explanatory variables are deleted iteratively through the use of indicators for
detecting multicollinearity, such as condition number of the correlation matrix
and variance inflation factor [42]. Chong and Jun [44] have compared the performance of subset selection methods with that of partial least squares regression
and suggested that a goodness-of-fit measure for evaluating a subset regression
model should be chosen carefully when multicollinearity is present. For subset
selection in the presence of multicollinearity, several researchers [49, 87] have proposed goodness-of-fit measures based on ridge regression. It is notable, however,that commonly used subset selection methods are heuristic algorithms, such as
stepwise regression [50]; hence, they do not necessarily find the best subset ofvariables in terms of a given goodness-of-fit measure.
The mixed-integer optimization (MIO) approach to subset selection has recently received much attention due to advances in algorithms and hardware [30,32, 92, 99, 101, 118, 119, 140, 141]. In contrast to heuristic algorithms, the MIO
approach has the potential to provide the best subset of variables under a given
goodness-of-fit measure. To deal with multicollinearity, Bertsimas and King [30]use a cutting-plane strategy, which iteratively adds constraints for deleting subsets of collinear variables. However, this strategy requires solving an enormous
number of mixed-integer quadratic optimization (MIQO) problems when multicollinearity exists in many different sets of variables.
The aim of this chapter is to devise a more sophisticated MIO approach tobest subset selection for eliminating multicollinearity. In particular, this chapter
addresses the following problem: Find a subset of variables that minimizes theresidual sum of squares under the constraint that the condition number of theassociated correlation matrix is bounded. We first formulate this subset selection problem as an MISDO problem. Furthermore, to increase computational
efficiency, we consider incorporating constraints based on the normal equations
into the MISDO problem. Regrettably, as we show later, current general-purpose
MISDO solvers do not deliver high computational performance and the problem
size to be handled is limited. Thus, we secondly propose to develop a specialized
cutting-plane algorithm for subset selection that works well in practice. Theoretically, this cutting-plane algorithm must solve an exponential number of MIQO
problems, which are NP-hard, in a worst-case situation. However, it has a practical advantage that we can use sophisticated MIO solvers such as Gurobi and
CPLEX to solve MIQO problems. In addition, we propose to use a backwardelimination method that searches a smaller subset of collinear variables to generate strong cutting planes, which improves the efficiency of the cutting-planealgorithm.
The effectiveness of our MIO approaches is assessed through computational
experiments using several datasets from the UCI Machine Learning Repository [109].
We succeeded in solving some of our MISDO problems for subset selection when
the number of candidate explanatory variables was less than 26. In addition, thecomputational results demonstrate that our specialized cutting-plane algorithm
yielded optimal solutions when the number of variables was 65, and it frequently
gave a better subset of variables than did the conventional local search algorithms
for large-sized datasets.
26

3.2

Linear regression and multicollinearity

This section contains a brief review of linear regression and subset selection foreliminating multicollinearity.
3.2.1

Linear regression

Let us suppose that we are given N samples, (xn , yn ) for n ∈ [N ]. Here, xn :=(xn1 , xn2 , . . . , xnP )> is a vector composed of P explanatory variables, and yn isan explained variable for each sample n ∈ [N ].
We focus on the following linear regression model:X
ap xnp + εn (∀n ∈ [N ]),yn = a> xn + εn =p∈[P ]

where a := (a1 , a2 , . . . , aP )> is a vector of regression coefficients to be estimatedand εi is a prediction residual for each sample n ∈ [N ]. We assume
P here that all
explanatory and explained variables are standardized so that ( n∈[N ] xnp )/N =PPP( n∈[N ] yn )/N = 0 and ( n∈[N ] (xnp )2 )/N = ( n∈[N ] (yn )2 )/N = 1 for all p ∈[P ]. Therefore, no intercept (constant term) is present in the regression model.
The above linear regression model can be rewritten as
y = Xa + ε,
where y := (y1 , y2 , . . . , yN )> , ε := (ε1 , ε2 , . . . , εN )> , and
x11 x12 · · · x1P
 x21 x22 · · · x2P X :=  .
..
..
..  .
.
 .
.
.
. 
xN 1 xN 2 · · · xN P
To estimate the regression coefficients, a, the ordinary least squares methodminimizes the residual sum of squares (RSS):ε> ε = (y − Xa)> (y − Xa).

(3.1)

After partial differentiation, the ordinary least squares method is equivalent tosolving a system of linear equations:
X > Xa = X > y.

(3.2)

This is the well-known normal equation.
3.2.2

Subset selection for eliminating multicollinearity

We use the condition number of the correlation matrix R := (rpq ) ∈ S P to
detect multicollinearity. Because the explanatory variables are standardized, thecorrelation matrix is calculated as
R=

X >X
,
N

and its condition number is defined asλmax (R)(λmin (R) > 0),λmin (R)cond(R) :=+∞(λmin (R) = 0),
27

(3.3)

where λmin (R) and λmax (R) are the minimum and maximum eigenvalues of matrix R, respectively. It follows from cond(R) = cond(X > X) that when cond(R)is very large, the coefficient matrix of equations (3.2) is ill-conditioned, which
implies that the regression coefficient estimates are subject to large numericalerrors.
To compute accurate estimates, we consider selecting a subset S ⊆ [P ] ofcandidate explanatory variables. Let us denote by RS the correlation sub-matrix
of a subset of variables, that is, RS := (rpq ; p, q ∈ S). To avoid multicollinearity,
the condition number of the sub-matrix should not exceed a user-defined parameter κ (> 1). The subset selection problem determines a subset S of explanatory
variables so that the residual sum of squares of a subset regression model isminimized:

2
X
X
yn −
(3.4a)
minimize
ap xnp a, Sn∈[N ]

subject top∈S

cond(RS ) ≤ κ,

(3.4b)S ⊆ [P ].

(3.4c)The subset selection problem (3.4) can be converted into an MIO problem. Letz = (z1 , z2 , . . . , zP )> be a vector of 0-1 decision variables for selecting explanatory
variables; that is, zp = 1 if p ∈ S; otherwise, zp = 0. The correlation sub-matrix
is then written as R(z) := (rpq ; zp = zq = 1). Consequently, the subset selection
problem is formulated as an MIO problem,
minimizea, z(y − Xa)> (y − Xa)subject to zp = 0 ⇒ ap = 0(∀p ∈ [P ]),(3.5a)(3.5b)

cond(R(z)) ≤ κ,

(3.5c)z ∈ {0, 1}P .

(3.5d)

Here, if zp = 0, then the pth explanatory variable is deleted from the regression
model because its coefficient is set to zero by the logical implications (3.5b). It
is known that these logical implications can be represented by using a big-M
method or a special ordered set type 1 (SOS1) constraint [14, 15].

3.3

Mixed-integer semidefinite optimization ApproachThis section presents an MISDO approach to best subset selection for eliminating
multicollinearity.
3.3.1

Formulation

A convex quadratic objective function (3.5a) is expressed as a linear objectivefunction with a PSD constraint [155, 161]. We begin by computing a decomposition of the form:
X >X = N R = V V >,where the square matrix V ∈ RP ×P can be created, e.g., by the Cholesky/eigenvalue
decomposition. Introducing a scalar decision variable f to be minimized, werewrite the associated constraint as a PSD constraint as follows:


IP
V >a
>
(y − Xa) (y − Xa) ≤ f ⇐⇒
 O,
a> V 2y > Xa − y > y + f28

where I is the identity matrix of size P , O is a zero matrix of appropriate size,
and A  B means that A − B is a PSD matrix.
We denote by Diag(x) the diagonal matrix whose diagonal entries are components of vector x. We also denote by A ◦ B the Hadamard product of matrices
A and B. The next theorem shows that the condition number constraint (3.5c)
is expressed as PSD constraints based on the following matrices:
1 − z1
0
···
0


.
 0
1 − z2 . .
0 
Diag(1 − z) =  .
..  ,
..
..
 ..
.
.
. 
0
0
· · · 1 − zP

r11 z1 z1
 r21 z2 z1R ◦ zz > = ..
.

r12 z1 z2
r22 z2 z2
..
.

···
···
..
.


r1P z1 zP
r2P z2 zP 

.
..

.

rP 1 zP z1 rp2 zP z2 · · · rP P zP zP
Theorem 3.1. Suppose that z ∈ {0, 1}P . Then, cond(R(z)) ≤ κ holds if andonly if there exists λ ∈ [1/κ, 1] such that
λIP − Diag(1 − z)  R ◦ zz >  κλIP .

(3.6)

Proof. Let Q be the number of nonzero elements of z. Without loss of generality,
we assume that(zp = 1 for p ∈ [Q],(3.7)
zp = 0 for p ∈ [P ] \ [Q]Since R(z) is a PSD matrix whose diagonal entries are all one, 0 ≤ λmin (R(z)) ≤1 ≤ λmax (R(z)) holds. It then follows from the definition (3.3) that cond(R(z)) ≤κ is equivalent toλmax (R(z)) ≤ κλmin (R(z)).
This also implies that 1/κ ≤ λmax (R(z))/κ ≤ λmin (R(z)). Therefore, it isnecessary to consider only λmin (R(z)) ∈ [1/κ, 1].
Using a PSD constraint for minimizing the maximal eigenvalue [155, 161], thecondition number constraint can be converted as follows:λmax (R(z)) ≤ κλmin (R(z))⇐⇒ ∃λ ∈ [1/κ, 1], λ ≤ λmin (R(z)) and λmax (R(z)) ≤ κλ⇐⇒ ∃λ ∈ [1/κ, 1], λIQ  R(z)  κλIQ

 



λIQ
O
R(z) OIQ
O⇐⇒ ∃λ ∈ [1/κ, 1], κλO (λ − 1)IP −Q
OOO IP −Q⇐⇒ ∃λ ∈ [1/κ, 1], λIP − Diag(1 − z)  R ◦ zz >  κλIP .

29

To linearize the bilinear term zz > in constraint (3.6), we introduce a symmetric matrix of decision variables:
w11 w21 · · · wP 1
 w21 w22 · · · wP 2 
W = .
..
..  .
.
.
.
 .
.
.
. 
wP 1 wP 2 · · · wP P
It is known that when z ∈ {0, 1}P , wpq = zp zq can be rewritten by means of its
convex and concave envelopes as follows [5, 116]:
wpq ≥ 0, wpq ≥ zp + zq − 1, wpq ≤ zp , wpq ≤ zq .
Consequently, the subset selection problem is cast into an MISDO problem,
minimizea, f, λ, W , z

subject tof

(3.8a)
IP
a> V

V >a
2y > Xa − y > y + f

zp = 0 ⇒ ap = 0 O,(∀p ∈ [p]),λIP − Diag(1 − z)  R ◦ W  κλIP ,wpp = zp(∀p ∈ [P ]),(3.8b)(3.8c)(3.8d)(3.8e)

wpq ≥ 0, wpq ≥ zp + zq − 1, wpq ≤ zp , wpq ≤ zq(∀(p, q) ∈ [P ] × [P ]),1/κ ≤ λ ≤ 1,
3.3.2

z ∈ {0, 1}P .

(3.8f)(3.8g)

Normal-equation-based constraints

MISDO problems can be handled by a branch-and-bound procedure, but it involves solving a large number of relaxed semidefinite optimization problems.
To improve computational efficiency, we consider including the normal equations (3.2) as the constraints of the MISDO problem. More precisely, when thepth explanatory variable is selected, the pth normal equation is placed as follows:
X > Xa + s = X > y,(3.9)

zp = 1 ⇒ s p = 0(3.10)(∀p ∈ [P ]),where s = (s1 , s2 , . . . , sP )> is a vector of auxiliary decision variables.
The next theorem shows that constraints (3.9) and (3.10) are necessary optimality conditions for problem (3.5).
Theorem 3.2. Let (a? , z ? ) be an optimal solution to problem (3.5). There exists
s? = (s?1 , s?2 , . . . , s?P )> such that
X > Xa? + s? = X > y,zp? = 1 ⇒ s?p = 0(∀p ∈ [P ]).

Proof. Without loss of generality, we can assume that
 
1
?
z =
, 1 ∈ RQ , 0 ∈ RP −Q .
030

According to z ? , we partition X and a? as follows:
X = X1 X2 , X1 ∈ RN ×Q , X2 ∈ RN ×(P −Q) ,
 ?
a1
?
a =
, a?1 ∈ RQ , a?2 ∈ RP −Q .
a?2
Because (a? , z ? ) is an optimal solution to problem (3.5), we have a?2 = 0. Moreover, a?1 minimizes RSS (3.1); that is, the following holds for a?1 and the associatednormal equation:
X1> X1 a?1 = X1> y.
Now we define s? as follows:
?



s =

0
>
X2 y − X2> X1 a?1.

It then follows that >   > 0
X1 X1 a?1
X1 y
>
?
?
X Xa + s =+== X > y,X2> y − X2> X1 a?1
X2> X1 a?1
X2> yzp? = 1 ⇒ s?p = 0

3.4(∀p ∈ [P ]).

Cutting-plane algorithm

This section presents a cutting-plane algorithm for solving Problem (3.5). The
motivation behind this approach is that, as we see later in numerical experiments,
the current general-purpose MISDO solver could solve only small-sized instancesof Problem (3.8) and suffer from numerical instability. In contrast to the MISDO
approach, where we can obtain an optimal solution by solving a single MIO
problem, the cutting-plane algorithm has a theoretical drawback: it must solve
an exponential number of MIQO problems, which are NP-hard, in a worst-casesituation. However, it has a practical advantage in using sophisticated MIO
solvers such as Gurobi and CPLEX to solve MIQO problems. In addition, we
can improve its efficiency by generating strong cutting planes by exploiting thestructure of the condition number constraint Eq. (3.5c).
The cutting-plane algorithm first removes the condition number constraint (3.5c)from the problem. Hence its feasible set is expressed as
F1 := {(a, z) ∈ RP × {0, 1}P | Eqs. (3.5b) and (3.5d)},(3.11)and the problem reduces to an MIQO problem, which can be handled by sophisticated MIO solvers such as Gurobi and CPLEX. The basic strategy of thealgorithm involves repeatedly solving such relaxed MIQO problems and iteratively adding valid inequalities, instead of imposing the condition number constraint (3.5c), to the MIQO problems.
At the tth iteration (t ≥ 1), our algorithm solves the follwing MIQO problem:

minimize(y − Xa)> (y − Xa)(3.12a)subject to(a, z) ∈ Ft ,(3.12b)a, z31

where Ft is a relaxed feasible set at the tth iteration such that Ft ⊆ F1 , which
will be updated according to Eq. (3.13) at the tth iteration. Let (at , zt ) bean optimal solution to Problem (3.12). If cond(R(zt )) ≤ κ, then (at , zt ) is an
optimal solution to Problem (3.5). In this case, we terminate the algorithm.
Otherwise, we update the feasible region to cut off the solution zt as follows:
Ft+1 ← Ft ∩ {(a, z) ∈ RP × {0, 1}P | zt> z ≤ zt> zt − 1}.

(3.13)

After that, we set t ← t + 1 and solve the MIQO problem (3.12) again. We
repeat this process until we obtain a solution that satisfies the condition number
constraint (3.5c).
We next show that the condition number is not increased by deleting explanatory variables.
Lemma 3.3. Suppose that z, z̄ ∈ {0, 1}P .
cond(R(z̄)).

If z ≥ z̄, then cond(R(z)) ≥Proof. It follows from the Cauchy’s interlace theorem (see, e.g., Theorem 4.3.17[84]) that0 ≤ λmin (R(z)) ≤ λmin (R(z̄)) ≤ λmax (R(z̄)) ≤ λmax (R(z)),
which completes the proof.
The next theorem states that the feasible set of the original problem (3.5) is
contained in Ft for all t.
Theorem 3.4. Suppose that z̄ ∈ {0, 1}P and cond(R(z̄)) > κ. If z ∈ {0, 1}Psatisfies cond(R(z)) ≤ κ, then it also satisfies z̄ > z ≤ z̄ > z̄ − 1.
Proof. We prove the proposition by contradiction:z̄ > z > z̄ > z̄ − 1 ⇒ z̄ > (z − z̄) ≥ 0∵ z, z̄ ∈ {0, 1}P⇒ z ≥ z̄⇒ cond(R(z)) ≥ cond(R(z̄))∵ Lemma 3.3⇒ cond(R(z)) > κ.

This cutting-plane algorithm is developed on the basis of Bertsimas andKing [30] and requires solving a large number of MIQO problems. To reduce
the number of MIQO problems to be solved, we develop stronger valid inequalities for approximating the condition number constraint (3.5c). To this end, we
employ a backward elimination method that searches a smaller subset of collinear
variables. Specifically, it starts with an incumbent solution (e.g., zt ) and deletes
explanatory variables one by one on the basis of the RSS (3.1). Finally, we obtainz̄ (≤ zt ) such that cond(R(z̄)) > κ. The feasible set is then updated:
Ft+1 ← Ft ∩ {(a, z) ∈ RP × {0, 1}P | z̄ > z ≤ z̄ > z̄ − 1}.

(3.14)

This cutting plane cuts off all z ∈ {0, 1}P satisfying z ≥ z̄; therefore, it is stronger
than the previous one (3.13) because zt ≥ z̄. Our cutting-plane algorithm issummarized by Algorithm 3.1
We next prove the finite convergence of the algorithm.

32

Algorithm 3.1 Cutting-plane algorithms for solving Problem (3.5)
Step 0 (Initialization) Set t ← 1. Let F1 be a feasible set (3.11).
Step 1 (Relaxed Problem) Solve Problem (3.12). Let (a, zt ) be an optimal
solution.
Step 2 (Muliticollinearity Detection) If cond(R(zt )) ≤ κ, terminate thealgorithm with the solution zt .
Step 3 (Backward Elimination) Find a solution z̄ ∈ {0, 1}P such thatz̄ ≤ zt and cond(R(z̄)) > κ by using a backward elimination method
starting with zt .
Step 4 (Cut Generation) Add cut (3.14) to update the feasible set Ft .
Step 5 Set t ← t + 1 and return to Step 1.

Theorem 3.5. Algorithm 3.1 provides an optimal solution to problem (3.5) in a
finite number of iterations.
Proof. Step 3 in Algorithm 3.1 removes zt from Ft in each iteration; therefore,
the algorithm terminates with a feasible solution (e.g., z = (1, 0, 0, . . . , 0)> ) after
at most 2P iterations.
Let (a? , z ? ) be an optimal solution to problem (3.5). Theorem 3.4 guarantees
that the feasible set of problem (3.5) is contained in Ft and hence that (y −Xat )> (y − Xat ) ≤ (y − Xa? )> (y − Xa? ) for all t. Therefore, the algorithm
provides an optimal solution to problem (3.5).

3.5

Numerical experiments

In this section, we assess the computational performance of our mixed-integer
optimization approaches to best subset selection for eliminating multicollinearity.
We downloaded six datasets for regression analysis from the UCI Machine
Learning Repository [109]. Tables 3.1 lists the instances used for computational
experiments, where N and P are the number of samples and number of candidate
explanatory variables, respectively. In the SolarFlareC instance, C-class flares
production was employed as an explained variable. In the ForestFires instance,interaction terms were created from the variables of the x-axis and y-axis spatial
coordinates. Each categorical variable was transformed into one or more dummy
variables. All explanatory and explained variables were standardized to a meanof zero and standard deviation of one, as mentioned in Section 2.1. Samples
containing missing values and redundant variables having the same value in all
samples were eliminated.
3.5.1

Computational performance of MISDO approaches

We first evaluate the computational performance of the following MISDO approaches:MISDO

MISDO formulation (3.8)

MISDONE MISDO formulation (3.8) with the normal-equation-based constraints (3.9) and (3.10)

33

Abbreviation
Servo
AutoMPG
SolarFlareC
BreastCancer
ForestFires
Automobile

Table 3.1: List of instances
N P Original dataset [109]
167
392
1066
194
517
159

19
25
26
32
63
65

Servo
Auto MPG
Solar Flare (C-class flares production)
Breast Cancer Wisconsin
Forest Fires
Automobile

Table 3.2: Values of big-M for constraints (3.16)
M
Instance

N

P

Method(κ = 100)(κ = 225)

Servo

167

19

MISDO
MISDONE

2.25
2.50

2.25
2.50

AutoMPG

392

25

MISDO
MISDONE

2.50
1.05

5.00
1.50

Here the logical implications (3.10) was represented by means of the big-M
method,− M (1 − zp ) ≤ sp ≤ M (1 − zp )(∀p ∈ [P ]),(3.15)

where M was set to 1000 in all instances of our experiments. Similarly, theimplications (3.8c) were represented with the big-M method,− M zp ≤ ap ≤ M zp(∀p ∈ [P ]).

(3.16)

However, we had difficulty in finding a unified value of M for constraints (3.16)
such that MISDO computation was not aborted due to numerical instability.
Hence, we tuned the values of M through preliminary experiments as shown in
Table 3.2.
These computations were performed on a Linux computer with an Intel Core2
Quad CPU (2.66 GHz) and 4 GB memory. MISDO problems were solved by using
SCIP-SDP1 2.0.0 [61] combined with SCIP2 3.2.0 [1] and SDPA3 7.3.8 [168].
The computation for solving the MISDO problem was terminated if it did notfinish by itself within 10000 s. In this case, the best feasible solution obtained
within 10000 s was taken as the result. Chatterjee and Hadi [42] mention thatwhen the value of the condition number exceeds 225, the deleterious effects ofmulticollinearity in the data become strong. Hence, the upper bound on thecondition number was set as κ = 100 and 225.
Remark 3.1. We can use the cutting-plane and branch-and-cut algorithms in
Chapter 2 for solving MISDO and MISDONE . However, the study in this chapter
was conducted chronologically prior to that in Chapter 2, and thus, we did notused them in this experiment.
1

http://www.opt.tu-darmstadt.de/scipsdp/
http://scip.zib.de/3
http://sdpa.sourceforge.net/
2

34

Table 3.3: Results of solving MISDO problems (κ = 100)
R2

Cond

|S|

Time (s)MISDOMISDONE

0.758770.75877

78.4
78.3

15
15

60.19
17.74

MISDO
MISDONE

0.874290.87430

76.9
92.8

19
21

>10000.00
5563.34

Instance

N

P

Method

Servo

167

19

AutoMPG

392

25

Table 3.4: Results of solving MISDO problems (κ = 225)
R2

Cond

|S|

Time (s)MISDOMISDONE

0.758770.75877

97.1
104.1

15
15

4.99
5.49

MISDO
MISDONE

0.874380.87438

142.0
184.3

21
22

>10000.00
336.14

Instance

N

P

Method

Servo

167

19

AutoMPG

392

25

Tables 3.3 and 3.4 show the computational results of solving MISDO problems with κ = 100 and 225, respectively. The results for only the small-sized
instances Servo and AutoMPG are shown in the tables because many of the largescale MISDO problems could not be solved because of numerical instability4 (i.e.,
violation of Slater’s condition). The column labeled “R2 ” shows the value of thecoefficient of determination of a subset regression model built by each method,
i.e.,
P
P
2
(y
−
n
n∈[N
]
p∈S ap xnp )
2
P
R =1−
,
2
n∈[N ] (yn − ȳ)Pwhere ȳ := ( n∈[N ] yn )/N . Note here that the largest R2 values for each instanceare indicated in bold. The column labeled “Cond” shows the condition number ofthe correlation sub-matrix RS , and the column labeled “|S|” shows the numberof selected explanatory variables. The column labeled “Time (s)” shows thecomputation time in seconds. Note that the computation of the cutting-planealgorithm was terminated if it did not finish by itself within 10000 s. In suchcases, the best feasible solution obtained within 10000 s was used as the result,and “N/A” means that no feasible solution was found.
Tables 3.3 and 3.4 show that all the MISDO problems for Servo were solved
within 61 s, and they were solved faster when κ = 225 than when κ = 100.
Moreover, the normal-equation-based constraints worked effectively in speedingup the computations. For instance, in Table ,3.3 the computation time of theMISDO formulation was reduced from 60.19 s to 17.74 s by incorporating thenormal-equation-based constraints into the problem.
In the case of AutoMPG, only the computations of MISDONE finished within
10000 s for both κ = 100 and 225. Furthermore, MISDONE attained the largest
R2 value for every instance in Tables 3.3 and 3.4. These results demonstrate theeffectiveness of the normal-equation-based constraints in the MISDO formulation
4
Using SCIP-SDP 2.1.0 instead of 2.0.0 and softening feasibility tolerances of SCIP and
SCIP-SDP solvers (feastol and sdpsolverfeastol, respectively) from 10−6 (default) to 10−4
could resolve numerical issues, but it missed true optimal solutions to some of the instances.
Hence, we used SCIP-SDP 2.0.0 with its default parameters.

35

for best subset selection to eliminate multicollinearity.
3.5.2

Computational performance of cutting-plane algorithms

Next, we compare the computational performance of the cutting-plane algorithms
with that of conventional local search algorithms for subset selection. The algorithms used in the comparison are listed below:
FwS

Forward selection method: Starts with S = ∅ and iteratively adds thevariable p (i.e., S ← S ∪{p}) that leads to the largest decrease in RSS (3.1);
this operation is repeated while cond(RS ) ≤ κ is satisfied.

BwE Backward elimination method: Starts with S = [P ] and iteratively eliminates the variable p (i.e., S ← S \ {p}) that leads to the smallest increase
in RSS (3.1); this operation is repeated until cond(RS ) ≤ κ holds.
CPA

Cutting-plane algorithm that omits Step 2 (Backward elimination) andappends cut (3.13).

CPABwE Cutting-plane algorithm that appends cut (3.14) strengthened bymeans of the backward elimination method.
These computations were performed on a Linux computer with an Intel Xeon
W3520 CPU (2.66 GHz) and 12 GB memory. The algorithms FwS and BwE
were implemented in MATLAB5 R2013a ; the algorithms CPA and CPABwE
were implemented in Python 2.7.3, with Gurobi Optimizer6 5.6.0 used to solve
relaxed MIQO problems (3.12a) and (3.12b). Here the logical implications (3.5b)
were incorporated in the form of SOS1 constraints, which imply that at most one
element in the set can have a nonzero value. Specifically, the SOS1 constraint
is imposed on {1 − zp , ap } (p ∈ [P ]) with the SOS type 1 function implemented
in Gurobi Optimizer. Therefore, if zp = 0, then 1 − zp has a nonzero value andap must be zero from the SOS1 constraints. The upper bound on the condition
number was set as κ = 100 or 225.
Tables 3.5 and 3.6 show the computational results obtained using the fouralgorithms with κ = 100 and 225, respectively. The column labeled “#Iter”
shows the number of iterations in the cutting-plane algorithms. Note that thecomputation of the cutting-plane algorithm was terminated if it did not finish
by itself within 10000 s. In such cases, the best feasible solution obtained within
10000 s was used as the result and “N/A” means that no feasible solution was
found.
We can see from Tables 3.5 and 3.6 that the local search algorithms FwS
and BwE finished their computations within 2 s for all the instances. Theirobtained solutions were, however, frequently inferior to those of the cuttingplane algorithms, especially when κ = 100 (see Table 3.5); for the BreastCancer
instance, CPABwE gave an R2 value of 0.29040, whereas FwS and BwE gave R2
values of 0.27580 and 0.26572, respectively.
We can also see that the computation finished much faster for CPABwE thanfor CPA. The main reason for this is that the number of iterations required for
CPABwE was significantly reduced by the strengthened cut (3.14). Indeed, in thecase of SolarFlareC in Table 3.5, CPA arrived at an optimal solution in 1246.97 s
after 4209 iterations, whereas CPABwE terminated with an optimal solution in
84.59 s after 331 iterations.
5
6

http://www.mathworks.com/products/matlab/
http://www.gurobi.com

36

Table 3.5: Results of local search algorithms and cutting-plane algorithms (κ =100)
Instance

N

P

Method

R2

Cond

|S|

#Iter

Time (s)

Servo

167

19

FwS
BwE
CPA
CPABwE

0.758620.758770.758770.75877

63.2
39.0
38.4
39.9

14
15
15
15

—
—
76
16

0.09
0.04
0.93
0.42

AutoMPG

392

25

FwS
BwE
CPA
CPABwE

0.87335
0.874290.87430
0.87430

87.6
86.2
93.2
92.8

21
19
21
21

—
—
1284
84

0.17
0.09
325.20
14.11

SolarFlareC

1066

26

FwS
BwE
CPA
CPABwE

0.197130.19705
0.19715
0.19715

4.3
2.0
18.7
34.3

19
15
19
19

—
—
4209
331

0.23
0.87
1246.97
84.59

BreastCancer

194

32

FwS
BwE
CPA
CPABwE

0.2758091.8
0.2657261.8
N/A
0.29040
95.8

15
9

—
—
>4364
>1044

0.17
0.27
>10000.00>10000.00

FwS
BwE
CPA
CPABwE

0.1639999.9
0.1648153.1
N/A
0.16537
96.5

52
50

—
—
>5901
53

1.68
1.32
>10000.00
301.07FwS
BwE
CPA
CPABwE

0.9644799.9
0.96571
67.3
N/A
0.96908
72.9

27
24

—
—
>6960
>277

0.76
1.89
>10000.00>10000.00

ForestFires

Automobile

517

159

63

65

14

59

25

Note that when the computation is terminated due to the time limit of 10000 s,
CPA does not provide a feasible solution because it is found only at the end of thealgorithm. In contrast, CPABwE can find a feasible solution of good quality in theearly stage of the algorithm by means of the backward elimination method. For
this reason, CPABwE always provided the best solution among the four methods
for all the instances in Tables 3.5 and 3.6.
Comparing the MISDO approaches and the cutting-plane algorithms, it needs
to be said that the computational performance of the MISDO approaches wasmuch lower than that of the cutting-plane algorithms. For instance in the case of
AutoMPG with κ = 225, MISDONE took 336.14 s to solve the problem (Table 3.4),
but CPABwE required only 1.76 s to solve the same problem (Table 3.6). These
results imply that our MISDO approach was not effective for the subset selection
problem at the current moment; however, since the computational performance ofMISDO algorithms is being improved, our MISDO formulation will work betterin the future.

37

Table 3.6: Results of local search algorithms and cutting-plane algorithms (κ =225)
Instance

N

P

Method

R2

Cond

|S|

#Iter

Time (s)

Servo

167

19

FwS
BwE
CPA
CPABwE

0.758770.758770.758770.75877

146.3
39.0
102.5
39.9

15
15
15
15

—
—
56
15

0.11
0.04
0.62
0.40

AutoMPG

392

25

FwS
BwE
CPA
CPABwE

0.874380.874380.874380.87438

185.3
181.1
157.1
173.2

22
22
22
22

—
—
60
18

0.18
0.05
2.04
1.76

SolarFlareC

1066

26

FwS
BwE
CPA
CPABwE

0.197130.19705
0.19715
0.19715

4.3
2.0
18.7
169.4

19
15
19
19

—
—
4209
750

0.22
0.48
1244.26
189.48

BreastCancer

194

32

FwS
BwE
CPA
CPABwE

0.30010 217.9
0.26572
61.8
N/A
0.30513 215.6

19
9
20

—
—
>4369
287

0.21
0.27
>10000.00
288.59ForestFires

517

63

FwS
BwE
CPA
CPABwE

0.16556
0.16556
0.16556
0.16556

214.2
212.6
209.5
209.0

60
60
60
60

—
—
59
12

1.74
0.42
7.60
27.37

Automobile

159

65

FwS
BwE
CPA
CPABwE

0.97124 224.4
0.97153 183.5
N/A
0.97391 224.4

39
29

—
—
>6973
>960

1.14
1.85
>10000.00>10000.00

3.6

36

Conclusion

This chapter addressed the problem of selecting the best subset of explanatory
variables subject to an upper bound on the condition number for eliminating
multicollinearity from linear regression models. The first contribution of thisstudy is a novel computational framework for eliminating multicollinearity basedon the MISDO formulation. This framework reformulates the subset selection
problem as a single MISDO problem. The second contribution is the establishment of a high-performance cutting-plane algorithm. In this algorithm, we
generate strong cutting-planes with backward elimination and reduce the numberof mixed-integer quadratic optimization problems to be solved. While numericalexperiments show that our MISDO formulation could only be applied to smallsized instances at present, this chapter provides a new statistical application of
mixed-integer semidefinite optimization formulation. Moreover, we found that
our cutting-plane algorithm frequently provided a better subset of variables than
did the common local search algorithms. This finding demonstrates the effectiveness of the specialized cutting-plane algorithm in eliminating multicollinearity
from a linear regression model.
38

Chapter 4
Cutting-plane Algorithm for
Cardinality-constrained Distributionally
Robust Portfolio Optimization
This chapter focuses on portfolio selection problems. In particular, we treat
a distributionally robust portfolio optimization model with a cardinality constraint, which limits the number of invested assets. This problem is formulated
as a MISDO problem, and thus, solving it exactly is computationally challenging
when the number of investable assets is large. To overcome this issue, we propose
a specialized cutting-plane algorithm to solve the cardinality-constrained distributionally robust optimization problem. We first reformulate the problem as a
bilevel optimization problem and design a cutting-plane algorithm for solving theupper-level problem. Then, to generate cutting planes efficiently, we apply thetechnique of positive semidefinite matrix completion to the lower-level problem
and reduce its problem size. The content of this chapter is included in Kobayashi
et al. [96].
We give an introduction to this study and summarize our contribution in
Section 4.1 In Section 4.2, we give a formulation of the cardinality-constrained
distributionally robust portfolio optimization problem. In Section 4.3, we explain
our cutting-plane algorithm for solving the problem. In Section 4.4, we describe
the reduction of the lower-level SDO problem. We report computational resultsin Section 4.5 and conclude in Section 4.6.

4.1

Introduction

Portfolio optimization models, originating from the mean-variance portfolio selection pioneered by Markowitz [114], have been actively studied by academic
researchers and institutional investors. In real-world portfolio optimization problems, we must depend on inaccurate estimates of asset returns, which can lead tolower investment performance. Accordingly, robust optimization [18, 20], whichcopes with uncertainty about input data, has played an important role in many
practical situations [48, 53, 72]. This chapter focuses on a distributionally robust
optimization model developed by Delage and Ye [47] for portfolio optimization.
In distributionally robust optimization, optimal solutions are evaluated under
the worst-case expectation with respect to a set of probability distributions of
uncertain parameters. This model was first introduced by Scarf [142] for an inventory problem. Distributionally robust optimization can be viewed as a framework
unifying traditional stochastic programming based on sample average approximation [93, 144, 145] and conventional robust optimization based on uncertainty setsof possible realizations of random variables [18, 20]. Moreover, some theoretical39

analyses support the effectiveness of distributionally robust optimization in datadriven decision-making under uncertainty [69, 126].
Moment-based uncertainty sets, which are sets of probability distributions ofasset returns whose moments satisfy certain conditions, are commonly adopted in
distributionally robust portfolio optimization [51, 66, 107, 131, 160]. El Ghaoui
et al. [51] considered minimizing the worst-case value-at-risk over an uncertainty
set of probability distributions such that the mean and covariance matrix are
elementwise bounded. Tütüncü and Koeing [160] used a similar uncertainty set
for solving robust mean-variance portfolio optimization problems. Goldfarb andIyengar [66] studied a distributionally robust portfolio optimization problem inwhich asset returns are formed by a linear factor model. Popescu [131] assumed
that the mean and covariance matrix are exactly known and employed an uncertainty set of probability distributions whose first and second moments are equal
to them. Similar assumptions were also made in other distributionally robust
optimization models [122, 175]. Goh and Sim [65] studied an uncertainty setof probability distributions whose mean belongs to a convex set, and covariancematrix is given. As an extension of the elementwise uncertainty sets [51, 160],
Lotfi and Zenios [107] considered an ellipsoidal uncertainty set of the tuple ofmean and covariance matrix when minimizing the worst-case value-at-risk andconditional value-at-risk.
Unlike these models that require prior knowledge of the region including
the mean and/or covariance matrix, Delage and Ye [47] proposed a data-driven
method for defining a moment-based uncertainty set and gave its probabilistic
guarantee. They also formulated the worst-case loss minimization problem overthis uncertainty set as an SDO problem, which can be solved in polynomial timeby interior-point methods [123]. Bertsimas et al. [29] proposed a practical framework based on a bootstrap hypothesis test to determine the uncertainty set ofprobability distributions from historical data. See Rahimian and Mehrotra [133]for a comprehensive survey on distributionally robust optimization models.
This chapter focuses on solving distributionally robust portfolio optimization
problems based on the moment-based uncertainty set [47] with a cardinality constraint to limit the number of invested assets. In real-world practice, when thenumber of invested assets is large, it is difficult for investors to keep track of each
asset, and substantial transaction costs are required [113, 128]. Thus, there is aneed to control the number of invested assets by introducing the cardinality constraint. However, solving such a cardinality-constrained distributionally robustmodel is computationally challenging. Due to the cardinality constraint and themoment-based uncertainty set, the proposed model is formulated as a MISDO
problem.
Recently, Bertsimas et al. [26] proposed a general framework of cutting-planealgorithms to exactly solve mixed-integer convex optimization problems with logical constraints. They reformulated this problem as a bilevel optimization problem
composed of lower- and upper-level problems. To solve the upper-level problem,
they devised a cutting-plane algorithm, which iteratively approximates the objective function by generating cutting planes from the solution to the lower-levelproblem with the strong duality theory. Bertsimas and Cory-Wright [24] applied
this algorithm to the cardinality-constrained mean-variance portfolio optimization. They demonstrated that their algorithm was much faster than state-of-theart MIO methods when solving large-scale problem instances. The cutting-planealgorithm was also applied to sparse principal component analysis [25] and sparse
inverse covariance estimation [33].
Motivated by these prior studies, we propose a specialized cutting-plane algo40

rithm to exactly solve the cardinality-constrained distributionally robust portfolio
optimization problem. At each iteration of this cutting-plane algorithm, we mustsolve an SDO problem, which is the dual of the lower-level problem for generating
cutting planes. However, the size of this SDO problem depends on the numberof investable assets, so solving this problem at each iteration is computationallyprohibitive when handling a large number of assets. To overcome this difficulty,
we reduce the size of the lower-level SDO problem through the application of
the matrix completion technique [60, 121]. Notably, the size of the reduced SDO
problem depends not on the number of investable assets but on the cardinality
parameter, which is usually set to a small number. To the best of our knowledge,
we are the first to develop an effective algorithm for exactly solving cardinalityconstrained distributionally robust portfolio optimization problems. Numericalexperiments using real-world datasets [16, 41, 59, 90] demonstrate the effectiveness of our method in terms of both computational and out-of-sample investment
performances.

4.2

Problem formulation

In this section, we formulate the cardinality-constrained distributionally robust
portfolio optimization problem that we consider in this chapter.
4.2.1

Cardinality constraint

Let x := (x1 , x2 , . . . , xN )> be a portfolio, where xn is the investment weight of
the nth asset. Throughout this paper, we consider the set of feasible portfolios:




X
N
X := x ∈ R
xn = 1, x ≥ 0 .
n∈[N ]The nonnegativity constraint on x prohibits short selling.
Let k ∈ [N ] be a user-defined parameter for limiting the cardinality (i.e.,
the number of assets to be held). We then impose the following cardinality
constraint [27, 34, 128] on portfolio x:
kxk0 ≤ k,(4.1)

where k · k0 is the `0 -norm (i.e., the number of nonzero entries). In practice,
this constraint is required by investors to reduce their portfolio monitoring andtransaction costs.
Let z := (z1 , z2 , . . . , zN )> be a vector composed of binary decision variables
for selecting assets; that is, zn = 1 if the nth asset is selected, and zn = 0
otherwise. We also introduce the feasible set corresponding to the cardinality
constraint (4.1):




X
k
ZN
:= z ∈ {0, 1}N
zn = k .
n∈[N ]

Then, the cardinality constraint (4.1) is represented by the following logical implication:zn = 0 ⇒ xn = 0 (∀n ∈ [N ]),(4.2)
kz ∈ ZN
.

(4.3)

41

4.2.2

Moment-based uncertainty set

We consider a measurable space (Ω, F). Let ξ̃ : Ω → RN be an F-measurable
function (N -dimensional random vector) representing the rate of random returnof each asset. Suppose that ξm ∈ RN is an observed historical return for eachperiod m ∈ [M ]. Then, the sample mean vector and the sample covariance matrix
are calculated as
1 X
1 X
µ̂ :=
ξm ,
Σ̂ :=
(ξm − µ̂)(ξm − µ̂)> .
(4.4)
M
M
m∈[M ]m∈[M ]Let M be the set of all probability measures in the measurable space (Ω, F),
and EF [ · ] be the expectation under the probability measure F ∈ M. Delage andYe [47] considered the moment-based uncertainty set of probability distributionsof ξ̃ with the assumption Σ̂  O. When the support of distributions is RN , thisuncertainty set is given by



>




EF [ξ̃] − µ̂ Σ̂−1 EF [ξ̃] − µ̂ ≤ κ1
, (4.5)D(µ̂, Σ̂, κ1 , κ2 ) := F ∈ ME [(ξ̃ − µ̂)(ξ̃ − µ̂)> ]  κ Σ̂2F

where κ1 ≥ 0 and κ2 ≥ 1 are user-defined uncertainty parameters about µ̂ andΣ̂, respectively. The first condition in Eq. (4.5) ensures that the expectation ofξ̃ lies in an ellipsoid of size κ1 centered at µ̂. The second condition in Eq. (4.5)
implies that the second central-moment matrix of ξ̃ is bounded above by κ2 Σ̂ inthe sense of the matrix inequality.
Throughout this chapter, we make the following mild assumption about the
moment-based uncertainty set (4.5).
Assumption 4.1. κ1 > 0 and Σ̂  O.
4.2.3Piecewise-linear utility and loss functions

Let ξ ∈ RN be a realization of the random return ξ̃. Following Delage and Ye [47],
we employ the following utility function, which is piecewise-linear and concave inthe portfolio net return ξ > x asnou(x, ξ) := min a(`) ξ > x + b(`) ,
(4.6)
`∈[L]

where a(`) ∈ R and b(`) ∈ R are respectively the slope and intercept of the `thlinear function for ` ∈ [L]. The loss function is then defined as the negative of
the utility function (4.6):nonoL(x, ξ) := − min a(`) ξ > x + b(`) = max −a(`) ξ > x − b(`) .
(4.7)
`∈[L]

4.2.4

`∈[L]

Portfolio optimization model

Our objective is to minimize the following worst-case expected loss with respect
to an underlying probability distribution F ∈ D(µ̂, Σ̂, κ1 , κ2 ):
n h
io
max
EF L(x, ξ̃) .
(4.8)
F ∈D(µ̂,Σ̂,κ1 ,κ2 )The `2 -regularization term is also incorporated into the objective from a perspective of robust optimization [24, 48, 70, 71].
42

We are now in a position to formulate our cardinality-constrained distributionally robust portfolio optimization model as follows:
n h
io
1 >
x x+
max
EF L(x, ξ̃)2γF ∈D(µ̂,Σ̂,κ1 ,κ2 )minimizex,zsubject to zn = 0 ⇒ xn = 0
x ∈ X,(∀n ∈ [N ]),kz ∈ ZN
,(4.9a)(4.9b)(4.9c)

where γ > 0 is a user-defined regularization parameter. By deriving the dual of
the inner maximization problem (4.8) as in Delage and Ye [47], Problem (4.9)
can equivalently be reformulated as the following MISDO problem:
1 >x x + (κ2 Σ̂ − µ̂µ̂> ) • Q + r + Σ̂ • P − 2µ̂> p + κ1 s (4.10a)x,P ,Q,p,q,r,s,z2γ
subject to p = −q/2 − Qµ̂,
(4.10b)(`)
Q
q/2 + a x/2
 O (∀` ∈ [L]), (4.10c)(`)>(q/2 + a x/2)
r + b(`)


P p
 O,
(4.10d)p> sminimizezn = 0 ⇒ x n = 0(∀n ∈ [N ]),(4.10e)kz ∈ ZN,x ∈ X,(4.10f)where P , Q ∈ S N , p, q ∈ RN , and r, s ∈ R are dual decision variables of Problem (4.8).

4.3

Cutting-plane algorithm

In this section, we give our cutting-plane algorithm for solving the cardinalityconstrained distributionally robust portfolio optimization problem (4.10).
4.3.1

Bilevel optimization reformulation

We extend the method of bilevel optimization reformulation [24] to the cardinalityconstrained distributionally robust portfolio optimization problem (4.10). Let us
denote by Z := Diag(z) a diagonal matrix whose diagonal entries are given by z.
The logical implication (4.10e) are then incorporated by replacing x with Zx in
Problem (4.10). Now, we can reformulate Problem (4.10) as a bilevel optimization problem. Specifically, the upper-level problem is written as the following
integer optimization problem:
minimizezf (z)subject to43kz ∈ ZN
,(4.11)

and the lower-level problem for calculating the objective function is expressed as
the following SDO problem:
f (z) = minimize

x,P ,Q,p,q,r,s

1 >x x + (κ2 Σ̂ − µ̂µ̂> ) • Q + r + Σ̂ • P − 2µ̂> p + κ1 s
2γ
(4.12a)subject to p = −q/2 − Qµ̂,Q
q/2 + a(`) Zx/2
O
(q/2 + a(`) Zx/2)>
r + b(`)(4.12b)(∀` ∈ [L]),(4.12c)




p
 O,
s

P
p>

(4.12d)Zx ∈ X .

(4.12e)Note that (Zx)> Zx has been replaced by x> x in Eq. (4.12a) because x = Zx
holds for x that minimizes Eq. (4.12a) as in Bertsimas and Cory-Wright [24].
We then derive cutting planes that underestimate f (z) in the upper-levelproblem (4.11) by exploiting the dual formulation of the lower-level problem(4.12). The following theorem yields the dual formulation of Problem (4.12),
where ω ◦ ω denotes the Hadamard product of the vector ω.
k , the strong duality holds for Problem (4.12), andTheorem 4.2. For all z ∈ ZN
the dual formulation of Problem (4.12) is represented as follows:f (z) = maximize

ω,B,β,η,λ,π

subject toXγ >z (ω ◦ ω) −η (`) b(`) + π
2
`∈[L]
X
(`) (`)
ω≥a β + π1,−(4.13a)(4.13b)

`∈[L]

>


X

B (`) = κ2 Σ̂ − µ̂µ̂> + µ̂ 

`∈[L]

X
`∈[L]



β (`)  + X

β (`)  µ̂> ,`∈[L](4.13c)Xβ(`)= λ + µ̂,(4.13d)η (`) = 0,

(4.13e)

`∈[L]

1−

X
`∈[L]

!
B (`)
β (`)
>
O
β (`)
η (`)


Σ̂ λ
 O,
λ> κ1

(∀` ∈ [L]),(4.13f)(4.13g)where ω ∈ RN , B := (B (`) ) ∈ RN ×N ×L , β := (β (`) ) ∈ RN ×L , η := (η (`) ) ∈ RL ,λ ∈ RN , and π ∈ R are dual decision variables.
Proof. See Appendix A.1.
According to Theorem 4.2, we can extend the definition of f (z) to the optimal
objective value of Problem (4.13) for real-valued z ∈ [0, 1]N . Then, we obtain two
key properties of f (z) as in Bertsimas and Cory-Wright [24].

44Lemma 4.3 (Convexity). The function f (z) is convex in z ∈ [0, 1]N .
Proof. Since f (z) is a pointwise maximum of linear functions in z ∈ [0, 1]N ,its epigraph, which is the intersection of epigraphs of the linear functions, isconvex.
Lemma 4.4 (Subgradient). Let ω ? (z) be an optimal solution of ω to Problem (4.13) with z ∈ {0, 1}N . Then, a subgradient of the function f (z) is givenbyγg(z) := − ω ? (z) ◦ ω ? (z) ∈ ∂f (z).
(4.14)
2
Proof. See, for example, Proposition 8.1.1 in Bertsekas et al. [23].
k , Problem (4.13) with z = ẑLemmas 4.3 and 4.4 verify that for each ẑ ∈ ZN
yields a linear underestimator of f (z) for z ∈ [0, 1]N as follows:f (z) ≥ f (ẑ) + g(ẑ)> (z − ẑ).
4.3.2

(4.15)

Algorithm description

We revise the cutting-plane algorithm [24] with the aim of solving the upperlevel problem (4.10). Let θLB be a lower bound of the optimal objective valueof Problem (4.11); this bound can easily be calculated by solving a continuous
relaxation version of Problem (4.11). Our cutting-plane algorithm starts withthe following initial feasible region:k
F1 := {(z, θ) ∈ ZN
× R | θ ≥ θLB },(4.16)

where θ is an auxiliary decision variable that corresponds to a lower estimate off (z).
At the tth iteration (t ≥ 1), our algorithm solves the surrogate upper-levelproblem:minimize θ subject to (z, θ) ∈ Ft ,(4.17)z,θ

where Ft is a feasible region at the tth iteration such that Ft+1 ⊆ Ft . Eq. (4.16)
ensures that the objective value of Problem (4.17) is bounded below, and thus,
there exists an optimal solution (zt , θt ) to Problem (4.17).
We next solve the dual lower-level problem (4.13) with z = zt . We thus
obtain the function value f (zt ) and its subgradient g(zt ) through Eq. (4.14). If
f (zt ) − θt ≤ ε with sufficiently small ε ≥ 0, then zt is an ε-optimal solution to
Problem (4.11), which means thatf ? ≤ f (zt ) ≤ f ? + ε,

(4.18)

where f ? is the optimal objective value of Problem (4.11). In this case, we
terminate the algorithm with the ε-optimal solution zt . Otherwise, we add theconstraint (4.15) with ẑ = zt to the feasible region as follows:k
Ft+1 ← Ft ∩ {(z, θ) ∈ ZN
× R | θ ≥ f (zt ) + g(zt )> (z − zt )}.

(4.19)

Note that this update cuts off the solution (zt , θt ) because θt < f (zt ).
We then set t ← t+1 and solve the surrogate upper-level problem (4.17) again
with the updated feasible region (4.19). We repeat this procedure until we findan ε-optimal solution ẑ. After termination of the algorithm, we can compute thecorresponding portfolio by solving the lower-level problem (4.12) with z = ẑ.
Our cutting-plane algorithm is summarized by Algorithm 4.1. we can prove
the finite convergence of Algorithm 4.1 as follows.
45

Algorithm 4.1 Cutting-plane algorithm for solving Problem (4.11)
Step 0 (Initialization) Let ε ≥ 0 be a tolerance for optimality. Define thefeasible region F1 as in Eq. (4.16). Set t ← 1 and UB0 ← ∞.
Step 1 (Surrogate Upper-level Problem) Solve Problem (4.17). Let (zt , θt )be an optimal solution, and set LBt ← θt .
Step 2 (Dual Lower-level Problem) Solve Problem (4.13) with z = zt tocalculate f (zt ) and ω ? (zt ). If f (zt ) < UBt−1 , set UBt ← f (zt ) andẑ ← zt ; otherwise, set UBt ← UBt−1 .
Step 3 (Termination Condition) If UBt − LBt ≤ ε, terminate the algorithmwith the ε-optimal solution ẑ.
Step 4 (Cut Generation) Calculate g(zt ) as in Eq. (4.14) and update thefeasible region as in Eq. (4.19). Set t ← t + 1 and return to Step 1.

Theorem 4.5. Algorithm 4.1 terminates in a finite number of iterations andoutputs an ε-optimal solution to Problem (4.11).
Proof. Let {(zt , θt ) | t ∈ [T ]} be a sequence of solutions generated by Algorithm 4.1. Suppose that there exists u < T such that zu = zT . Since(zT , θT ) ∈ Fu+1 , it follows from Eq. (4.19) that
LBT = θT ≥ f (zu ) + g(zu )> (zT − zu ) = f (zT ),which verifies that zT is an optimal solution to Problem (4.11). Since there arek , the algorithm terminates
at most a finite number of possible solutions z ∈ ZN
with an ε-optimal solution after a finite number of iterations.

4.4

Problem reduction of the dual lower-level problem

Recall that Step 2 of Algorithm 4.1 solves Problem (4.13), which is an SDO
problem including positive semidefinite constraints (4.13f) and (4.13g) on (N +1) × (N + 1) symmetric matrices. It is clearly difficult to directly solve Problem(4.13) when N is very large. To remedy this situation, we reduce its problem size
by applying the technique of positive semidefinite matrix completion [60, 121] tothe lower-level SDO problem (4.13).
4.4.1

Reduced problem formulation

For any vector v := (vn ) ∈ RN and matrix M := (mnn0 ) ∈ RN ×N , we write thesubvector and submatrix corresponding to z, z 0 ∈ {0, 1}N as
vz := (vn )n∈N (z) = (vn | zn = 1) ∈ R|N (z)| ,0

Mz,z0 := (Mnn0 )(n,n0 )∈N (z)×N (z0 ) = (Mnn0 | zn = zn0 = 1) ∈ R|N (z)|×|N (z )| ,where N (z) := {n ∈ [N ] | zn = 1}.
k , we have |N (z)| = k. Then ω ∈ RN can be replaced with its
When z ∈ ZN
subvector ωz ∈ Rk in the objective (4.13a) as follows:γ >γz (ω ◦ ω) = ωz> ωz .
22

46

(4.20)

The cardinality parameter k is usually much smaller than N . We exploit thisproblem structure to reduce the size of Problem (4.13).
k is formulated as
A reduced version of Problem (4.13) for z ∈ ZN
f 0 (z) =Xγ >η (`) b(`) + π
ωz ωz −
ωz ,Bz,z ,βz ,η,λz ,π
2
`∈[L]
X
(`)
a(`) βz + π1,
subject to ωz ≥maximize−(4.21a)(4.21b)

`∈[L]

>


(`)

X


Bz,z = κ2 Σ̂z,z − µ̂z µ̂>z + µ̂zX(`)

βz  + 

`∈[L]

`∈[L]




X

`∈[L]

(4.21c)βz = λz + µ̂z ,(`)(4.21d)X(4.21e)

X
`∈[L]

1−

η (`) = 0,

`∈[L]




(`)
(`)
Bz,z
βz


O
(`) >
βz
η (`)


Σ̂z,z λz
 O,
λ>
κ1
z
(`)

(` ∈ [L]),(4.21f)(4.21g)(`)

where ωz ∈ Rk , Bz,z := (Bz,z ) ∈ Rk×k×L , βz := (βz ) ∈ Rk×L , and λz ∈ Rk
are reduced versions of decision variables.
In the next subsection, we verify that the reduced problem (4.21) is equivalent to the original problem (4.13) in the sense that an optimal solution to theoriginal problem (4.13) can be recovered from an optimal solution to the reduced
problem (4.21). We also prove that f 0 (z) defined by the reduced problem (4.21)
is equal to f (z) defined by the original problem (4.13).
4.4.2

Equivalence of the original and reduced problems

We first prove the following lemma.
Lemma 4.6. Let (ωz , Bz,z , βz , η, λz , π) be a feasible solution to Problem (4.21)
k . We also define β̄ := (β̄ (`) ) ∈ RN ×L as
for z ∈ ZN
( (`)
(`)
β̄z := βz ,
(`)

(`)

β̄1−z := Σ̂1−z,z (Σ̂z,z )−1 (βz − η (`) µ̂z ) + η (`) µ̂1−z(∀` ∈ [L]).

(4.22)

It then follows that
κ2 Σ̂ X`∈[L]
η (`) >01η (`)(β̄ (`) − η (`) µ̂)(β̄ (`) − η (`) µ̂)> .

(4.23)

Proof. See Appendix A.2.
We next prove that feasible B and β for the original problem (4.13) can becompleted from a feasible solution to the reduced problem (4.21).

47

(`)

βz  µ̂>z,Lemma 4.7. Let (ωz , Bz,z , βz , η, λz , π) be a feasible solution to Problem (4.21)
k , and β̄ is defined by Eq. (4.22). Then, there exists B̄ := (B̄ (`) ) ∈for z ∈ ZN
N
×N
×L
R
such that (B, β) = (B̄, β̄) satisfies Eqs. (4.13c) and (4.13f).
(`)

Proof. For ` ∈ [L] such that η (`) = 0, we have βz = 0 from the constraint(4.21f), thus β̄ (`) = 0 from the definition (4.22). In this case, the solution(B̄ (`) , β̄ (`) , η (`) ) = (O, 0, 0) satisfies the constraint (4.13f) and can be left outof Eq. (4.13c). Therefore, we can assume without loss of generality that η (`) > 0
for all ` ∈ [L] in Eqs. (4.13c) and (4.13f).
From Eqs. (4.21e) and (4.23), we can see that

> 

X
X
X
e (`) ,
κ2 Σ̂ − µ̂µ̂> + µ̂ β̄ (`)  + β̄ (`)  µ̂> B(4.24)
`∈[L]

where

`∈[L]

e (`) := 1 β̄ (`) (β̄ (`) )>
B
η (`)

`∈[L]

(∀` ∈ [L]).

Then we have (`)   (`) >

e (`)
1
β̄
β̄
B
β̄ (`)
= (`)
O
η (`)
η (`)
(β̄ (`) )> η (`)η(∀` ∈ [L]).

(4.25)

From Eqs. (4.24) and (4.25), we can see that B̄ satisfying Eqs. (4.13c) and (4.13f)
e (`) for ` ∈is obtained by adding appropriate positive semidefinite matrices to B
[L].
Now, we can prove our main theorem.
Theorem 4.8. Let (ωz , Bz,z , βz , η, λz , π) be an optimal solution to Problem
k , and (ω̄, β̄, λ̄) be defined by Eq. (4.22) and(4.21) for z ∈ ZN
ω̄z := ωz ,hPi
(4.26)(`) β̄ (`) + π1ω̄1−z :=
a
,
1−z
`∈[L]
+
(
λ̄z := λz ,(4.27)
λ̄1−z := Σ̂1−z,z (Σ̂z,z )−1 λz ,where [v]+ := (max{0, vn }) ∈ RN for v := (vn ) ∈ RN . Then, there exists B̄ suchthat (ω̄, B̄, β̄, η, λ̄, π) is an optimal solution to Problem (4.13). In addition, wek.
have f (z) = f 0 (z) for all z ∈ ZN
Proof. We assume without loss of generality that 1z=
, 1 ∈ Rk , 0 ∈ RN −k .
0(4.28)

Accordingly, we partition vector v ∈ RN and symmetric matrix M ∈ S N fornotational simplicity as follows:
 
v1
v=
, v1 = vz , v2 = v1−z ,(4.29)v2

M11 M12>
M=, M11 = Mz,z , M12 = M21= Mz,1−z , M22 = M1−z,1−z .
M21 M22(4.30)
48

Suppose that (ω, B, β, η, λ, π) is an optimal solution to Problem (4.13). Then,(ω1 , B11 , β1 , η, λ1 , π) is a feasible solution to Problem (4.21), and its objective
value is equal to f (z) due to Eq. (4.20). This proves that f (z) ≤ f 0 (z).
Next, we show that f (z) ≥ f 0 (z). Let (ω1 , B11 , β1 , η, λ1 , π) be an optimal
solution to Problem (4.21). From the definition (4.26), the constraint (4.13b) issatisfied by (ω̄, β̄, π). Since the constraints (4.21d) and (4.21e) are satisfied, theconstraint (4.13d) is satisfied by (β̄, λ̄) as


P
(`)
X
β
P`∈[L] 1
 ∵ Eqs. (4.21e), (4.22)
β̄ (`) = 
(`)
−1
Σ̂
(
Σ̂
)
β
−
µ̂
+
µ̂
21
11
1
2
1
`∈[L]
`∈[L]


λ1 + µ̂1∵ Eq. (4.21d)=Σ̂21 (Σ̂11 )−1 λ1 + µ̂2= λ̄ + µ̂.

∵ Eq. (4.27)

Since the constraint (4.21g) is satisfied, we can use the following Schur complement property (see, e.g., Section A.5.5 [36]):Σ̂11 λ1−1
 O ⇒ κ1 − λ>(4.31)
1 Σ̂11 λ1  O.
λ>κ1
1
From Eq. (4.27), we have

 



Σ̂11
Σ̂12
λ1
Σ̂11 Σ̂12 λ̄1Σ̂ λ̄= Σ̂21 Σ̂22 λ̄2  = Σ̂21
Σ̂22
Σ̂21 Σ̂−1
11 λ1 ,λ̄> κ1
−1
>
>
>
>
λ̄1
λ̄2 κ1λ1 (Σ̂21 Σ̂11 λ1 )
κ1which is positive semidefinite because




Σ̂11
Σ̂12
λ1
I
O 0
I O −Σ̂−1
11 λ1


O I
O
I 0 Σ̂21
0
Σ̂22
Σ̂21 Σ̂−1
11 λ1
−1
>
>
>
>
−1
>
>
1
−λ1 Σ̂11 0
00
1
λ1 (Σ̂21 Σ̂11 λ1 )κ1Σ̂11 Σ̂120  O. ∵ Eq. (4.31)= Σ̂21 Σ̂22
0
−1
>
>
>
0
0
κ1 − λ1 Σ̂11 λ1
Therefore, the constraint (4.13g) is satisfied by λ̄.
From Lemma 4.7, the constraints (4.13c) and (4.13f) are satisfied by (B̄, β̄).
Therefore, (ω̄, B̄, β̄, η, λ̄, π) is a feasible solution to Problem (4.13), and its objective value is equal to f 0 (z) due to Eq. (4.20). This implies that f (z) ≥ f 0 (z). As
a result, we have f (z) = f 0 (z), and thus, (ω̄, B̄, β̄, η, λ̄, π) is an optimal solution
to Problem (4.13).
Remark 4.1. Note that Eq. (4.26) defines a minimum-norm solution to Problem(4.13). This solution is known to generate strong cutting planes [24].
Remark 4.2. When κ1 = 0, unlike in Assumption 4.1, we can set EF [ξ̃] = µ̂ anddelete the first condition in Eq. (4.5). In this case, we can also prove the theorem
corresponding to Theorem 4.8.
From Theorem 4.8, we can revise Step 2 of Algorithm 4.1 as follows:
Step 2 (Reduced Dual Lower-level Problem) Solve Problem (4.21) with z =zt . Let (ωz , Bz,z , βz , η, λz , π) be an optimal solution. Calculate f (zt ) =f 0 (zt ) and ω ? (zt ) = ω̄ from Eqs. (4.22) and (4.26). If f (zt ) < UBt−1 , set
UBt ← f (zt ) and ẑ ← zt ; otherwise, set UBt ← UBt−1 .
49Remark 4.3. Our cutting-plane algorithm can be extended to a feasible set of
portfolios with additional linear constraints:




X
X 0 := x ∈ RN
xn = 1, x ≥ 0, Cx ≤ d ,n∈[N ]

where C ∈ RJ×N and d ∈ RJ are given constants (see, e.g., Bertsimas andCory-Wright [24] and Kobayashi et al. [95] for details).

4.5

Numerical experiments

In this section, we report numerical results to evaluate the efficacy of our methodfor distributionally robust portfolio optimization with the cardinality constraint.
We first examine the computational performance of our cutting-plane algorithm
and then demonstrate the out-of-sample investment performance of our portfolio
optimization model. All experiments were performed on a Windows 10 PC withan Intel Xeon E-2274G CPU (4.00 GHz) and 32 GB of memory.
4.5.1

Computational performance

In this subsection, we evaluate the computational performance of our cuttingplane algorithm by comparison with other MISDO algorithms.
Experimental design
Table 4.1 lists the datasets used in our experiments, where N is the numberof assets. From the data library on the website of Kenneth R. French [59], we
downloaded two historical datasets (i.e., ind49 and sbm100) of US stock returns.
We used monthly data from January 2010 to December 2019 to calculate thesample estimates (4.4), µ̂ and Σ̂ of asset returns. From the OR-Library [16, 41],
we downloaded two datasets (i.e., port2 and port5) of the sample estimates µ̂and Σ̂, which were multiplied by 100 to be consistent with the other datasets.
From Kaggle datasets [90], we downloaded a historical dataset (i.e., s&p500) of
US stock returns. To compute the sample estimates µ̂ and Σ̂, we used daily datafrom February 8th, 2013 to February 7th, 2018, where 32 companies including
missing values were omitted from the S&P500 index.
Table 4.1: Dataset description
Abbr.
ind49
port2
sbm100
port5
s&p500

N
49
89
100
225
468

Original dataset
49 Industry Portfolios [59]port2 (Portfolio optimization: Single period) [16, 41]
100 Portfolios Formed on Size and Book-to-Market [59]port5 (Portfolio optimization: Single period) [16, 41]468 companies in the S&P 500 index [90]As the utility function (4.6), we used a piecewise-linear approximation of thefollowing normalized exponential utility function:ũ(y) :=

µmax (1 − exp(−αy/µmax )),α(4.32)

where µmax is the maximum entry of the sample mean vector µ̂, and α > 0 is arisk-aversion parameter. As Figure 4.1 illustrates, we set α = 10 and employed50

three tangent lines (i.e., L = 3) at y ∈ {0, µmax /2, µmax } for piecewise-linear
approximation.
a(1) y + b(1)ua(2) y + b(2)a(3) y + b(3)ũ(y)µmax (1−exp(−10))10µmax (1−exp(−5))
10

O

µmax
2

µmax

y

Figure 4.1: Piecewise-linear approximation of the utility function (4.32) with
L=3
We compare the computational performance of the following methods forsolving the MISDO problem (4.10):
SCIP:

MISDO solver SCIP-SDP1 [62],CPA:

cutting-plane algorithm (Algorithm 4.1),

CPA+ : cutting-plane algorithm (Algorithm 4.1) using the matrix-completionbased problem reduction (i.e., Step 2 in Section 4.4.2).
These methods used Mosek2 9.2.40 to solve SDO problems. CPA and CPA+ were
implemented in Python 3.7 with Gurobi Optimizer3 8.1.11 to solve the surrogate
upper-level problem (4.17). We also used the lazy constraint callback to add
cutting planes during the branch-and-bound procedure. We set ε = 10−5 as thetolerance for optimality. The computation of each method was terminated if it
did not finish by itself within 3,600 s. In these cases, the results obtained within
3,600 s were taken as the final outcome.
We also applied the branch-and-cut algorithm in Chapter 2 for solving someinstances. The preliminary experiments showed that it was very slow and failed tocomplete the computations even for small-sized instances. Therefore, we did not
employ the branch-and-cut algorithm as a baseline in the numerical experiments.
The following column labels are used in Tables 4.2–4.4:
Obj: objective value of the obtained best feasible solution,
GAP(%): absolute difference between lower and upper bounds on the optimal
objective value divided by the upper bound,Time:#Cuts:

computation time in seconds,
number of cutting planes generated by the cutting-plane algorithms,#Nodes: number of nodes explored in the branch-and-bound procedure.
1

http://www.opt.tu-darmstadt.de/scipsdp/
https://www.mosek.com/
3https://www.gurobi.com/
251

Note that the best values of “Time” are indicated in bold for each problem
instance, and those of “Obj” are also indicated in bold for the port2 and s&p500
datasets. Also, “OM” in the Obj column indicates that the computation did notstart due to the out-of-memory condition.
Results for different values of the cardinality parameter k
Table 4.2 gives the numerical results of each method for√ the cardinality parameterk ∈ {5, 10, 15}. Here, we set the parameters γ = 10/ N for the `2 -regularization
term and (κ1 , κ2 ) = (1, 4) for the uncertainty set.
First, we focus on our cutting-plane algorithms (i.e., CPA and CPA+ ). BothCPA and CPA+ failed to finish solving problem instances for the port2 ands&p500 datasets; however, CPA+ solved other problem instances much fasterthan did CPA, especially for large N . In the case of the sbm100 dataset withk = 5, although CPA was terminated due to the time limit, CPA+ solved the
problem instance completely in 1.9 s. These results verify the effectiveness of ourmatrix-completion-based problem reduction of the dual lower-level problem.
Next, we compare our cutting-plane algorithms with the MISDO solver SCIP.
CPA+ was the fastest when it finished computations within the time limit. SCIP
returned an incorrect optimal objective value for the dataset ind49 with k ∈{10, 15} due to the numerical instability. For the port2 dataset, all three methods
failed to complete the computations within the time limit, but CPA+ foundsolutions of better quality than did SCIP and CPA for k ∈ {10, 15}.
For the port5 dataset, SCIP did not finish solving even a first continuous
relaxation problem within the time limit, and thus, failed to provide a feasible
solution for all k ∈ {5, 10, 15}. CPA did not start computations due to the outof-memory condition. In contrast, CPA+ succeeded in computing solutions withguaranteed optimality for all k ∈ {5, 10, 15}.
For the s&p500 dataset involving 468 investable assets, only CPA+ worked
normally to find feasible solutions. Optimality gaps attained by CPA+ for k ∈{5, 10, 15} are 21.0%, 19.1%, and 9.9%, respectively, and this gap is expected tobe smaller if longer time can be spent on the computation. These results support
the potential of CPA+ to give good-quality solutions to large problem instanceswith a limited memory capacity.
The computation time of CPA+ tended to be shorter for larger k. For the
port5 dataset, CPA+ finished the computations in 316.2 s, 84.6 s, and 16.16 s fork ∈ {5, 10, 15}, respectively. Here, both #Cuts and #Nodes were much smaller
for k = 15 than for k = 5, which is part of the reason why CPA+ got to be fasterwith larger k.
Results for different values of the uncertainty parameters (κ1 , κ2 )Table 4.3 gives the numerical results of each method for the pair of uncertainty
parameters (κ1 , κ2 ) ∈ {(0.5, 2), (1, 4), (2, 8)}. √Here, we set the parameters k = 10
for the cardinality constraint and γ = 10/ N for the `2 -reguralization term.
Table 4.3 shows that CPA+ consistently outperformed other methods regardless of values of the uncertainty parameters (κ1 , κ2 ). For the port2 dataset,all methods failed to complete the computations within the time limit for all(κ1 , κ2 ), but CPA+ found solutions of better quality than did other methods for(κ1 , κ2 ) ∈ {(1, 4), (2, 8)}. Moreover, the computation time of CPA+ was stableagainst the change in (κ1 , κ2 ). Indeed, for the port5 dataset, CPA+ finishedcomputations in 60.9 s, 84.6 s, and 62.6 s for (κ1 , κ2 ) ∈ {(0.5, 2), (1, 4), (2, 8)},52

respectively. In addition, for the s&p500 dataset, optimality gaps given by CPA+
were 17.9%, 19.1%, and 20.1% for (κ1 , κ2 ) ∈ {(0.5, 2), (1, 4), (2, 8)}, respectively.
These results suggest that the computational performance of the CPA+ is not
greatly affected by the uncertainty parameters (κ1 , κ2 ).
Results for different values of the `2 -regularization parameter γ
Table 4.4 gives the √numerical√ results√of each method for the `2 -regularization
parameter γ ∈ {1/ N , 10/ N , 100/ N }. Here, we set the parameters k =10 for the cardinality constraint and (κ1 , κ2 ) = (1, 4) for the uncertainty set.
Table 4.4 shows that CPA+ performed very well compared to other methods.
CPA+ was always faster than other methods for the ind49, sbm100, and port5
datasets. For the port2 dataset, all methods failed to solve problem instances tooptimality, but CPA√
+ found solutions of better quality than did other methods
except for γ = 100/ N . For the s&p500 dataset, only CPA+ provided feasible
solutions without causing an out-of-memory condition. The computation time of
CPA+ tended to be shorter for smaller γ. For the port5 √dataset,√CPA+ finished√the computations in 66.0 s, 84.6 s, and 417.1 s for γ ∈ {1/ N , 10/ N , 100/ N },respectively. A similar tendency was also shown by the results for the s&p500
dataset, wheregaps√ optimality√√ obtained by CPA+ were 4.2%, 19.1%, and 79.4%
for γ ∈ {1/ N , 10/ N , 100/ N }, respectively.
4.5.2

Out-of-sample investment performance

Finally, we investigate the investment performance of our cardinality-constrained
distributionally robust portfolio optimization model in a practical situation.
Experimental design
For the purpose of comparison, we consider the following cardinality-constrained
mean-variance portfolio optimization model:
minimizexx> Σ̂x(4.33a)subject to µ̂> x ≥ r̄,(4.33b)zn = 0 ⇒ xn = 0
x ∈ X,kz ∈ ZN,(∀n ∈ [N ]),(4.33c)(4.33d)where r̄ is a user-defined parameter of the required return level.
We compare the out-of-sample investment performance of portfolios determined by the following optimization models:MV:

cardinality-constrained mean-variance optimization model (4.33), which
was solved using Gurobi Optimizer 8.1.11;

DR:

cardinality-constrained distributionally robust optimization model (4.10),
which was solved by our cutting-plane algorithm CPA+ .

For the MV model, we set the required return level r̄ to the first quartile of entriesof the sample mean µ̂. For the DR model, we tuned the uncertainty parameters(κ1 , κ2 ) based on the 80% confidence region estimated by the tailored bootstrap√
method [29] with the bootstrap sample size of 100,000. We set γ = 10/ N forthe `2 -regularization term, and the utility function (4.6) in the same way as in
Section 4.5.1.
53

Cumulative Return

4.0

MV (k = 5)
MV (k = 30)

DR (k = 5)
DR (k = 30)

3.0
2.0
1.0
2001

2005

2010
Year

2015

2020

Figure 4.2: Evolution of cumulative total returns over the years 2000–2020
We evaluated the out-of-sample investment performance based on a rolling
horizon strategy. From Yahoo! Finance [167], we downloaded weekly data ofJapanese stock returns from January 2000 to December 2020, where top 30 companies were selected from the Nikkei 225 index according to market capitalization
as of December 2020. We solved portfolio optimization models by calculating
sample estimates µ̂ and Σ̂ in the first training period (156 weeks). We next
calculated weekly returns by applying the obtained portfolios to the subsequenttesting period (52 weeks). We repeated this process at intervals of 52 weeks until
the end of the entire data period.
We evaluate the out-of-sample investment performance based on the following
cumulative total return:Y(1 + Rm ),m∈[M ]

where Rm is the rate of portfolio return in the mth week during the testing period
M.
Results of the investment performance
Figure 4.2 shows the evolution of (out-of-sample) cumulative total returns produced by each optimization model with the cardinality parameter k ∈ {5, 30}.
The highest investment performance was achieved by the DR model with k = 5,
which suggests that the out-of-sample investment performance can be improvedby the interplay between the distributionally robust optimization and the cardinality constraint. In contrast, the investment performance of the MV model was
lower with k = 5 than with k = 30 probably because portfolio diversification isprevented by the cardinality constraint.
A notable result is the investment performance in 2020, when the stock market was significantly affected by the spread of COVID-19. In this year, while
the cumulative return was not increased by the MV model with k ∈ {5, 30},that was improved by the DR model with k ∈ {5, 30}. This implies that distributionally robust optimization models are helpful in coping with an unexpected
situation. These results demonstrate the practical effectiveness of our cardinalityconstrained distributionally robust portfolio optimization model in terms of theout-of-sample investment performance.

54

4.6

Conclusion

In this chapter, we consider moment-based distributionally robust portfolio optimization problems [47] with the cardinality constraint. Due to the discreteness, this problem is formulated mixed-integer semidefinite optimization problem, which is very hard to exactly solve when the number of investable assets islarge. We reformulated the problem as a bilevel optimization problem and devised a specialized cutting-plane algorithm for solving the upper-level problem.
In addition, we apply the matrix completion [60, 121] to the lower-level problem
to efficiently generate cutting planes.
The computational results indicate that our cutting-plane algorithm was very
effective, especially when the number of investable assets is large. Our algorithm
succeeded in obtaining an optimal solution within 3,600 seconds, to a problem
including 225 assets. For a large-sized problem instance with 468 assets, while theexisting method did not start its computation due to lack of memory, our algorithm attained a solution of good quality. In addition, the out-of-sample investment performances given by the cardinality-constrained distributionally robustmodel outperformed the cardinality-constrained mean-variance model.

55

√Table 4.2: Numerical results with γ = 10/ N and (κ1 , κ2 ) = (1, 4) for k ∈{5, 10, 15}
Data

N

k

MethodObj

Gap(%)

Time

#Cuts

#Nodes

ind49

49

5

SCIP
CPA
CPA+
SCIP
CPA
CPA+
SCIP
CPA
CPA+
SCIP
CPA
CPA+
SCIP
CPA
CPA+
SCIP
CPA
CPA+
SCIP
CPA
CPA+
SCIP
CPA
CPA+
SCIP
CPA
CPA+
SCIP
CPA
CPA+
SCIP
CPA
CPA+
SCIP
CPA
CPA+
SCIP
CPA
CPA+
SCIP
CPA
CPA+
SCIP
CPA
CPA+

3.108
3.108
3.108
3.036
3.034
3.034
3.037
3.033
3.033
1.951
2.181
2.071
1.761
1.924
1.727
1.749
2.009
1.635
4.494
3.940
3.940
3.955
3.935
3.935
3.952
3.935
3.935
∞
OM
2.812
∞
OM
2.687
∞
OM
2.677
OM
OM
1.051
OM
OM
0.910
OM
OM
0.830

0.0
2.2
0.0
0.0
0.0
0.0
0.0
0.0
0.0
15.1
374.1
68.9
8.7
188.6
51.6
8.0
98.1
17.5
12.1
0.3
0.0
0.4
0.0
0.0
0.3
0.0
0.0
100.0
—
0.0
100.0
—
0.0
100.0
—
0.0
—
—
21.0
—
—
19.1
—
—
9.9

627.9
>3600
26.5
303.2
400.6
9.3
408.6
107.8
5.2
>3600
>3600
>3600
>3600
>3600
>3600
>3600
>3600
>3600
>3600
>3600
1.9
>3600
1753.42.0>3600
1910.1
4.4
>3600
—
316.2
>3600
—
84.6
>3600
—
16.6
—
—
>3600
—
—
>3600
—
—
>3600

—
>195
182
—
22
22
—
6
6
—
>23
>20640
—
>20
>7937
—
>23
>4014
—
>12
13
—
5
5
—
5
5
—
—
770
—
—
191
—
—
18
—
—
>9794
—
—
>1996
—
—
>3435

61
>940
1047
17
20
20
43
0
0
>27
>71
>196390>20>39>147890>20>87>76327
>5
>1
8
>10
0
0
>9
0
0
>1
—
9402
>1
—
2077
>1
—
37
—
—
>477699——>118070——>270889

10
15

port2

86

5
10
15

sbm100

100

5
10
15

port5

225

5
10
15

s&p500

468

5
10
15

56

√
Table 4.3: Numerical results with k = 10 and γ = 10/ N for(κ1 , κ2 ) ∈ {(0.5, 2), (1, 4), (2, 8)}.
DataN

(κ1 , κ2 )MethodObj

Gap(%)

Time

#Cuts

#Nodesind49

49(0.5,2)

SCIP
CPA
CPA+
SCIP
CPA
CPA+
SCIP
CPA
CPA+
SCIP
CPA
CPA+
SCIP
CPA
CPA+
SCIP
CPA
CPA+
SCIP
CPA
CPA+
SCIP
CPA
CPA+
SCIP
CPA
CPA+
SCIP
CPA
CPA+
SCIP
CPA
CPA+
SCIP
CPA
CPA+
SCIP
CPA
CPA+
SCIP
CPA
CPA+
SCIP
CPA
CPA+

1.948
1.946
1.946
3.036
3.034
3.034
4.599
4.588
4.588
1.156
1.336
1.177
1.761
1.924
1.727
2.549
2.766
2.484
2.612
2.575
2.575
3.955
3.935
3.935
5.910
5.874
5.874
∞
OM
1.915
∞
OM
2.687
∞
OM
3.776
OM
OM
0.659
OM
OM
0.910
OM
OM
1.238

0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
5.9
135.1
34.7
8.7
188.6
51.6
7.8
255.2
61.7
1.3
0.0
0.0
0.4
0.0
0.0
0.3
0.0
0.0
100.0
—
0.0
100.0
—
0.0
100.0
—
0.0
—
—
17.9
—
—
19.1
—
—
20.1

476.6
409.2
9.3
303.2
400.6
9.3
289.6
323.3
7.5
>3600
>3600
>3600
>3600
>3600
>3600
>3600
>3600
>3600
>3600
1929.42.0>3600
1753.42.0>3600
2169.4
2.0
>3600
—
60.9
>3600
—
84.6
>3600
—
62.6
—
—
>3600
—
—
>3600
—
—
>3600

—
22
22
—
22
22
—
18
18
—
>22
>8448
—
>20
>7937
—
>20
>8415
—
5
5
—
5
5
—
6
5
—
—
144
—
—
191
—
—
142
—
—
>648
—
—
>1996
—
—
>5944

29
20
20
17
20
20
19
8
8
>22
>11
>158415>20>39>147890>21>59>139514
>8
0
0
>10
0
0
>10
0
0
>1
—
2177
>1
—
2077
>1
—
1486
—
—
>2558
—
—
>118070——>364115

(1,4)
(2,8)

port2

86

(0.5,2)
(1,4)
(2,8)

sbm100100(0.5,2)
(1,4)
(2,8)

port5

225

(0.5,2)
(1,4)
(2,8)

s&p500468(0.5,2)
(1,4)
(2,8)

57

Table 4.4:
results√ Numerical√√ with k = 10 and (κ1 , κ2 ) = (1, 4) forγ ∈ {1/ N , 10/ N , 100/ N }.
Data

N

ind49

49

γ
√
1/ N
√
10/ N
√
100/ N

port2

86

√
1/ N
√
10/ N
√
100/ N

sbm100

100

√
1/ N
√
10/ N
√
100/ N

port5

225

√
1/ N
√
10/ N
√
100/ N

s&p500

468

√
1/ N
√
10/ N
√
100/ N

MethodObj

Gap(%)

Time

#Cuts

#Nodes

SCIP
CPA
CPA+
SCIP
CPA
CPA+
SCIP
CPA
CPA+
SCIP
CPA
CPA+
SCIP
CPA
CPA+
SCIP
CPA
CPA+
SCIP
CPA
CPA+
SCIP
CPA
CPA+
SCIP
CPA
CPA+
SCIP
CPA
CPA+
SCIP
CPA
CPA+
SCIP
CPA
CPA+
SCIP
CPA
CPA+
SCIP
CPA
CPA+
SCIP
CPA
CPA+

3.415
3.415
3.415
3.036
3.034
3.034
2.984
2.984
2.984
2.154
2.251
2.154
1.761
1.924
1.727
1.699
2.143
1.711
5.029
4.594
4.594
3.955
3.935
3.935
3.806
3.801
3.801
∞
OM
3.380
∞
OM
2.687
∞
OM
2.611
OM
OM
1.861
OM
OM
0.910
OM
OM
0.789

0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
15.3
43.2
6.2
8.7
188.6
51.6
6.8
2805.9
459.1
9.2
0.0
0.0
0.4
0.0
0.0
0.0
0.0
0.0
100.0
—
0.0
100.0
—
0.0
100.0
—
0.0
—
—
4.2
—
—
19.1
—
—
79.4

433.2
339.8
8.2
303.2
400.6
9.3
323.9
398.4
10.0
>3600
>3600
>3600
>3600
>3600
>3600
>3600
>3600
>3600
>3600
2599.22.9>3600
1753.42.0
2094.42183.0
2.1
>3600
—
66.0
>3600
—
84.6
>3600
—
417.1
—
—
>3600
—
—
>3600
—
—
>3600

—
19
19
—
22
22
—
22
22
—
>23
>7591
—
>20
>7937
—
>21
>7487
—
7
7
—
5
5
—
6
5
—
—
144
—
—
191
—
—
968
—
—
>5452
—
—
>1996
—
—
>6580

30
18
18
17
20
20
20
37
42
>22
>159
>208109>20>39>147890>18>131>82148
>7
0
0
>10
0
0
3
0
0
>1
—
2007
>1
—
2077
>1
—
8718
—
—
>456468——>118070——>23941858

Chapter 5
Bilevel Cutting-plane Algorithm for
Cardinality-constrained Mean-CVaR
Portfolio Optimization
In this chapter, we examine cardinality-constrained conditional value-at-risk (CVaR)
minimization problems and extend the cutting-plane algorithm for the distributionallyrobust portfolio optimization problem discussed in Chapter 4. While the cardinalityconstrained CVaR minimization problem is formulated as an MILO problem,
which can be handled by state-of-the-art MIO solvers, the problem size of CVaR
minimization problems depends not only on the number of investable assets but
also on the number of asset return scenarios. Thus, the computational efficiency
decreases when the number of scenarios is large. To overcome this challenge,
we propose a specialized cutting-plane algorithm named the bilevel cutting-planealgorithm for exactly solving the cardinality-constrained mean-CVaR portfolio
optimization problem. Our proposal extends the cutting-plane algorithm discussed in Chapter 4 so that it incorporates another cutting-plane algorithm that
efficiently minimizes CVaR. The content of this chapter is included in Kobayashi
et al. [94].
We give an introduction to this study and summarize our contribution in
Section 5.1. In Section 5.2, we show the definition of CVaR and formulate thecardinality-constrained mean-CVaR portfolio optimization problem as an MIO
problem. In Section 5.3, we describe our bilevel cutting-plane algorithm for solving the problem. We report computational results in Section 5.4 and conclude in
Section 5.5.

5.1

Introduction

This chapter focuses on mean-risk portfolio optimization models using the conditional value-at-risk (CVaR) [137, 138] as a risk measure. CVaR is a downsiderisk measure for evaluating a potential heavy loss. It is known as a coherent risk
measure that has the desirable properties of monotonicity, translation invariance,
positive homogeneity, and subadditivity [12, 129]. Additionally, CVaR is monotonic with respect to second-order stochastic dominance [129], which means that
CVaR minimization is consistent with the preference of any rational risk-averse
decision-maker. In the standard formulation based on sample average approximation [137, 138, 145], many scenarios are required to approximate CVaR accurately [91, 150]. To resolve this computational difficulty, various efficient algorithms have been proposed, including the dual solution method [125], nonsmooth
optimization algorithms [17, 85, 104, 137], the factor model [100], cutting-planealgorithms [4, 80, 103, 149], the level method [52], smoothing methods [6, 157],
59

and successive regression approximations [3].
We consider solving mean-CVaR portfolio optimization problems with a cardinality constraint for limiting the number of invested assets. Because of the
cardinality constraint, this model is formulated as an MILO problem, which isa special case of MISDO problems. While state-of-the-art MIO solvers suchas Gurobi and CPLEX can be applied to the problem, the problem size of the
cardinality-constrained mean-CVaR model depends not only on the number of investable assets but also on the number of asset return scenarios. When we aim toapproximate CVaR accurately [91, 150], sufficiently many scenarios are required,
which decreases computational efficiency even if we use the state-of-the-art MIO
solvers. For this problem, some heuristic optimization algorithms have been developed, including continuous-relaxation-based heuristics [9], the `1 -norm-based
approximation of the cardinality constraint [43], and a Scholtes-type regularization method for complementarity constraints [37]. However, these algorithms
cannot guarantee global optimality of obtained solutions.
In this chapter, we extend the cutting-plane algorithm for distributionallyrobust portfolio optimization in Chapter 4 and propose a high-performance algorithm named the bilevel cutting-plane algorithm for exactly solving the cardinalityconstrained mean-CVaR portfolio optimization problem. Unlike the distributionallyrobust portfolio optimization problems discussed in Chapter 4, the size of thelower-level problem of the mean-CVaR model depends on not only the numberof investable assets but also that of scenarios, which degrades the computational
performance of the algorithm. To speed up the computations of the lower-levelproblem, we apply the cutting-plane algorithm [4, 80, 103, 149] that was developed for efficiently minimizing CVaR. As a result, these two types of cuttingplane algorithms are integrated into our bilevel cutting-plane algorithm to achieve
faster computations. To our knowledge, we are the first to develop an effective
algorithm for exactly solving cardinality-constrained mean-CVaR portfolio optimization problems.
We conducted numerical experiments using some benchmark datasets [16, 41,
59, 90] to evaluate the efficiency of our bilevel cutting-plane algorithm. Numericalresults demonstrate that our algorithm is faster than other MIO approaches,
especially for large problem instances. Remarkably, our algorithm attained an
optimal solution within 3600 s to a problem involving 225 assets and 100,000
scenarios.

5.2

Problem formulation

In this section, we formulate the cardinality-constrained mean-CVaR portfolio
optimization problem that we consider in this chapter.
5.2.1

Conditional value-at-risk

As we did in Chapter 4, we define x := (x1 , x2 , . . . , xN )> to be a portfolio satisfyingX
xn = 1, x ≥ 0,n∈[N ]

where xn is the investment weight of the nth asset. We consider a probability
space (Ω, F, P). Let ξ̃ : Ω → RN be an F-measurable function (N -dimensionalrandom vector) representing the rate of random return of each asset, and G(r) :=P({ω ∈ Ω | ξ̃(ω) ≤ r}) be the corresponding cumulative distribution function of60

the asset return vector r ∈ RN . In this chapter, we define the loss function as
the negative of the portfolio net return:L(x, r) := −r > x.
Let β ∈ (0, 1) be a probability level parameter, which is frequently set close toone. Then, β-CVaR can be regarded as the approximate conditional expectationof a random loss exceeding the β-value-at-risk (β-VaR). To calculate CVaR, we
use the following function:Z
1
Fβ (a, x) := a +[L(x, r) − a]+ dG(r),(5.1)
1 − β RN
where a is an auxiliary decision variable, and [v]+ is a positive part of v (i.e.,[v]+ := max{0, v}). Then, β-CVaR minimization is posed as follows [137, 138]:
subject to x ∈ X ,minimize Fβ (a, x)
a,xwhere X is a set of feasible portfolios. Note that an optimal solution of a corresponds to β-VaR, which is a β-quantile of the probability distribution of theportfolio loss L(x, ξ̃).
Because multiple integration in Eq. (5.1) is computationally expensive, thesample average approximation [137, 138, 145] is commonly used. Note also thateven when the number of scenarios is small, the cardinality constraint helps asolution based on the sample average approximation to nearly converge to theexact solution with high probability [2, 105].
Accordingly, we consider a finite sample space Ω := {ωs | s ∈ [S]}, where
S is the number of scenarios. Here, r (s) := ξ̃(ωs ) is the sth scenario of asset
returns for s ∈ [S], and p := (p1 , p2 , . . . , pS )> is a vector composed of occurrence
probabilities ps := P({ωs }). We then haveFβ (a, x) ≈ a +i
1 X hps −(r (s) )> x − a .
1−β
+
s∈[S]

5.2.2

Portfolio optimization modelThroughout this chapter, we consider the following set of feasible portfolios:
X := {x ∈ RN | Cx ≤ d, 1> x = 1, x ≥ 0},
where C ∈ RM ×N and d ∈ RM are given. It is supposed that the linear constraint
Cx ≤ d contains the expected return constraint:µ̂> x ≥ µ̄,(5.2)

where µ̂ := (µ̂1 , µ̂2 , . . . , µ̂N )> is a vector composed of expected returns of assets,
and µ̄ is a parameter of the required return level.
As we did in Chapter 4, we also impose a cardinality constraint on the portfoliox:kxk0 ≤ k,(5.3)
where k ≥ 1 be a user-defined parameter for limiting the cardinality and k · k0
is the `0 -norm (i.e., the number of nonzero entries) [27, 34, 128]. Let z :=(z1 , z2 , . . . , zN )> be a vector of binary decision variables for selecting assets; that61

is, zn = 1 if the nth asset is selected, and zn = 0 otherwise. We also introduce
the feasible set corresponding to the cardinality constraint (5.3):




X
k
zn = k .
ZN:= z ∈ {0, 1}Nn∈[N ]

In line with previous studies [24, 48, 70, 71], we incorporate the `2 -regularization
term into the objective from a perspective of robust optimization. The cardinalityconstrained mean-CVaR portfolio optimization model is formulated as follows:
1 >x x+a+v
a,v,x,z2γi
1 X hps −(r (s) )> x − a ,subject to v ≥1−β+minimize(5.4a)(5.4b)s∈[S]

z n = 0 ⇒ xn = 0
x ∈ X,(∀n ∈ [N ]),kz ∈ ZN
,(5.4c)(5.4d)

where v is an auxiliary decision variable, and γ > 0 is a regularization parameter.
Throughout this chapter, we assume that Problem (5.4) is feasible.

5.3

Cutting-plane algorithms

In this section, we present our bilevel cutting-plane algorithm for solving the
cardinality-constrained mean-CVaR portfolio optimization problem (5.4).
5.3.1

Bilevel optimization reformulation

Here, we extend the method of bilevel optimization reformulation in Chapter 4 tothe mean-CVaR portfolio optimization problem (5.4). Let us define Z := Diag(z)as a diagonal matrix whose diagonal entries are given by z. We first eliminate
the logical implication (5.4c) by replacing x with Zx in Problem (5.4). We next
reformulate Problem (5.4) as a bilevel optimization problem. Specifically, theupper-level problem is posed as the following integer optimization problem:
minimize f (z)zksubject to z ∈ ZN
,

(5.5)

and its objective function is defined by the following lower-level problem:
1 >x x+a+v
a,v,x
2γ
i
1 X h
subject to v ≥ps −(r (s) )> Zx − a ,1−β+f (z) = minimize(5.6a)(5.6b)s∈[S]

Zx ∈ X .

(5.6c)

Note that (Zx)> Zx has been replaced with x> x in the objective (5.6a); thisis because x = Zx holds after minimization (5.6a) as in Bertsimas and CoryWright [24].
We then derive cutting planes for approximating f (z) in the upper-level problem (5.5) by exploiting the dual formulation of the lower-level problem (5.6). For
this purpose, Constraint (5.6b), which is nonlinear and nondifferentiable, needs
to be transformed into tractable constraints. A commonly used method is the62lifting representation [52, 137], which converts Problem (5.6) into the following
optimization problem with linear constraints:
1 >x x+a+v
2γ
1 X
subject to v ≥p s qs ,1−βf (z) = minimize(5.7a)

a,q,v,x(5.7b)s∈[S]

qs ≥ −(r(s) >) Zx − a (∀s ∈ [S]),(5.7c)(∀s ∈ [S]),(5.7d)

qs ≥ 0Zx ∈ X ,(5.7e)

where q := (q1 , q2 , . . . , qS )> is a vector of auxiliary decision variables. The following theorem gives the dual formulation of Problem (5.7).
Theorem 5.1. Suppose that Problem (5.7) is feasible with z ∈ {0, 1}N . Then,the strong duality holds, and the dual formulation of Problem (5.7) is represented
as follows:f (z) = maximize
α,ζ,λ,ω

subject toγ >z (ω ◦ ω) − d> ζ + λ2Xω≥
αs r (s) − C > ζ + λ1,−(5.8a)(5.8b)s∈[S]X

αs = 1,(5.8c)s∈[S]

ps(∀s ∈ [S]),1−βα ≥ 0, ζ ≥ 0,αs ≤(5.8d)(5.8e)

where α := (α1 , α2 , . . . , αS )> , ζ ∈ RM , λ ∈ R, and ω ∈ RN are dual decisionvariables, and ω ◦ ω denotes the Hadamard product of the vector ω.
Proof. See Appendix A.3.
Remark 5.1. Note that the(5.8) is always feasible with α = p,
P dual problem(s)ζ = 0, λ = 0, and ω = s∈[S] ps r . From the strong duality, the dual problem (5.8) is unbounded (i.e., f (z) = ∞) if and only if the corresponding primal
problem (5.7) is infeasible. Our objective is to minimize f (z), so we can assume
without loss of generality that the primal problem (5.7) is feasible.
Remark 5.2. When zn = 0, we can eliminate ωn from the dual problem (5.8).
Suppose that after the elimination, (α? , ζ ? , λ? ) is an optimal solution of (α, ζ, λ)to Problem (5.8). We then recover ωn as


X
?
?
ωn? = αs? rn(s) − c>(n ∈ [N ] with zn = 0),
nζ +λ
s∈[S]

+

where cn is the nth column vector of C. This solution satisfies Constraint (5.8b)
and thus optimal because its coefficient is zero in the objective (5.8a).
Theorem 5.1 allows us to redefine f (z) as the optimal objective value of
Problem (5.8) for z ∈ [0, 1]N . Then, as we showed Lemmas 4.3 and 4.4 in
Chapter 4, we obtain two properties of f (z).
Lemma 5.2 (Convexity). The function f (z) is convex in z ∈ [0, 1]N .
63

Lemma 5.3 (Subgradient). Suppose that Problem (5.7) is feasible with z ∈{0, 1}N , and that ω ? (z) is an optimal solution of ω to Problem (5.8). Then, a
subgradient of the function f (z) is given byγg(z) := − ω ? (z) ◦ ω ? (z) ∈ ∂f (z).
(5.9)
2k , Problem (5.8) with
From Lemmas 5.2 and 5.3, we see that for each ẑ ∈ ZN
z = ẑ yields a linear underestimator of f (z) for z ∈ [0, 1]N as follows:f (z) ≥ f (ẑ) + g(ẑ)> (z − ẑ).
5.3.2

(5.10)

Upper-level cutting-plane algorithm

Here, we extend the cutting-plane algorithm in Chapter 4 to the upper-levelproblem (5.5) for mean-CVaR portfolio optimization. We refer to this algorithm
as the upper-level cutting-plane algorithm, which finds a sequence of solutions torelaxed versions of Problem (5.5).
We first define the initial feasible region ask
F1 := {(z, θ) ∈ ZN
× R | θ ≥ θLB },(5.11)

where θ is an auxiliary decision variable that serves as a lower bound on f (z),and θLB is a lower bound of the optimal objective value of Problem (5.5). At thetth iteration (t ≥ 1), our algorithm solves the following optimization problem:
minimize θz,θ

subject to (z, θ) ∈ Ft ,(5.12)

where Ft is a relaxed feasible region at the tth iteration such that Ft+1 ⊆ Ft .
According to Eq. (5.11), the objective value of Problem (5.12) is bounded below. Therefore, unless Ft = ∅, Problem (5.12) has an optimal solution, which isdenoted by (zt , θt ).
After obtaining (zt , θt ), we solve the dual lower-level problem (5.8) with z =zt . If Problem (5.8) is unbounded, Problem (5.7) is infeasible because of thestrong duality. In this case, we update the feasible region to cut off the solution
zt as follows:k
Ft+1 ← Ft ∩ {(z, θ) ∈ ZN
× R | zt> (1 − z) + (1 − zt )> z ≥ 1}.

(5.13)

If Problem (5.8) is bounded, we obtain the function value f (zt ) and its subgradient g(zt ) as in Lemma 5.3.
If f (zt ) − θt ≤ ε with sufficiently small ε ≥ 0, then zt is an ε-optimal solution
to Problem (5.5), which means thatf ? ≤ f (zt ) ≤ f ? + ε,

(5.14)

where f ? is the optimal objective value of Problem (5.5). In this case, we terminate the algorithm with the ε-optimal solution zt . Otherwise, we add theconstraint (5.10) to the feasible region:k
Ft+1 ← Ft ∩ {(z, θ) ∈ ZN
× R | θ ≥ f (zt ) + g(zt )> (z − zt )}.

(5.15)

Note that this update cuts off the solution (zt , θt ) because θt < f (zt ).
After updating the feasible region, we set t ← t + 1 and solve Problem (5.12)
again. This procedure is repeated until an ε-optimal solution ẑ is found. Aftertermination of the algorithm, we can compute the corresponding portfolio bysolving Problem (5.7) with z = ẑ. Our upper-level cutting-plane algorithm issummarized by Algorithm 5.1.
64

Algorithm 5.1 Upper-level cutting-plane algorithm for solving Problem (5.5)
Step 0 (Initialization) Let ε ≥ 0 be a tolerance for optimality. Define thefeasible region F1 as in Eq. (5.11). Set t ← 1 and UB0 ← ∞.
Step 1 (Relaxed Problem) Solve Problem (5.12). Let (zt , θt ) be an optimal
solution, and set LBt ← θt .
Step 2 (Cut Generation) Solve Problem (5.8) with z = zt to calculate f (zt )and ω ? (zt ).
(a) If Problem (5.8) is unbounded, set UBt ← UBt−1 and update thefeasible region as in Eq. (5.13).
(b) Otherwise, perform the following procedures:
i. Calculate g(zt ) as in Eq. (5.9).
ii. Update the feasible region as in Eq. (5.15).
iii. If f (zt ) < UBt−1 , set UBt ← f (zt ) and (ẑ, θ̂) ← (zt , θt );otherwise, set UBt ← UBt−1 .
Step 3 (Termination Condition) If UBt − LBt ≤ ε, terminate the algorithmwith the ε-optimal solution ẑ.
Step 4 Set t ← t + 1 and return to Step 1.

5.3.3

Lower-level cutting-plane algorithm

Note that Step 2 of Algorithm 5.1 solves Problem (5.8), which contains O(S)decision variables and O(S) constraints; so, solving it directly is computationally expensive especially when the number of scenarios is very large. To avoid
this difficulty, we propose solving the primal lower-level problem (5.6) by instead
using the other cutting-plane algorithm [4, 80, 103, 149]. Hereinafter, we refer to this algorithm as the lower-level cutting-plane algorithm to distinguish itfrom the upper-level cutting-plane algorithm (Algorithm 5.1). We will explain in
Section 5.3.4 how to derive dual solutions for cut generation from Problem (5.6).
We first use the cutting-plane representation [52, 103] to reformulate Problem (5.6) as follows:
1 >x x+a+v
2γ
1 X
ps (−(r (s) )> Zx − a)subject to v ≥1−βf (z) = minimize(5.16a)a,v,x(∀J ⊆ [S]),(5.16b)s∈JZx ∈ X ,(5.16c)where J denotes a subset of the scenario set [S]. We pick out the constraint (5.16b)with J = [S] to define the initial relaxed feasible region:


1 X

v≥
ps (−(r (s) )> Zx − a),

1−βG1 (z) := (a, v, x) ∈ R × R × RN
.
s∈[S]




v ≥ 0, Zx ∈ X(5.17)
65

It is clear from Eq. (5.4b) that v ≥ 0 is a valid constraint. At the tth iteration (t ≥1), our algorithm solves the following optimization problem:
minimize
a,v,x

subject to1 >x x+a+v
2γ(a, v, x) ∈ Gt (z),(5.18a)(5.18b)where Gt (z) is a relaxed feasible region at the tth iteration such that Gt (z) ⊆G1 (z). The objective value of Problem (5.18) is bounded below by the initial
feasible region (5.17). Therefore, unless Gt (z) = ∅, Problem (5.18) has an optimal
solution, which is denoted by (at , vt , xt ). Note that Problem (5.18) contains onlyN + 2 decision variables, which are independent of the number of scenarios.
After obtaining (at , vt , xt ), we define a subset of scenarios
Jt ← {s ∈ [S] | −(r (s) )> Zxt − at > 0}

(5.19)1 X
ps (−(r (s) )> Zxt − at ).
1−β(5.20)and calculate
vt0 ←

s∈Jt

Because the revised solution (at , vt0 , xt ) satisfies all the constraints of Problem (5.16),
this solution gives an upper bound of f (z).
If vt0 − vt ≤ δ with sufficiently small δ ≥ 0, we terminate the algorithm withthe δ-optimal solution (at , vt0 , xt ) to Problem (5.16). Otherwise, we update thefeasible region as follows:
(
)
1 X
N
(s) >
Gt+1 (z) ← Gt (z)∩ (a, v, x) ∈ R × R × Rv≥
ps (−(r ) Zx − a) .
1−β
s∈Jt

(5.21)
Because vt < vt0 , this update cuts off the solution (at , vt , xt ).
Next, we set t ← t + 1 and solve Problem (5.18) again. The lower-level
cutting-plane algorithm is summarized by Algorithm 5.2. As shown in previous
studies [4, 80, 103, 149], Algorithm 5.2 terminates in a finite number of iterations
and outputs a δ-optimal solution (â, v̂, x̂) to the lower-level problem (5.6). Thisalgorithm is empirically much faster than solving Problem (5.7) based on thelifting representation when the number of scenarios is very large.
Suppose that Algorithm 5.2 terminates at the T th iteration and outputs a
δ-optimal solution (â, v̂, x̂), where (aT , vT , xT ) is an optimal solution to Problem (5.18) with t = T . The resultant family of scenario subsets Kδ (z) ={[S], J1 , J2 , . . . , JT } is also provided by Algorithm 5.2. We then obtain lower
and upper bounds on f (z) as follows:f (z) − δ ≤ fδ (z) ≤ f (z) ≤ fˆδ (z) ≤ f (z) + δ,

(5.22)

wherefδ (z) :=

5.3.4

1 >
x xT + aT + vT ,
2γ T

1 >
fˆδ (z) :=x̂ x̂ + â + v̂.
2γ

Efficient cut generation for the upper-level problem

Now, we turn to the computation of dual solutions from primal solutions provided
by the lower-level cutting-plane algorithm (Algorithm 5.2). Such dual solutions

66

Algorithm 5.2 Lower-level cutting-plane algorithm for solving Problem (5.16)
Step 0 (Initialization) Let δ ≥ 0 be a tolerance for optimality. Define thefeasible region G1 (z) as in Eq. (5.17). Set t ← 1.
Step 1 (Relaxed Problem) Solve Problem (5.18). Let (at , vt , xt ) be an optimal solution.
Step 2 (Upper Bound) Calculate Jt and vt0 as in Eqs. (5.19) and (5.20).
Step 3 (Termination Condition) If vt0 −vt ≤ δ, terminate the algorithm withthe δ-optimal solution (â, v̂, x̂) := (at , vt0 , xt ) and the resultant family of
scenario subsets Kδ (z) := {[S], J1 , J2 , . . . , Jt }.
Step 4 (Cut Generation) Update the feasible region as in Eq. (5.21).
Step 5 Set t ← t + 1 and return to Step 1.

are required for cut generation in the upper-level cutting-plane algorithm (Algorithm 5.1).
Suppose that Algorithm 5.2 results in the following relaxed problem with
K = Kδ (z):1 >x x+a+v
2γ
1 X
subject to v ≥ps (−(r (s) )> Zx − a)1−βfK (z) := minimize(5.23a)a,v,x(∀J ∈ K),

(5.23b)s∈Jv ≥ 0,

Zx ∈ X .

(5.23c)The following theorem derives the dual formulation of Problem (5.23).
Theorem 5.4. Suppose that Problem (5.23) is feasible with z ∈ {0, 1}N . Then,the strong duality holds, and the dual formulation of Problem (5.23) can be writtenasfK (z) = maximize
α,ζ,λ,ω

subject toγ >z (ω ◦ ω) − d> ζ + λ
2
X
X
ω≥
αJ
ps r (s) − C > ζ + λ1,
−

J ∈K

X

(5.24a)(5.24b)s∈JαJ ≤ 1,(5.24c)X(5.24d)

J ∈K

X

αJ

J ∈K

α ≥ 0,

ps = 1 − β,s∈Jζ ≥ 0,(5.24e)where α := (αJ )J ∈K , ζ ∈ RM , λ ∈ R, and ω ∈ RN are dual decision variables.
Proof. See Appendix A.4.
Note that |Kδ (z)| is usually very small even when there are many scenarios.
Therefore, Problem (5.24) can be solved efficiently regardless of the number ofscenarios.
From Theorem 5.4, we can redefine fK (z) to be the optimal objective valueof Problem (5.24) for z ∈ [0, 1]N as we did for f (z). Let ωδ? (z) be the optimal67

solution of ω to Problem (5.24) with K = Kδ (z). We then defineγgδ (z) := − ωδ? (z) ◦ ωδ? (z).
2(5.25)

Note that when δ > 0, fδ (z) and gδ (z) are not exactly the same as f (z) andg(z), respectively. However, the following theorem verifies that fδ (z) and gδ (z)still provide a linear underestimator of f (z).
Theorem 5.5. Suppose that ẑ ∈ {0, 1}N and Problem (5.23) is feasible with(z, K) = (ẑ, Kδ (ẑ)). Then, it holds that for all z ∈ [0, 1]N ,f (z) ≥ fδ (ẑ) + gδ (ẑ)> (z − ẑ).

(5.26)

Proof. Because K = Kδ (ẑ), it follows from Lemma 5.3 that gδ (ẑ) is a subgradient of fK (z) at z = ẑ. It then holds that for all z ∈ [0, 1]N ,f (z) ≥ fK (z) ≥ fK (ẑ) + gδ (ẑ)> (z − ẑ).

(5.27)

This proof is completed because fδ (ẑ) = fK (ẑ) with K = Kδ (ẑ).
5.3.5

Bilevel cutting-plane algorithm

We are now ready to describe our bilevel cutting-plane algorithm for solving
the upper-level problem (5.5). We use the lower-level cutting-plane algorithm(Algorithm 5.2) to accelerate the cut generation at Step 2 of Algorithm 5.1.
Suppose that (zt , θt ) is an optimal solution to the relaxed problem (5.12) atthe tth iteration (t ≥ 1). We then solve the primal lower-level problem (5.16) withz = zt by means of the lower-level cutting-plane algorithm (Algorithm 5.2), which
provides fδ (zt ), fˆδ (zt ), and Kδ (zt ). We next solve the reduced version (5.24) of
the dual lower-level problem with (z, K) = (zt , Kδ (zt )), thereby yielding the
optimal solution ωδ? (zt ). After that, we calculate gδ (zt ) as in Eq. (5.25) andupdate the feasible region as follows:k
Ft+1 ← Ft ∩ {(z, θ) ∈ ZN
× R | θ ≥ fδ (zt ) + gδ (zt )> (z − zt )}.

(5.28)

Our bilevel cutting-plane algorithm is summarized by Algorithm 5.3. Notethat at Step 2 (Algorithm 5.3), the portfolio x̂ is generated by Algorithm 5.2 ateach iteration.
Remark 5.3. Several algorithms (e.g., primal-dual interior-point methods [166]and alternating direction methods of multipliers [35]) solve primal and dual problems simultaneously. By applying such algorithms in Step 1 of Algorithm 5.2, we
can skip solving the dual lower-level problem (5.24) at Step 2 of Algorithm 5.3.
Remark 5.4. Note that Step 1 of Algorithm 5.3 solves the MIO problem (5.12) ateach iteration; this amounts to a “multi-tree” implementation, where a branchand-bound algorithm is repeatedly executed from scratch. To improve computational efficiency, we can use the function of lazy constraint callback, which
is offered by modern optimization software (e.g., CPLEX or Gurobi). This function enables a “single-tree” implementation [132]; that is, cutting planes (5.28)
are dynamically generated during the process of the branch-and-bound algorithm.
To prove the convergence properties of Algorithm 5.3, we show the following
lemma.

68

Algorithm 5.3 Bilevel cutting-plane algorithm for solving Problem (5.5)
Step 0 (Initialization) Let ε ≥ 0 and δ ≥ 0 be tolerances for optimality.
Define the feasible region F1 as in Eq. (5.11). Set t ← 1 and UB0 ← ∞.
Step 1 (Relaxed Problem) Solve Problem (5.12). Let (zt , θt ) be an optimal
solution, and set LBt ← θt .
Step 2 (Cut Generation) Solve Problem (5.16) with z = zt by means of
Algorithm 5.2 to calculate fδ (zt ), fˆδ (zt ), and Kδ (zt ).
(a) If Problem (5.16) is infeasible, set UBt ← UBt−1 and update thefeasible region as in Eq. (5.13).
(b) If Problem (5.16) is feasible, perform the following procedures:
i. Solve Problem (5.24) with (z, K) = (zt , Kδ (zt )) to calculateωδ? (zt ).
ii. Calculate gδ (zt ) as in Eq. (5.25).
iii. Update the feasible region as in Eq. (5.28).
iv. If fˆδ (zt ) < UBt−1 , set UBt ← fˆδ (zt ) and (ẑ, θ̂) ← (zt , θt );otherwise, set UBt ← UBt−1 .
Step 3 (Termination Condition) If zu = zt for some u < t or UBt −LBt ≤ ε,then terminate the algorithm with the max{δ, ε}-optimal solution ẑ.
Step 4 Set t ← t + 1 and return to Step 1.
Lemma 5.6. Suppose that {(zt , θt ) | t = 1, 2, . . . , T } is a sequence of solutions
generated by Algorithm 5.3. If there exists u < T such that zu = zT , then zT isa δ-optimal solution to Problem (5.5), meaning thatf ? ≤ f (zT ) ≤ f ? + δ.
Proof. Note that (zT , θT ) is contained in the feasible region (5.28) at t = u.
Because zu = zT , it follows that
θT ≥ fδ (zu ) + gδ (zu )> (zT − zu ) = fδ (zT ).
From Eqs. (5.14) and (5.22), we havef ? ≤ f (zT ),f (zT ) − δ ≤ fδ (zT ).

Because θT ≤ f ? , it follows thatf ? ≤ f (zT ) ≤ fδ (zT ) + δ ≤ θT + δ ≤ f ? + δ,
which completes the proof.
From Lemma 5.6, we can establish the convergence properties of Algorithm
5.3.
Theorem 5.7. Algorithm 5.3 terminates in a finite number of iterations andoutputs a max{δ, ε}-optimal solution.
69

k , Algorithm
Proof. Because there are at most a finite number of solutions z ∈ ZN
5.3 terminates in a finite number of iterations with zu = zt for some u < t. In
this case, a δ-optimal solution is found according to Lemma 5.6. If UBt −LBt ≤ ε
is fulfilled first, Algorithm 5.3 terminates with an ε-optimal solution.

5.4

Numerical experiments

In this section, we report numerical results to evaluate the efficiency of our algorithms for solving cardinality-constrained mean-CVaR portfolio optimization
problems.
5.4.1

Problem instances

Table 5.1 lists the datasets used in our experiments, where N is the number ofassets. From the data library on the website of Kenneth R. French [59], we downloaded two historical datasets (i.e., ind49 and sbm100) of US stock returns. We
used monthly data from January 2010 to December 2019 to compute the mean
vector µ̂ and covariance matrix Σ̂ of asset returns. From Kaggle datasets [90], we
downloaded a dataset (i.e., sp200) of historical stock returns, where top 200 companies were selected from the S&P 500 index according to market capitalization
as of December 2020. To compute µ̂ and Σ̂, we used daily data from February
8th, 2013 to February 7th, 2018. From the OR-Library [16, 41], we downloaded
three datasets (i.e., port1, port2, and port5) of µ̂ and Σ̂, which were multipliedby 100 to be consistent with the other datasets.
Table 5.1: Dataset description
Abbr.
ind49
sbm100
sp200
port1
port2
port5

N
49
100
200
31
89
225

Original dataset
49 Industry Portfolios [59]100 Portfolios Formed on Size and Book-to-Market [59]Top 200 companies in the S&P 500 index [90]port1 (Portfolio optimization: Single period) [16, 41]
port2 (Portfolio Foptimization: Single period) [16, 41]
port5 (Portfolio optimization: Single period) [16, 41]

For each dataset, we randomly generated S scenarios of asset returns {r (s) |s ∈ [S]} from N(µ̂, Σ̂), a normal distribution with the parameters (µ̂, Σ̂). We set
the occurrence probability as ps = 1/S for all s ∈ [S]. The required return level
in Eq. (5.2) was set as µ̄ = 0.3µmin + 0.7µmax , where µmin and µmax were theaverage returns of the bottom- and top-k assets, respectively; this setting ensures
feasibility of Problem (5.4). The probability level of CVaR was set as β = 0.9.
5.4.2

Methods for comparison

A standard MIO formulation of Problem (5.4) uses the big-M method, which
replaces the logical implication (5.4c) with0 ≤ x n ≤ zn(∀n ∈ [N ]),which is valid because of the constraint x ∈ X in Problem (5.4).

70(5.29)

Another state-of-the-art MIO formulation uses the perspective reformulation [74, 75]:
1 X
yn + a + v
a,v,x,y,z2γn∈[N ]i
1 X hsubject to v ≥ps −(r (s) )> x − a ,1−β+minimize(5.30a)(5.30b)s∈[S]x2n ≤ yn zn ,x ∈ X,yn ≥ 0(∀n ∈ [N ]),kz ∈ ZN,(5.30c)(5.30d)where y := (y1 , y2 , . . . , yN )> is a vector of auxiliary decision variables. In recent
years, the perspective reformulation has been actively studied [56, 57, 79, 173].
We compare the computational performances of the following methods:
BigM solves Problem (5.4) after replacing Eq. (5.4c) with Eq. (5.29);
Persp solves Problem (5.30) as a mixed-integer second-order cone optimization
problem;
Lift replaces Eqs. (5.4b) and (5.30b) with the lifting representation (5.7b)–(5.7d);
Cut applies the cutting-plane algorithm [4, 80, 103, 149] to Problems (5.4)
and (5.30);
CPsolves Problem (5.5) by means of the upper-level cutting-plane algorithm(Algorithm 5.1);

BCP solves Problem (5.5) by means of the bilevel cutting-plane algorithm (Algorithm 5.3);
BCPc solves Problem (5.5) by means of the bilevel cutting-plane algorithm(Algorithm 5.3), where the callback function is employed in implementation.
All experiments were performed on a Windows 10 PC with an Intel Core i74790 CPU (3.6.0GHz) and 16 GB of memory. All methods were implemented in
Python 3.7 with Gurobi Optimizer 8.1.11 . We employed the primal-dual interiorpoint method offered by Gurobi at Step 1 of Algorithm 5.2. We also tested
the performance of an alternating direction method of multipliers, but it madeAlgorithm 5.2 slower because of numerical inaccuracies.
We set ε = δ = 10−5 for tolerances for optimality. Except for BCP, we usedlazy constraint callback for adding cutting planes during the branch-andbound procedure.
The computation of each method was terminated if it did not finish by itself
within 3600 s. In these cases, the results obtained within 3600 s were taken as
the final outcome.
5.4.3

Evaluation metrics

The row labels used in the tables of experimental results are defined as follows:
Solvednumber of problem instances solved to optimality within 3600 s;

1

The source code is available athttps://github.com/KenKoba2119/cardinality-constrained_cvar_optimization.

71

Time

computation time in seconds;

Obj objective value of the obtained best feasible solution;
Gap(%) absolute difference between lower and upper bounds on the optimal
objective value divided by the upper bound;#Nodes number of nodes explored in the branch-and-bound algorithm;#Cuts number of upper-level cutting planes generated.
Note that the best values of Time are indicated in bold for each instance, andthose of Obj and Gap(%) are also indicated in bold for only the sp200, port2,
and port5 datasets.
To aggregate results of “Time,” “#Nodes,” and “#Cuts” in Tables 5.2 and 5.5
and Figure 5.1, we use the shifted geometric mean, which is used in Chapter 2.
Recall the shifted geometric mean of values x1 , x2 , . . . , xN is defined as
1/N


Y


(xn + s)− s,n∈[N ]

where s is the shift parameter. Following Achterberg [1] and Mittlemann [77],we set s = 10 for “Time.” We respectively set s to 100 and 10 for “#Nodes” and“#Cuts” in proportion to their scale. We use the arithmetic mean for “Gap(%).”
5.4.4

Results for various numbers of scenarios

We evaluate the computational performance of each method with the numberof scenarios S ∈ {103 , 104 , 105 }. Here,
we set the parameters k = 10 for the√
cardinality constraint and γ = 10/ N for the `2 -regularization term. Table 5.2
summarizes numerical results for the six datasets (Table 5.1), and Tables 5.3 and5.4 give the numerical results for each dataset.
First, we focus on the results of our cutting-plane algorithms (i.e., CP, BCP,
and BCPc). BCP and BCPc were faster than CP for S ≥ 104 (Table 5.2). In
the case of the ind49 dataset (Table 5.3), for example, although CP was thefastest among the three methods for S = 103 , BCP and BCPc were much fasterthan CP for S ≥ 104 . These results suggest the effectiveness of the lower-level
cutting-plane algorithm, which is used by BCP and BCPc to solve the lower-levelproblem efficiently regardless of the number of scenarios.
The differences in computation time between BCP and BCPs were relatively
small (Table 5.2). For the port2 dataset (Table 5.4), BCP and BCPc failed tocomplete the computations within 3600 s even when S = 103 ; however, BCPc
found solutions of better quality than did BCP. BCP provides at most one feasible solution at every iteration t, whereas BCPc explores feasible solutions more
frequently through the callback procedure. For this reason, BCPc is capable ofyielding solutions of good quality even if the computation is terminated due tothe time limit.
Next, we compare our bilevel cutting-plane algorithms (i.e., BCP and BCPc)
with the MIO formulations (i.e., BigM and Persp). Persp+Lift was the fastest
among the six methods for S ≤ 104 , whereas BCP and BCPc were faster thanthe other methods for S = 105 (Table 5.2). Indeed, for the sbm100 dataset with
S = 105 (Table 5.3), BCP was much faster than the four methods related to
BigM and Persp. Moreover, for the port5 dataset with S = 105 (Table 5.4), only
72

BCP and BCPc finished solving the problem within 3600 s, which is a remarkable
result.
When S = 105 , Persp+Lift finished solving only one problem instance, andits computation time was much longer than those of BCP and BCPc (Table 5.2).
A main reason for this is that Persp+Lift solves at each node a second-ordercone optimization problem whose problem size depends on S. In the case of
the port5 dataset with S = 103 (Table 5.4), Persp+Cut returned an incorrect
optimal objective value because of numerical instability.
According to the convergence properties of the sample average approximation,a huge number of scenarios are required for calculating CVaR accurately [150].
Even with only 15 investable assets, it was necessary to have at least 5000 scenarios to ensure the stability of the CVaR optimization model [91]. These facts
support the practicality of our bilevel cutting-plane algorithm, which has thepotential to deal with many scenarios as shown in Tables 5.2–5.4.
5.4.5

Sensitivity to hyperparameter values

We examine the sensitivity of the computational performance to the hyperparameters k for the cardinality constraint and γ for the `2 -regularization term.
Here, we focus on the four datasets, namely, ind49, sbm100, port1, and port5.
Sensitivity to the cardinality parameter k
Table 5.5 summarizes numerical results for the √four datasets with the cardinality
parameter k ∈ {5, 10, 15}. Here, we set γ = 10/ N for the `2 -regularization termand S = 105 as the number of scenarios. The numerical results for each dataset
are given in Appendix B (Table B.1).
Table 5.5 shows that when k = 5, although only Persp+Cut finished solving
all the problem instances, the relative gaps of BCP and BCPc were very small(i.e., 0.36 and 1.68), which means that these methods found solutions of goodquality. When k ∈ {10, 15}, both BCP and BCPc were faster than the othermethods and solved all the problem instances to optimality within 3600 s.
BCP and BCPc tended to be faster with larger k. Here, #Cuts was also muchsmaller for k = 15 than for k = 5, which is part of the reason why BCP and BCPc
performed better with larger k.
Sensitivity to the regularization parameter γ
Figure 5.1 summarizes results of Time, GAP(%), #Nodes,√and #Cuts for the four
datasets with the regularization parameter γ ∈ {100.2i / N | i = 0, 1, . . . , 20}.
Here, we set k = 10 for the cardinality constraint and S = 105 as the numberof scenarios. The numerical results for each dataset are shown in Appendix B(Figures B.1–B.4).
We can see from Figure 5.1 that the performances of BCP and BCPc were
not greatly affected by γ, and their computation times were consistently shorter
than the other methods for every γ. Additionally, BCP and BCPc tended to befaster with smaller γ.
We next analyze the impact of the `2 -regularization term on the objective
value. Suppose that (a? (γ), v ? (γ), x? (γ), z ? (γ)) is an optimal solution computedby the bilelevel cutting-plane algorithm to Problem (5.4). Figure 5.2 shows73

BigM+Lift
BigM+Cut

Persp+Lift
Persp+Cut

CP
BCP

BCPc

50

BigM+Lift
BigM+Cut

Persp+Lift
Persp+Cutp

10 2 / N

CP
BCP

BCPc

40

Relative gap (%)

Computation time (s)

104

103

102

1001 p10 / Np

10 1 / Np

p

10 2 / N

γ

20
10
0p

p

10 3 / N30

10 4 / N

10 0 / N(a) Time
BigM+Lift
BigM+Cut

105p

10 1 / Nγp

10 3 / Np

10 4 / N(b) Gap(%)

Persp+Lift
Persp+Cut

CP

BCPc

CP

103

BCP

BCPc

p

10 3 / N

Number of cuts

Number of nodes

104
103
102

102

101
1000 p10 / Np

10 1 / Np

p

10 2 / N

10 3 / Nγ1001 pp

10 4 / N

10 / N(c) #Nodesp

10 1 / N

10 2 / Nγpp

10 4 / N(d) #Cuts

Figure 5.1: Summary of numerical results for the three√ datasets (port1, ind49,
50.2i
and sbm100) with (S, k) = (10 , 10) for γ ∈ {10 / N | i = 0, 1, . . . , 20}“Reg + CVaR” and “CVaR” defined as
1 ? > ?
x (γ) x (γ) + a? (γ) + v ? (γ)
|
{z
}
2γ
|
{z
}
CVaR

(5.31)

Reg√with the regularization parameter γ ∈ {100.2i / N | i = 0, 1, . . . , 20}. These
values were averaged over the three datasets (port1, ind49, sbm100), where
port5 was omitted because most of the methods failed to find an optimal solution
due to the time limit. The numerical results for each dataset are shown in B
(Figure B.5).
We find from Figure √5.2 that the objective value (i.e., Reg + CVaR) was close2
to CVaR when γ ≥ 10
√ / N . In addition, CVaR almost converged to its minimalvalue when γ ≥ 10/ N . This result suggests√that we can obtain a sufficiently
good solution minimizing CVaR with γ ≥ 10/ N .

5.5

Conclusion

In this chapter, we dealt with mean-CVaR portfolio optimization problems basedon the sample average approximation with the cardinality constraint. It is very
hard to exactly solve the problem when it involves a large number of investableassets. In addition, accurate approximation of CVaR based on the sample average approximation requires sufficiently many scenarios, which decreases computational efficiency. To overcome these challenges, we proposed a specialized
74

4.6

Reg+CVaR

CVaR

p

10 3 / N

Objective value

4.5
4.4
4.3
4.2
4.1
4.0
3.9 p

10 0 / Np

10 1 / N

10 2 / Nγpp

10 4 / N

Figure 5.2: Summary of the objective value (5.31) for the three√ datasets (port1,
ind49, and sbm100) with (S, k) = (105 , 10) for γ ∈ {100.2i / N | i = 0, 1, . . . , 20}
cutting-plane algorithm named the bilevel cutting-plane algorithm for solving
the cardinality-constrained mean-CVaR portfolio optimization problem. We extended the cutting-plane algorithm discussed in Chapter 4 so that it incorporates
another cutting-plane algorithm that efficiently minimizes CVaR. We also provedthat our algorithms give a solution with guaranteed global optimality in a finite
number of iterations.
The computational results indicate that our cutting-plane algorithms werevery effective especially when the number of scenarios was large. Remarkably,
our bilevel cutting-plane algorithm attained an optimal solution within 3600 s to a
problem involving 225 assets and 100,000 scenarios. Furthermore, our algorithms
performed well for most of the hyperparameter values.

75Table 5.2: Summaryof numerical results for the six datasets (Table 5.1) with√(k, γ) = (10, 10/ N ) for S ∈ {103 , 104 , 105 }
BigM

S

Persp

Lift

Cut

Lift

Cut

CP

BCP

BCPc

103

SolvedTime
Gap(%)#Nodes
#Cuts

6
6.8
0.00
334.1
—

5
107.6
3.19
16351.5
—

6
5.5
0.00
106.1
—

4
109.6
3.58
12296.4
—

5
79.2
0.73
2234.6
322.3

5
116.7
1.80
—
147.4

5
71.0
1.15
1305.6166.1

104

SolvedTime
Gap(%)#Nodes
#Cuts

6
106.6
0.00
546.3
—

4
207.7
2.61
20097.8
—

6
89.1
0.00
228.3
—

4
185.4
2.02
13973.2
—

4
338.1
2.49
1507.4
216.8

4
154.6
1.97
—
96.8

4
172.0
1.34
1378.5165.2

105SolvedTime
Gap(%)#Nodes#Cuts

31610.4
21.09
53.4
—

3
650.0
5.30
9265.2—

13245.9
2.95
9.2
—

3
535.7
5.09
6763.1—31265.0
4.77
493.9
91.4

4
370.8
1.88
—
57.5

4
497.3
3.05
476.9
76.3

76

√
Table 5.3: Numerical results for the datasets [59, 90] with (k, γ) = (10, 10/ N )for S ∈ {103 , 104 , 105 }
Data
ind49

N
49

103

104

105

sbm100

100

103

104

105

sp200

200

BigM

S

103

104

105

Time
Obj
Gap(%)#Nodes
#Cuts
Time
Obj
Gap(%)#Nodes
#Cuts
Time
Obj
Gap(%)#Nodes
#Cuts
Time
Obj
Gap(%)#Nodes
#Cuts
Time
Obj
Gap(%)#Nodes
#Cuts
Time
Obj
Gap(%)#Nodes
#Cuts
Time
Obj
Gap(%)#Nodes
#Cuts
Time
Obj
Gap(%)#Nodes
#Cuts
Time
Obj
Gap(%)#Nodes
#Cuts

Persp

Lift

Cut

Lift

Cut

1.8
3.095
0.00
16
—
18.0
3.343
0.00
1
—
709.8
3.379
0.00
31
—
3.3
4.337
0.00
1
—
39.1
4.397
0.00
1
—
1107.1
4.364
0.00
1
—
30.9
1.190
0.00
7790
—
1010.7
1.267
0.00
22932
—
>3600
1.499
17.48
>68
—

4.8
3.095
0.00
3792
—
21.2
3.343
0.00
3357
—
125.0
3.379
0.00
4101
—
1.1
4.337
0.00
396
—
8.1
4.397
0.00
1820
—
65.1
4.364
0.00
412
—
1379.01.190
0.00323767
—
>3600
1.275
2.51
>968714
—
>3600
1.322
9.20
>59214
—

2.3
3.095
0.00
23
—
26.0
3.343
0.00
1
—
>3600
3.380
0.00
>1
—
3.7
4.337
0.00
1
—
45.0
4.397
0.00
1
—
>3600
4.364
0.00
>1
—
12.5
1.190
0.00
645
—
250.8
1.267
0.00
1052
—
>3600
1.326
8.12
>1
—

11.8
3.095
0.00
4506
—
26.2
3.343
0.00
3815
—
220.0
3.379
0.00
3489
—
4.2
4.337
0.00
4626
—
9.3
4.397
0.00
3274
—
84.6
4.364
0.00
4676
—
>3600
1.274
6.88
>87982
—
>3600
1.321
5.20
>100939
—
>3600
1.302
8.47
>24933
—

77

CP

BCP

BCPc

6.1
3.095
0.00
211
76
47.1
3.343
0.00
89
47
662.2
3.379
0.00
250
68
4.9
4.337
0.00
14
21
51.1
4.397
0.00
58
20
482.7
4.364
0.00
10
17
277.1
1.190
0.00
20663
899
>3600
1.285
5.10
>21183
>1122
>3600
1.331
18.21
>1873
>145

7.3
3.095
0.00
—
40
5.2
3.343
0.00
—
1
43.1
3.379
0.00
—
1
0.6
4.337
0.00
—
1
3.3
4.397
0.00
—
1
27.6
4.364
0.00
—
1
3441.1
1.190
0.00
—
757
>3600
1.273
1.20
>—
>771
>3600
1.282
1.97
—
>533

8.6
3.095
0.00
33
45
8.3
3.343
0.00
0
7
103.6
3.379
0.00
0
12
0.8
4.337
0.00
0
7
4.0
4.397
0.00
0
7
36.6
4.364
0.00
0
7
226.6
1.190
0.00
20670
953
>3600
1.267
1.53
>77070
>2941
>3600
1.282
4.73
>6095
>406

√
Table 5.4: Numerical results for the datasets [16, 41] with (k, γ) = (10, 10/ N )for S ∈ {103 , 104 , 105 }
Data
port1

N
31

103

104

105

port2

89

103

104

105

port5

225

103

104

105

†

BigM

S
Time
Obj
Gap(%)#Nodes
#Cuts
Time
Obj
Gap(%)#Nodes
#Cuts
Time
Obj
Gap(%)#Nodes
#Cuts
Time
Obj
Gap(%)#Nodes
#Cuts
Time
Obj
Gap(%)#Nodes
#Cuts
Time
Obj
Gap(%)#Nodes
#Cuts
Time
Obj
Gap(%)#Nodes
#Cuts
Time
Obj
Gap(%)#Nodes
#Cuts
Time
Obj
Gap(%)#Nodes
#Cuts

Persp

Lift

Cut

Lift

Cut

1.1
4.404
0.00
1
—
11.7
4.374
0.00
1
—
468.6
4.264
0.00
1
—
7.4
1.773
0.00
1432
—
311.2
1.938
0.00
5715
—
>3600
2.030
9.03
>481
—
8.3
2.933
0.00
368
—
246.1
3.147
0.00
428
—
>3600
∞
100.00
>0
—

2.3
4.404
0.00
2301
—
11.4
4.374
0.00
850
—
163.3
4.264
0.00
2433
—
>3600
2.119
19.13
>287,645
—
>3600
2.054
13.17
>261,628
—
>3600
2.140
18.61
>99,640
—
250.0
2.933
0.00
45,800
—
668.1
3.147
0.00
42,363
—
>3600
3.246
3.98
>20,833
—

1.5
4.404
0.00
1
—
24.0
4.374
0.00
1
—
1933.2
4.264
0.00
1
—
7.5
1.773
0.00
577
—
429.0
1.938
0.00
2429
—
>3600
2.039
9.17
>1
—
7.9
2.933
0.00
21
—
113.0
3.147
0.00
317
—
>3600
3.143
0.44
>61
—

4.3
4.404
0.00
2494
—
17.7
4.374
0.00
3005
—
15.8
4.264
0.00
253
—
>3600
2.002
14.61
>152,730
—
>3600
1.981
6.90
>204,649
—
>3600
2.074
17.76
>40,749—40.6
2.650†
0.00
4674
—
211.4
3.147
0.00
9056
—
>3600
3.276
4.29
>16,790
—

CP

BCP

BCPc

3.0
4.404
0.00
58
36
18.0
4.374
0.00
46
17
265.8
4.264
0.00
39
25
>3600
1.833
8.52
>216,174>15,007
>3600
2.009
9.83
>22,111
>3654
>3600
1.988
8.35
>1756
>406
144.8
2.933
0.00
7982
775
1387.1
3.147
0.00
8269
701
>3600
3.173
2.05
>2140
>219

0.8
4.404
0.00
—
1
4.8
4.374
0.00
—
1
39.5
4.264
0.00
—
1
>3600
2.022
10.80
—
>7190
>3600
2.181
10.61
—
>2695
>3600
2.136
9.30
—
>350
157.8
2.933
0.00
—
446
500.1
3.147
0.00
—
517
2351.6
3.138
0.00
—
352

1.7
4.404
0.00
0
7
8.2
4.374
0.00
0
8
64.6
4.264
0.00
0
7
>3600
1.773
6.90
>36,644
>3899
>3600
1.943
6.53
>18,077
>2548
>3600
1.950
13.55
>927
>285
130.3
2.933
0.00
7499
489
587.5
3.147
0.00
7347
727
3300.8
3.138
0.00
5697
519

Gurobi returned an incorrect optimal objective value because of numerical instability.

78

Table 5.5: Summary of numerical results √for the four datasets (ind49, sbm100,
port1, and port5) with (S, γ) = (105 , 10/ N ) for k ∈ {5, 10, 15}
BigM

k

Persp

Lift

Cut

Lift

Cut

CP

BCP

BCPc

5

SolvedTime
Gap(%)#Nodes#Cuts

21863.3
50.00
36.9
—

2
795.3
4.36
8430.2—

13422.7
7.12
35.2
—

4
766.0
0.00
7914.4—2
1902.8
4.26
968.1
167.7

2
789.4
0.36
—
166.8

2
905.2
1.68
1817.6224.8

10SolvedTime
Gap(%)#Nodes#Cuts

21789.2
50.00
7.3
—

2
592.9
1.99
5440.8—03600.0
0.22
27.5
—

2
719.7
2.14
8262.3—2
1431.3
1.02
563.0
92.5

4
315.0
0.00
—
53.1

4
480.8
0.00
661.4
91.1

15

SolvedTime
Gap(%)#Nodes#Cuts

21662.2
50.00
0.5
—

4
305.9
0.00
2539.6—13558.8
0.00
1.0
—

4
268.0
0.00
3769.4—2
1268.8
0.98
494.1
76.5

4
53.5
0.00
—
1.0

4
88.1
0.00
0.0
8.9

79

Chapter 6
Conclusion and Prospects
Now let us conclude this thesis. Specifically, we will summarize the contributions
described herein and offer some possible future directions for research.

6.1

SummaryThis thesis focused on solving mixed-integer semidefinite optimization problems.
The mixed-integer semidefinite optimization problem involves minimizing or maximizing a linear objective function subject to constraints whereby a given matrix formed from the decision variables is positive semidefinite, and some of thevariables are integer-valued. Since this problem includes nonlinearity and discreteness, various practical optimization problems can be formulated as a mixedinteger semidefinite optimization problem.
In this thesis, we devised efficient cutting-plane algorithms for solving mixedinteger semidefinite optimization problems. First, we proposed a general-purpose
cutting-plane algorithm for solving the standard form of mixed-integer semidefinite optimization problems. This method allows us to handle the general mixedinteger semidefinite optimization problem with existing state-of-the-art mixedinteger optimization solvers such as Gurobi and CPLEX. Second, we devised
specialized cutting-plane algorithms for some important applications of mixedinteger semidefinite optimization. In statistics and financial engineering, wesometimes encounter large-sized mixed-integer semidefinite problem instances,and such problems are difficult to solve even with general-purpose algorithms.
We focused on a variable selection and portfolio selection problem formulated
as mixed-integer semidefinite problems and proposed efficient cutting-plane algorithms that exploit the structures of these problems.
In Chapter 2, we described a general-purpose framework for solving the standard form of mixed-integer semidefinite optimization problems. First, we formulated a cutting-plane algorithm to solve mixed-integer semidefinite optimization problems and proved its convergence properties. Then, we developed abranch-and-cut algorithm where cutting planes are added dynamically to therelaxed mixed-integer linear optimization problem during a branch-and-boundprocedure. The experimental results confirmed that our branch-and-cut algorithm could solve three different mixed-integer semidefinite optimization problems: random instances, computing restricted isometry constants, and robust
truss topology design. Our algorithm is based on the branch-and-bound procedure for solving mixed-integer semidefinite optimization problems, so it can use
a warm-starting strategy to solve a series of continuous relaxation problems efficiently. Moreover, our algorithm does not solve any semidefinite optimization
problems, which makes the mixed-integer semidefinite optimization computation80

very stable and practical.
In Chapter 3, we addressed the problem of selecting the best subset of explanatory variables subject to an upper bound on the condition number for eliminating multicollinearity from linear regression models. The first contribution ofthis study is a novel computational framework for eliminating multicollinearity
based on the mixed-integer semidefinite optimization formulation. This framework reformulates the subset selection problem as a single mixed-integer semidefinite optimization problem. The second contribution is the establishment of ahigh-performance cutting-plane algorithm. In this algorithm, we generate strong
cutting planes with a heuristic search and reduce the number of relaxation problems to be solved. While numerical experiments show that our mixed-integer
semidefinite optimization formulation can only be applied to small-sized instances
at present, this chapter provides a new statistical application of mixed-integer
semidefinite optimization formulation. Moreover, we found that our cuttingplane algorithm frequently provided a better subset of variables than did thecommon local search algorithms.
In Chapter 4, we studied moment-based distributionally robust portfolio optimization problems with a cardinality constraint. Because of the discreteness ofthe cardinality constraint, this problem is formulated as a mixed-integer semidefinite optimization problem, which is hard to solve exactly when the number ofinvestable assets is large. We reformulated the problem as a bilevel optimization
problem and devised a cutting-plane algorithm for solving the upper-level problem. In addition, we applied the technique of positive semidefinite matrix completion to the lower-level problem in order to efficiently generate cutting planes.
The computational results indicate that our cutting-plane algorithm is more effective than the existing general-purpose mixed-integer semidefinite optimization
solver, especially when the number of investable assets is large. In addition, theout-of-sample investment performances given by the cardinality-constrained distributionally robust model were better than those of the cardinality-constrained
mean-variance model.
Chapter 5 dealt with cardinality-constrained mean-CVaR portfolio optimization problems, and we extended the cutting-plane algorithm for the distributionallyrobust portfolio optimization problem discussed in Chapter 4. This model isformulated as a mixed-integer linear optimization problem, which is a special
case of mixed-integer semidefinite optimization problems. Since its problem size
depends not only on the number of investable assets but also on the numberof asset return scenarios, the computational efficiency decreases when the number of scenarios is large. To overcome this challenge, we propose a specialized
cutting-plane algorithm named the bilevel cutting-plane algorithm for solving
the cardinality-constrained mean-CVaR portfolio optimization problem. We extended the cutting-plane algorithm discussed in Chapter 4 so that it incorporates
another cutting-plane algorithm that efficiently minimizes CVaR. The computational results indicate that our cutting-plane algorithms are very effective, especially when the number of scenarios is large.

6.2

Future directions

The results in this thesis suggest that the cutting-plane algorithms could be a
powerful tool for solving mixed-integer semidefinite optimization problems. However, there is still much room for research on mixed-integer semidefinite optimization.

81

While our general-purpose algorithms delivered better computational performances than the existing methods for several problem instances in Chapter 2, ouralgorithms are not always faster for general problem instances. Thus, we haveto continue improving the efficiency of our algorithms to make mixed-integer
semidefinite optimization more practical. In Chapters 3 to 5, we devised specialized cutting-plane algorithms that exploit the individual problems’ structure
to improve the computational efficiency. In particular, we proposed to generate strong cutting planes with a heuristic search for subset selection problems in
Chapter 3. Also, we used the technique of positive semidefinite matrix completion
to calculate cutting planes efficiently in Chapter 4. To develop more efficient andpractical general-purpose algorithms, these techniques should be incorporated
into our general-purpose algorithms for mixed-integer semidefinite optimization.
Thus, we would be interested in characterizing whether these specialized techniques can be applied to a given mixed-integer semidefinite optimization problem.
Of course, it also makes sense to extend the specialized cutting-plane algorithms not only to mixed-integer semidefinite optimization problems but also toother problems appearing in the corresponding fields. In the case of subset selection, it is interesting to extend the cutting-plane algorithm so that it can deal
with other variable selection problems for statistics and machine learning fields.
For portfolio optimization, various practical and distributionally robust risk measures have been proposed. We expect that our cutting-plane algorithm can beextended to handle such risk measures.
While mixed-integer semidefinite optimization has the potential to be applied
to a wide range of practical decision-making problems, vigorous research has juststarted. We will continue refining both general-purpose and specialized algorithms for mixed-integer semidefinite optimization, and expand the horizons of
its applications. As our modern society is becoming more and more complex, we
hope that this work will be a new avenue for solving some of its fundamental
problems.

82

References[1] T. Achterberg. SCIP: solving constraint integer programs. Mathematical
Programming Computation, 1(1):1–41, 2009.
[2] A. Agarwal, S. N. Negahban, and M. J. Wainwright. Stochastic optimization and sparse statistical recovery: Optimal algorithms for high dimensions. In Proceedings of the 25th International Conference on Neural Information Processing Systems, pages 1538–1546, 2012.
[3] K. C. Ágoston. CVaR minimization by the SRA algorithm. Central European Journal of Operations Research, 20(4):623–632, 2012.
[4] S. Ahmed. Convexity and decomposition of mean-risk stochastic programs.
Mathematical Programming, 106(3):433–446, 2006.
[5] F. A. Al-Khayyal and J. E. Falk. Jointly constrained biconvex programming. Mathematics of Operations Research, 8(2):273–286, 1983.
[6] S. Alexander, T. F. Coleman, and Y. Li. Minimizing CVaR and VaR fora portfolio of derivatives. Journal of Banking and Finance, 30(2):583–605,
2006.
[7] F. Alizadeh. Interior point methods in semidefinite programming with applications to combinatorial optimization. SIAM Journal on Optimization,
5(1):13–51, 1995.
[8] D. Aloise and P. Hansen. A branch-and-cut SDP-based algorithm for minimum sum-of-squares clustering. Pesquisa Operacional, 29(3):503–516, 2009.
[9] E. Angelelli, R. Mansini, and M. G. Speranza. A comparison of MAD
and CVaR models with real features. Journal of Banking and Finance,
32(7):1188–1197, 2008.
[10] M. F. Anjos, B. Ghaddar, L. Hupp, F. Liers, and A. Wiegele. Solving kway graph partitioning problems to optimality: The impact of semidefinite
relaxations and the bundle method. In M. Jünger and G. Reinelt, editors,
Facets of Combinatorial Optimization, pages 355–386. Springer, 2013.
[11] M. Armbruster, M. Fügenschuh, C. Helmberg, and A. Martin. LP and
SDP branch-and-cut algorithms for the minimum graph bisection problem:
a computational comparison. Mathematical Programming Computation,
4(3):275–306, 2012.
[12] P. Artzner, F. Delbaen, J. M. Eber, and H. David. Coherent measures ofrisk. Mathematical Finance, 9(3):203–228, 1999.
[13] R. Baraniuk. Compressive sensing [lecture notes]. IEEE Signal Processing
Magazine, 24(4):118–121, 2007.
83

[14] E. M. L. Beale. Two transportation problems. In Proceedings of the Third
International Conference on Operational Research, pages 780–788, 1963.
[15] E. M. L. Beale and J. Tomlin. Special facilities in a general mathematicalprogramming system for non-convex problems using ordered sets of variables. In Proceedings of the Fifth International Conference on Operational
Research, pages 447–454, 1970.
[16] J. E. Beasley. OR-library: Distributing test problems by electronic mail.
The Journal of the Operational Research Society, 41(11):1069–1072, 1990.
[17] G. Beliakov and A. Bagirov. Non-smooth optimization methods for computation of the conditional value-at-risk and portfolio optimization. Optimization, 55(5-6):459–479, 2006.
[18] A. Ben-Tal, L. El Ghaoui, and A. Nemirovski. Robust Optimization. Princeton Series in Applied Mathematics. Princeton University Press, 2009.
[19] A. Ben-Tal and A. Nemirovski. Robust truss topology design via semidefinite programming. SIAM Journal on Optimization, 7(4):991–1016, 1997.
[20] A. Ben-Tal and A. Nemirovski. Robust optimization-methodology and applications. Mathematical Programming, 92(3):453–480, 2002.
[21] S. J. Benson, Y. Ye, and X. Zhang. Solving large-scale sparse semidefinite
programs for combinatorial optimization. SIAM Journal on Optimization,
10(2):443–461, 2000.
[22] L. Berk and D. Bertsimas. Certifiably optimal sparse principal component
analysis. Mathematical Programming Computation, 11(3):381–420, 2019.
[23] D. Bertsekas, A. Nedić, and A. Ozdaglar. Convex Analysis and Optimization. Athena Scientific, 2003.
[24] D. Bertsimas and R. Cory-Wright. A scalable algorithm for sparse portfolio
selection. INFORMS Journal on Computing, 2022. https://doi.org/10.
1287/ijoc.2021.1127.
[25] D. Bertsimas, R. Cory-Wright, and J. Pauphilet. Solving large-scale sparse
PCA to certifiable (near) optimality. arXiv preprint arXiv:2005.05195,
2020.
[26] D. Bertsimas, R. Cory-Wright, and J. Pauphilet. A unified approach tomixed-integer optimization problems with logical constraints. SIAM Journal on Optimization, 31(3):2340–2367, 2021.
[27] D. Bertsimas, C. Darnell, and R. Soucy. Portfolio construction throughmixed-integer programming at Grantham, Mayo, Van Otterloo and Company. Interfaces, 29(1):49–66, 1999.
[28] D. Bertsimas, I. Dunning, and M. Lubin. Reformulation versus cuttingplanes for robust optimization. Computational Management Science,13(2):195–217, 2016.
[29] D. Bertsimas, V. Gupta, and N. Kallus. Data-driven robust optimization.
Mathematical Programming, 167(2):235–292, 2017.

84

[30] D. Bertsimas and A. King. OR forum—an algorithmic approach to linear
regression. Operations Research, 64(1):2–16, 2016.
[31] D. Bertsimas and A. King. Logistic regression: From art to science. Statistical Science, 32(3):367–384, 2017.
[32] D. Bertsimas, A. King, and R. Mazumder. Best subset selection via amodern optimization lens. The Annals of Statistics, 44(2):813–852, 2016.
[33] D. Bertsimas, J. Lamperski, and J. Pauphilet. Certifiably optimal sparse
inverse covariance estimation. Mathematical Programming, 184(1-2):491–
530, 2019.
[34] D. Bienstock. Computational study of a family of mixed-integer quadratic
programming problems. Mathematical Programming, 74(2):121–140, 1996.
[35] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein. Distributed
optimization and statistical learning via the alternating direction methodof multipliers. Foundations and Trends in Machine Learning, 3(1):1–122,
2011.
[36] S. Boyd and L. Vandenberghe. Convex Optimization. Cambridge University
Press, 2004.
[37] M. Branda, M. Bucher, M. Červinka, and A. Schwartz. Convergence ofa scholtes-type regularization method for cardinality-constrained optimization problems with an application in sparse robust portfolio optimization.
Computational Optimization and Applications, 70(2):503–530, 2018.
[38] G. Braun, S. Fiorini, S. Pokutta, and D. Steurer. Approximation limits oflinear programs (beyond hierarchies). Mathematics of Operations Research,
40(3):756–772, 2015.
[39] S. B. Çay, I. Pólik, and T. Terlaky. Warm-start of interior point methods forsecond order cone optimization via rounding over optimal jordan frames.
Technical Report 17T-00, Lehigh University, 2017.
[40] A. Cerveira, A. Agra, F. Bastos, and J. Gromicho. A new branch andbound method for a discrete truss topology design problem. Computational
Optimization and Applications, 54(1):163–187, 2013.
[41] T. J. Chang, N. Meade, E. B. John, and Y. M. Sharaiha. Heuristics for
cardinality constrained portfolio optimisation. Computers and Operations
Research, 27(13):1271–1302, 2000.
[42] S. Chatterjee and A. S. Hadi. Regression Analysis by Example, Fifth Edition. John Wiley & Sons, 2012.
[43] R. Cheng and J. Gao. On cardinality constrained mean-cvar portfolio optimization. In Proceedings of the 27th Chinese Control and Decision Conference, pages 1074–1079, 2015.
[44] I. G. Chong and C. H. Jun. Performance of some variable selection methods
when multicollinearity is present. Chemometrics and Intelligent Laboratory
Systems, 78(1-2):103–112, 2005.

85

[45] C. Coey, M. Lubin, and J. P. Vielma. Outer approximation with conic
certificates for mixed-integer convex problems. Mathematical Programming
Computation, 12(2):249–293, 2020.
[46] J. Czyzyk, M. Mesnier, and J. More. The NEOS server. IEEE Computational Science and Engineering, 5(3):68–75, 1998.
[47] E. Delage and Y. Ye. Distributionally robust optimization under moment
uncertainty with application to data-driven problems. Operations Research,
58(3):595–612, 2010.
[48] V. DeMiguel, L. Garlappi, F. J. Nogales, and R. Uppal. A generalized
approach to portfolio optimization: Improving performance by constraining
portfolio norms. Management Science, 55(5):798–812, 2009.
[49] A. V. Dorugade and D. N. Kashid. Variable selection in linear regression
based on ridge estimator. Journal of Statistical Computation and Simulation, 80(11):1211–1224, 2010.
[50] M. Efroymson. Multiple regression analysis. In A. Ralston and H. Wilf,
editors, Mathematical Methods for Digital Computers, pages 191–203. John
Wiley & Sons, 1960.
[51] L. El Ghaoui, M. Oks, and F. Oustry. Worst-case value-at-risk and robust portfolio optimization: A conic programming approach. Operations
Research, 51(4):543–556, 2003.
[52] C. I. Fábián. Handling CVaR objectives and constraints in two-stage
stochastic models. European Journal of Operational Research, 191(3):888–
911, 2008.
[53] F. J. Fabozzi, D. Huang, and G. Zhou. Robust portfolios: contributions
from operations research and finance. Annals of Operations Research,
176(1):191–220, 2010.
[54] D. E. Farrar and R. R. Glauber. Multicollinearity in regression analysis:
The problem revisited. The Review of Economics and Statistics, 49(1):92,
1967.
[55] S. Foucart and M. J. Lai. Sparsest solutions of underdetermined linear
systems via `q -minimization for 0 < q ≤ 1. Applied and Computational
Harmonic Analysis, 26(3):395–407, 2009.
[56] A. Frangioni, F. Furini, and C. Gentile. Approximated perspective relaxations: a project and lift approach. Computational Optimization andApplications, 63(3):705–735, 2016.
[57] A. Frangioni and C. Gentile. SDP diagonalizations and perspective cuts fora class of nonseparable MIQP. Operations Research Letters, 35(2):181–185,
2007.
[58] l. E. Frank and J. H. Friedman. A statistical view of some chemometrics
regression tools. Technometrics, 35(2):109–135, 1993.
[59] K. R. French. Kenneth R. French - data library. https://mba.
tuck.dartmouth.edu/pages/faculty/ken.french/data_library.html,
Accessed 17 July 2020.
86

[60] M. Fukuda, M. Kojima, K. Murota, and K. Nakata. Exploiting sparsity
in semidefinite programming via matrix completion I: General framework.
SIAM Journal on Optimization, 11(3):647–674, 2001.
[61] T. Gally and M. E. Pfetsch. Computing restricted isometry constants via
mixed-integer semidefinite programming. Oprimization Online, 2016.
[62] T. Gally, M. E. Pfetsch, and S. Ulbrich. A framework for solving
mixed-integer semidefinite programs. Optimization Methods and Software,
33(3):594–632, 2018.
[63] A. Gleixner, L. Eifler, T. Gally, G. Gamrath, P. Gemander, R. L. Gottwald,
G. Hendel, C. Hojny, T. Koch, M. Miltenberger, B. Müller, M. Pfetsch,
C. Puchert, D. Rehfeldt, F. Schlösser, F. Serrano, Y. Shinano, J. M. Viernickel, S. Vigerske, D. Weninger, J. T. Witt, and J. Witzig. The SCIP
optimization suite 5.0. Technical Report 17–61, Zuse Institute Berlin, 2017.
[64] J.-L. Goffin and J.-P. Vial. Convex nondifferentiable optimization: A survey
focused on the analytic center cutting plane method. Optimization Methods
and Software, 17(5):805–867, 2002.
[65] J. Goh and M. Sim. Distributionally robust optimization and its tractableapproximations. Operations Research, 58(4-part-1):902–917, 2010.
[66] D. Goldfarb and G. Iyengar. Robust portfolio selection problems. Mathematics of Operations Research, 28(1):1–38, 2003.
[67] R. E. Gomory. Outline of an algorithm for integer solutions to linear
programs. Bulletin of the American Mathematical Society, 64(5):275–279,
1958.
[68] J. Gondzio. Warm start of the primal-dual method applied in the cuttingplane scheme. Mathematical Programming, 83(1-3):125–143, 1998.
[69] J. Gotoh, M. J. Kim, and A. E. B. Lim. Calibration of distributionallyrobust empirical optimization models. Operations Research, 69(5):1630–
1650, 2021.
[70] J. Gotoh, K. Shinozaki, and A. Takeda. Robust portfolio techniques for
mitigating the fragility of CVaR minimization and generalization to coherent risk measures. Quantitative Finance, 13(10):1621–1635, 2013.
[71] J. Gotoh and A. Takeda. On the role of norm constraints in portfolio
selection. Computational Management Science, 8(4):323–353, 2011.
[72] C. Gregory, K. Darby-Dowman, and G. Mitra. Robust optimization andportfolio selection: The cost of robustness. European Journal of Operational
Research, 212(2):417–428, 2011.
[73] I. E. Grossmann. Review of nonlinear mixed-integer and disjunctive programming techniques. Optimization and Engineering, 3(3):227–252, 2002.
[74] O. Günlük and J. Linderoth. Perspective reformulations of mixed integer
nonlinear programs with indicator variables. Mathematical Programming,
124(1-2):183–205, 2010.
[75] O. Günlük and J. Linderoth. Perspective reformulation and applications.
In Mixed Integer Nonlinear Programming, pages 61–89. Springer, 2011.
87

[76] R. Gunst and J. Webster. Regression analysis and problems of multicollinearity. Communications in Statistics, 4(3):277–292, 1975.
[77] H. Mittelmann. Benchmarks for optimization software. http://plato.
asu.edu/bench.html, Accessed 6 Junuary 2021.
[78] A. S. Hadi and R. F. Ling. Some cautionary notes on the use of principal
components regression. The American Statistician, 52(1):15, 1998.
[79] S. Han, A. Gómez, and A. Atamtürk. 2×2 convexifications for convex quadratic optimization with indicator variables. arXiv preprint
arXiv:2004.07448, 2020.
[80] W. K. Haneveld and M. H. van der Vlerk. Integrated chance constraints:
Reduced forms and an algorithm. Computational Management Science,3(4):245–269, 2006.
[81] C. Helmberg and F. Rendl. Solving quadratic (0, 1)-problems by semidefinite programs and cutting planes. Mathematical Programming, 82(3):291–
315, 1998.
[82] J.-B. Hiriart-Urruty and C. Lemaréchal. Convex Analysis and Minimization
Algorithms I. Springer, 1993.
[83] A. E. Hoerl and R. W. Kennard. Ridge regression: Biased estimation for
nonorthogonal problems. Technometrics, 12(1):55–67, 1970.
[84] R. A. Horn and C. R. Johnson. Matrix Analysis. Cambridge University
Press, 2nd edition, 2012.
[85] G. Iyengar and A. K. C. Ma. Fast gradient descent method for Mean-CVaR
optimization. Annals of Operations Research, 205(1):203–212, 2013.
[86] J. J. E. Kelley. The cutting-plane method for solving convex programs.
Journal of the Society for Industrial and Applied Mathematics, 8(4):703–
712, 1960.
[87] N. H. Jadhav, D. N. Kashid, and S. R. Kulkarni. Subset selection in multiple
linear regression in the presence of outlier and multicollinearity. Statistical
Methodology, 19:44–59, 2014.
[88] I. T. Jolliffe. A note on the use of principal components in regression.
Applied Statistics, 31(3):300, 1982.
[89] S. Joshi and S. Boyd. Sensor selection via convex optimization. IEEE
Transactions on Signal Processing, 57(2):451–462, 2009.
[90] Kaggle. S&P 500 stock data. https://www.kaggle.com/camnugent/
sandp500, Accessed 23 December 2020.
[91] M. Kaut, H. Vladimirou, S. W. Wallace, and S. A. Zenios. Stability analysis of portfolio management with conditional value-at-risk. Quantitative
Finance, 7(4):397–409, 2007.
[92] K. Kimura and H. Waki. Minimization of akaike's information criterion in
linear regression analysis via mixed integer nonlinear program. Optimization Methods and Software, 33(3):633–649, 2017.
88

[93] A. J. Kleywegt, A. Shapiro, and T. H. de Mello. The sample averageapproximation method for stochastic discrete optimization. SIAM Journal
on Optimization, 12(2):479–502, 2002.
[94] K. Kobayashi and Y. Takano. A branch-and-cut algorithm for solving
mixed-integer semidefinite optimization problems. Computational Optimization and Applications, 75(2):493–513, 2019.
[95] K. Kobayashi, Y. Takano, and K. Nakata. Bilevel cutting-plane algorithm
for cardinality-constrained mean-CVaR portfolio optimization. Journal of
Global Optimization, 81(2):493–528, 2021.
[96] K. Kobayashi, Y. Takano, and K. Nakata. Cardinality-constrained distributionally robust portfolio optimization. arXiv preprint arXiv:2112.12454,
2021.
[97] H. Konno, J. Gotoh, T. Uno, and A. Yuki. A cutting plane algorithm
for semi-definite programming problems with applications to failure discriminant analysis. Journal of Computational and Applied Mathematics,
146(1):141–154, 2002.
[98] H. Konno, N. Kawadai, and H. Tuy. Cutting plane algorithms for nonlinear
semi-definite programming problems with applications. Journal of Global
Optimization, 25(2):141–155, 2003.
[99] H. Konno and Y. Takaya. Multi-step methods for choosing the best setof variables in regression analysis. Computational Optimization and Applications, 46(3):417–426, 2010.
[100] H. Konno, H. Waki, and A. Yuuki. Portfolio optimization under lower
partial risk measures. Asia-Pacific Financial Markets, 9(2):127–140, 2002.
[101] H. Konno and R. Yamamoto. Choosing the best set of variables in regression analysis using integer programming. Journal of Global Optimization,
44(2):273–282, 2009.
[102] K. Krishnan and J. E. Mitchell. A unifying framework for several cutting
plane methods for semidefinite programming. Optimization Methods andSoftware, 21(1):57–74, 2006.
[103] A. Künzi-Bay and J. Mayer. Computational aspects of minimizing conditional value-at-risk. Computational Management Science, 3(1):3–27, 2006.
[104] C. Lim, H. D. Sherali, and S. Uryasev. Portfolio optimization by minimizing Conditional Value-at-Risk via nondifferentiable optimization. Computational Optimization and Applications, 46(3):391–415, 2010.
[105] H. Liu, X. Wang, T. Yao, R. Li, and Y. Ye. Sample average approximation
with sparsity-inducing penalty for high-dimensional stochastic programming. Mathematical Programming, 178(1-2):69–108, 2018.
[106] J. Lofberg. YALMIP : a toolbox for modeling and optimization in MATLAB. In 2004 IEEE International Conference on Robotics and Automation,
pages 284–289.
[107] S. Lotfi and S. A. Zenios. Robust VaR and CVaR optimization under joint
ambiguity in distributions, means, and covariances. European Journal ofOperational Research, 269(2):556–576, 2018.
89

[108] M. Lubin, E. Yamangil, R. Bent, and J. P. Vielma. Polyhedral approximation in mixed-integer convex optimization. Mathematical Programming,
172(1-2):139–168, 2017.
[109] M. Lichman. UCI machine learning repository. http://archive.ics.uci.
edu/ml, Accessed 5 September 2021.
[110] M. Mäkelä. Survey of bundle methods for nonsmooth optimization. Optimization Methods and Software, 17(1):1–29, 2002.
[111] N. M. Manousakis and G. N. Korres. Semidefinite programming for optimal
placement of PMUs with channel limits considering pre-existing SCADA
and PMU measurements. In 2016 Power Systems Computation Conference,
pages 1–7, 2016.
[112] E. R. Mansfield and B. P. Helms. Detecting multicollinearity. The American
Statistician, 36(3):158, 1982.
[113] R. Mansini, W. Ogryczak, and M. G. Speranza. Twenty years of linear programming based portfolio optimization. European Journal of Operational
Research, 234(2):518–535, 2014.
[114] H. Markowitz. Portfolio Selection. The Journal of Finance, 7(1):77–91,
1952.
[115] W. F. Massy. Principal components regression in exploratory statisticalresearch. Journal of the American Statistical Association, 60(309):234–256,1965.
[116] G. P. McCormick. Computability of global solutions to factorable nonconvex programs: Part I — convex underestimating problems. Mathematical
Programming, 10(1):147–175, 1976.
[117] H. Mittelmann. An independent benchmarking of SDP and SOCP solvers.
Mathematical Programming, 95(2):407–430, 2003.
[118] R. Miyashiro and Y. Takano. Mixed integer second-order cone programming
formulations for variable selection in linear regression. European Journalof Operational Research, 247(3):721–731, 2015.
[119] R. Miyashiro and Y. Takano. Subset selection by mallows’ cp : A mixed integer programming approach. Expert Systems with Applications, 42(1):325–
331, 2015.
[120] MOSEK ApS. MOSEK. https://www.mosek.com/, Accessed 20 October
2021.
[121] K. Nakata, K. Fujisawa, M. Fukuda, M. Kojima, and K. Murota. Exploiting
sparsity in semidefinite programming via matrix completion II: implementation and numerical results. Mathematical Programming, 95(2):303–327,
2003.
[122] K. Natarajan, M. Sim, and J. Uichanco. Tractable robust expected utility
and risk models for portfolio optimization. Mathematical Finance: An
International Journal of Mathematics, Statistics and Financial Economics,
20(4):695–731, 2010.

90

[123] Y. Nesterov and A. Nemirovskii. Interior-Point Polynomial Algorithms inConvex Programming. Society for Industrial and Applied Mathematics,
1994.
[124] N. Noyan, B. Balcik, and S. Atakan. A stochastic optimization model for
designing last mile relief networks. Transportation Science, 50(3):1092–
1113, 2015.
[125] W. Ogryczak and T. Śliwiński. On solving the dual for portfolio selection
by optimizing Conditional Value at Risk. Computational Optimization andApplications, 50(3):591–595, 2011.
[126] B. P. G. V. Parys, P. M. Esfahani, and D. Kuhn. From data to decisions: Distributionally robust optimization is optimal. Management Science, 67(6):3387–3402, 2020.
[127] J. Peng and Y. Xia. A new theoretical framework for k-means-type clustering. In W. Chu and L. T. Young, editors, Foundations and Advances in
Data Mining, pages 79–96. Springer, 2005.
[128] A. F. Perold. Large-scale portfolio optimization. Management Science,30(10):1143–1160, 1984.
[129] G. C. Pflug. Some remarks on the value-at-risk and the conditional valueat-risk. In S. Uryasev, editor, Nonconvex Optimization and Its Applications,pages 272–281. Springer, 2000.
[130] A. Philipp, S. Ulbrich, Y. Cheng, and M. Pesavento. Multiuser downlink
beamforming with interference cancellation using a SDP-based branch-andbound algorithm. In 2014 IEEE International Conference on Acoustics,
Speech and Signal Processing, pages 7724–7728, 2014.
[131] I. Popescu. Robust mean-covariance solutions for stochastic optimization.
Operations Research, 55(1):98–112, 2007.
[132] I. Quesada and I. Grossmann. An LP/NLP based branch and bound algorithm for convex MINLP optimization problems. Computers and Chemical
Engineering, 16(10-11):937–947, 1992.
[133] H. Rahimian and S. Mehrotra. Distributionally robust optimization: A
review. arXiv preprint arXiv:1908.05659, 2019.
[134] F. Rendl. Semidefinite relaxations for integer programming. In M. Jünger,
T. M. Liebling, D. Naddef, G. L. Nemhauser, W. R. Pulleyblank,
G. Reinelt, G. Rinaldi, and L. A. Wolsey, editors, 50 Years of Integer
Programming 1958-2008, pages 687–726. Springer, 2010.
[135] F. Rendl, G. Rinaldi, and A. Wiegele. Solving max-cut to optimality byintersecting semidefinite and polyhedral relaxations. Mathematical Programming, 121(2):307–335, 2010.
[136] R. T. Rockafellar. Convex analysis. Princeton Mathematical Series. Princeton University Press, 1970.
[137] R. T. Rockafellar and S. Uryasev. Optimization of conditional value-at-risk.
The Journal of Risk, 2(3):21–41, 2000.

91

[138] R. T. Rockafellar and S. Uryasev. Conditional value-at-risk for general loss
distributions. Journal of Banking and Finance, 26(7):1443–1471, 2002.
[139] C. Rowe and J. Maciejowski. An efficient algorithm for mixed integer
semidefinite optimisation. In Proceedings of the 2003 American Control
Conference, pages 4730–4735, 2003.
[140] T. Sato, Y. Takano, and R. Miyashiro. Piecewise-linear approximation
for feature subset selection in a sequential logit model. Journal of the
Operations Research Society of Japan, 60(1):1–14, 2017.
[141] T. Sato, Y. Takano, R. Miyashiro, and A. Yoshise. Feature subset selection for logistic regression via mixed integer optimization. Computational
Optimization and Applications, 64(3):865–880, 2016.
[142] H. Scarf. A min-max solution of an inventory problem. Studies in The
Mathematical Theory of Inventory and Production, pages 201–209, 1958.
[143] N. Schwertman and D. Allen. Smoothing an indefinite variance-covariancematrix. Journal of Statistical Computation and Simulation, 9(3):183–194,
1979.
[144] A. Shapiro. Monte carlo sampling methods. In Handbooks in Operations
Research and Management Science, pages 353–425. Elsevier, 2003.
[145] A. Shapiro, D. Dentcheva, and A. Ruszczyński. Lectures on Stochastic
Programming: Modeling and Theory. Society for Industrial and Applied
Mathematics, 2009.
[146] A. Skajaa, E. D. Andersen, and Y. Ye. Warmstarting the homogeneousand self-dual interior point method for linear and conic quadratic problems.
Mathematical Programming Computation, 5(1):1–25, 2012.
[147] R. Sotirov. SDP relaxations for some combinatorial optimization problems.
In M. F. Anjos and J. B. Lasserre, editors, Handbook on Semidefinite, Conic
and Polynomial Optimization, pages 795–819. Springer, 2012.
[148] J. F. Sturm. Using SeDuMi 1.02, a Matlab toolbox for optimization oversymmetric cones. Optimization Methods and Software, 11(1-4):625–653,
1999.
[149] Y. Takano, K. Nanjo, N. Sukegawa, and S. Mizuno. Cutting plane algorithms for mean-CVaR portfolio optimization with nonconvex transaction
costs. Computational Management Science, 12(2):319–340, 2014.
[150] A. Takeda and T. Kanamori. A robust approach based on conditional
value-at-risk measure to statistical learning problems. European Journal ofOperational Research, 198(1):287–296, 2009.
[151] R. Tamura, K. Kobayashi, Y. Takano, R. Miyashiro, K. Nakata, andT. Matsui. Best subset selection for eliminating multicollinearity. Journal of the Operations Research Society of Japan, 60(3):321–336, 2017.
[152] K. Tanaka and R. Yamamoto. Identification of best discrimination surface
by mixed-integer semi-definite programming for support vector machine.
International Journal of Financial Engineering, 2021. https://doi.org/
10.1142/S2424786321500420.
92[153] J. A. Taylor, N. Luangsomboon, and D. Fooladivanda. Allocating sensors
and actuators via optimal estimation and control. IEEE Transactions on
Control Systems Technology, 25(3):1060–1067, 2017.
[154] R. Tibshirani. Regression shrinkage and selection via the lasso. Journalof the Royal Statistical Society: Series B (Methodological), 58(1):267–288,
1996.
[155] M. J. Todd. Semidefinite optimization. Acta Numerica, 10:515–560, 2001.
[156] K. C. Toh, M. J. Todd, and R. H. Tütüncü. SDPT3 — A Matlab software
package for semidefinite programming, version 1.3. Optimization Methods
and Software, 11(1-4):545–581, 1999.
[157] X. Tong, L. Qi, F. Wu, and H. Zhou. A smoothing method for solving portfolio optimization with CVaR and applications in allocation of generation
asset. Applied Mathematics and Computation, 216(6):1723–1740, 2010.
[158] M. Torchio, L. Magni, and D. M. Raimondo. A mixed integer SDP approach for the optimal placement of energy storage devices in power gridswith renewable penetration. In Proceedings of the 2015 American Control
Conference, pages 3892–3897, 2015.
[159] S. F. Tóth, M. E. McDill, N. Könnyü, and S. George. Testing the useof lazy constraints in solving area-based adjacency formulations of harvest
scheduling models. Forest Science, 59(2):157–176, 2013.
[160] R. Tütüncü and M. Koenig. Robust asset allocation. Annals of Operations
Research, 132(1-4):157–187, 2004.
[161] L. Vandenberghe and S. Boyd. Semidefinite programming. SIAM Review,38(1):49–95, 1996.
[162] T. Westerlund and F. Pettersson. An extended cutting plane method forsolving convex MINLP problems. Computers and Chemical Engineering,
19:131–136, 1995.
[163] H. P. Williams. Model Building in Mathematical Programming. John Wiley& Sons, 2013.
[164] H. Wold. Estimation of principal components and related models by iterative least squares. In P. Krishnaiaah, editor, Multivariate Analysis, pages
391–420. Academic Press, 1966.
[165] S. Wold, A. Ruhe, H. Wold, and I. W. J. Dunn. The collinearity problem
in linear regression. the partial least squares (PLS) approach to generalized
inverses. SIAM Journal on Scientific and Statistical Computing, 5(3):735–
743, 1984.
[166] S. J. Wright. Primal-Dual Interior-Point Methods. Society for Industrial
and Applied Mathematics, 1997.
[167] Yahoo Japan Corporation. Yahoo! Finance. https://finance.yahoo.co.
jp, Accessed 7 October 2021.

93[168] M. Yamashita, K. Fujisawa, M. Fukuda, K. Kobayashi, K. Nakata, and
M. Nakata. Latest developments in the sdpa family for solving large-scale
sdps. In M. Anjos and J. Lasserre, editors, Handbook on Semidefinite, Conic
and Polynomial Optimization, pages 687–713. Springer, 2012.
[169] M. Yamashita, K. Fujisawa, and M. Kojima. Implementation and evaluation of SDPA 6.0 (semidefinite programming algorithm 6.0). Optimization
Methods and Software, 18(4):491–505, 2003.
[170] R. Yokoyama, Y. Shinano, S. Taniguchi, M. Ohkura, and T. Wakui. Optimization of energy supply systems by MILP branch and bound method
in consideration of hierarchical relationship between design and operation.
Energy Conversion and Management, 92:92–104, 2015.
[171] K. Yonekura and Y. Kanno. Global optimization of robust truss topology
via mixed integer semidefinite programming. Optimization and Engineering, 11(3):355–379, 2010.
[172] Y. Zhang, S. Shen, and S. A. Erdogan. Solving 0–1 semidefinite programsfor distributionally robust allocation of surgery blocks. Optimization Letters, 12(7):1503–1521, 2018.
[173] X. Zheng, X. Sun, and D. Li. Improving the performance of MIQP
solvers for quadratic programs with cardinality and minimum thresholdconstraints: A semidefinite program approach. INFORMS Journal on Computing, 26(4):690–703, 2014.
[174] H. Zou and T. Hastie. Regularization and variable selection via the elastic
net. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 67(2):301–320, 2005.
[175] S. Zymler, D. Kuhn, and B. Rustem. Distributionally robust joint chanceconstraints with second-order moment information. Mathematical Programming, 137(1-2):167–198, 2011.

94

Appendix A
ProofsIn this chapter, we give the complete proofs of Theorem 4.2, Lemma 4.6, Theorem 5.1, and Theorem 5.4.

A.1

Proof of Theorem 4.2

The Lagrange function of Problem (4.12) is expressed as
L(x, P , Q, p, q, r, s; α, B, β, η, Λ, λ, ν, π, ρ)
1 >:=x x + (κ2 Σ̂ − µ̂µ̂> ) • Q + r + Σ̂ • P − 2µ̂> p + κ1 s
2γ
− α> (p + q/2 + Qµ̂)
! 

X
B (`)
β (`)
Q
q/2 + a(`) Zx/2
−•>(q/2 + a(`) Zx/2)>
r + b(`)
η (`)
β (`)
`∈[L]

 

Λ λ
P p
−
•
− π(1> Zx − 1) − ρ> Zx,
λ> ν
p> s
!


B (`)
β (`)
Λ λ
>
 O (` ∈ [L]), O, π ∈ R, and
λ> ν
β (`)
η (`)
ρ ≥ 0 are Lagrange multipliers. The Lagrange dual of Problem (4.12) is thenposed aswhere α ∈ RN ,

max

min

α,B,β,η,Λ,λ,ν,π,ρ x,P ,Q,p,q,r,s

L(x, P , Q, p, q, r, s; α, B, β, η, Λ, λ, ν, π, ρ). (A.1)

Now, let us focus on the inner minimization problem:minx,P ,Q,p,q,r,s

L(x, P , Q, p, q, r, s; α, B, β, η, Λ, λ, ν, π, ρ).

(A.2)

Note that Problem (A.2) is an unconstrained convex quadratic optimization problem and its objective function is linear in (P , Q, p, q, r, s). Since Problem (A.2)
must be bounded, the Lagrange multipliers are required to satisfy the following

95

conditions:∇P L = Σ̂ − Λ = O,

(A.3)

X
1
∇Q L = κ2 Σ̂ − µ̂µ̂> − (µ̂α> + αµ̂> ) −B (`) = O,
2

(A.4)

`∈[L]

∇p L = −2µ̂ − α − 2λ = 0,X1β (`) = 0,
∇q L = − α −
2
`∈[L]
X
η (`) = 0,
∇r L = 1 −

(A.5)
(A.6)
(A.7)

`∈[L]

∇s L = κ1 − ν = 0.

(A.8)

Also, the following optimality condition should be satisfied:X
1∇x L = x − Z a(`) β (`) + π1 + ρ = 0.
γ

(A.9)

`∈[L]

According to the conditions (A.3)–(A.9), the optimal objective value of Problem(A.2) is calculated asXγ− ω>Z 2ω −
η (`) b(`) + π,
(A.10)
2
`∈[L]

where
ω=

X

a(`) β (`) + π1 + ρ.

`∈[L]Since z ∈ {0, 1}N , it holds that ω > Z 2 ω = z > (ω ◦ ω). Thus, the Lagrange
dual (A.1) of Problem (4.12) is formulated as follows:Xγ
maximize− z > (ω ◦ ω) −η (`) b(`) + π(A.11a)ω,α,B,β,η,Λ,λ,ν,π,ρ
2
`∈[L]
X
(`) (`)
subject to ω =a β + π1 + ρ,
(A.11b)`∈[L]Σ̂ − Λ = O,

(A.11c)1κ2 Σ̂ − µ̂µ̂> − (µ̂α> + αµ̂> ) −2

XB (`) = O,

(A.11d)`∈[L]− 2µ̂ − α − 2λ = 0,X1β (`) = 0,
− α−
2
`∈[L]
X
1−
η (`) = 0,

(A.11e)(A.11f)(A.11g)`∈[L]κ1 − ν = 0,
β (`)

>

η (`)

β (`)


(A.11h)

B (`)
Λ λ
λ> ν

!
O

(` ∈ [L]),(A.11i) O,(A.11j)ρ ≥ 0.

(A.11k)96

We then delete ρ, Λ α, and ν from Problem (A.11) by substituting Eqs. (A.11b),(A.11c), (A.11f), and (A.11h) into other constraints. We now obtain the desired
formulation (4.13).
To prove the strong duality, we next show that both the primal problem (4.12)
and the dual problem (4.13) are strictly feasible [7]. Let us set x = z/(1> z),
q = 0, r > − min`∈[L] {b(`) }, and s > 0 in the primal problem (4.12). Then, there
exists Q satisfyingQ
q/2 + a(`) Zx/2
 O (∀` ∈ [L]).
(q/2 + a(`) Zx/2)>
r + b(`)
We next set p = −Qµ̂, and then there exists P satisfying


P p
 O,
p> s
which implies that the primal problem (4.12) is strictly feasible.
For the dual problem (4.13), we set λ = 0 and (B, β, η) as
B (`) =

1(κ2 Σ̂ + µ̂µ̂> ),L

β (`) =1µ̂,Lη (`) =

1
L(∀` ∈ [L])such that the constraints (4.13c)–(4.13e) are satisfied. We next set ω =π1. As for the constraint (4.13f), we have
!


B (`)
β (`)
1 κ2 Σ̂ + µ̂µ̂> µ̂>= O (∀` ∈ [L]),
µ̂>
1
L
β (`)
η (`)

P

`∈[L] a(`) β (`) +(A.12)

because the Schur complement is given from Assumption 4.1 by κ
12κ2 Σ̂ + µ̂µ̂> − µ̂µ̂> =Σ̂  O.
L
L
Also for the constraint (4.13g), Assumption 4.1 ensures that

 

Σ̂ λ
Σ̂ 0
=
 O,
0 κ1
λ> κ1

(A.13)(A.14)

which shows that the dual problem (4.13) is strictly feasible.

A.2

Proof of Lemma 4.6

For notational simplicity, we assume η (`) > 0 for all ` ∈ [L] without loss ofgenerality. As in the proof of Theorem 4.8, we use the notations (4.29) and(4.30) under the assumption (4.28). We define φ(`) := β̄ (`) − η (`) µ̂ for ` ∈ [L],and due to the definition (4.22), we have
(`)

(`)

φ1 = β1 − η (`) µ̂1 ,


(`)
(`)
(`)
φ2 = Σ̂21 (Σ̂11 )−1 (β1 − η (`) µ̂1 ) + η (`) µ̂2 − η (`) µ̂2 = Ψφ1 ,(A.15)(A.16)

whereΨ := Σ̂21 (Σ̂11 )−1 .

97(A.17)

Substituting them into the right-hand side of Eq. (4.23), we obtain
X 1(β̄ (`) − η (`) µ̂)(β̄ (`) − η (`) µ̂)>
(`)
η
`∈[L]
! 

(`)
(`)
(`)
(`)
X 1
φ1 (φ1 )> φ1 (φ2 )>
Φ
ΦΨ>
=
=
,
(`) >
(`)
(`)
ΨΦ ΨΦΨ>η (`) φ(`)
φ2 (φ2 )>
2 (φ1 )
`∈[L]
where
Φ :=

X 1 (`) (`)φ (φ1 )> .
η (`) 1

(A.18)(A.19)`∈[L]

From the constraint (4.21f), we can use the following Schur complement property (see, e.g., Section A.5.5 [36]):


(`)
(`)
B11
β1
1 (`)  (`) >
(`)
>

  O ⇒ B11
ββ1(∀` ∈ [L]).
(A.20)(`)η (`) 1β1
η (`)
It then follows that
>
X 1 (`) (`)
X (`)
X (`)

Φ=
β1 (β1 )> − 
β1  µ̂>
β1  + µ̂1 µ̂>
1 − µ̂1
1
(`)
η
`∈[L]

`∈[L]

`∈[L]

∵ Eqs. (4.21e), (A.15), (A.19)

X (`)
X (`)
X (`)


B11 + µ̂1 µ̂>
β1  − β1  µ̂>∵ Eq. (A.20)
1 − µ̂1
1
>



`∈[L]

= κ2 Σ̂11 .

`∈[L]



`∈[L]

∵ Eq. (4.21c)(A.21)

From Assumption 4.1, we can also use the following Schur complement property (see, e.g., Section A.5.5 [36]):Σ̂11 Σ̂>21  O ⇒ Σ̂ − Σ̂ (Σ̂ )−1 Σ̂>  O.
22
21
11
21
Σ̂21 Σ̂22(A.22)

We now obtain the desired result:


X 1
Φ
ΦΨ>
(`)
(`)
(`)
(`)
>
(
β̄
−
η
µ̂)(
β̄
−
η
µ̂)
=
∵ Eq. (A.18)
ΨΦ ΨΦΨ>
η (`)
`∈[L]


Σ̂11
Σ̂>
21
 κ2
∵ Eqs. (A.17), (A.21)
Σ̂21 Σ̂21 (Σ̂11 )−1 Σ̂>21Σ̂11 Σ̂>21
 κ2. ∵ Eq. (A.22)
Σ̂21 Σ̂22

98

A.3

Proof of Theorem 5.1

Problem (5.7) is formulated as follows:
1 >x x+a+v
a,q,v,x2γ
1 X
subject to v ≥p s qs ,1−βf (z) = minimize

qs ≥ −(r(A.23a)(A.23b)s∈S(s) >) Zx − a (∀s ∈ S),(A.23c)CZx ≤ d,(A.23d)1> Zx = 1,

(A.23e)Zx ≥ 0,(A.23f)q ≥ 0.

(A.23g)

The Lagrange function of Problem (A.23) is expressed as
1 >1 X
x x+a+v−η v−
p s qsL(a, q, v, x; η, α, ζ, λ, π, θ) :=
2γ
1−β
s∈S

X 
(s) >
−
αs qs + (r ) Zx + a

!

s∈S



− ζ > (d − CZx) − λ 1> Zx − 1 − π > Zx − θ > q,
where η ≥ 0, α := (αs )s∈S ≥ 0, ζ ≥ 0, λ ∈ R, π ≥ 0 and θ ≥ 0 are Lagrange
multipliers. Then, the Lagrange dual problem of Problem (A.23) is posed as:
max

min

η≥0, α≥0, ζ≥0, a∈R, q∈RS ,λ∈R, π≥0, θ≥0 v∈R, x∈RN

L(a, q, v, x; η, α, ζ, λ, π, θ).

(A.24)

Recall that Problem (A.23) is feasible. Also, the objective function is proper
convex, and all the constraints are linear in Problem (A.23). Then, the strong
duality holds; see, for example, Section 5.2.3 in Boyd and Vandenberghe [36]. As
a result, f (z) is equal to the optimal objective value of Problem (A.24). Now,
let us focus on the inner minimization problem:mina∈R, q∈RS ,v∈R, x∈RN

L(a, q, v, x; η, α, ζ, λ, π, θ).

(A.25)

Note that Problem (A.25) is an unconstrained convex quadratic optimization
problem and its objective function is linear in (a, q, v). Because Problem (A.25)
must be bounded, the Lagrange multipliers are required to satisfy the following
conditions:X∇a L = 1 −αs = 0,(A.26)s∈Sηp − α − θ = 0,
∇q L =
1−β∇v L = 1 − η = 0.
Also, the following optimality condition should be satisfied:
!
X
1
(s)
>
∇x L = x − Z
αs r − C ζ + λ1 + π = 0.
γ
s∈S

99

(A.27)(A.28)(A.29)

According to Eqs. (A.26), (A.27), (A.28), and (A.29), the optimal objective
value of Problem (A.25) is calculated asγ− ω > Z 2 ω − d> ζ + λ,2where ω ∈ RN is a vector of auxiliary decision variables satisfyingXω=
αs r (s) − C > ζ + λ1 + π.
s∈SBecause z ∈ {0, 1}N , it holds that ω > Z 2 ω = z > (ω ◦ ω). Therefore, theLagrange dual problem (A.24) is formulated as follows:γ >z (ω ◦ ω) − d> ζ + λ
α,ζ,λ,ω
2
X
subject to ω ≥αs r (s) − C > ζ + λ1,−

f (z) = maximizes∈SX

αs = 1,s∈S

ps(∀s ∈ S),1−βα ≥ 0, ζ ≥ 0,
αs ≤where we substitute η = 1, and eliminate the nonnegative variables π and θ.

A.4

Proof of Theorem 5.4

Problem (5.23) is formulated as follows:
1 >x x+a+v
2γ
1 X
subject to v ≥ps (−(r (s) )> Zx − a)1−βfK (z) = minimize(A.30a)a,v,x(∀J ∈ K),

(A.30b)s∈Jv ≥ 0,(A.30c)CZx ≤ d,(A.30d)>1 Zx = 1,(A.30e)Zx ≥ 0.

(A.30f)

The Lagrange function of Problem (A.30) is expressed as
1 X  (s) >v+
ps (r ) Zx + a
1−β
J ∈K
s∈J


− ξv − ζ > (d − CZx) − λ 1> Zx − 1 − π > Zx,

X
1 >L(a, v, x; α, ξ, ζ, λ, π) :=x x+a+v−
αJ
2γ

where α := (αJ )J ∈K ≥ 0, ξ ≥ 0, ζ ≥ 0, λ ∈ R and π ≥ 0 are Lagrange
multipliers. Recall that Problem (A.30) is feasible. Then, the strong duality
holds as well as the proof of Theorem 5.1, and fK (z) is equal to the optimal
objective value of the following Lagrange dual problem:
max

min

α≥0, ξ≥0, ζ≥0, a∈R, v∈R, x∈RN
λ∈R, π≥0100

L(a, v, x; α, ξ, ζ, λ, π).

(A.31)!

Now, let us consider the inner minimization problem:
L(a, v, x; α, ξ, ζ, λ, π).

mina∈R, v∈R, x∈RN

(A.32)

Similarly to the proof of Theorem 5.1, the following conditions must be satisfied:X
1 X
αJ
ps = 0,
1−β
J ∈K
s∈J
X
∇v L = 1 −αJ − ξ = 0,
∇a L = 1 −(A.33)(A.34)J ∈K

1∇x L = x − Z
γ

X
1 X
αJ
ps r (s) − C > ζ + λ1 + π
1−β
J ∈K

!
= 0.

(A.35)s∈J

According to Eqs. (A.33), (A.34), and (A.35), the optimal objective value of
Problem (A.32) is calculated asγ− ω > Z 2 ω − d> ζ + λ,2where ω ∈ RN is a vector of auxiliary decision variables satisfying
ω=

X
1 X
ps r (s) − C > ζ + λ1 + π.
αJ
1−β
J ∈K

s∈J

Because z ∈ {0, 1}N , it holds that ω > Z 2 ω = z > (ω ◦ ω). Therefore, theLagrange dual problem (A.31) is formulated as follows:γ >z (ω ◦ ω) − d> ζ + λ
α,ζ,λ,ω
2X
1 X
ps r (s) − C > ζ + λ1,
subject to ω ≥
αJ
1−β
s∈J
J ∈K
X
αJ ≤ 1,fK (z) = maximize

−

J ∈K

X
J ∈K

αJ

X

ps = 1 − β,s∈Jα ≥ 0, ζ ≥ 0,
where we eliminate the nonnegative variables ξ and π.

101

Appendix B
Detailed Experimental Results
In this chapter, we show detailed experimental results for sensitivity analysis
reported in Section 5.4.5. Table B.1 gives the numerical results for the four
datasets (ind49, sbm100, port1, √and port5) with the cardinality parameter k ∈{5, 10, 15}. Here, we set γ = 10/ N for the `2 -regularization term and S = 105
as the number of scenarios. Table B.1 shows that BCP and BCPc were almost
always faster than the other methods when k ∈ {10, 15}. Also, BCP and BCPc
tended to be faster with larger k. In the case of the port5 dataset, for example,
BCP solved the problem with k = 15 in 73.7 s, whereas it failed to complete thecomputation within 3600 s with k = 5.
Figures B.1–B.4 show the numerical results for the four datasets (ind49,√
sbm100, port1, and port5) with the regularization parameter γ ∈ {100.2i / N |i = 0, 1, . . . , 20}. Here, we set k = 10 for the cardinality constraint and S = 105 as
the number of scenarios. We can see from Figures B.1–B.4 that the performancesof BCP and BCPc were not greatly affected by γ, and they were generally fasterthan the other methods when they completed the computations within 3600 s.
In addition, BCP and BCPc tended to be faster with smaller γ for each dataset.
Figure B.5 shows “Reg+CVaR” and “CVaR” defined in Eq. (5.31) for thethree datasets√ (ind49, sbm100, and port1) with the regularization parameterγ ∈ {100.2i / N | i = 0, 1, . . . , 20}. For each dataset,√ the objective value (i.e.,“Reg+CVaR”) was close to CVaR when γ√≥ 102 / N , and CVaR almost converged to its minimum value when γ ≥ 10/ N .

102

Table B.1: Numerical results √for the four datasets (port1, ind49, sbm100, andport5) with (S, γ) = (105 , 10/ N ) for k ∈ {5, 10, 15}
Data
port1

N
30

5

10

15

ind49

49

5

10

15

sbm100

100

5

10

15

port5

225

BigM

k

5

10

15

Time
Obj
Gap(%)#Nodes
#Cuts
Time
Obj
Gap(%)#Nodes
#Cuts
Time
Obj
Gap(%)#Nodes
#Cuts
Time
Obj
Gap(%)#Nodes
#Cuts
Time
Obj
Gap(%)#Nodes
#Cuts
Time
Obj
Gap(%)#Nodes
#Cuts
Time
Obj
Gap(%)#Nodes
#Cuts
Time
Obj
Gap(%)#Nodes
#Cuts
Time
Obj
Gap(%)#Nodes
#Cuts
Time
Obj
Gap(%)#Nodes
#Cuts
Time
Obj
Gap(%)#Nodes
#Cuts
Time
Obj
Gap(%)#Nodes
#Cuts

Persp

Lift

Cut

Lift

Cut

643.6
4.436
0.00
32
—
468.6
4.264
0.00
1
—
469.4
4.222
0.00
1
—
733.7
3.443
0.00
83
—
709.8
3.379
0.00
31
—
555.2
3.366
0.00
1
—
1260.7
4.367
0.00
92
—
1107.1
4.364
0.00
1
—
1051.6
4.364
0.00
1
—
>3600
∞
100.00>0—>3600∞
100.00>0—>3600∞
100.00
>0
—

221.0
4.436
0.00
2636
—
163.3
4.264
0.00
2433
—
120.1
4.222
0.00
879
—
160.0
3.443
0.00
2475
—
125.0
3.379
0.00
4101
—
128.2
3.366
0.00
3392
—
179.8
4.367
0.00
1092
—
65.1
4.364
0.00
412
—
43.2
4.364
0.00
90
—
>3600
3.491
8.72
>41,433
—
>3600
3.246
3.98
>20,833—
1154.4
3.109
0.00
8454
—

1496.1
4.436
0.00
13
—
1933.2
4.264
0.00
1
—
3220.8
4.222
0.00
1
—
>3600
3.443
0.00
>61
—
>3600
3.380
0.00
>1
—
>3600
3.366
0.00
>1
—
2941.5
4.367
0.00
3
—
>3600
4.364
0.00
>1
—
3438.1
4.364
0.00
1
—
>3600
3.755
14.25
>42
—
>3600
3.143
0.44
>61
—
>3600
3.109
0.01
>1
—

200.0
4.436
0.00
3542
—
15.8
4.264
0.00
253
—
217.1
4.222
0.00
4290
—
434.3
3.443
0.00
5172
—
220.0
3.379
0.00
3489
—
76.7
3.366
0.00
1043
—
147.8
4.367
0.00
1866
—
84.6
4.364
0.00
4676
—
44.4
4.364
0.00
5738
—
2263.93.295
0.00
19,851
—
>3600
3.276
4.29
>16,790—1116.1
3.109
0.00
5696
—

103

CP

BCP

BCPc

644.0
4.436
0.00
445
105
265.8
4.264
0.00
39
25
157.6
4.222
0.00
3
8
1455.9
3.443
0.00
2614
226
662.2
3.379
0.00
250
68
448.6
3.366
0.00
67
39
690.8
4.367
0.00
116
53
482.7
4.364
0.00
10
17
437.5
4.364
0.00
9
12
>3600
3.397
8.52
>1390
>249
>3600
3.173
2.05
>2140
>219
>3600
3.170
1.97
>2516
>218

279.6
4.436
0.00
—
59
39.5
4.264
0.00
—
1
39.1
4.222
0.00
—
1
638.5
3.443
0.00
—
132
43.1
3.379
0.00
—
1
48.7
3.366
0.00
—
1
38.3
4.367
0.00
—
5
27.6
4.364
0.00
—
1
29.5
4.364
0.00
—
1
>3600
3.295
0.72
—
>667
2351.6
3.138
0.00
—
352
73.7
3.109
0.00
—
1

319.8
4.436
0.00
76
62
64.6
4.264
0.00
0
7
89.8
4.222
0.00
0
7
907.5
3.443
0.00
970
172
103.6
3.379
0.00
0
12
82.3
3.366
0.00
0
7
48.7
4.367
0.00
7
12
36.6
4.364
0.00
0
7
48.0
4.364
0.00
0
7
>3600
3.295
3.37
>10,767>861
3300.8
3.138
0.00
5697
519
121.4
3.109
0.00
0
11

BigM+Lift
BigM+Cut

Persp+Lift
Persp+Cut

CP
BCP

BCPc

0.012

103

102

Persp+Lift
Persp+Cut

CP
BCP

BCPc

0.008
0.006
0.004
0.002

101 p

10 0 / Np

10 1 / Np

10 2 / Nγp

10 3 / N

0.0000 p

p

10 4 / N

10 / N(a) Time

104

BigM+Lift
BigM+Cut

Persp+Lift
Persp+Cut

CP

BCPc

102
101

10 / Np

10 1 / Np

10 2 / Nγp

10 2 / Nγp

10 3 / N

CP

102

103

1000 pp

10 1 / Np

10 3 / Np

10 4 / N(b) Gap(%)

Number of cuts

Number of nodes

BigM+Lift
BigM+Cut

0.010

Relative gap (%)

Computation time (s)

104

BCPc

p

10 3 / N

101

1000 pp

10 4 / N

BCP

10 / Np

10 1 / N(c) #Nodes

10 2 / Nγpp

10 4 / N(d) #Cuts√
Figure B.1: Results for port1 with (S, k) = (105 , 10) for γ ∈ {100.2i / N | i =0, 1, . . . , 20}

BigM+Lift
BigM+Cut

Persp+Lift
Persp+Cutp

10 2 / N

CP
BCP

BCPc

Relative gap (%)

Computation time (s)

104

103

102

1001 p10 / N

10 1 / Npγp

10 3 / Np

10 4 / N

0.016
0.014
0.012
0.010
0.008
0.006
0.004
0.002
0.0000 p

10 / N(a) Time

105

BigM+Lift
BigM+Cut

Persp+Lift
Persp+Cut

Persp+Lift
Persp+Cutpp

10 1 / N

10 2 / N

γ

CP
BCP

BCPc

p

10 3 / Np

10 4 / N(b) Gap(%)
CP

BCPc

CP

103

Number of cuts

104

Number of nodes

BigM+Lift
BigM+Cut

103
102

BCP

BCPc

p

10 3 / N

102

101

101
1000 p10 / Np

10 1 / Np

10 2 / Nγp

10 3 / Np

10 4 / N1000 p10 / N(c) #Nodesp

10 1 / N

10 2 / Nγpp

10 4 / N(d) #Cuts√
Figure B.2: Results for ind49 with (S, k) = (105 , 10) for γ ∈ {100.2i / N | i =0, 1, . . . , 20}

104

104

BigM+Lift
BigM+Cut

Persp+Lift
Persp+Cut

CP
BCP

BCPc

100

BigM+Lift
BigM+Cut

Persp+Lift
Persp+Cutp

10 2 / N

BCPc

Relative gap (%)

Computation time (s)

80

103

102

1001 p10 / Np

10 1 / Np

10 2 / Nγp

10 3 / N

60
40
20
0p

p

10 0 / N

10 4 / N

105

BigM+Lift
BigM+Cut

Persp+Lift
Persp+Cutp

10 1 / N(a) Timeγp

10 3 / Np

10 4 / N(b) Gap(%)
CP

BCPc

CP

103

Number of cuts

104

Number of nodes

CP
BCP

103
102

BCP

BCPc

p

10 3 / N

102

101

101
1000 p10 / Np

10 1 / Np

10 2 / Nγp

10 3 / N1000 pp

10 4 / N

10 / Np

10 1 / N(c) #Nodes

10 2 / Nγpp

10 4 / N(d) #Cuts√
Figure B.3: Results for sbm100 with (S, k) = (105 , 10) for γ ∈ {100.2i / N | i =0, 1, . . . , 20}

BigM+Lift
BigM+Cut

Persp+Lift
Persp+Cut

CP
BCP

BCPc

100

BigM+Lift
BigM+Cut

Persp+Lift
Persp+Cutp

10 2 / N

CP
BCP

BCPc

80

Relative gap (%)

Computation time (s)

104

103

102

1001 p10 / Np

10 1 / Np

10 2 / Nγp

10 3 / Np

10 4 / N

60
40
20
0p

10 0 / N(a) Time

106

BigM+Lift
BigM+Cut

Persp+Lift
Persp+Cutp

10 2 / Np

10 1 / Nγp

10 3 / Np

10 4 / N(b) Gap(%)
CP

BCPc

CP

103

BCP

BCPc

p

10 3 / N

Number of cuts

Number of nodes

105
104
103
102
101
1000 p10 / N

10 1 / Npγp

10 3 / Np

10 4 / N

1002 p10 / N(c) #Nodesp

10 1 / N

10 2 / Nγpp

10 4 / N(d) #Cuts√
Figure B.4: Results for port5 with (S, k) = (105 , 10) for γ ∈ {100.2i / N | i =0, 1, . . . , 20}105Reg+CVaR

CVaR

p

10 3 / N

Objective value

4.6
4.5
4.4
4.3
4.2 p

10 0 / Np

10 1 / N

10 2 / Nγpp

10 4 / N(a) port1
3.9

Reg+CVaR

CVaR

p

10 3 / N

Objective value

3.8
3.7
3.6
3.5
3.4
3.3 p

10 0 / Np

10 1 / N

10 2 / Nγpp

10 4 / N(b) ind49

Reg+CVaR

CVaR

p

10 3 / N

Objective value

5.0
4.8
4.6
4.4
4.2 p

10 0 / Np

10 1 / N

10 2 / Nγpp

10 4 / N(c) sbm100

Figure B.5: Objective value (5.31) for the three√ datasets (port1, ind49, andsbm100) with (S, k) = (105 , 10) for γ ∈ {100.2i / N | i = 0, 1, . . . , 20}106