# Prescriptive PCA: Dimensionality Reduction for Two-stage Stochastic
Optimization

Long He
School of Business
George Washington University
Washington, DC 20052
<EMAIL_ADDRESS>
&Ho-Yin Mak
McDonough School of Business
Georgetown University
Washington, DC 20057
<EMAIL_ADDRESS>

###### Abstract

In this paper, we consider the alignment between an upstream dimensionality
reduction task of learning a low-dimensional representation of a set of high-
dimensional data and a downstream optimization task of solving a stochastic
program parameterized by said representation. In this case, standard
dimensionality reduction methods (e.g., principal component analysis) may not
perform well, as they aim to maximize the _amount_ of information retained in
the representation and do not generally reflect the _importance_ of such
information in the downstream optimization problem. To address this problem,
we develop a _prescriptive_ dimensionality reduction framework that aims to
minimize the degree of suboptimality in the optimization phase. For the case
where the downstream stochastic optimization problem has an expected value
objective, we show that prescriptive dimensionality reduction can be performed
via solving a distributionally-robust optimization problem, which admits a
semidefinite programming relaxation. Computational experiments based on a
warehouse transshipment problem and a vehicle repositioning problem show that
our approach significantly outperforms principal component analysis with real
and synthetic data sets.

##  1 Introduction

Common practice of data analytics to business planning and operations usually
involves a _pipeline_ structure, consisting of a sequential set of processes
converting collected raw data to business decisions to be implemented. For
prescribing operational decisions, mathematical programming models (e.g., for
production planning) often fit in the downstream stages of the pipeline, with
inputs fed from upstream processes, e.g., learning statistical models for
product demand with machine learning methods. This approach where learning and
optimization are performed separately and sequentially, while intuitive to
implement, can be suboptimal as the learning phase often does not account for
how its outputs are used as inputs in the downstream optimization phase [1].

In this paper, we consider a prescriptive analytics pipeline that aims to
prescribe a set of optimal decisions in a stochastic optimization problem
given high-dimensional input data. Specifically, we consider two sequential
phases: (1) a dimensionality reduction phase that learns a low-dimensional
representation of the probability distribution that generates the high-
dimensional data; and (2) a stochastic optimization phase that prescribes an
optimal solution to a stochastic program whose parameter uncertainty is
governed by said distribution. Such a scenario has a variety of practical
applications. For example, an urban mobility platform may optimize
repositioning of its fleet based on a probability distribution of travel
demand learned from data. In such applications, while the raw origin-
destination travel data could be high dimensional, the variations in the
underlying demand distribution are often governed by a low-dimensional
structure involving a smaller number of factors.

Our research objectives are as follows. First, we investigate the limitation
of using standard dimensionality reduction methods (such as principal
component analysis, PCA) in such a two-phase scenario. Particularly, we
demonstrate that standard PCA can fail to identify factors that are relevant
to the downstream stochastic program. Second, observing the shortcoming of PCA
in aligning the two phases of the problem, we propose a prescriptive
alternative that learns a low-dimensional representation of data that
minimizes downstream suboptimality. Using distributionally-robust optimization
techniques, we show that the problem of minimizing a proxy (upper bound) for
the degree of suboptimality in the downstream stochastic program can be
formulated as a bi-convex problem, to which a local optimal solution can be
found via an alternating algorithm that involves solving a semidefinite
program in each iteration. Third, using synthetic and real data, we
investigate the effectiveness of our approach based on a joint production and
inventory allocation problem for a supply chain network, and a vehicle
repositioning problem for an urban mobility platform.

###  1.1 Related Literature

Our work is related to the literature on mathematical programming and machine
learning.

_Interface of Machine Learning and Optimization._ Aligned with the prevalent
application of machine learning, numerous researchers have studied end-to-end
techniques that integrate machine learning and mathematical programming in
prescriptive settings. A prevalent approach that integrates machine learning
and mathematical programming in prescriptive settings is “predict-then-
optimize”, which involves first making accurate predictions from data, using
machine learning tools, and then solving optimization problems taking such
predictions as input [2, 3]. Noting that the criteria for improving
predictions and improving decisions are often not aligned, a growing stream of
recent work looks into more integrated, end-to-end approaches to guide
decisions directly from historical data that leverage contextual information
[4, 5, 6]. To address the potential misalignment between the loss function for
a predictive model and the objective function in the downstream optimization
model, [1] define suitable loss functions that take the downstream
optimization problem into account when measuring errors in predictions. They
further elaborate on how to train predictive models using these loss functions
under this smart predict-then-optimize framework. [7] further provide
generalization bounds for the framework. These mentioned papers focus on
integrating predictive models, trained through supervised learning, with a
downstream optimization phase. Our paper, on the other hand, considers the
upstream phase of dimensionality reduction, a class of unsupervised learning.
As opposed to making accurate predictions to inform the downstream
optimization task, we aim to identify a low-dimensional space of features that
is informative for the downstream optimization problem.

_Dimensionality Reduction and Optimization._ PCA has been adopted as the
standard approach for linear dimensionality reduction given a sample
covariance matrix [8]. Casting the dimensionality reduction problem as one of
minimizing the error in approximating a matrix subject to structural
(including rank) constraints, researchers have proposed mathematical-
programming-based approaches to different variants of the problem [9, 10, 11,
12]. Although the studies mentioned above take a mathematical programming
approach as in our paper, they consider the classic objective of minimizing
reconstruction error. As we argue later, when the low-dimensional model is fed
into a subsequent stochastic programming problem, this objective of error
minimization is not necessarily aligned with one of identifying good solutions
in the sense of the downstream objective function. To address this issue, [13]
propose the directed PCA approach that estimates the covariance matrix by
balancing between the PCA objective and empirical optimization. While its aim
is similar to ours, both our setting and methodological approach are
significantly different. In particular, [13] consider a downstream single-
stage stochastic convex optimization problem with uncertain objective
coefficients where the decision variables are unconstrained; in contrast, we
consider a downstream two-stage stochastic program with recourse decisions
that depend on the realization of uncertainty. Therefore, the low-dimensional
representation governs the space of recourse decisions in the downstream
problem. Further, the presence of first-stage decisions in our setting
requires an approach that yields a covariance matrix that is without using the
mean (i.e., invariant to the first-stage decisions) of the data. This is one
of the motivations for us to use distributionally-robust optimization rather
than Bayesian optimization as in [13].

_Decision Rules in Two-Stage Optimization under Uncertainty._ Our setting
concerns learning a low-dimensional representation to characterize the
uncertainty associated with a two-stage stochastic program with recourse.
Conceptually, this objective closely links with the literature on decision
rules in multi-stage optimization with recourse. In this literature, instead
of allowing recourse decisions to be optimized to the specific realizations of
uncertainty, they are confined to the space of parametric functions, known as
decision rules, of the realized uncertain parameters. This enables the problem
to be (heuristically) solved by optimizing over the space of parameters of
these decision rule functions. The decision rule approach has been adopted in
both robust optimization and stochastic programming settings [14, 15, 16, 17,
18, 19]. Various studies in the literature have pointed out the importance of
careful parameterization of the primitive uncertainties of the problem in
enhancing the performances of LDRs and their generalizations [20, 21, 22].
Similarly to these works, our analysis suggests that learning an appropriate
initial representation from data can help significantly improve the
performance of linear decision rules in a data-driven stochastic programming
setting.

##  2 Dimensionality Reduction for Stochastic Optimization

We consider a two-stage stochastic program with recourse:

$\displaystyle\min_{\mathbf{x}\in\mathbf{X}}$
$\displaystyle\mathbf{c}^{T}\mathbf{x}+E[h(\tilde{\bm{z}}-\mathbf{D}^{T}\mathbf{x})]$
(1) where $\displaystyle
h(\bm{z})=\min_{\mathbf{y}}\mathbf{b}^{T}\mathbf{y},\mbox{ s.t.
}\mathbf{A}\mathbf{y}\geq\bm{z}.$ (2)

In (1), the first-stage decision variables $\mathbf{x}$ are chosen under
uncertainty, as characterized by the random variable $\tilde{\bm{z}}$ (in
$\mathbb{R}^{n}$); then, once the values for $\tilde{\bm{z}}$ are realized,
the second-stage decisions $\mathbf{y}$ are chosen to optimize the recourse
problem (2). We assume that the problem has complete recourse, i.e., (2) is
feasible for any value of $\bm{z}$. This implies strong duality:

$h(\bm{z})=\max_{\mathbf{w}\geq 0}\mathbf{w}^{T}\bm{z},\mbox{ s.t.
}\mathbf{A}^{T}\mathbf{w}=\mathbf{b}.$ (3)

We consider the case where $n$ is large, i.e., $\tilde{\bm{z}}$ resides in a
high-dimensional space. In a data-driven setting, the distribution of
$\tilde{\bm{z}}$ is unknown; Instead, a set of training data is available. Let
$\boldsymbol{\mu}$ and $\boldsymbol{\Sigma}$ denote estimates of the mean and
covariance matrix from data (e.g., the sample mean and covariance). With high
dimensionality of $\tilde{\bm{z}}$, it is common practice to model it with a
low-dimensional factor-based representation. It is known that evaluating the
expectation of the random objective function for a stochastic program is
$\\#P$-hard [23]. As a computationally-efficient approach, decision rule
models consider the $n$-dimensional uncertain problem parameters to be
linearly dependent on a set of $k$ _primitive uncertainties_ or _factors_ [16,
18], where $k<<n$, and aim to optimize decision rules defined on said factors.
Thus, the effectiveness of decision rule approaches critically depends on
identifying such a low-dimensional factor model that closely represents the
uncertainties pertaining to the original problem.

An intuitive approach would be to apply a standard dimensionality reduction
algorithm on $\tilde{\bm{z}}$ (such as PCA) and then feeding the resulting
model to the stochastic program (1). However, this naïve sequential approach
would not perform well generally, because the dimensionality reduction
algorithm does not take into account the downstream optimization task. For
example, one may apply PCA to identify the rank-$k$ projection that captures
the maximal amount of variance in the data. Intuitively, this corresponds to
finding the $k$ basis directions along which the data exhibits the largest
variation; However, these are not generally the most _relevant_ directions of
variation for the downstream stochastic program (e.g., for defining effective
decision rules). To address this limitation, we propose a _prescriptive_
dimensionality reduction framework that identifies a low-dimensional
projection of the data that minimizes a measure of suboptimality in the
downstream stochastic program.

###  2.1 The Limitation of PCA

To illustrate, we consider the evaluation of the downstream stochastic
program’s objective (or more specifically, the component that depends on the
model of uncertainty, i.e., the recourse objective $h(\cdot)$) based on a
projection of the data onto some lower-dimensional subspace. In particular,
suppose each data point $\bm{z}$ is projected onto a $k$-dimensional subspace
as $\hat{\bm{z}}=\mathbf{V}\mathbf{V}^{T}\bm{z}$ where
$\mathbf{V}\in\mathbb{R}^{n\times k}$ and $rank(\mathbf{V})=k$. Note that
$\mathbf{V}\mathbf{V}^{T}$ is a symmetric $n\times n$ matrix with rank $k$.
For example, in the case of PCA, we have $\mathbf{V}=\mathbf{V}_{[k]}$, the
$n\times k$ matrix whose columns correspond to the eigenvectors associated
with the $k$ largest eigenvectors. The second-stage objective value under the
projected data is

$\displaystyle h(\hat{\bm{z}})=h(\mathbf{V}\mathbf{V}^{T}\bm{z})$
$\displaystyle=$ $\displaystyle\max_{\mathbf{w}\geq
0}\mathbf{w}^{T}\mathbf{V}\mathbf{V}^{T}\bm{z},\mbox{ s.t.
}\mathbf{w}\in\boldsymbol{P},$ (4)

where the polyhedron $\boldsymbol{P}=\\{\mathbf{w}\geq
0|\mathbf{A}^{T}\mathbf{w}=b\\}$. Then, the following suggests that the
second-stage objective value evaluated under the projected data is equivalent
to the optimal objective value of a counterpart problem defined over the
projected feasible region:
$h(\hat{\bm{z}})=\\{\max\hat{\mathbf{w}}^{T}\bm{z},\mbox{ s.t.
}\hat{\mathbf{w}}\in\hat{\mathbf{P}}\\}$ where
$\hat{\boldsymbol{P}}=\\{\hat{\mathbf{w}}|\hat{\mathbf{w}}=\mathbf{V}\mathbf{V}^{T}\mathbf{w},\mathbf{w}\in\boldsymbol{P}\\}$.

Under PCA, the data is projected onto the $k$-eigenspace of the covariance
matrix. Thus, when the (dual) problem is not _aligned_ with said eigenspace,
the PCA solution could perform badly. In particular, if the projected
polyhedron $\hat{\mathbf{P}}$ is orthogonal to the first $k$ eigenvectors of
$\boldsymbol{\Sigma}$ (i.e., the columns of $\mathbf{V}_{[k]}$), the recourse
objective under the PCA projection will have $h(\hat{\bm{z}})\equiv 0$ for all
$\bm{z}$. That is, the PCA solution may yield a projection that, while
capturing the maximum _amount_ of variation in the data, fails to capture any
_relevant_ variation in terms of optimizing the second-stage problem. This
occurs if the data is projected onto a subspace (the $k$-eigenspace) that is
orthogonal to the _dual_ feasible region of the recourse problem.

###  2.2 Prescriptive PCA

To address the above limitation of PCA, we propose an alternative to PCA,
which we refer to as prescriptive PCA (PPCA) that aligns with the downstream
stochastic program. Formulating the prescriptive PCA problem as a mathematical
program, we will show that a distributionally-robust bound on the expected
reconstruction error can be computed by solving semidefinite programs.

Following the previous discussion, we seek a projection
$\mathbf{V}\mathbf{V}^{T}\tilde{\bm{z}}$ that yields a small expected
reconstruction error (or loss) in terms of the second-stage objective value,
i.e.,

$L(\mathbf{V})=\left|E[h(\tilde{\bm{z}})]-E[h(\mathbf{V}\mathbf{V}^{T}\tilde{\bm{z}})]\right|.$

To this end, we derive an upper bound on $L(\mathbf{V})$ that can be computed
efficiently.

Recall that we seek an approximation independent of the mean of
$\tilde{\bm{z}}$, denoted $\boldsymbol{\mu}$. Let
$\tilde{\bm{z}}_{0}=\tilde{\bm{z}}-\boldsymbol{\mu}$ be the centered random
variable. We seek an upper bound on $L(\mathbf{V})$ that only depends on
$\tilde{\bm{z}}_{0}$, but not $\boldsymbol{\mu}$. Following (4), we have:

######  Proposition 1.

Suppose the linear program (2) has complete recourse. Let $\bm{z}_{0}$ be a
realization of $\tilde{\bm{z}}_{0}$. Then, for any $\bm{z}_{1},\bm{z}_{e}$
such that $\bm{z}_{1}+\bm{z}_{e}=\bm{z}_{0}$, it holds that:

$h(\boldsymbol{\mu}+\bm{z}_{0})\leq
h(\boldsymbol{\mu}+\bm{z}_{1})+h(\bm{z}_{e})\mbox{, and
}h(\boldsymbol{\mu}+\bm{z}_{1})\leq
h(\boldsymbol{\mu}+\bm{z}_{0})+h(-\bm{z}_{e}).$

###### Proof.

We note that

$\displaystyle h(\boldsymbol{\mu}+\bm{z}_{0})$
$\displaystyle=\max_{j\in\\{1,\cdots,J\\}}\left[\mathbf{w}_{j}^{T}(\mathbf{\mu}+\bm{z}_{1})+\mathbf{w}_{j}^{T}\bm{z}_{e}\right]$
$\displaystyle\leq\max_{j\in\\{1,\cdots,J\\}}\mathbf{w}_{j}^{T}(\mathbf{\mu}+\bm{z}_{1})+\max_{j\in\\{1,\cdots,J\\}}\mathbf{w}_{j}^{T}\bm{z}_{e}$
$\displaystyle=h(\boldsymbol{\mu}+\bm{z}_{1})+h(\bm{z}_{e}).$

Because $\bm{z}_{1}=\bm{z}_{0}-\bm{z}_{e}$, it follows similarly that
$h(\boldsymbol{\mu}+\bm{z}_{1})\leq
h(\boldsymbol{\mu}+\bm{z}_{0})+h(-\bm{z}_{e})$ ∎

Proposition 1 implies that the error in evaluating the objective value under
the approximation is bounded by the objective value evaluated under the error
term. More specifically, the expected approximation error is bounded as
follows.

######  Proposition 2.

Consider an approximation $\tilde{\bm{z}}_{0}\approx\tilde{\bm{z}}_{1}$, with
error $\tilde{\bm{z}}_{e}=\tilde{\bm{z}}_{0}-\tilde{\bm{z}}_{1}$ (with
probability one). The absolute error on the evaluated expectation of the
recourse problem is bounded above by:

$\displaystyle\left|E[h(\boldsymbol{\mu}+\tilde{\bm{z}}_{1})]-E[h(\boldsymbol{\mu}+\tilde{\bm{z}}_{0})]\right|\leq\max\left\\{E[h(\tilde{\bm{z}}_{e})],E[h(-\tilde{\bm{z}}_{e})]\right\\}.$
(5)

###### Proof.

Following Proposition 1, it holds that:
$E[h(\boldsymbol{\mu}+\tilde{\bm{z}}_{0})-h(\boldsymbol{\mu}+\tilde{\bm{z}}_{1})]\leq
E[h(\tilde{\bm{z}}_{2})]$ and
$E[h(\boldsymbol{\mu}+\tilde{\bm{z}}_{1})-h(\boldsymbol{\mu}+\tilde{\bm{z}}_{0})]\leq
E[h(-\tilde{\bm{z}}_{2})]$,  which imply (5). ∎

The inequality (5) bounds the error of the approximating
$E[h(\boldsymbol{\mu}+\tilde{\bm{z}}_{0})]$ by
$E[h(\boldsymbol{\mu}+\tilde{\bm{z}}_{1})]$. This bound can serve as a
surrogate of the degree of suboptimality of a chosen approximation
$\tilde{\bm{z}}_{1}$ to $\tilde{\bm{z}}_{0}$ that satisfies desired structural
properties (e.g., has a low-rank covariance matrix). A key feature of this
bound is that it is independent of the mean $\boldsymbol{\mu}$, and thus is
also invariant up to the addition of any constant to $\tilde{\bm{z}}$. This is
a critical property, as approximating the stochastic program (1) requires an
approximation for $E[h(\tilde{\bm{z}}-\mathbf{D}^{T}\mathbf{x})]$ that holds
for any first-stage decision $\mathbf{x}$ (where the deterministic term
$\mathbf{D}^{T}\mathbf{x}$ can be absorbed in the mean).

The bound (5) is not necessarily tight in general. In particular, if the
random variable $\boldsymbol{\mu}+\tilde{\bm{z}}_{0}$ is split into two
components, $\boldsymbol{\mu}+\tilde{\bm{z}}_{1}$ and $\tilde{\bm{z}}_{e}$
with comparable weights, it is conceivable that the bound is loose. On the
other hand, the bound becomes tighter if the approximation
$\boldsymbol{\mu}+\tilde{\bm{z}}_{1}$ carries dominant weight compared with
the residual term $\tilde{\bm{z}}_{e}$, which tends to be the case as the
right-hand side of (5) are to be minimized. Next, we show that a tight bound
on this surrogate can be computed efficiently, and propose an algorithm to
minimize it.

###  2.3 Distributionally-Robust Bound

Proposition 2 suggests a surrogate for the reconstruction error in
approximating the expected second-stage objective value under a given
projection $\tilde{\bm{z}}_{1}=\mathbf{V}\mathbf{V}^{T}\tilde{\bm{z}}_{0}$, by
evaluating the expected value of the second-stage objective under the residual
$\tilde{\bm{z}}_{0}$. Ideally, evaluating the bound in Proposition 2 requires
knowledge of the distribution of $\tilde{\bm{z}}_{0}$. However, in practice,
it is desirable to evaluate this bound without making distributional
assumptions. Invoking results from the distributionally-robust optimization
literature, we show that a bound on this expected value can be computed by
solving a semidefinite program. This bound enables us to formulate a
parsimonious PPCA procedure that uses only the covariance matrix as a
sufficient statistic. We use the following as a direct result of Theorem 2 in
[24].

######  Proposition 3.

Suppose the mean and covariance matrix of $\tilde{\bm{z}}$ are given by
$\boldsymbol{\mu}$ and $\boldsymbol{\Sigma}$, respectively. Then, under any
distribution of $\tilde{\bm{z}}$,
$E[h(\tilde{\bm{z}})]\leq\bar{h}(\boldsymbol{\mu},\boldsymbol{\Sigma})$,
where:

$\begin{array}[]{r@{}l}\bar{h}(\boldsymbol{\mu},\boldsymbol{\Sigma})\equiv\max_{\mathbf{p,Y,X}}&tr(\mathbf{Y})\\\
\mbox{s.t.}&\mathbf{A}^{T}\mathbf{p}=\mathbf{b}\\\
&diag(\mathbf{A}^{T}\mathbf{X}\mathbf{A})=\mathbf{b}^{2}\\\
&\begin{pmatrix}1&\boldsymbol{\mu}^{T}&\mathbf{p}^{T}\\\
\boldsymbol{\mu}&\boldsymbol{\Sigma}&\mathbf{Y}^{T}\\\
\mathbf{p}&\mathbf{Y}&\mathbf{X}\end{pmatrix}\succeq 0\\\ &\mathbf{p,X}\geq
0\end{array}$ (6)

In Proposition 3, problem (6) is a relaxation of the tight upper bound proved
in Theorem 2 of [24]. Using this result, we can bound the expected value on
the right hand side of (5) by $\bar{h}(\mathbf{0},\boldsymbol{\Sigma}_{e})$,
where $\boldsymbol{\Sigma}_{e}$ is the covariance matrix of the residuals
$\tilde{\bm{z}}_{e}$ (which has zero mean by construction). Note that both
$\tilde{\bm{z}}_{e}$ and $-\tilde{\bm{z}}_{e}$ have the same mean (zero) and
covariance matrix $\boldsymbol{\Sigma}_{e}$, thus applying the bound gets rid
of the max operator in (5). This yields a _distribution-free_ performance
bound on any approximation $\tilde{\bm{z}}_{0}\approx\tilde{\bm{z}}_{1}$.

Then, the best (e.g., low-rank) approximation can be obtained by minimizing
the performance bound. Specifically, we look for a rank-$k$ projection
$\tilde{\bm{z}}_{1}=\mathbf{V}\mathbf{V}^{T}\tilde{\bm{z}}_{0}$ that minimizes
the distribution-free performance bound. In the spirit of parsimony, we
operate over the space of covariance matrices. Thus we focus on the covariance
matrix of the projected data, $\boldsymbol{\Sigma}_{1}$. Note that its
eigenvalue decomposition is given by
$\boldsymbol{\Sigma}_{1}=\mathbf{V}\mathbf{E}_{k}\mathbf{V}^{T}$, where
$\mathbf{E}_{k}$ is a $k\times k$ matrix with positive values in the diagonal
entries and zero elsewhere. Thus, the problem of optimizing over the $n\times
k$ matrix $\mathbf{V}$ is equivalent to optimizing over
$\boldsymbol{\Sigma}_{1}$ subject to $rank(\boldsymbol{\Sigma}_{1})\leq k$.
Then, we can formulate the PPCA problem for dimensionality reduction as:

$\displaystyle\min_{\mathbf{\Sigma}_{1},\mathbf{\Sigma}_{e}}\quad$
$\displaystyle\bar{h}(\mathbf{0},\mathbf{\Sigma}_{e})+\theta\langle\mathbf{\Sigma}_{1},\mathbf{\Sigma}_{e}\rangle$
(7) s.t.
$\displaystyle\mathbf{\Sigma}_{1}+\mathbf{\Sigma}_{e}=\mathbf{\Sigma}_{0}$
(10) $\displaystyle\mathbf{\Sigma}_{1},\mathbf{\Sigma}_{e}\succeq 0$
$\displaystyle\mathbf{\Sigma}_{1}\in\mathbb{W},$

where $\mathbb{W}=\left\\{\mathbf{S}\in\mathbb{R}^{n\times
n}:\mbox{rank}(\mathbf{S})\leq k\right\\}$ denotes the set of $n\times n$
symmetric matrices with rank not exceeding $k$ and $\langle\cdot,\cdot\rangle$
denotes the matrix inner product, i.e.,
$\langle\mathbf{A},\mathbf{B}\rangle=tr(\mathbf{A}^{T}\mathbf{B})$. Instead of
minimizing the Frobenius norm of the reconstruction error (for the covariance
matrix) in PCA, PPCA minimizes the distributionally-robust bound on expected
optimality loss regularized with
$\langle\mathbf{\Sigma}_{1},\mathbf{\Sigma}_{e}\rangle$, which is a relaxation
of the requirement on the columns of $\mathbf{\Sigma}_{1}$ and
$\mathbf{\Sigma}_{e}$ be orthogonal (i.e.,
$\langle\mathbf{\Sigma}_{1},\mathbf{\Sigma}_{e}\rangle=0$). The first two
constraints in the above formulation require the second moments of both
$\tilde{\bm{z}}_{1}$ and $\tilde{\bm{z}}_{e}$ to be valid, i.e., there exist
valid multivariate distributions with the corresponding covariance matrices
that sum up to $\mathbf{\Sigma}_{0}$. The constraint (10) is written in a
generic form such that it could enforce any desired structural properties on
the projection, e.g., low rank and sparsity.

The subproblem (6) to compute $\bar{h}(\mathbf{0},\mathbf{\Sigma}_{e})$ is in
the maximization form. Thus, problem (7) is a min-max problem; in fact, it can
be interpreted as a robust optimization problem. In the literature, min-max
robust optimization formulations are typically reformulated as minimization
problems using the duality of the inner problem. Following this standard
approach, we have the following result.

######  Proposition 4.

The PPCA problem can be reformulated as:

$\displaystyle\min\quad$
$\displaystyle\boldsymbol{\alpha}^{T}\mathbf{b}+\boldsymbol{\beta}^{T}\mathbf{b}^{2}+g_{1}+\langle\mathbf{\Lambda}+\theta\mathbf{\Sigma}_{1},\mathbf{\Sigma}_{e}\rangle$
(11) s.t. $\displaystyle\eqref{eq:Sigma}-\eqref{eq:rank},$ (16)
$\displaystyle\mathbf{G}=\begin{pmatrix}g_{1}&\boldsymbol{g}_{2}^{T}&\mathbf{g}_{3}^{T}\\\
\boldsymbol{g}_{2}&\mathbf{G}_{22}&\mathbf{G}_{32}^{T}\\\
\mathbf{g}_{3}&\mathbf{G}_{32}&\mathbf{G}_{33}\end{pmatrix}\succeq 0$
$\displaystyle\begin{array}[]{r@{}ll}&g_{1}=\nu&\quad\mathbf{G}_{22}=\mathbf{\Lambda}\\\
&\mathbf{g}_{2}=\frac{1}{2}\boldsymbol{\gamma}&\quad\mathbf{G}_{32}=-\frac{1}{2}\mathbf{I}\\\
&\mathbf{g}_{3}\leq\frac{1}{2}\mathbf{A}\boldsymbol{\alpha}&\quad\mathbf{G}_{33}=\sum_{i=1}^{m}\beta_{i}\mathbf{a}_{i}\mathbf{a}_{i}^{T},\end{array}$

provided there exists a feasible solution satisfying (16) with strict positive
definiteness.

The problem (11) is not a convex optimization problem and is difficult to
solve in general. In particular, the objective (11) is non-convex due to the
bilinear inner product
$\langle\mathbf{\Lambda}+\theta\mathbf{\Sigma}_{1},\mathbf{\Sigma}_{e}\rangle$.
Furthermore, fixing any feasible value of
$\mathbf{\Lambda}+\theta\mathbf{\Sigma}_{1}$ and solving the remaining problem
yields an upper bound on the original problem. This suggests the alternating
algorithm 1.

Algorithm 1 Alternating Algorithm

Input: problem parameters $\mathbf{A}$, $\mathbf{b}$,
$\boldsymbol{\Sigma}_{0}$, and penalty $\theta$

Initialize $\hat{\boldsymbol{\Sigma}}_{e}$ such that
$\hat{\boldsymbol{\Sigma}}_{e}\succeq 0$ and
$\hat{\boldsymbol{\Sigma}}_{1}=\boldsymbol{\Sigma}_{0}-\hat{\boldsymbol{\Sigma}}_{e}\succeq
0$.

repeat

1\. Given $\mathbf{\Sigma}_{e}=\hat{\boldsymbol{\Sigma}}_{e}$ and
$\mathbf{\Sigma}_{1}=\hat{\boldsymbol{\Sigma}}_{1}$, solve problem (11) over
$(\boldsymbol{\nu},\boldsymbol{\gamma},\alpha,\boldsymbol{\beta},\boldsymbol{\rho},\mathbf{\Lambda},\mathbf{G})$.
Save the optimal value of $\boldsymbol{\Lambda}$ as
$\hat{\boldsymbol{\Lambda}}$.

2\. Given
$\mathbf{\Lambda}+\theta\mathbf{\Sigma}_{1}=\hat{\mathbf{\Lambda}}+\theta\hat{\boldsymbol{\Sigma}}_{1}$,
solve problem (11) over
$(\boldsymbol{\nu},\boldsymbol{\gamma},\alpha,\boldsymbol{\beta},\boldsymbol{\rho},\boldsymbol{\Lambda},\mathbf{G},\mathbf{\Sigma}_{1},\mathbf{\Sigma}_{e})$.
Save the optimal value of $\mathbf{\Sigma}_{e}$ as
$\hat{\boldsymbol{\Sigma}}_{e}$ and $\mathbf{\Sigma}_{1}$ as
$\hat{\boldsymbol{\Sigma}}_{1}$.

until no improvement in the objective value

Output: $\mathbf{\Sigma}_{1}$

It is easy to see that, the objective value (weakly) improves in every
iteration of the algorithm. Therefore, though the algorithm may not
necessarily converge to the global optimal solution, it will converge to a
local minimum. We also remark that it is possible to initialize the algorithm
with the PCA solution. Thus the algorithm, even if terminated early, can
guarantee to produce a solution better than PCA in terms of the worst-case
expected performance bound.

We further remark that it is straightforward to include any convex
regularization term in problem (11). In particular, penalizing the Frobenius
norm of the reconstruction error $\boldsymbol{\Sigma}_{e}$ works well in our
computational studies. For example, (11) can be replaced with:

$\min\eta\left(\boldsymbol{\alpha}^{T}\mathbf{b}+\boldsymbol{\beta}^{T}\mathbf{b}^{2}+g_{1}+\langle\mathbf{\Lambda}+\theta\mathbf{\Sigma}_{1},\mathbf{\Sigma}_{e}\rangle\right)+(1-\eta)\|\boldsymbol{\Sigma}_{e}\|_{F}+\rho\|\Sigma_{1}\|_{*},$

for some $\eta\in[0,1]$. The cases with $\eta=0$ and $\eta=1$ reduce to the
conventional PCA and the unregularized PPCA (11), respectively. Our
computational experiments suggest that a lightly regularized objective (e.g.,
$\eta=0.95$) tends to work better than the unregularized version.

###  2.4 Tight Distributionally Robust Bound

Algorithm 1 yields a low-rank covariance matrix and computes an upper bound on
the worst-case error in approximating the expected second-stage objective
value. In general, this bound is not necessarily tight for three reasons.
First, the performance bound (5) makes use of a subadditivity relation and is
not tight in general. Second, the distributionally-robust formulation (6) is
not exactly tight, in that there is no guarantee that a feasible distribution
(or a sequence thereof) exactly (or asymptotically) achieves this bound.
Third, the alternative algorithm converges to a local minimum to problem (11)
which is not necessarily globally optimal.

Our computational experiments show that the algorithm is effective in
identifying low-rank data projections that perform well in downstream
optimization problems in practice. Yet, it is of theoretical interest to close
these three gaps. Below, we provide a formulation that computes a tight upper
bound on the worst-case expected approximation error given a projection vector
$\mathbf{V}$.

Let the projected polyhedron be
$\hat{\boldsymbol{P}}=\\{\hat{\mathbf{w}}|\hat{\mathbf{w}}=\mathbf{V}\mathbf{V}^{T}\mathbf{w},\mathbf{w}\in\boldsymbol{P}\\}$.
For a data point $\tilde{\bm{z}}$ and a projection
$\mathbf{V}\mathbf{V}^{T}\tilde{\bm{z}}$, the approximation error is:

$\displaystyle
H(\tilde{\bm{z}})=h(\mathbf{V}\mathbf{V}^{T}\tilde{\bm{z}})-h(\tilde{\bm{z}})=\max_{\hat{\mathbf{w}}\in\hat{\mathbf{P}}}\quad$
$\displaystyle\hat{\mathbf{w}}^{T}\tilde{\bm{z}}-\max_{\mathbf{w}\in\mathbf{P}}\mathbf{w}^{T}\tilde{\bm{z}}$
$\displaystyle=\max_{\mathbf{w,y,s}}\quad$
$\displaystyle\mathbf{V}\mathbf{V}^{T}\mathbf{w}\tilde{\bm{z}}-\mathbf{b}^{T}\mathbf{y}$
s.t. $\displaystyle\mathbf{A}\mathbf{y}-\mathbf{s}=\tilde{\bm{z}}$
$\displaystyle\mathbf{A}^{T}\mathbf{w}=\mathbf{b}$
$\displaystyle\mathbf{w,s}\geq 0.$

Given the mean $\boldsymbol{\mu}$ and covariance matrix $\boldsymbol{\Sigma}$
of $\tilde{\bm{z}}$ respectively, the tight distributionally robust bound on
the mean absolute error (MAE) can be evaluated as:

$\displaystyle
Z_{P}=\sup_{\tilde{\bm{z}}\sim(\bm{\mu},\bm{\Sigma})}\left|\mathbb{E}[H(\tilde{\bm{z}})]\right|.$
(17)

######  Theorem 1.

The distributionally robust bound on the MAE in low-rank approximation with
$\mathbf{V}$ is given by $Z_{P}=\max\\{Z_{C}^{+},Z_{C}^{-}\\}$, where
$Z_{C}^{+}$ and $Z_{C}^{-}$ can be evaluated by the following convex
optimization problems:

$\displaystyle\begin{array}[]{ll}Z_{C}^{+}=\max&\langle\mathbf{VV}^{T},\mathbf{Y}_{w}\rangle-\mathbf{b}^{T}\mathbf{p}_{y}\\\
\qquad\text{ s.t. }&\mathbf{a}_{i}^{T}\mathbf{p}_{w}=b_{i},\quad\forall
i=1,...,m\\\
&\mathbf{a}_{i}^{T}\mathbf{X}_{w}\mathbf{a}_{i}=b_{i}^{2},\quad\forall
i=1,...,m\\\
&\langle\mathbf{A}^{T}\mathbf{A},\mathbf{X}_{y}\rangle=\langle\mathbf{A}^{T},\mathbf{Y}_{y}+\mathbf{Z}_{sy}^{T}\rangle\\\
&\langle\mathbf{A}^{T}\mathbf{A},\mathbf{X}_{y}\rangle=\langle\mathbf{I},\mathbf{X}_{s}+2\mathbf{Y}_{s}+\bm{\Sigma}\rangle\\\
&\left(\begin{array}[]{ccccc}1&\bm{\mu}^{T}&\mathbf{p}_{w}^{T}&\mathbf{p}_{y}^{T}&\mathbf{p}_{s}^{T}\\\
\bm{\mu}&\bm{\Sigma}&\mathbf{Y}_{w}^{T}&\mathbf{Y}_{y}^{T}&\mathbf{Y}_{s}^{T}\\\
\mathbf{p}_{w}&\mathbf{Y}_{w}&\mathbf{X}_{w}&\mathbf{Z}_{yw}^{T}&\mathbf{Z}_{sw}^{T}\\\
\mathbf{p}_{y}&\mathbf{Y}_{y}&\mathbf{Z}_{yw}&\mathbf{X}_{y}&\mathbf{Z}_{sy}^{T}\\\
\mathbf{p}_{s}&\mathbf{Y}_{s}&\mathbf{Z}_{sw}&\mathbf{Z}_{sy}&\mathbf{X}_{s}\\\
\end{array}\right)\in\left\\{\mathbf{M}\left|\exists\begin{array}[]{l}V_{1}\in\mathbb{R}_{+}^{1\times
l}\\\ \mathbf{V}_{2}\in\mathbb{R}^{n\times l}\\\
\mathbf{V}_{3}\in\mathbb{R}_{+}^{n\times l}\\\
\mathbf{V}_{4}\in\mathbb{R}^{m\times l}\\\
\mathbf{V}_{5}\in\mathbb{R}_{+}^{n\times l}\end{array}\right.\right.\text{
s.t. }\left.\mathbf{M}=\left(\begin{array}[]{l}V_{1}\\\ \mathbf{V}_{2}\\\
\mathbf{V}_{3}\\\ \mathbf{V}_{4}\\\
\mathbf{V}_{5}\end{array}\right)\left(\begin{array}[]{l}V_{1}\\\
\mathbf{V}_{2}\\\ \mathbf{V}_{3}\\\ \mathbf{V}_{4}\\\
\mathbf{V}_{5}\end{array}\right)^{T}\right\\},\end{array}$

and

$\displaystyle\begin{array}[]{ll}Z_{C}^{-}=\max&\mathbf{b}^{T}\mathbf{p}_{y}-\langle\mathbf{VV}^{T},\mathbf{Y}_{w}\rangle\\\
\qquad\text{ s.t. }&\mbox{Constraints in }Z_{C}^{+}.\end{array}$

##  3 Computational Experiments with Synthetic Data

We first apply the PPCA approach to a stochastic programming problem with a
set of synthetic, simulated data. This allows for the evaluation of the
effectiveness of our proposed approach under a controlled setting. In the next
section, we further illustrate its application based on real-life data.

###  3.1 Problem Setting

We consider a joint production and inventory allocation problem with
transshipment in a network consisting of a set of demand nodes $I$, and a set
of production nodes $I^{\prime}$. Each production node $i\in I^{\prime}$ has a
production capacity of $S_{i}$, and each demand node $j\in I$ faces stochastic
demand $\tilde{z}_{j}$ (where $\tilde{\bm{z}}\in\mathbb{R}^{|I|}$ denotes the
vector of demand). The problem consists of two stages. In the first stage, the
firm determines how much to produce in each site $i\in I^{\prime}$ (with unit
production cost $f_{i}$), as well as the shipment quantity $x_{ij}$ from each
$i\in I^{\prime}$ to each demand location $j\in I$ (with unit shipment cost
$c_{ij}$. Then, in the second stage, demand is realized, and the firm can
fulfill demand with on-hand inventory with possible transshipments: in
particular, it determines the transshipment quantity $y_{ij}$ for demand nodes
$i,j\in I$ (at unit transshipment cost $\bar{c}_{ij}$). Unmet demand at $i\in
I$, denoted by $w_{i}$, will be penalized with a unit shortage cost of
$p_{i}$. This problem can be formulated as a two-stage stochastic program as
below:

$\displaystyle\min_{\mathbf{X}\geq 0}$ $\displaystyle\sum_{i\in
I^{\prime},j\in I}(f_{i}+c_{ij})x_{ij}+\mathbb{E}[h(\mathbf{X,\tilde{z}})]$
(21) s.t. $\displaystyle\sum_{j\in I}x_{ij}\leq S_{i},\forall i\in
I^{\prime},$

where $\mathbf{X}=(x_{ij},\forall i\in I^{\prime},j\in I)$ and $h(\cdot)$
denotes the second stage cost, given by:

$\displaystyle h(\mathbf{X,z})=$ $\displaystyle\min_{\mathbf{Y,w}\geq
0}\sum_{i\in I,j\in I}\bar{c}_{ij}y_{ij}+\sum_{i\in I}w_{i}p_{i}$ s.t.
$\displaystyle-\sum_{j\in I}y_{ij}+\sum_{j\in I}y_{ji}+w_{i}\geq
z_{i}-\sum_{j\in I^{\prime}}x_{ji},\forall i\in I.$

Here, $z_{i}$ is the realized demand at node $i$ and
$\mathbf{p}=(p_{i},\forall i\in I)$ is the vector of penalty costs. Note that
this problem has complete recourse.

###  3.2 Synthetic Data Generation

We generate a synthetic data set based on the 49-node problem instance for
facility location problems from [25], where the 49 demand nodes (set $I$) are
the 48 continental U.S. state capitals and Washington, DC., and the shipping
costs ($c_{ij}$ and $\bar{c}_{ij}$) are proportional to the great circle
distances between any pair of locations $i$ and $j$. We generate stochastic
demand at node $i$ as follows based on the notion of primitive uncertainties
[16]:

$\tilde{z}_{i}=\phi_{i}(\xi_{i1}\tilde{\zeta}_{1}+\xi_{i2}\tilde{\zeta}_{2}+\cdots+\xi_{iK}\tilde{\zeta}_{K})^{+},$

where $(\zeta_{1},\cdots,\zeta_{K})$ are the primitive uncertainties following
some joint distribution (e.g., independent Gaussian distributions), $\xi_{ik}$
denotes fixed coefficients sampled from Uniform($-0.8,1$), and $\psi_{i}>0$ is
a scaling vector proportional to the corresponding demand nodes $i$’
population. Moreover, the operator $\times$ denotes elementwise multiplication
and $(\cdot)^{+}$ denotes $\max(0,\cdot)$. In our experiments, we set $K=25$
and $\mathcal{F}$ to be componentwise independent, which implies that the
covariance matrix of $\tilde{\bm{z}}$ has a rank up to 25. We then sample 100
and 1000 observations as the training and test sets in each experiment
instance.

Furthermore, to evaluate the potential impact of data perturbation or
contamination common in real-life applications, we run a set of experiments
where the training data includes a random noise $\varepsilon_{i}$ following a
Gaussian distribution. In this case, the training data is generated as
follows:

$\tilde{z}_{i}=\phi_{i}(\xi_{i1}\tilde{\zeta}_{1}+\xi_{i2}\tilde{\zeta}_{2}+\cdots+\xi_{iK}\tilde{\zeta}_{K}+\varepsilon_{i})^{+}.$

We also randomly select five of the 49 nodes as the production sites (set
$I^{\prime}$), each with production capacity $S_{i}$ set to $40\%$ of the sum
of the mean demand over the 49 nodes. The production cost $f_{i}$ at site
$i\in I^{\prime}$ is sampled from the Uniform($10,20$) distribution. The
shipping cost in the first stage $c_{ij}=0.015\times$ the distance from node
$i$ to $j$; and the transshipment cost in the second stage is set to
$\bar{c}_{ij}=0.02\times$ the distance from node $i$ to $j$. Finally, the
penalty cost per unit of lost sales is $p_{i}=100$ for all $i\in I$.

###  3.3 Performance Evaluation

We first obtain a low-dimensional representation by solving the (regularized)
PPCA problem with Algorithm 1. This yields a low-rank covariance matrix
$\boldsymbol{\Sigma}_{1}$ that approximates $\boldsymbol{\Sigma}_{0}$,
estimated from the training sample after centering the data. By re-solving the
problem with varying weights $\rho$ on the nuclear norm regularization term,
we obtain $\boldsymbol{\Sigma}_{1}$ with different ranks (values of $k$). For
each $\Sigma_{1}$, we can recover the associated low-dimensional data
projection $\mathbf{V}$ via eigenvalue decomposition.

Given $\mathbf{V}$, the $\mathbf{V}\mathbf{V}^{T}\tilde{\bm{z}}$ gives the
projection of the $n$-dimensional random vector $\tilde{\bm{z}}$ onto a
$k$-dimensional subspace of $\mathbb{R}_{n}$. We test the performance of
approximating the stochastic program based on the alternative low-dimensional
projections identified with PPCA (our proposed approach) and PCA (as a
benchmark) with the same $k$. In particular, we obtain first-stage production
and shipping decisions ($\mathbf{X}$) by solving the LDR-based approximation
for stochastic programs discussed in [16]. Importantly, projecting the demand
vector onto a $k$-dimensional subspace implies modeling the $n$-dimensional
demand vector with $k$ features, i.e., the (prescriptive) principal
components. Following the LDR approach in [16], we restrict the recourse
decisions (transshipments, $\mathbf{y}$) to be affine functions of the
principal components (i.e., with $k+1$ degrees of freedom), instead of the
original demand ($n+1$ degrees of freedom). Thus, dimensionality reduction
effectively reduces the complexity of the LDR formulation. For example, as
illustrated in Figure 1, the computational times can be reduced by as much as
98.7% when the LDR is defined based on two-dimensional primitive uncertainties
(identified with PPCA) than the original 49-dimensional demand vector.

Figure 1: Computational time for solving stochastic program using LDR against
the number of dimensions ($k$) when $\zeta_{i}\sim$Normal(2,1).

To evaluate the performances of the low-rank projections, we evaluate the
first-stage decisions obtained by solving the resulting LDR formulations [16],
via sample average approximation (SAA) approach with the test data. For each
solution, we compute the optimality gap, i.e., the relative difference from
the “true” optimal cost assuming knowledge of the test data (evaluated via
SAA). We first report the experiment where the primitive uncertainties
$\tilde{\zeta}_{k}$’s are independent and identically distributed following
univariate normal distributions, and the training data is sampled without
noise. Figure 2 shows the optimality gaps using PPCA and PCA at varying values
of $k$. We find that PPCA is very effective in identifying low-dimensional
representations for the stochastic program. In particular, it allows for
projecting the 49-dimensional demand data onto subspaces of below 5 dimensions
($k<5$), while maintaining very small (2-3%) optimality gaps in all cases. For
regular PCA, projecting demand data onto $k=5$ or below can lead to poor
performance (e.g., optimality gap exceeding 20%). The performance of PCA is
only able to match PPCA when $k$ is sufficiently large ($k\geq 10$), i.e.,
leading to less efficient representations of the data.

Figure 2: $\zeta_{i}\sim$Normal(2,1) without noise.

A prime motivation for dimensionality reduction is to filter out noises in
training data. In practice, the training samples for many stochastic
optimization problems may be subject to noises, e.g., from data collection or
demand censoring. To test the effectiveness of PPCA under noisy data, we
consider the noise terms $\varepsilon_{i}$’s to be drawn independently from a
univariate Normal distribution with zero mean and a standard deviation that is
proportional to the mean demand. We repeat the computational experiment given
this set of noisy training data and compare the performances of PPCA and
regular PCA. Compared with Figure 2, in Figure 3, we see that the performance
of prescriptive PCA is robust with respect to noise; however, PCA performs
very poorly under noisy training data, even at high values of $k$.
Furthermore, the performance of low-dimensional representation deteriorates
more significantly when the variances of random demands are smaller. It is
because, as the true variances of demands decrease, a larger proportion of the
variation in the training data comes from noise. Importantly, we find that the
performance of PCA deteriorates substantially more than PPCA, showcasing the
robustness of PPCA under noisy data. Unlike the case of noise-free data
(Figure 2), PCA no longer achieves similarly low optimality gaps than does
PPCA unless $k$ is very high ($k\geq 20$; recall that the rank of the true
demand distribution is 25).

Figure 3: $\zeta_{i}\sim$Normal(2,1) with noise.

##  4 Case Study: NYC Taxi Pre-Allocation

Having illustrated the effectiveness of PPCA with a synthetic data set, we
further examine its performance with the New York City taxi data set [26]. We
consider a mobility platform, e.g., a ridesharing platform with autonomous
vehicles, serving the 59 taxi zones on Manhattan island and assume that the
platform faces passenger demand as recorded in the taxi trip data. To ensure
high availability of cars in close vicinity of passengers, the platform needs
to reposition idle cars to prospective demand (pick-up) locations ahead of
demand realization especially during the morning or afternoon peak hours when
demand is highest.

Following [27] that considers the vehicle pre-allocation problem for a single
bottleneck demand period, we formulate this problem as a two-stage stochastic
program: In the first stage (before the peak hour), the platform reposition
idle cars at certain (location-dependent) costs; Then in the second stage,
peak-hour trip demand is realized, and the platform has to match realized
passenger demand with the available cars across the city, to minimize travel
distances for cars or waiting times for passengers. Considering a risk-neutral
objective, this problem can be expressed as a special case of problem (21) in
Section 3.1 by re-interpreting the notation as follows. The set $I=I^{\prime}$
is the set of taxi zones, where each zone $i$ is endowed with $S_{i}$ vehicles
at the beginning of the first stage and faces random demand $\tilde{z}_{i}$ to
be realized in the second stage. The first-stage decisions involve
repositioning $x_{ij}$ vehicles between each pair of zones $i,j\in I$, at a
unit cost $c_{ij}$ per vehicle. For each unit of realized demand at zone $j\in
I$ in the second stage, a matching cost (e.g., customer waiting cost) of
$c_{ij}$ is incurred if it is met by a vehicle positioned at zone $i$, and a
penalty cost $w_{j}$ is incurred if it is unmet. Unlike the inventory
transshipment setting in Section 6, there is no production cost (i.e., $f_{i}$
= 0).

We use the trip records of Yellow Taxis in Manhattan from 8:00am – 8:59am
daily between June to August 2020 provided by [26]. For each day, we count the
number of trips in each taxi zone, resulting in 92 observations of a
59-dimensional random demand vector. Moreover, the costs of repositioning cars
(in the first stage) and traveling to meet customer demand (in the second)
between zones are proportional to the average trip fare between those zones
observed in the data. The penalty cost for unmet demand in each zone is
estimated based on the average fare for all trips originating from that zone.
Finally, we assume the initial total supply of cars to be equal to the mean
demand, and randomly distributed in 10 randomly selected zones. We then evenly
split the 92 observations of trip demand into training and test sets.

Figure 4: The optimality gaps using PCA and PPCA for NYC Taxi Pre-allocation.

Following similar procedures as in the previous section, we solve for the
first-stage vehicle repositioning decisions in the low-dimensional subspace
identified by PPCA and PCA respectively based on the training data. We then
evaluate the out-of-sample performance of the solution with test data. The
optimality gaps are shown in Figure 4. First, we observe that both PPCA and
regular PCA achieve lower optimality gaps as $k$ increases, as expected. In
addition, we can see that the optimality gaps tend to be higher than observed
for synthetic data (without noise) across the range of $k$. This is
inevitable, as the assumption that training and test data are identical and
independently distributed only holds approximately in real-life data. Thus,
the case with real-life data is more comparable with the case of synthetic
data with (some) noise in the training data.

Similar to the case with synthetic data, we find that the performance of PPCA
dominates that of regular PCA: the former is able to achieve lower optimality
gaps (better solution quality) with equal or lower values of $k$ (i.e., the
computational burden in solving the stochastic program). For example, the
performances of PPCA with $k<10$ is similar to that of PCA with $k\approx 20$.
This reaffirms the insight that PPCA projects the data along dimensions that
retains more relevant information with respect to the stochastic program than
does regular PCA.

Moreover, we visualize the top two principal components (PCs) identified by
PPCA and regular PCA in Figure 5. The first PCs, as identified by both
methods, are almost identical. The correlation between the two sets of
loadings is over 0.99. This indicates that both methods are consistent in
finding the first dominant factor of variation in the data. However, the
second PCs differ significantly. From PCA, the second PC highlights an axis of
strong demand variation, where zones 42 (Upper Manhattan/ Harlem) and 43
(Central Park) exhibit clear negative correlation and all other zones carry
relatively uniform weights. From PPCA, however, zone 43 not only exhibits a
clear negative correlation with zone 42, but also a cluster of zones in
Midtown. In our stochastic program, such a factor highlights the importance of
moving vehicles between the Central Park and Midtown areas in the recourse
problem. While this pattern does not necessarily carry the largest variation
as opposed to the PCA solution, it contains relevant information for the
optimization problem that the PCA does not capture.

(a) First PC (PPCA).

(b) Second PC (PPCA)

(c) First PC (PCA)

(d) Second PC (PCA)

Figure 5: The top two principal components (PCs) from PPCA and PCA.

##  5 Conclusion

As the integration of data science and machine learning tools with
mathematical programming becomes more prevalent in practice, there is a
stronger need for the development of methods to align these tools. Our paper
contributes to this growing literature by proposing a prescriptive
dimensionality reduction method that uncovers a low-dimensional linear
projection of data preserving optimality in the sense of solving a downstream
stochastic program, as opposed to preserving variation in the data as in
standard approaches.

Our work can be extended in several directions. First, while our analysis has
focused on rank reduction of the covariance matrix, i.e., the case of PCA, a
similar approach can be adopted for different dimensionality reduction
settings such as sparse PCA and factor analysis. This would involve imposing
the applicable structural constraint in the optimization problem (7), and the
main challenge would lie in developing reformulations and/or approximations
amenable to computationally-efficient solution algorithms.

Second, while our focus was on the unsupervised learning context of linear
dimensionality reduction in this paper, it is possible to follow a similar
approach for supervised learning tasks. A direct extension of our approach
could be applied to the case of low-rank linear regression. Specifically,
instead of seeking a low-rank projection of given data points, one would seek
a low-rank regression model that predicts the dependent variable to be used in
a downstream optimization problem given input features. It would be
interesting to explore similar directions with other machine learning methods,
such as ensemble methods (e.g., random forests).

## References

  * [1] Adam N Elmachtoub and Paul Grigas.  Smart “predict, then optimize”.  Management Science, 68(1):9–26, 2022.
  * [2] Kris Johnson Ferreira, Bin Hong Alex Lee, and David Simchi-Levi.  Analytics for an online retailer: Demand forecasting and price optimization.  Manufacturing & service operations management, 18(1):69–88, 2016\.
  * [3] Chloe Kim Glaeser, Marshall Fisher, and Xuanming Su.  Optimal retail location: Empirical methodology and application to practice: Finalist–2017 M&SOM practice-based research competition.  Manufacturing & Service Operations Management, 21(1):86–102, 2019\.
  * [4] Gah-Yi Ban and Cynthia Rudin.  The big data newsvendor: Practical insights from machine learning.  Operations Research, 67(1):90–108, 2019.
  * [5] Gah-Yi Ban, Jérémie Gallien, and Adam J Mersereau.  Dynamic procurement of new products with covariate information: The residual tree method.  Manufacturing & Service Operations Management, 21(4):798–815, 2019\.
  * [6] Dimitris Bertsimas and Nathan Kallus.  From predictive to prescriptive analytics.  Management Science, 66(3):1025–1044, 2020.
  * [7] Othman El Balghiti, Adam N Elmachtoub, Paul Grigas, and Ambuj Tewari.  Generalization bounds in the predict-then-optimize framework.  Advances in neural information processing systems, 32, 2019.
  * [8] Carl Eckart and Gale Young.  The approximation of one matrix by another of lower rank.  Psychometrika, 1(3):211–218, 1936.
  * [9] Alexandre d’Aspremont, Laurent El Ghaoui, Michael Jordan, and Gert Lanckriet.  A direct formulation for sparse pca using semidefinite programming.  Advances in neural information processing systems, 17, 2004.
  * [10] Alexandre d’Aspremont, Francis Bach, and Laurent El Ghaoui.  Optimal solutions for sparse principal component analysis.  Journal of Machine Learning Research, 9(7), 2008.
  * [11] Dimitris Bertsimas, Ryan Cory-Wright, and Jean Pauphilet.  Solving large-scale sparse pca to certifiable (near) optimality.  Journal of Machine Learning Research, 23(13):1–35, 2022.
  * [12] Dimitris Bertsimas, Martin S Copenhaver, and Rahul Mazumder.  Certifiably optimal low rank factor analysis.  The Journal of Machine Learning Research, 18(1):907–959, 2017.
  * [13] Yi-Hao Kao and Benjamin Van Roy.  Learning a factor model via regularized pca.  Machine learning, 91(3):279–303, 2013.
  * [14] Aharon Ben-Tal, Alexander Goryashko, Elana Guslitzer, and Arkadi Nemirovski.  Adjustable robust solutions of uncertain linear programs.  Mathematical programming, 99(2):351–376, 2004.
  * [15] Dimitris Bertsimas, Dan A Iancu, and Pablo A Parrilo.  Optimality of affine policies in multistage robust optimization.  Mathematics of Operations Research, 35(2):363–394, 2010.
  * [16] Xin Chen, Melvyn Sim, Peng Sun, and Jiawei Zhang.  A linear decision-based approximation approach to stochastic programming.  Operations Research, 56(2):344–357, 2008.
  * [17] Chuen-Teck See and Melvyn Sim.  Robust approximation to multiperiod inventory management.  Operations research, 58(3):583–594, 2010.
  * [18] Joel Goh and Melvyn Sim.  Distributionally robust optimization and its tractable approximations.  Operations research, 58(4-part-1):902–917, 2010.
  * [19] Daniel Kuhn, Wolfram Wiesemann, and Angelos Georghiou.  Primal and dual linear decision rules in stochastic and robust optimization.  Mathematical Programming, 130(1):177–209, 2011.
  * [20] Xin Chen and Yuhan Zhang.  Uncertain linear programs: Extended affinely adjustable robust counterparts.  Operations Research, 57(6):1469–1482, 2009.
  * [21] Jianzhe Zhen, Dick Den Hertog, and Melvyn Sim.  Adjustable robust optimization via fourier–motzkin elimination.  Operations Research, 66(4):1086–1100, 2018.
  * [22] Dimitris Bertsimas, Melvyn Sim, and Meilin Zhang.  Adaptive distributionally robust optimization.  Management Science, 65(2):604–618, 2019.
  * [23] Grani A Hanasusanto, Daniel Kuhn, and Wolfram Wiesemann.  A comment on “computational complexity of stochastic programming problems”.  Mathematical Programming, 159(1):557–569, 2016.
  * [24] Karthik Natarajan and Chung-Piaw Teo.  On reduced semidefinite programs for second order moment bounds with applications.  Mathematical Programming, 161(1-2):487–518, 2017.
  * [25] Mark S Daskin.  Network and discrete location, john willey & sons.  Inc., New York, 1995.
  * [26] NYC Taxi and Limousine Commission.  Tlc trip record data.
  * [27] Zhaowei Hao, Long He, Zhenyu Hu, and Jun Jiang.  Robust vehicle pre-allocation with uncertain covariates.  Production and Operations Management, 29(4):955–972, 2020.
  * [28] Karthik Natarajan, Chung Piaw Teo, and Zhichao Zheng.  Mixed 0-1 linear programs under objective uncertainty: A completely positive representation.  Operations research, 59(3):713–728, 2011.
  * [29] Qingxia Kong, Shan Li, Nan Liu, Chung-Piaw Teo, and Zhenzhen Yan.  Appointment scheduling under time-dependent patient no-show behavior.  Management Science, 66(8):3480–3500, 2020.

##  Appendix A Proofs

###  A.1 Proof of Proposition 4

We can rewrite the inner maximization problem as follows (with dual variables
associated with each constraint in parentheses):

$\begin{array}[]{r@{}l}\bar{h}(\mathbf{0},\mathbf{\Sigma})=\max_{u,\mathbf{v},\mathbf{p},\mathbf{Y},\mathbf{X}}&tr(\mathbf{Y})\\\
\text{s.t.}&\mathbf{a}_{i}^{T}\mathbf{p}=b_{i},\quad\mbox{ for
}i=1,\cdots,m\quad(\alpha_{i})\\\
&\mathbf{a}_{i}^{T}\mathbf{X}\mathbf{a}_{i}=b_{i}^{2},\quad\mbox{ for
}i=1,\cdots,m\quad(\beta_{i})\\\ &u=1\quad(\nu)\\\ &v_{i}=0,\quad\mbox{ for
}i=1,\cdots,n\quad(\gamma_{i})\\\ &S_{ij}=\Sigma_{ij},\quad\mbox{ for
}i,j=1,\cdots,n\quad(\lambda_{ij})\\\ &p_{i}\geq 0,\quad\mbox{ for
}i=1,\cdots,n\quad(\rho_{i})\\\
&\begin{pmatrix}u&\boldsymbol{v}^{T}&\mathbf{p}^{T}\\\
\boldsymbol{v}&\mathbf{S}&\mathbf{Y}^{T}\\\
\mathbf{p}&\mathbf{Y}&\mathbf{X}\end{pmatrix}\succeq
0\quad(\mathbf{G}).\end{array}$

The strict feasibility condition implies strong duality holds by Slater’s
theorem. With the covariance matrix $\mathbf{\Sigma}_{e}$, the Lagrangian
$L(\boldsymbol{\alpha},\boldsymbol{\beta},\nu,\boldsymbol{\gamma},\boldsymbol{\Lambda},\boldsymbol{\rho},\mathbf{G},u,\mathbf{v},\mathbf{p},\mathbf{Y},\mathbf{X})$
is given by

$\displaystyle\langle\mathbf{I},\mathbf{Y}\rangle+\sum_{i=1}^{m}\alpha_{i}(b_{i}-\mathbf{a}_{i}^{T}\mathbf{p})+\sum_{i=1}^{m}\beta_{i}(b_{i}^{2}-\mathbf{a}_{i}^{T}\mathbf{X}\mathbf{a}_{i})$
$\displaystyle\qquad+\nu(1-u)+\sum_{i=1}^{n}\gamma_{i}(-v_{i})+\langle\mathbf{\Lambda},\mathbf{\Sigma}_{e}-\mathbf{S}\rangle+\sum_{i=1}^{n}\rho_{i}p_{i}$
$\displaystyle\qquad+\left\langle\begin{pmatrix}g_{1}&\boldsymbol{g}_{2}^{T}&\mathbf{g}_{3}^{T}\\\
\boldsymbol{g}_{2}&\mathbf{G}_{22}&\mathbf{G}_{32}^{T}\\\
\mathbf{g}_{3}&\mathbf{G}_{32}&\mathbf{G}_{33}\end{pmatrix},\begin{pmatrix}u&\boldsymbol{v}^{T}&\mathbf{p}^{T}\\\
\boldsymbol{v}&\mathbf{S}&\mathbf{Y}^{T}\\\
\mathbf{p}&\mathbf{Y}&\mathbf{X}\end{pmatrix}\right\rangle.$

To minimize the Lagrangian, note that it is always possible to find
$\mathbf{G}$ that is positive semidefinite such that the inner product goes to
negative infinity, unless the primal decision matrix is itself positive
semidefinite. The Lagrangian dual can then be written as

$\min_{\boldsymbol{\alpha},\boldsymbol{\beta},\nu,\boldsymbol{\gamma},\mathbf{\Lambda}}\max_{u,\mathbf{v},\mathbf{p},\mathbf{Y},\mathbf{X}}L(\cdot),$

subject to $\mathbf{G}\succeq 0$. Then, the inner maximization is bounded only
if constraints (16) are satisfied. Thus, the dual can be rewritten as:

$\displaystyle\min_{\alpha,\beta,\nu,\gamma,\Lambda,\rho,\mathbf{G}}$
$\displaystyle\boldsymbol{\alpha}^{T}\mathbf{b}+\boldsymbol{\beta}^{T}\mathbf{b}^{2}+g_{1}+\langle\mathbf{\Lambda},\mathbf{\Sigma}_{e}\rangle,\mbox{s.t.}\mathbf{G}\succeq
0,\mbox{ and }\eqref{eq:dual-constraints}.$

Embedding this into the outer minimization problem yields (11).

###  A.2 Proof of Theorem 1

We first rewrite $Z_{P}$ as below:

$\displaystyle Z_{P}$
$\displaystyle=\sup_{\tilde{\bm{z}}\sim(\bm{\mu},\bm{\Sigma})}\max\left\\{\mathbb{E}[H(\tilde{\bm{z}})],\mathbb{E}[-H(\tilde{\bm{z}})]\right\\}$
$\displaystyle=\max\left\\{\sup_{\tilde{\bm{z}}\sim(\bm{\mu},\bm{\Sigma})}\mathbb{E}[H(\tilde{\bm{z}})],\sup_{\tilde{\bm{z}}\sim(\bm{\mu},\bm{\Sigma})}\mathbb{E}[-H(\tilde{\bm{z}})]\right\\}.$

Define
$Z_{P}^{+}=\sup_{\tilde{\bm{z}}\sim(\bm{\mu},\bm{\Sigma})}\mathbb{E}[H(\tilde{\bm{z}})]$
and
$Z_{P}^{-}=\sup_{\tilde{\bm{z}}\sim(\bm{\mu},\bm{\Sigma})}\mathbb{E}[-H(\tilde{\bm{z}})]$.
Then, $Z_{P}=\max\\{Z_{P}^{+},Z_{P}^{-}\\}$. In the following, we prove
$Z_{P}^{+}=Z_{C}^{+}$ following similar results in [28] and [29]; the proof
for $Z_{P}^{-}=Z_{C}^{-}$ will be similar.

######  Lemma 1.

$Z_{C}^{+}\geq Z_{P}^{+}$.

###### Proof.

Any feasible solution
$(w(\tilde{\bm{z}}),y(\tilde{\bm{z}}),s(\tilde{\bm{z}}))$ to $Z_{P}^{+}$
satisfies

$\displaystyle\mathbf{A}\mathbf{y}(\tilde{\bm{z}})-\mathbf{s}(\tilde{\bm{z}})=\tilde{\bm{z}},\mathbf{A}^{T}\mathbf{w}(\tilde{\bm{z}})=\mathbf{b}.$

The first constraint, i.e.,
$\mathbf{a}_{i}^{T}\mathbf{y}_{i}(\tilde{\bm{z}})=\tilde{z}_{i}+s_{i}(\tilde{\bm{z}})\forall
i$, leads to

$\mathbb{E}[(\mathbf{a}_{i}^{T}\mathbf{y}_{i}(\tilde{\bm{z}}))^{2}]=\mathbb{E}[(\tilde{z}_{i}+s_{i}(\tilde{\bm{z}}))^{2}],\mbox{
and }$
$\mathbb{E}[(\mathbf{a}_{i}^{T}\mathbf{y}_{i}(\tilde{\bm{z}}))^{2}]=\mathbb{E}[(\mathbf{a}_{i}^{T}\mathbf{y}_{i}(\tilde{\bm{z}}))(\tilde{z}_{i}+s_{i}(\tilde{\bm{z}}))],\forall
i.$

Rearranging the terms, we have

$\displaystyle\langle\mathbf{A}^{T}\mathbf{A},\mathbb{E}[\mathbf{y}(\tilde{\bm{z}})\mathbf{y}(\tilde{\bm{z}})^{T}]\rangle=\langle\mathbf{A}^{T},\mathbb{E}[\mathbf{y}(\tilde{\bm{z}})(\mathbf{s}(\tilde{\bm{z}})+\tilde{\bm{z}})^{T}]\rangle$
$\displaystyle\langle\mathbf{A}^{T}\mathbf{A},\mathbb{E}[\mathbf{y}(\tilde{\bm{z}})\mathbf{y}(\tilde{\bm{z}})^{T}]\rangle=\langle\mathbf{I},\mathbb{E}[(\mathbf{s}(\tilde{\bm{z}})+\tilde{\bm{z}})(\mathbf{s}(\tilde{\bm{z}})+\tilde{\bm{z}})^{T}]\rangle.$

From the second constraint, we have

$\mathbb{E}[\mathbf{a}_{i}^{T}\mathbf{w}(\tilde{\bm{z}})]=b_{i},\mbox{ and
}\mathbb{E}[(\mathbf{a}_{i}^{T}\mathbf{w}(\tilde{\bm{z}}))^{2}]=\mathbf{a}_{i}^{T}\mathbb{E}[\mathbf{w}(\tilde{\bm{z}})\mathbf{w}(\tilde{\bm{z}})^{T}]\mathbf{a}_{i}=b_{i}^{2},\forall
i.$

Let $\mathbf{p}_{w}=\mathbb{E}[\mathbf{w}(\tilde{\bm{z}})]$,
$\mathbf{p}_{y}=\mathbb{E}[\mathbf{y}(\tilde{\bm{z}})]$,
$\mathbf{p}_{s}=\mathbb{E}[\mathbf{s}(\tilde{\bm{z}})]$,
$\mathbf{X}_{w}=\mathbb{E}[\mathbf{w}(\tilde{\bm{z}})\mathbf{w}(\tilde{\bm{z}})^{T}]$,
$\mathbf{X}_{y}=\mathbb{E}[\mathbf{y}(\tilde{\bm{z}})\mathbf{y}(\tilde{\bm{z}})^{T}]$,
$\mathbf{X}_{s}=\mathbb{E}[\mathbf{s}(\tilde{\bm{z}})\mathbf{s}(\tilde{\bm{z}})^{T}]$,
$\mathbf{Y}_{w}=\mathbb{E}[\mathbf{w}(\tilde{\bm{z}})\tilde{\bm{z}}^{T}]$,
$\mathbf{Y}_{y}=\mathbb{E}[\mathbf{y}(\tilde{\bm{z}})\tilde{\bm{z}}^{T}]$,
$\mathbf{Y}_{s}=\mathbb{E}[\mathbf{s}(\tilde{\bm{z}})\tilde{\bm{z}}^{T}]$,
$\mathbf{Z}_{yw}=\mathbb{E}[\mathbf{w}(\tilde{\bm{z}})\mathbf{w}(\tilde{\bm{z}})^{T}]$,
$\mathbf{Z}_{sw}=\mathbb{E}[\mathbf{s}(\tilde{\bm{z}})\mathbf{w}(\tilde{\bm{z}})^{T}]$,
$\mathbf{Z}_{sy}=\mathbb{E}[\mathbf{s}(\tilde{\bm{z}})\mathbf{y}(\tilde{\bm{z}})^{T}]$.
They satisfy

$\displaystyle\left(\begin{array}[]{ccccc}1&\bm{\mu}^{T}&\mathbf{p}_{w}^{T}&\mathbf{p}_{y}^{T}&\mathbf{p}_{s}^{T}\\\
\bm{\mu}&\bm{\Sigma}&\mathbf{Y}_{w}^{T}&\mathbf{Y}_{y}^{T}&\mathbf{Y}_{s}^{T}\\\
\mathbf{p}_{w}&\mathbf{Y}_{w}&\mathbf{X}_{w}&\mathbf{Z}_{yw}^{T}&\mathbf{Z}_{sw}^{T}\\\
\mathbf{p}_{y}&\mathbf{Y}_{y}&\mathbf{Z}_{yw}&\mathbf{X}_{y}&\mathbf{Z}_{sy}^{T}\\\
\mathbf{p}_{s}&\mathbf{Y}_{s}&\mathbf{Z}_{sw}&\mathbf{Z}_{sy}&\mathbf{X}_{s}\\\
\end{array}\right)=\mathbb{E}\left[\left(\begin{array}[]{c}1\\\
\tilde{\bm{z}}\\\ \mathbf{w}(\tilde{\bm{z}})\\\ \mathbf{y}(\tilde{\bm{z}})\\\
\mathbf{s}(\tilde{\bm{z}})\end{array}\right)\left(\begin{array}[]{c}1\\\
\tilde{\bm{z}}\\\ \mathbf{w}(\tilde{\bm{z}})\\\ \mathbf{y}(\tilde{\bm{z}})\\\
\mathbf{s}(\tilde{\bm{z}})\end{array}\right)^{T}\right].$

Thus,
$(\mathbf{p}_{w},\mathbf{p}_{y},\mathbf{p}_{s},\mathbf{X}_{w},\mathbf{X}_{y},\mathbf{X}_{s},\mathbf{Y}_{w},\mathbf{Y}_{y},\mathbf{Y}_{s},\mathbf{Z}_{yw},\mathbf{Z}_{sw},\mathbf{Z}_{sy})$
is a feasible solution to $Z_{C}^{+}$ with the objective value
$\langle\mathbf{VV}^{T},\mathbf{Y}_{w}\rangle-\mathbf{b}^{T}\mathbf{p}_{y}=\langle\mathbf{VV}^{T},\mathbb{E}[\mathbf{w}(\tilde{\bm{z}})\tilde{\bm{z}}^{T}]\rangle-\mathbf{b}^{T}\mathbb{E}[\mathbf{y}(\tilde{\bm{z}})]=\mathbb{E}[H(\tilde{\bm{z}})]$.
That is, $Z_{C}^{+}$ is a relaxation of $Z_{P}^{+}$ and we have $Z_{C}^{+}\geq
Z_{P}^{+}$. ∎

We then consider a feasible solution to $Z_{C}^{+}$. By definition, it
satisfies

$\displaystyle\left(\begin{array}[]{ccccc}1&\bm{\mu}^{T}&\mathbf{p}_{w}^{T}&\mathbf{p}_{y}^{T}&\mathbf{p}_{s}^{T}\\\
\bm{\mu}&\bm{\Sigma}&\mathbf{Y}_{w}^{T}&\mathbf{Y}_{y}^{T}&\mathbf{Y}_{s}^{T}\\\
\mathbf{p}_{w}&\mathbf{Y}_{w}&\mathbf{X}_{w}&\mathbf{Z}_{yw}^{T}&\mathbf{Z}_{sw}^{T}\\\
\mathbf{p}_{y}&\mathbf{Y}_{y}&\mathbf{Z}_{yw}&\mathbf{X}_{y}&\mathbf{Z}_{sy}^{T}\\\
\mathbf{p}_{s}&\mathbf{Y}_{s}&\mathbf{Z}_{sw}&\mathbf{Z}_{sy}&\mathbf{X}_{s}\\\
\end{array}\right)=\sum_{l\in L}\left(\begin{array}[]{l}\alpha_{l}\\\
\bm{\beta}_{l}\\\ \bm{\gamma}_{l}^{w}\\\ \bm{\gamma}_{l}^{y}\\\
\bm{\gamma}_{l}^{s}\end{array}\right)\left(\begin{array}[]{l}\alpha_{l}\\\
\bm{\beta}_{l}\\\ \bm{\gamma}_{l}^{w}\\\ \bm{\gamma}_{l}^{y}\\\
\bm{\gamma}_{l}^{s}\end{array}\right)^{T}.$

Define $L_{+}=\\{l\in L|\alpha_{l}>0\\}$ and $L_{0}=\\{l\in
L|\alpha_{l}=0\\}$.

######  Lemma 2.

The above decomposition satisfies
$A\bm{\gamma}_{l}^{y}-\bm{\gamma}_{l}^{s}=\bm{\beta}_{l},\forall l\in L$.

###### Proof.

By the decomposition, we have

$\displaystyle\bm{\Sigma}=\sum_{l\in
L}\bm{\beta}_{l}\bm{\beta}_{l}^{T},\mathbf{X}_{y}=\sum_{l\in
L}\bm{\gamma}_{l}^{y}{\bm{\gamma}_{l}^{y}}^{T},\mathbf{X}_{s}=\sum_{l\in
L}\bm{\gamma}_{l}^{s}{\bm{\gamma}_{l}^{s}}^{T},$
$\displaystyle\mathbf{Y}_{y}=\sum_{l\in
L}\bm{\gamma}_{l}^{y}{\bm{\beta}_{l}}^{T},\mathbf{Y}_{s}=\sum_{l\in
L}\bm{\gamma}_{l}^{s}{\bm{\beta}_{l}}^{T},\mathbf{Z}_{sy}^{T}=\sum_{l\in
L}\bm{\gamma}_{l}^{y}{\bm{\gamma}_{l}^{s}}^{T}.$

The constraint
$\langle\mathbf{A}^{T}\mathbf{A},\mathbf{X}_{y}\rangle=\langle\mathbf{A}^{T},\mathbf{Y}_{y}+\mathbf{Z}_{sy}^{T}\rangle$
in $Z_{C}^{+}$ is thus equivalent to
$\langle\mathbf{A}^{T}\mathbf{A},\sum_{l\in
L}\bm{\gamma}_{l}^{y}{\bm{\gamma}_{l}^{y}}^{T}\rangle=\langle\mathbf{A}^{T},\sum_{l\in
L}\bm{\gamma}_{l}^{y}{\bm{\beta}_{l}}^{T}+\sum_{l\in
L}\bm{\gamma}_{l}^{y}{\bm{\gamma}_{l}^{s}}^{T}\rangle$, which can be further
rearranged into

$\sum_{l\in
L}(\mathbf{A}\bm{\gamma}_{l}^{y})^{T}(\mathbf{A}\bm{\gamma}_{l}^{y})=\sum_{l\in
L}(\bm{\beta}_{l}+\bm{\gamma}_{l}^{s})^{T}(\mathbf{A}\bm{\gamma}_{l}^{y}).$

Moreover, the constraint
$\langle\mathbf{A}^{T}\mathbf{A},\mathbf{X}_{y}\rangle=\langle\mathbf{I},\mathbf{X}_{s}+2\mathbf{Y}_{s}+\bm{\Sigma}\rangle$
is equivalent to $\langle\mathbf{A}^{T}\mathbf{A},\sum_{l\in
L}\bm{\gamma}_{l}^{y}{\bm{\gamma}_{l}^{y}}^{T}\rangle=\langle\mathbf{I},\sum_{l\in
L}\bm{\gamma}_{l}^{s}{\bm{\gamma}_{l}^{s}}^{T}+2\sum_{l\in
L}\bm{\gamma}_{l}^{s}{\bm{\beta}_{l}}^{T}+\sum_{l\in
L}\bm{\beta}_{l}\bm{\beta}_{l}^{T}\rangle$, which can be further rearranged
into

$\sum_{l\in
L}(\mathbf{A}\bm{\gamma}_{l}^{y})^{T}(\mathbf{A}\bm{\gamma}_{l}^{y})=\sum_{l\in
L}(\bm{\beta}_{l}+\bm{\gamma}_{l}^{s})^{T}(\bm{\beta}_{l}+\bm{\gamma}_{l}^{s}).$

Combining both equations above, we have

$\displaystyle\sum_{l\in
L}(\mathbf{A}\bm{\gamma}_{l}^{y}-(\bm{\beta}_{l}+\bm{\gamma}_{l}^{s}))^{T}(\mathbf{A}\bm{\gamma}_{l}^{y}-(\bm{\beta}_{l}+\bm{\gamma}_{l}^{s}))$
$\displaystyle=\sum_{l\in
L}(\mathbf{A}\bm{\gamma}_{l}^{y})^{T}(\mathbf{A}\bm{\gamma}_{l}^{y})-2\sum_{l\in
L}(\bm{\beta}_{l}+\bm{\gamma}_{l}^{s})^{T}(\mathbf{A}\bm{\gamma}_{l}^{y})+\sum_{l\in
L}(\bm{\beta}_{l}+\bm{\gamma}_{l}^{s})^{T}(\bm{\beta}_{l}+\bm{\gamma}_{l}^{s})$
$\displaystyle=0.$

Note that
$(\mathbf{A}\bm{\gamma}_{l}^{y}-(\bm{\beta}_{l}+\bm{\gamma}_{l}^{s}))^{T}(\mathbf{A}\bm{\gamma}_{l}^{y}-(\bm{\beta}_{l}+\bm{\gamma}_{l}^{s}))\geq
0$. We conclude
$(\mathbf{A}\bm{\gamma}_{l}^{y}-(\bm{\beta}_{l}+\bm{\gamma}_{l}^{s}))^{T}(\mathbf{A}\bm{\gamma}_{l}^{y}-(\bm{\beta}_{l}+\bm{\gamma}_{l}^{s}))=0$,
i.e.,
$\mathbf{A}\bm{\gamma}_{l}^{y}-\bm{\gamma}_{l}^{s}=\bm{\beta}_{l},\forall l\in
L$. ∎

######  Lemma 3.

The decomposition satisfies
$\mathbf{A}^{T}\frac{\bm{\gamma}_{l}^{w}}{\alpha_{l}}=\mathbf{b},\forall l\in
L_{+}$ and $\bm{\gamma}_{l}^{w}=\mathbf{0},\forall l\in L_{0}$.

###### Proof.

By the decomposition, we have $\mathbf{p}_{w}=\sum_{l\in
L}\alpha_{l}\gamma_{l}^{w}$ and $\mathbf{X}_{w}=\sum_{l\in
L}\bm{\gamma}_{l}^{w}{\bm{\gamma}_{l}^{w}}^{T}$. The first two constraints
$\mathbf{a}_{i}^{T}\mathbf{p}_{w}=b_{i}$ and
$\mathbf{a}_{i}^{T}\mathbf{X}_{w}\mathbf{a}_{i}=b_{i}^{2}$ in $Z_{C}^{+}$ are
equivalent to

$\mathbf{a}_{i}^{T}\sum_{l\in L}\alpha_{l}\gamma_{l}^{w}=b_{i}\mbox{ and
}\mathbf{a}_{i}^{T}\sum_{l\in
L}\bm{\gamma}_{l}^{w}{\bm{\gamma}_{l}^{w}}^{T}\mathbf{a}_{i}=b_{i}^{2}$

Note that $\sum_{l\in L}\alpha_{l}^{2}=1$. The above equations lead to

$\left(\sum_{l\in L}\alpha_{l}^{2}\right)\left(\sum_{l\in
L}(\mathbf{a}_{i}^{T}\bm{\gamma}_{l}^{w})^{2}\right)=b_{i}^{2}=\left(\sum_{l\in
L}\alpha_{l}\mathbf{a}_{i}^{T}\gamma_{l}^{w}\right)^{2}.$

By the equality conditions of Cauchy-Schwartz inequality, for each $i$,
$\exists\zeta_{i}$ such that
$\zeta_{i}\alpha_{l}=\mathbf{a}_{i}^{T}\bm{\gamma}_{l}^{w},\forall l\in L$.
Recall that $b_{i}=\sum_{l\in
L}\alpha_{l}\mathbf{a}_{i}^{T}\gamma_{l}^{w}=\zeta_{i}\sum_{l\in
L}\alpha_{l}^{2}=\zeta_{i}$. Then, for all $l\in L_{+}$, we have
$\mathbf{a}_{i}^{T}\frac{\bm{\gamma}_{l}^{w}}{\alpha_{l}}=b_{i}$, which leads
to $\mathbf{A}^{T}\frac{\bm{\gamma}_{l}^{w}}{\alpha_{l}}=\mathbf{b},\forall
l\in L_{+}$; for all $l\in L_{0}$, we have
$\mathbf{a}_{i}^{T}\bm{\gamma}_{l}^{w}=b_{i}\alpha_{k}=0$, which leads to
$\mathbf{A}^{T}\bm{\gamma}_{l}^{w}=\mathbf{0},\forall l\in L_{0}$. Suppose the
feasible region $\boldsymbol{P}=\\{\mathbf{w}\geq
0|\mathbf{A}^{T}\mathbf{w}=b\\}$ for the second stage problem is bounded and
nonempty, then we conclude $\bm{\gamma}_{l}^{w}=\mathbf{0},\forall l\in
L_{0}$.∎

By the above lemma, we can decompose the optimal solution to $Z_{C}^{+}$ as

$\displaystyle\left(\begin{array}[]{ccccc}1&\bm{\mu}^{T}&\mathbf{p}_{w}^{*T}&\mathbf{p}_{y}^{*T}&\mathbf{p}_{s}^{*T}\\\
\bm{\mu}&\bm{\Sigma}&\mathbf{Y}_{w}^{*T}&\mathbf{Y}_{y}^{*T}&\mathbf{Y}_{s}^{*T}\\\
\mathbf{p}_{w}^{*}&\mathbf{Y}_{w}^{*}&\mathbf{X}_{w}^{*}&\mathbf{Z}_{yw}^{*T}&\mathbf{Z}_{sw}^{*T}\\\
\mathbf{p}_{y}^{*}&\mathbf{Y}_{y}^{*}&\mathbf{Z}_{yw}^{*}&\mathbf{X}_{y}^{*}&\mathbf{Z}_{sy}^{*T}\\\
\mathbf{p}_{s}^{*}&\mathbf{Y}_{s}^{*}&\mathbf{Z}_{sw}^{*}&\mathbf{Z}_{sy}^{*}&\mathbf{X}_{s}^{*}\\\
\end{array}\right)$ $\displaystyle=\sum_{l\in
L_{+}}\alpha_{l}^{2}\left(\begin{array}[]{l}1\\\
\frac{\bm{\beta}_{l}}{\alpha_{l}}\\\ \frac{\bm{\gamma}_{l}^{w}}{\alpha_{l}}\\\
\frac{\bm{\gamma}_{l}^{y}}{\alpha_{l}}\\\
\frac{\bm{\gamma}_{l}^{s}}{\alpha_{l}}\end{array}\right)\left(\begin{array}[]{l}1\\\
\frac{\bm{\beta}_{l}}{\alpha_{l}}\\\ \frac{\bm{\gamma}_{l}^{w}}{\alpha_{l}}\\\
\frac{\bm{\gamma}_{l}^{y}}{\alpha_{l}}\\\
\frac{\bm{\gamma}_{l}^{s}}{\alpha_{l}}\end{array}\right)^{T}+\sum_{l\in
L_{0}}\left(\begin{array}[]{l}0\\\ \bm{\beta}_{l}\\\ \bm{0}\\\
\bm{\gamma}_{l}^{y}\\\
\bm{\gamma}_{l}^{s}\end{array}\right)\left(\begin{array}[]{l}0\\\
\bm{\beta}_{l}\\\ \bm{0}\\\ \bm{\gamma}_{l}^{y}\\\
\bm{\gamma}_{l}^{s}\end{array}\right)^{T}.$

Let $\epsilon\in(0,1)$. We define a sequence of random variables
$\tilde{\bm{z}}_{\epsilon}$ and their corresponding feasible solutions
$\mathbf{w}^{*}\left(\tilde{\bm{z}}_{\epsilon}\right),\mathbf{y}^{*}\left(\tilde{\bm{z}}_{\epsilon}\right),\mathbf{s}^{*}\left(\tilde{\bm{z}}_{\epsilon}\right)$
(by Lemmas 2 and 3) with the following distribution

$\displaystyle\mathbb{P}\left(\left(\tilde{\bm{z}}_{\epsilon},\mathbf{w}^{*}\left(\tilde{\bm{z}}_{\epsilon}\right),\mathbf{y}^{*}\left(\tilde{\bm{z}}_{\epsilon}\right),\mathbf{s}^{*}\left(\tilde{\bm{z}}_{\epsilon}\right)\right)=\left(\frac{\boldsymbol{\beta}_{l}}{\alpha_{l}},\frac{\boldsymbol{\gamma}_{l}^{w}}{\alpha_{l}},\frac{\boldsymbol{\gamma}_{l}^{y}}{\alpha_{l}},\frac{\boldsymbol{\gamma}_{l}^{s}}{\alpha_{l}}\right)\right)$
$\displaystyle=\left(1-\epsilon^{2}\right)\alpha_{l}^{2},\forall l\in L_{+},$
$\displaystyle\mathbb{P}\left(\left(\tilde{\bm{z}}_{\epsilon},\mathbf{w}^{*}\left(\tilde{\bm{z}}_{\epsilon}\right),\mathbf{y}^{*}\left(\tilde{\bm{z}}_{\epsilon}\right),\mathbf{s}^{*}\left(\tilde{\bm{z}}_{\epsilon}\right)\right)=\left(\frac{\sqrt{|L_{0}|}\boldsymbol{\beta}_{l}}{\epsilon},\mathbf{w}^{\dagger},\frac{\sqrt{|L_{0}|}\boldsymbol{\gamma}_{l}^{y}}{\epsilon},\frac{\sqrt{|L_{0}|}\boldsymbol{\gamma}_{l}^{s}}{\epsilon}\right)\right)$
$\displaystyle=\frac{\epsilon^{2}}{\left|L_{0}\right|},\forall l\in L_{0},$

where $\mathbf{w}^{\dagger}\in\mathbf{P}$ is any feasible solution
$\mathbf{w}$ to $H(\tilde{\bm{z}})$.

We can verify that

$\displaystyle\mathbb{E}\left[\left(\begin{array}[]{c}1\\\
\tilde{\bm{z}}_{\epsilon}\\\
\mathbf{w}^{*}\left(\tilde{\bm{z}}_{\epsilon}\right)\\\
\mathbf{y}^{*}\left(\tilde{\bm{z}}_{\epsilon}\right)\\\
\mathbf{s}^{*}\left(\tilde{\bm{z}}_{\epsilon}\right)\end{array}\right)\left(\begin{array}[]{c}1\\\
\tilde{\bm{z}}_{\epsilon}\\\
\mathbf{w}^{*}\left(\tilde{\bm{z}}_{\epsilon}\right)\\\
\mathbf{y}^{*}\left(\tilde{\bm{z}}_{\epsilon}\right)\\\
\mathbf{s}^{*}\left(\tilde{\bm{z}}_{\epsilon}\right)\end{array}\right)^{T}\right]$
$\displaystyle=\sum_{l\in
L_{+}}(1-\epsilon^{2})\alpha_{l}^{2}\left(\begin{array}[]{c}1\\\
\frac{\bm{\beta}_{l}}{\alpha_{l}}\\\ \frac{\bm{\gamma}_{l}^{w}}{\alpha_{l}}\\\
\frac{\bm{\gamma}_{l}^{y}}{\alpha_{l}}\\\
\frac{\bm{\gamma}_{l}^{s}}{\alpha_{l}}\end{array}\right)\left(\begin{array}[]{c}1\\\
\frac{\bm{\beta}_{l}}{\alpha_{l}}\\\ \frac{\bm{\gamma}_{l}^{w}}{\alpha_{l}}\\\
\frac{\bm{\gamma}_{l}^{y}}{\alpha_{l}}\\\
\frac{\bm{\gamma}_{l}^{s}}{\alpha_{l}}\end{array}\right)^{T}+\sum_{l\in
L_{0}}\frac{\epsilon^{2}}{\left|L_{0}\right|}\left(\begin{array}[]{c}1\\\
\frac{\sqrt{|L_{0}|}\bm{\beta}_{l}}{\epsilon}\\\ \mathbf{w}^{\dagger}\\\
\frac{\sqrt{|L_{0}|}\bm{\gamma}_{l}^{y}}{\epsilon}\\\
\frac{\sqrt{|L_{0}|}\bm{\gamma}_{l}^{s}}{\epsilon}\end{array}\right)\left(\begin{array}[]{c}1\\\
\frac{\sqrt{|L_{0}|}\bm{\beta}_{l}}{\epsilon}\\\ \mathbf{w}^{\dagger}\\\
\frac{\sqrt{|L_{0}|}\bm{\gamma}_{l}^{y}}{\epsilon}\\\
\frac{\sqrt{|L_{0}|}\bm{\gamma}_{l}^{s}}{\epsilon}\end{array}\right)^{T}$
$\displaystyle\overrightarrow{\epsilon\downarrow
0}\left(\begin{array}[]{ccccc}1&\bm{\mu}^{T}&\mathbf{p}_{w}^{*T}&\mathbf{p}_{y}^{*T}&\mathbf{p}_{s}^{*T}\\\
\bm{\mu}&\bm{\Sigma}&\mathbf{Y}_{w}^{*T}&\mathbf{Y}_{y}^{*T}&\mathbf{Y}_{s}^{*T}\\\
\mathbf{p}_{w}^{*}&\mathbf{Y}_{w}^{*}&\mathbf{X}_{w}^{*}&\mathbf{Z}_{yw}^{*T}&\mathbf{Z}_{sw}^{*T}\\\
\mathbf{p}_{y}^{*}&\mathbf{Y}_{y}^{*}&\mathbf{Z}_{yw}^{*}&\mathbf{X}_{y}^{*}&\mathbf{Z}_{sy}^{*T}\\\
\mathbf{p}_{s}^{*}&\mathbf{Y}_{s}^{*}&\mathbf{Z}_{sw}^{*}&\mathbf{Z}_{sy}^{*}&\mathbf{X}_{s}^{*}\\\
\end{array}\right)$

As $\epsilon\downarrow 0$, the random vectors
$\left(\tilde{\bm{z}}_{\epsilon},\mathbf{w}^{*}\left(\tilde{\bm{z}}_{\epsilon}\right),\mathbf{y}^{*}\left(\tilde{\bm{z}}_{\epsilon}\right),\mathbf{s}^{*}\left(\tilde{\bm{z}}_{\epsilon}\right)\right)$
converge almost surely (a.s.) to
$\left(\tilde{\bm{z}}^{*},\mathbf{w}^{*}\left(\tilde{\bm{z}}^{*}\right),\mathbf{y}^{*}\left(\tilde{\bm{z}}^{*}\right),\mathbf{s}^{*}\left(\tilde{\bm{z}}^{*}\right)\right)$
with distribution below:

$\mathbb{P}\left(\left(\tilde{\bm{z}}^{*},\mathbf{w}^{*}\left(\tilde{\bm{z}}^{*}\right),\mathbf{y}^{*}\left(\tilde{\bm{z}}^{*}\right),\mathbf{s}^{*}\left(\tilde{\bm{z}}^{*}\right)\right)=\left(\frac{\boldsymbol{\beta}_{l}}{\alpha_{l}},\frac{\boldsymbol{\gamma}_{l}^{w}}{\alpha_{l}},\frac{\boldsymbol{\gamma}_{l}^{y}}{\alpha_{l}},\frac{\boldsymbol{\gamma}_{l}^{s}}{\alpha_{l}}\right)\right)=\alpha_{l}^{2},\forall
l\in L_{+}.$

We then have

$\displaystyle\lim_{\epsilon\downarrow
0}\mathbb{E}\left[H(\tilde{\bm{z}}_{\epsilon})\right]$
$\displaystyle=\lim_{\epsilon\downarrow 0}\left(\sum_{l\in
L_{+}}(1-\epsilon^{2})\alpha_{l}^{2}\left(\mathbf{V}\mathbf{V}^{T}\frac{\bm{\gamma}_{l}^{w}}{\alpha_{l}}\frac{\bm{\beta}_{l}}{\alpha_{l}}-\mathbf{b}^{T}\frac{\bm{\gamma}_{l}^{y}}{\alpha_{l}}\right)\right.$
$\displaystyle\quad\left.+\sum_{l\in
L_{0}}\frac{\epsilon^{2}}{\left|L_{0}\right|}\left(\mathbf{V}\mathbf{V}^{T}\mathbf{w}^{\dagger}\frac{\sqrt{|L_{0}|}\boldsymbol{\beta}_{l}}{\epsilon}-\mathbf{b}^{T}\frac{\sqrt{|L_{0}|}\boldsymbol{\gamma}_{l}^{y}}{\epsilon}\right)\right)$
$\displaystyle=\sum_{l\in
L_{+}}\left(\mathbf{V}\mathbf{V}^{T}\bm{\gamma}_{l}^{w}\bm{\beta}_{l}-\mathbf{b}^{T}\bm{\gamma}_{l}^{y}\alpha_{l}\right)=\mathbb{E}\left[H(\tilde{\bm{z}}^{*})\right].$

By the similar limiting argument on page 719 of [28], we have
$Z_{P}^{+}=\sup_{\tilde{\bm{z}}\sim(\bm{\mu},\bm{\Sigma})}\mathbb{E}[H(\tilde{\bm{z}})]\geq\mathbb{E}\left[H(\tilde{\bm{z}}^{*})\right]$
and thus

$\displaystyle Z_{P}^{+}$
$\displaystyle=\sup_{\tilde{\bm{z}}\sim(\bm{\mu},\bm{\Sigma})}\mathbb{E}[H(\tilde{\bm{z}})]$
$\displaystyle\geq\mathbb{E}\left[H(\tilde{\bm{z}}^{*})\right]=\mathbf{V}\mathbf{V}^{T}\sum_{l\in
L_{+}}\bm{\gamma}_{l}^{w}\bm{\beta}_{l}-\mathbf{b}^{T}\sum_{l\in
L_{+}}\bm{\gamma}_{l}^{y}\alpha_{l}$
$\displaystyle\quad=\langle\mathbf{VV}^{T},\mathbf{Y}_{w}^{*}\rangle-\mathbf{b}^{T}\mathbf{p}_{y}^{*}=Z_{C}^{+}.$

Combining with Lemma 1, we conclude $Z_{C}^{+}=Z_{P}^{+}$. To prove
$Z_{C}^{-}=Z_{P}^{-}$, we only need to change the sign of the objective and
the proof follows as $Z_{C}^{-}$ and $Z_{C}^{+}$ share the same constraints.