\section{Introduction}\label{sec:intro}

Vector and multiobjective optimization is the field of maximizing or minimizing multiple objective functions simultaneously.  As with traditional (scalar) optimization, these problems can be classified as, e.g., linear, convex, or nonconvex.  Within this work we focus on convex vector optimization problems.  Such problems are widely applicable in practice.  Specifically, multiobjective optimization has been applied to fields as diverse as finance (for the mean-risk problem of~\cite{markowitz1952} in~\cite{KR18} as well as for set-valued risk measures~\cite{FR14-alg}), electric and power systems \cite{cassidy2016risk}, structural design \cite{marler2004survey} and references therein, and chemical engineering \cite{bhaskar2000applications}.

Convex vector optimization problems are those with two or more convex objectives and convex constraints.  Already multiple algorithms exist to approximate the solutions for these problems (or for specific subclasses of such problems).  We refer the interested reader to \cite{RW05}
for an overview of some algorithms for convex vector optimization problems.  For numerical tractability, these methods result in (polyhedral) approximations of the efficient frontier.  That is, rather than finding the exact solution, only an (inner and/or outer) approximation of the efficient frontier can be found.  We wish to highlight algorithms based on Benson's outer approximation algorithm~\cite{benson98} for linear vector optimization problems, that have been generalized to the convex case in~\cite{ESS11, LRU14, DLSW21, KU22}.
While powerful, these algorithms suffer from curse of dimensionality as the computational complexity of such problems grows exponentially in the dimension of the problem as, e.g., the vertex enumeration subproblem is NP-hard~\cite{khachiyan2008generating}.

Within this work we construct a new approximation method for convex vector optimization problems using neural networks.  The primary contributions and innovations of this work are two-fold.
\begin{enumerate}
\item We propose a machine learning methodology to (approximately) solve convex vector optimization problems by considering the weighted-sum scalarizations of the multiobjective problem.  Due to the use of neural networks, the proposed algorithm can tractably be applied to large vector optimization problems.  That is, we are able to overcome the curse of dimensionality in solving vector optimization problems.
\item By proposing two neural networks that are jointly trained, the proposed machine learning methodology provides both inner and outer approximations of the weakly efficient frontier of the relevant convex vector optimization problem.  In constructing the algorithm in this way, we are able to quantify (an upper bound to) the error of these approximations at each point on the weakly efficient frontier.  By having inner and outer approximations with a known error, this machine learning methodology can be applied in practice without the typical risks of machine learning -- the uncertain approximation error.
\end{enumerate}

The organization of this paper is as follows.
In Section~\ref{sec:background}, the background, notation, and structure of both convex vector optimization problems and neural networks is provided.
In Section~\ref{sec:nn}, the proposed neural network architecture for approximating the weakly efficient frontier of a strictly convex vector optimization problem is provided.  Specifically, a primal neural network is proposed in Section~\ref{sec:nn-primal}, a dual neural network is provided in Section~\ref{sec:nn-dual}, and the joint loss function for these neural networks is given in Section~\ref{sec:nn-loss}. Then, in Section~\ref{sec:approx}, inner and outer approximations of the upper image of the problem are constructed from these neural networks.  Within this section, we additionally provide a function which provides an upper bound to the errors for the inner approximation at each point on the weakly efficient frontier.  We implement the proposed machine learning approach for three strictly convex optimization problems in Section~\ref{sec:cs}.  We use these problems both to validate the proposed approach and to demonstrate its applicability to large-scale problems thus overcoming the curse of dimensionality.
Finally, in Section~\ref{sec:convex}, we extend the machine learning framework of Section~\ref{sec:nn} to allow for convex, but not strictly convex, vector optimization problems.  This is demonstrated on two small problems -- a linear vector optimization problem and the mean-risk portfolio optimization problem -- that permit simple visualizations.

\section{Background}\label{sec:background}
Within this section, we provide the background material necessary for this work.  First, in Section~\ref{sec:background-vop}, the vector optimization setting that is considered throughout this work is presented along with basic definitions and results.  Second, in Section~\ref{sec:background-nn}, some background on neural networks and their usage for approximating, i.e.\ learning, a target function, are given.

\subsection{Vector optimization problems}\label{sec:background-vop}

Consider a decision maker attempting to simultaneously minimize $P$ objective functions $f_i: \mathbb{X} \to \mathbb{R}$ for $i=1,...,P$ over the set of decisions $\mathbb{X} \subseteq \mathbb{R}^N$.    Throughout, we assume that the decision space $\mathbb{X} \neq \emptyset$ is nonempty and that the ordering cone in the objective space is $\mathbb{R}^P_+$.  To simplify notation, we denote $f(x) := (f_1(x),...,f_P(x))^\top$ for any $x \in \mathbb{X}$.  That is, throughout this work we consider a vector optimization problem (VOP) of the form
\begin{equation}\label{eq:vop}
\min\{f(x) \; | \; x \in \mathbb{X}\}.
\end{equation}

Before continuing to the main results of this work, we provide some basic definitions and results on vector optimization.  We refer the interested reader to, e.g.,~\cite{jahn2011vector} for more details on these problems and definitions. We denote by $f[\mathbb{X}]=  \{f(x) \; | \; x \in \mathbb{X}\}$ the image of the feasible set.
\begin{definition}\label{defn:bounded}
The VOP \eqref{eq:vop} is \textbf{\emph{bounded}} if there exists some $y \in \mathbb{R}^P$ such that $f[\mathbb{X}] \subseteq \{y\} + \mathbb{R}^P_+$.
\end{definition}

\begin{definition}\label{defn:pareto}
Consider the VOP \eqref{eq:vop}.
$x^* \in \mathbb{X}$ is a \textbf{\emph{(weak) minimizer}} if $f(x^*) - f(x) \in \mathbb{R}^P_+$ for some $x \in \mathbb{X}$ implies $f(x) = f(x^*)$ (resp.\ $f(x^*) - f(x) \not\in \mathbb{R}^P_{++}$ for every $x \in \mathbb{X}$).
The set of (weak) minimizers is denoted by $\mathcal{X}^* \subseteq \mathbb{X}$ (resp.\ $\mathcal{X} \subseteq \mathbb{X}$).
The \textbf{\emph{(weak) efficient frontier}} is given by $f[\mathcal{X}^*] := \{f(x^*) \; | \; x^* \in \mathcal{X}^*\}$ (resp.\ $f[\mathcal{X}]$).
The \textbf{\emph{upper image}} is provided by $\mathcal{P} := \operatorname{cl}\left(f[\mathcal{X}] + \mathbb{R}^P_+\right)$.
\end{definition}

Within this work we are specifically interested in problems with feasible space
\begin{equation}\label{eq:constraints}
\mathbb{X} := \{x \in \mathbb{R}^N \; | \; g(x) \in -\mathbb{R}^M_+\}
\end{equation}
for $g: \mathbb{R}^N \to \mathbb{R}^M$.
In particular, we are interested in \emph{convex} vector optimization problems (CVOP), i.e., problems of the form
\begin{equation}\label{eq:cvop}\tag{P}
\min\{f(x) \; | \; g_j(x) \leq 0,~j = 1,...,M\}
\end{equation}
such that $f_i,g_j$ are convex functions for every $i = 1,...,P$ and $j = 1,...,M$.

It is well known that (weak) minimizers of~\eqref{eq:cvop} are related to solutions of the weighted-sum scalarization problem.
\begin{proposition}\label{prop:scalar}
Consider the convex vector optimization problem~\eqref{eq:cvop}.
\begin{enumerate}
\item $x^*$ is a weak minimizer if and only if there exists some $w \in \mathcal{W} := \{w \in \mathbb{R}^P_+ \; | \; {\bf 1}^\top w = 1\}$, where we denote ${\bf 1} := (1,1,...,1)^\top \in \mathbb{R}^P$, such that $x^*$ solves the weighted-sum scalarization problem
    \begin{equation}\label{eq:scalarization}\tag{wP}
    \min\Big\{\sum_{i = 1}^P w_i f_i(x) \; | \; g_j(x) \leq 0,~j = 1,...,M\Big\}.
    \end{equation}
\item $x^*$ is a minimizer if it solves the weighted-sum scalarization problem~\eqref{eq:scalarization} for some $w \in \mathcal{W} \cap \mathbb{R}^P_{++}$.
\end{enumerate}
\end{proposition}
\begin{proof}
The first statement follows from Corollary 5.29 of~\cite{jahn2011vector}.  The second statement follows from Theorem 5.18(b) of~\cite{jahn2011vector}.
\end{proof}

\begin{assumption}\label{ass:slater}
Throughout this work we will make the following assumptions:
\begin{itemize}
\item $f_i$ are strictly convex and continuously differentiable;
\item $g_j$ are convex and continuously differentiable;
\item \eqref{eq:cvop} is bounded;
\item Slater's condition holds, i.e., there exists some $\bar x \in \mathbb{R}^N$ such that $g_j(\bar x) < 0$ for every $j = 1,...,M$.
\end{itemize}
\end{assumption}

\begin{remark}
\label{rem:convex}
Let us remark on some possible relaxations of these assumptions.
\begin{enumerate}
\item\label{rem:convex-convex}
Note that the strict convexity of the objective functions $f_i$ as encoded in Assumption~\ref{ass:slater} could be relaxed for all results presented within this work so long as the weighted-sum scalarization problem~\eqref{eq:scalarization} is strictly convex for all scalarizations $w \in \mathcal{W}$ of interest. This will be used in the example presented in Section~\ref{sec:mean-var}. It will also prove useful in Section~\ref{sec:convex} to allow for general convex, but not strictly convex, objective functions $f_i$ by augmenting the problem with an additional strictly convex objective function. This augmentation ensures that the new objective has strictly convex scalarizations for all $w \in \mathcal{W}$ of interest.
\item\label{rem:convex-equality}
As presented above, we focus solely on inequality constrained CVOPs within this work.
We wish to note that linear equality constraints can be introduced within this setting as well.
Consider the CVOP
\begin{equation*}
\min\{f(\hat{x}) \; | \; g_j(\hat{x}) \leq 0,~j = 1,...,M,\; A\hat{x} = b\}
\end{equation*}
for $A \in \mathbb{R}^{O \times N}$ and $b \in \mathbb{R}^O$.  Without loss of generality, we will assume $A$ is of full rank.
Let $A_\perp \in \mathbb{R}^{N \times (N-O)}$ provide a basis for the null space of $A$; that is, $A A_\perp z = {\bf 0}$ for any $z \in \mathbb{R}^{N-O}$.
Furthermore, let $\tilde{x} \in \mathbb{R}^N$ define a particular solution of the linear constraints, i.e, $A\tilde{x} = b$.
The equality constrained CVOP can then be reformulated in the form~\eqref{eq:cvop} as
\begin{equation*}
\min\{f(A_{\perp}x + \tilde{x}) \; | \; g_j(A_{\perp}x + \tilde{x}) \leq 0,~j = 1,...,M\}
\end{equation*}
with $N-O$ variables.  Provided this reformulated problem satisfies~\ref{ass:slater}, the methodology considered within this work can be applied directly to this reformulation. This will be used in the examples considered in Sections~\ref{sec:mean-var} and~\ref{sec:convex}.
\ite
Herein we take the ordering cone $\mathbb{R}^P_+$.  All results can be considered in which the minimization is taken with respect to an ordering cone $C \subseteq \mathbb{R}^p$ with nonempty interior instead.  The only modification necessary is that the scalarizations are taken with respect to $w \in C^+ := \{w \in \mathbb{R}^p \; | \; c^\top w \geq 0 \; \forall c \in C\}$ such that $w^\top c = 1$ for some fixed element $c \in \int C$ provided the objective $f$ is $C$-convex.
\end{enumerate}
\end{remark}

The weakly efficient frontier can be traced through considerations of $w\in\mathcal{W}\mapsto f(x^*(w))$ where $x^*(w)$ is a solution to the scalarization problem~\eqref{eq:scalarization} in Proposition~\ref{prop:scalar} for every $w \in \mathcal{W}$.

In addition to the primal scalarization problem~\eqref{eq:scalarization} we will also consider its Lagrange dual problem.  For that, consider the dual function $d: \mathbb{R}^M \times \mathcal{W} \to \mathbb{R}$ defined by
\begin{equation}\label{eq:dual}
d(\lambda,w) := \inf_{x \in \mathbb{R}^N} \Big[\sum_{i = 1}^P w_i f_i(x) + \sum_{j = 1}^M \lambda_j g_j(x)\Big].
\end{equation}
The dual problem is then to maximize this dual function subject to the dual feasibility constraints, i.e.,
\begin{equation}\label{eq:scalarization-dual}\tag{wD}
\max \Big\{\inf_{x \in \mathbb{R}^N} \Big[\sum_{i = 1}^P w_i f_i(x) + \sum_{j = 1}^M \lambda_j g_j(x)\Big] \; | \; \lambda \in \mathbb{R}^M_+\Big\}.
\end{equation}
By Assumption~\ref{ass:slater}, strong duality is satisfied for these scalar problems, i.e.\ the optimal values of~\eqref{eq:scalarization} and~\eqref{eq:scalarization-dual} are equal, and a dual solution exists.

The method presented in this paper to approximate the weakly efficient frontier by machine learning methods will rely heavily on the following well-known KKT conditions. Note that the KKT conditions~\eqref{eq:kkt}-\eqref{eq:kkt-cs} below are the usual KKT conditions of the scalarization problem~\eqref{eq:scalarization}, but can also be obtained using the multiobjective KKT conditions for the CVOP~\eqref{eq:cvop} directly (without the need to go through the scalarization first), see~\cite{eichfelder2021proximity}.

\begin{theorem}\label{thm:kkt}
$(x^*(w),\lambda^*(w)) \in \mathbb{R}^N \times \mathbb{R}^M$ are a pair of primal ($x^*(w)$) and dual ($\lambda^*(w)$) solutions of the scalarization problems~\eqref{eq:scalarization}, respectively~\eqref{eq:scalarization-dual}, w.r.t.\ $w \in \mathcal{W}$ if and only if they jointly satisfy the following KKT conditions
\begin{align}
\label{eq:kkt}  {\bf 0} &= Jf(x^*(w))^\top w + Jg(x^*(w))^\top \lambda^*(w) \\
\label{eq:kkt-primal}  {\bf 0} &\geq g(x^*(w)) \\
\label{eq:kkt-dual}  {\bf 0} &\leq \lambda^*(w) \\
\label{eq:kkt-cs}  0 &= \lambda_i^*(w) g_i(x^*(w))  \quad i=1,...,m,
\end{align}
where $J$ denotes the Jacobian of the corresponding function.
\end{theorem}
\begin{proof}
Due to the assumption that Slater's condition holds (Assumption~\ref{ass:slater}), this result can be found in, e.g.,~\cite[Chapter 5.5.3]{boyd2004convex}.
\end{proof}

\subsection{Neural networks}\label{sec:background-nn}
Consider, now, the task of approximating a function $y^*: \mathcal{T} \to \mathbb{R}^m$ for some \emph{compact} set of input variables $\mathcal{T} \subseteq \mathbb{R}^n$.  Such an approximation problem is, fundamentally, a regression problem.  Within this work, we will focus on feed-forward neural networks to regress the function $y^*$. These functions can be viewed as a multi-stage regression model.  We refer the interested reader to, e.g.,~\cite[Chapter 6]{goodfellow2016deep} for more details on feed-forward neural networks.
\begin{definition}\label{defn:nn}
An $\ell$-hidden layer \textbf{\emph{neural network}} with $h_l$ nodes in layer $l \in \{0,1,...,\ell+1\}$ (such that $h_0 := n$ and $h_{\ell+1} := m$ for the input and output layers respectively) is a mapping $y_{\ell+1}: \mathcal{T} \times \Theta \to \mathbb{R}^m$ that can be decomposed as
\begin{align*}
y_l(t,\theta) &:= \Phi_l(\theta_{l,b} + \theta_{l,w} y_{l-1}(t,\theta)) \quad \forall l \in \{1,...,\ell+1\}
\end{align*}
with $y_0(t,\theta) := t$
for \textbf{\emph{activation functions}} $\Phi_l: \mathbb{R}^{h_l} \to \mathbb{R}^{h_l}$ for every $l = 1,...,\ell+1$ and parameter space $\Theta \subseteq \{(\theta_{l,b},\theta_{l,w})_{l = 1}^{\ell+1} \; | \; \theta_{l,b} \in \mathbb{R}^{h_l}, \; \theta_{l,w} \in \mathbb{R}^{h_l \times h_{l-1}} \; \forall l = 1,...,\ell+1\}$.  The parameter $\theta_{l,b}$ is often called the \textbf{\emph{bias}} and $\theta_{l,w}$ is often called the \textbf{\emph{weights}} of the $l^{th}$ hidden layer.
Jointly, the number of hidden layers $\ell$, the number of nodes within each hidden layer $h_l$ for $l \in \{1,...,\ell\}$, the activation functions $\Phi_l$ for $l \in \{1,...,\ell+1\}$, and the restrictions on the parameter space $\Theta$ are called the \textbf{\emph{hyperparameters}} of the neural network.
\end{definition}
\begin{remark}\label{rem:density}
The connection structure of a neural network in layer $l$ is defined through the sparsity structure of the weights $\theta_{l,w}$ imposed within the parameter space $\Theta$.  A dense neural network, which we consider throughout the remainder of this work, is one in which $\Theta = \{(\theta_{l,b},\theta_{l,w})_{l = 1}^{\ell+1} \; | \; \theta_{l,b} \in \mathbb{R}^{h_l}, \; \theta_{l,w} \in \mathbb{R}^{h_l \times h_{l-1}} \; \forall l = 1,...,\ell+1\}$.
\end{remark}
Before continuing, we wish to consider some common choices of activation functions used within the neural network literature.
\begin{example}\label{ex:activation}
Each of the following activation functions $\phi: \mathbb{R} \to \mathbb{R}$ are often applied component-wise in practice (i.e., such that $\Phi( \textbf{z}) = (\phi(z_1),...,\phi(z_k))^\top$ for $\textbf{z} \in \mathbb{R}^k$).  Let $z \in \mathbb{R}$ be arbitrary.
\begin{enumerate}
\item \textbf{Linear}: Consider the identity mapping $\phi(z) := z$.  This activation function is often denoted as the linear activation function as applying it along with the bias and weights creates a linear regression within that layer of the neural network.  Because the composition of linear functions is again linear, this activation function is typically only used for the output layer $\ell+1$ in practice.
\item \textbf{Rectified linear unit [ReLU]}: Consider the positive mapping $\phi(z) := z^+$.  This activation function takes just the positive part of the input, setting all negative values to $0$.
\item \textbf{Smooth ReLU}: Consider the mapping $\phi(z) := \log(1+\exp(z))$.  By construction, this activation function provides a positive value that limits to the identity mapping for positive inputs and 0 for negative inputs.  However, unlike the ReLU activation function, this activation is smooth.
\item \textbf{Hyperbolic tangent}: Consider the mapping $\phi(z) := \operatorname{tanh}(z)$.  This activation function provides an S-shaped curve that tempers extreme values as it limits to plus or negative 1 for positive of negative inputs respectively.
\end{enumerate}
\end{example}

One fundamental reason for the prevalence of neural networks within the machine learning community is the so-called \emph{universal approximation theorem}.  This result says that for any \emph{continuous} target function $y^*: \mathcal{T} \to \mathbb{R}^m$, there exists a \emph{single-layer} ($\ell = 1$) neural network $y: \mathcal{T} \to \mathbb{R}^m$ that can uniformly approximate $y^*$.  This result has been expanded to deep neural networks ($\ell \geq 2$) in, e.g.,~\cite{kidger2020universal}.
\begin{theorem}\label{thm:uat}\cite[Theorem 3.1]{pinkus1999approximation}
Let $y^*: \mathcal{T} \to \mathbb{R}^m$ be a continuous function with compact domain and $\phi: \mathbb{R} \to \mathbb{R}$ be a continuous activation function (applied component-wise to vectors as in Example~\ref{ex:activation}).  For every $\epsilon > 0$ there exists a 1-layer neural network $y: \mathcal{T} \times \Theta \to \mathbb{R}^m$ with $h_1$ hidden nodes such that
\[\sup_{t \in \mathcal{T}} \inf_{\theta \in \Theta} \|y^*(t) - y(t,\theta)\|_{2} < \epsilon\]
if and only if $\phi$ is not a polynomial.
\end{theorem}

As noted within the statement of the universal approximation theorem (Theorem~\ref{thm:uat}), there only exists parameters $\theta \in \Theta$ such that this approximation holds.  To attempt to find these parameters, we undertake a \emph{training} regiment.
This is accomplished by minimizing some \textbf{loss function} $y(\cdot,\theta) \mapsto L(y(\cdot,\theta)) \in \mathbb{R}_+$ over the neural network (inputted as the function $t \mapsto y(t,\theta)$), i.e., we seek the parameters $\theta^* \in \Theta$ such that
\begin{align}
\label{eq:MLopt}
\theta^* \in \argmin_{\theta \in \Theta} L(y(\cdot,\theta)).
\end{align}
However, this minimizing argument need not exist in general and, furthermore, this may be a nonconvex optimization problem; we remark on these issues further within Remark~\ref{rem:training} below.
This process of minimizing the loss function is often called \textbf{training} the neural network and it is how the network \textbf{learns}.
Typically the loss function is evaluated at some \emph{finite} \textbf{training points} $\{t_1,...,t_K\} \subseteq \mathcal{T}$ so that it is defined via a decomposition $L(y(\cdot,\theta)) := \sum_{k = 1}^K \bar L(t_k , y(t_k,\theta))$ for some function $\bar L: \mathcal{T} \times \mathbb{R}^m \to \mathbb{R}_+$.  The loss functions $L$ and $\bar L$ differ insofar as $L$ takes the entire neural network as an input whereas $\bar L$ only considers the value of the neural network at (pre-specified) training points.  The loss $\bar L(t,y(t,\theta))$ can be thought of as a distance the neural network $y(t,\theta)$ is from its target value $y^*(t)$. In particular, we use the convention that $\bar L(t,y(t,\theta)) = 0$ if and only if $y(t,\theta) = y^*(t)$.
This approach of decomposing the loss function over training points is the method we will use within Section~\ref{sec:nn-loss}.

\begin{remark}\label{rem:training}
As a note of caution, the universal approximation theorem only guarantees there exists some neural network which can uniformly approximate the target.  However, there is no guarantee that we actually constructed that neural network in practice when setting hyperparameters and training the parameters.  In particular, optimizing the parameters $\theta$ by minimizing the loss function often requires solving a highly nonlinear optimization problems. Though this problem is typically high-dimensional and nonlinear, multiple algorithms exist to efficiently \emph{approximate} $\theta^*$.  Within this work, we use the Adam optimizer~\cite{kingma2014adam} which is a gradient descent-based approach developed for large-scale problems. As a brief preview for the results of this work, we are essentially substituting the hard problem of convex vector optimization (Problem~\eqref{eq:cvop}) with the previous scalar optimization approach (Problem~\eqref{eq:MLopt}) encoded within these neural network optimization methodologies.
\end{remark}

For the remainder of this work, to ease notation, we will drop the dependence of neural networks on their parameters $\theta$.

\section{A primal-dual neural network}\label{sec:nn}
Within this section our goal is to construct two neural networks which approximate the primal and dual solutions ($x^*(\cdot),\lambda^*(\cdot)$) to the scalarization problems~\eqref{eq:scalarization}.  These neural networks will be constructed to guarantee primal and dual feasibility as encoded in \eqref{eq:kkt-primal} and \eqref{eq:kkt-dual}, i.e., a primal feasible neural network $x: \mathcal{W} \to \mathbb{X}$ and a dual neural network $\lambda: \mathcal{W} \to \mathbb{R}^M_+$. In addition, we propose a loss function to \emph{jointly} train these neural networks based on the KKT conditions of Theorem~\ref{thm:kkt}.

\subsection{Primal feasible neural network}\label{sec:nn-primal}
Herein we define a neural network with $P$ inputs (given by $w \in \mathcal{W}$) and $N$ outputs that guarantees primal feasibility~\eqref{eq:kkt-primal}.  In order to accomplish this, we want to construct the final layer of the neural network $x:\mathcal{W} \to \mathbb{R}^N$ so that the constraints encoded within $\mathbb{X}$ are satisfied.

Within the following proposition, we construct a projection mapping which, if used as the final activation function for an arbitrary neural network, guarantees primal feasibility.  As such, we use the general recursive structure of a neural network to guarantee that the primal neural network $x:\mathcal{W} \to \mathbb{R}^N$ is feasible for every $w \in \mathcal{W}$.
\begin{proposition}\label{prop:primal}
Let $z: \mathcal{W} \to \mathbb{R}^N$ define an arbitrary neural network.  Let $x: \mathcal{W} \to \mathbb{R}^N$ such that
\begin{align}
\label{eq:x} x(w) &:= (1-t^*(z(w)))z(w) + t^*(z(w))\bar x, \\
\label{eq:t} t^*(z) &:= \max_{j: g_j(z) > 0} \left\{\frac{g_j(z)}{g_j(z) - g_j(\bar x)}\right\},
\end{align}
where $t^*(z) = 0$ if $g(z) \in -\mathbb{R}^M_+$ and $\bar x$ is as defined in Assumption~\ref{ass:slater}.
Then $x(w) \in \mathbb{X}$ for every $w \in \mathcal{W}$.
\end{proposition}
\begin{proof}
Recall by construction of the feasible set that $x(w) \in \mathbb{X}$ if and only if $g_j(x(w)) \leq 0$ for every $j = 1,...,M$.
Consider the $j^{th}$ inequality constraint.  Note that $t^*(z) \in [0,1)$ for any $z$ by construction of $\bar x$ strictly feasible.  Therefore, by convexity and the definition of $t^*$,
\begin{align*}
g_j(x(w)) &= g_j((1-t^*(z(w))) z(w) + t^*(z(w))\bar x)\\
&\leq (1-t^*(z(w))) g_j(z(w)) + t^*(z(w)) g_j(\bar x) \leq 0.
\end{align*}
\end{proof}
We often refer to the final activation function $z \mapsto (1-t^*(z)) z + t^*(z) \bar x$ provided in Proposition~\ref{prop:primal} as a projection operator as it projects the neural network into the feasible region $\mathbb{X}$.

\begin{remark}\label{rem:tolerance}
Due to rounding errors in practice, it is recommended to add a tolerance level to the mapping $t^*$ defined in~\eqref{eq:t}.  That is, define the neural network $x: \mathcal{W} \to \mathbb{X}$ using $t^*(\cdot;\epsilon): \mathbb{R}^N \to [0,1)$ provided by
\[t^*(z;\epsilon) := \max_{j: g_j(z) \geq -\epsilon} \left\{\frac{g_j(z) + \epsilon}{g_j(z) - g_j(\bar x)}\right\}\]
for fixed tolerance $\epsilon \in (0 , -\max_j g_j(\bar x))$.
\end{remark}

For the remainder of this work we assume the primal neural network is constructed so as to be feasible, i.e., as in Proposition~\ref{prop:primal} with the final activation function of $z \mapsto (1-t^*(z)) z + \bar x$ to a generic neural network $z: \mathcal{W} \to \mathbb{R}^N$.
\begin{remark}\label{rem:primal-feasible}
We propose the projection operator as a general formulation for a final activation function which guarantees primal feasibility.  If the feasible set $\mathbb{X}$ takes a standard form, then a more direct approach can be utilized instead.  For instance, within Section~\ref{sec:mean-var} below, $\mathbb{X}$ is the unit simplex; therein one could instead also guarantee primal feasibility by taking the normalized rectified linear units (i.e., $x_i(w) := z_i(w)^+ / \sum_{j = 1}^N z_i(w)^+$ for every $i$ if $\sum_{j = 1}^N z_i(w)^+ > 0$ and $x(w) = \bar{x}$ otherwise) without needing to directly eliminate the equality constraint through the use of $A_{\perp}$ and changing the input space as described in Remark~\ref{rem:convex}\eqref{rem:convex-equality}.
\end{remark}

\subsection{Dual feasible neural networks}\label{sec:nn-dual}
Herein we define the dual neural network with $P$ inputs (given by $w \in \mathcal{W}$) and $M$ outputs which guarantees dual feasibility~\eqref{eq:kkt-dual}.  Similar to the approach taken for the primal neural network, we want to construct the final layer of this neural network so that feasibility is guaranteed.

Specifically, feasibility for the dual variables associated with the inequality constraints requires non-negativity~\eqref{eq:kkt-dual}.  This can be accomplished by simply applying the rectified linear unit as introduced in Example~\ref{ex:activation} (or the smooth approximation of the positive function) as a final activation function to a neural network with $M$ outputs already, i.e., $\lambda(w) := \tilde \lambda(w)^+$ for neural network $\tilde \lambda: \mathcal{W} \to \mathbb{R}^M$.

\subsection{Loss function}\label{sec:nn-loss}

In order to train the parameters of our two neural networks, we need to also construct a loss function which we seek to minimize.  In particular, we design a single loss function to jointly train both neural networks simultaneously.  To accomplish this, we will use the KKT conditions presented in Theorem~\ref{thm:kkt}.

Specifically, to find the parameters for our neural networks, we select $K$ (representative) choices of the scalarization weights $w_k \in \mathcal{W}$ for $k = 1,...,K$ and seek to minimize the deviation our neural networks have from the KKT conditions.  Due to the design of our neural networks $(x(\cdot),\lambda(\cdot))$ as presented above, the primal and dual feasibility conditions~\eqref{eq:kkt-primal} and~\eqref{eq:kkt-dual} are automatically satisfied and, as such, do not need to be considered in our loss function.  Therefore, the deviation from the KKT conditions can be fully characterized by the norms of the right-hand sides of the first order condition~\eqref{eq:kkt} and complimentary slackness condition~\eqref{eq:kkt-cs}. As these errors may be of different orders of magnitude, we choose to weight the error in the complimentary slackness condition by $\eta > 0$.  That is, our loss function $L$ for jointly training the two neural networks is provided by
\begin{equation}\label{eq:loss}
L(x(\cdot),\lambda(\cdot);\eta) := \sum_{k = 1}^K \left(\begin{array}{l} \|Jf(x(w_k))^\top w_k + Jg(x(w_k))^\top \lambda(w_k)\|_2^2 \\ \quad + \eta \|\operatorname{diag}(\lambda(w_k)) g(x(w_k))\|_2^2\end{array}\right).
\end{equation}
\begin{remark}\label{rem:tolerance-loss}
The loss associated with the complimentary slackness conditions may be modified to account for the rounding errors as discussed in Remark~\ref{rem:tolerance}.  For instance, we can consider the complimentary slackness error to be $\|\operatorname{diag}(\lambda(w_k)) \operatorname{diag}(\mathbbm{1}_{g_j(x(w_k)) \leq -\epsilon}) g(x(w_k))\|_2^2$ for tolerance $\epsilon > 0$ with indicator function $\mathbbm{1}_{(\cdot)}$ taking value 1 if the argument is true and 0 otherwise.  In this way, we treat any solution within $\epsilon$ of the boundary to be sitting on the boundary when considering complimentary slackness.
\end{remark}

\begin{remark}\label{rem:dualgap}
We wish to note that the KKT condition based loss function~\eqref{eq:loss} introduced above, and used throughout the numerical examples herein, is not the only loss function that can be taken for this problem.  As Assumption~\ref{ass:slater} implies strong duality, the duality gap can be used as a loss function instead via
\[L(x(\cdot),\lambda(\cdot)) := \sum_{k = 1}^K \left[w_k^\top f(x(w_k)) - d(\lambda(w_k),w_k)\right],\]
where $d: \mathbb{R}^M \times \mathcal{W} \to \mathbb{R}$ is the dual function provided in~\eqref{eq:dual}.
Notably, this construction removes the direct need to jointly train the primal and dual neural networks as the primal neural network could be trained with loss function $L_p(x(\cdot)) := \sum_{k = 1}^K w_k^\top f(x(w_k))$ and the dual neural network could be trained with loss function $L_d(\lambda(\cdot)) := -\sum_{k = 1}^K d(\lambda(w_k),w_k)$ to achieve the same effect.
\end{remark}

\section{Approximating the weakly efficient frontier}\label{sec:approx}
Within this section we present a method to formally use the neural networks constructed in Section~\ref{sec:nn} to construct inner and outer approximations of the weakly efficient frontier.

Before formally providing the results on the inner and outer approximations, we wish to recall the (Lagrange) dual function~\eqref{eq:dual} to the scalarization problem~\eqref{eq:scalarization}.  That is, we consider the mapping $d: \mathbb{R}^M \times \mathcal{W} \to \mathbb{R}$ defined by
\begin{equation*}
d(\lambda,w) := \inf_{x \in \mathbb{R}^N} \left[w^\top f(x) + \lambda^\top g(x)\right].
\end{equation*}
The dual problem is then to maximize this dual function subject to the dual feasibility constraints~\eqref{eq:kkt-dual} as provided in~\eqref{eq:scalarization-dual}.

\begin{lemma}\label{lemma:approx}
The primal neural network $x: \mathcal{W} \to \mathbb{X}$ provides an inner approximation of the upper set $\mathcal{P}$ and the dual neural network $\lambda: \mathcal{W} \to \mathbb{R}^M_+$ provide an outer approximation of the upper set $\mathcal{P}$ via the relation
\[\operatorname{cl}\operatorname{co}\bigcup_{w \in \mathcal{W}}\left[f(x(w)) + \mathbb{R}^P_+\right] \subseteq \mathcal{P} \subseteq \bigcap_{w \in \mathcal{W}} \left\{y \in \mathbb{R}^P \; | \; w^\top y \geq d(\lambda(w);w)\right\}.\]
\end{lemma}
\begin{proof}
We will prove this result by showing, first, the inner approximation and, second, the outer approximation.
Let $y \in \bigcup_{w \in \mathcal{W}} [f(x(w)) + \mathbb{R}^P_+]$.  There exists some $\tilde w \in \mathcal{W}$ such that $y - f(x(\tilde w)) \in \mathbb{R}^P_+$.  Therefore, as $x^*(\tilde w)\in\mathbb{X}$ by Proposition~\ref{prop:primal}, for any $w \in \mathcal{W}$, it holds
\[w^\top y \geq w^\top f(x(\tilde w)) \geq \inf\{w^\top f(x) \; | \; x \in \mathbb{X}\}.\]
By Proposition~\ref{prop:scalar}, and noting that the upper set $\mathcal{P}$ is closed and convex by assumption, the inner approximation holds.
Now, let $y \in \mathcal{P}$.  There exists some $\tilde w \in \mathcal{W}$ such that $y - f(x^*(\tilde w)) \in \mathbb{R}^P_+$ for solution $x^*: \mathcal{W} \to \mathbb{X}$ to \eqref{eq:scalarization}.  Therefore, for any $w \in \mathcal{W}$
\begin{align*}
w^\top y &\geq w^\top f(x^*(\tilde w)) \geq \inf\{w^\top f(x) \; | \; x \in \mathbb{X}\} \\
&= \sup\{d(\lambda,w) \; | \; \lambda \in \mathbb{R}^M_+\} \geq d(\lambda(w),w),
\end{align*}
where the equality holds by strong duality.
(Though we utilize strong duality here, in fact only weak duality is required.)
As a direct consequence, the outer approximation holds.
\end{proof}

Beyond defining the inner and outer approximations as done in Lemma~\ref{lemma:approx}, we want to also quantify the error from these approximations. To do this, we first wish to define a variable approximation error.  This is accomplished by quantifying the approximation error in any scalarization direction.  As opposed to the standard, constant, approximation error (which we discuss more in Proposition~\ref{prop:defn-approx}) this variable approximation error allows us to more accurately assess the local errors coming from our machine learning approach.  As far as the authors are aware, this definition is novel to the literature.
\begin{definition}\label{defn:approx}
$\tilde\mathcal{P} \subseteq \mathbb{R}^P$ is an \textbf{\emph{$\epsilon(\cdot)$-inner approximation}} of~\eqref{eq:cvop} for the function $\epsilon: \mathcal{W} \to \mathbb{R}_+$ if
\[\tilde\mathcal{P} \subseteq \mathcal{P} \subseteq \bigcap_{w \in \mathcal{W}} \operatorname{cl}\left[\tilde\mathcal{P} - \epsilon(w){\bf 1} + G(w)\right]\]
for the upper set $\mathcal{P}$, and where $G(w) := \{y \in \mathbb{R}^P \; | \; w^\top y \geq 0\}$ for any $w \in \mathcal{W}$.
\end{definition}
Notably, the quantification of the errors in the $\epsilon(\cdot)$-inner approximation are variable.  In contrast, typically (see, e.g., \cite{LRU14}) a constant error $\epsilon \in \mathbb{R}_{++}$ is taken instead; in such a setting an $\epsilon$-inner approximation $\tilde\mathcal{P} \subseteq \mathbb{R}^P$ of $\mathcal{P}$ is defined by the relations
\[\tilde\mathcal{P} \subseteq \mathcal{P} \subseteq \tilde\mathcal{P} - \epsilon{\bf 1}.\]
In the following proposition we show that when the error function is constant, Definition~\ref{defn:approx} reduces to the classical definition of an approximation.
\begin{proposition}\label{prop:defn-approx}
Let $\tilde\mathcal{P} \subseteq \mathbb{R}^P$ be a closed and convex upper set and let $\epsilon:\mathcal{W} \to \mathbb{R}_+$ be a constant function, i.e., $\epsilon(w) \equiv \epsilon_0 > 0$ for every $w \in \mathcal{W}$.  $\tilde\mathcal{P}$ is an $\epsilon(\cdot)$-inner approximation of $\mathcal{P} \subseteq \mathbb{R}^P$ if and only if $\tilde\mathcal{P} \subseteq \mathcal{P} \subseteq \tilde\mathcal{P} - \epsilon_0 {\bf 1}$.
\end{proposition}
\begin{proof}
As $\tilde\mathcal{P}$ is a closed and convex upper set, it trivially follows (by the separating hyperplane theorem) that
\[\tilde\mathcal{P} - \epsilon_0 {\bf 1} = \bigcap_{w \in \mathcal{W}}[\tilde\mathcal{P} - \epsilon_0{\bf 1} + G(w)]\]
and the result is proven.
\end{proof}

Within Corollary~\ref{cor:approx} below, we provide an analytical form for the functional approximation error $\epsilon(\cdot)$ for the inner approximation defined within Lemma~\ref{lemma:approx}.  Specifically, this function depends on both the primal and dual neural networks to quantify an upper bound to the ``distance'' between the inner and outer approximations at any point on the weakly efficient frontier.
\begin{corollary}\label{cor:approx}
The primal neural network $x: \mathcal{W} \to \mathbb{X}$ provides an $\epsilon(\cdot)$-inner approximation $\operatorname{cl}\operatorname{co}\bigcup_{w \in \mathcal{W}} [f(x(w)) + \mathbb{R}^P_+]$ of the upper set $\mathcal{P}$ with $\epsilon: \mathcal{W} \to \mathbb{R}_+$ defined by
\[\epsilon(w) := w^\top f(x(w)) - d(\lambda(w),w)\]
for any $w \in \mathcal{W}$ based, additionally, on the dual neural network $\lambda: \mathcal{W} \to \mathbb{R}^M_+$.
\end{corollary}
\begin{proof}
Let $\tilde\mathcal{P} := \operatorname{cl}\operatorname{co}\bigcup_{w \in \mathcal{W}} [f(x(w)) + \mathbb{R}^P_+]$.  It immediately follows that
\begin{align*}
\operatorname{cl}[\tilde\mathcal{P} - \epsilon(w) {\bf 1} + G(w)] &= \left\{y \in \mathbb{R}^P \; | \; w^\top y \geq \inf_{\tilde y \in \tilde\mathcal{P}} w^\top [\tilde y - \epsilon(w){\bf 1}]\right\}\\
&= \left\{y \in \mathbb{R}^P \; | \; w^\top y \geq \inf_{\tilde y \in \tilde\mathcal{P}} w^\top \tilde y - \epsilon(w)\right\}\\
&\supseteq \left\{y \in \mathbb{R}^P \; | \; w^\top y \geq w^\top f(x(w)) - \epsilon(w)\right\}\\
&= \left\{y \in \mathbb{R}^P \; | \; w^\top y \geq d(\lambda(w),w)\right\}.
\end{align*}
The result trivially follows from Lemma~\ref{lemma:approx}.
\end{proof}

\section{Numerical case studies}\label{sec:cs}
Within this section we will consider three multiobjective problems to numerically study the neural network approach presented above.  The first problem (Section~\ref{sec:2obj}) has 2 objective functions but many decision variables and constraints. In the second problem (Section~\ref{sec:many}), we investigate the performance of the proposed methodology as the number of objectives grow.  The third problem (Section~\ref{sec:mean-var}) presents the mean-variance optimization problem~\cite{markowitz1952} over a large number of assets.  All computations of the neural networks were completed with PyTorch on a local machine using the Adam optimizer~\cite{kingma2014adam} with learning rate $10^{-4}$.
Throughout these examples, to generate the loss function, we weight the complimentary slackness condition~\eqref{eq:kkt-cs} with $\eta = 10$.
We wish to note that all hyperparameters -- except the terminal activation functions as presented in Section~\ref{sec:nn} -- are chosen arbitrarily and were not found via, e.g., a grid search.

\subsection{Two objective problem}\label{sec:2obj}
Consider the following problem with $P = 2$ objectives:
\begin{align*}
\min\{(\|x\|_2^2 , \|x-{\bf 2}\|_2^2)^\top/N \; | \; x \in [0,1]^N\}.
\end{align*}
Note that this problem was considered as ``Test Instance 2'' in~\cite{eichfelder2021approximation}.
Due to the symmetric structure of this problem, the true unique minimizer for any scalarization $w \in \mathcal{W}$ can trivially be deduced as $x^*: \mathcal{W} \to [0,1]^N$ defined by
\[x^*(w) := \begin{cases} 2w_2{\bf 1} &\text{if } w_2 \leq \frac{1}{2} \\ {\bf 1} &\text{if } w_2 > \frac{1}{2} \end{cases}.\]
As this problem has a closed form set of minimizers, the efficient frontier can be exactly provided as $\{(\|x^*(w)\|_2^2,\|x^*(w)-2\|_2^2)^\top/N \; | \; w \in \mathcal{W}\}$.  This permits us to present an exact comparison between the machine learning methodology to the ground truth.

In order to fully consider our outer approximation as presented in Lemma~\ref{lemma:approx}, we need to also discuss the Lagrange dual problem of the weighted-sum scalarizations.
Let us encode the box constraints through the linear inequalities $Ax \leq b$ for $A := (I , -I)^\top \in \mathbb{R}^{N \times 2N}$ and $b := ({\bf 1}^\top , {\bf 0}^\top)^\top \in \mathbb{R}^{2N}$.  In such a way we consider $M = 2N$ linear inequality constraints.
Due to the quadratic structure of the scalarizations, the Lagrange dual function $d: \mathbb{R}^M \times \mathcal{W} \to \mathbb{R}$ can be directly computed by the quadratic structure
\[d(\lambda,w) := -\frac{N}{4}\lambda^\top A A^\top \lambda + (2w_2 A{\bf 1} - b)^\top \lambda + 4w_1w_2.\]

With this setup, we can consider a specific instantiation of the problem.  In particular, for test purposes, we will simply consider this bi-objective problem with $N = 40$ decision variables.  As a direct consequence of the box constraints, there are $M = 80$ inequality constraints.\footnote{The neural network was considered with objective function $x \mapsto N f(x)$ as, numerically, the training process was found to discount the loss associated with the objective too heavily otherwise.}
Furthermore, to demonstrate the power of the neural network methodology proposed above, we train the primal and dual neural networks with only 4 scalarizations: $w \in \{(0,1)^\top,(\frac{1}{3},\frac{2}{3})^\top,(\frac{2}{3},\frac{1}{3})^\top,(1,0)^\top\}$.
The primal neural network $x: \mathcal{W} \to \mathbb{X}$ is designed with three hidden layers, each with 800 hidden nodes; the terminal (projection) activation function is considered with a tolerance of $5 \times 10^{-5}$ to guarantee primal feasibility as discussed in Remark~\ref{rem:tolerance} and with strictly feasible point $\bar x := \frac{1}{2}\times{\bf 1}$.  The dual neural network $\lambda: \mathcal{W} \to \mathbb{R}^{80}_+$ is designed with three hidden layers, each with 1600 hidden nodes.
All intermediate activation functions (i.e., all but the terminal activation functions) are chosen to be the hyperbolic tangent function with linear linking between layers.
Finally, we train these networks over 1000 epochs.

\begin{figure}[h]
\centering
\begin{subfigure}[t]{0.45\textwidth}
\centering
\includegraphics[width=\textwidth]{efficient.eps}
\caption{Plot of weak efficient frontier (black) and approximations through neural networks (blue and red) or direct computation (black dashed).}
\label{fig:2obj-efficient}
\end{subfigure}
~
\begin{subfigure}[t]{0.45\textwidth}
\centering
\includegraphics[width=\textwidth]{efficient-zoom.eps}
\caption{Zoomed in plot of the weak efficient frontier (black) and approximations through neural networks (blue and red) or direct computation (black dashed).}
\label{fig:2obj-efficient-zoom}
\end{subfigure}
~
\begin{subfigure}[t]{0.45\textwidth}
\centering
\includegraphics[width=\textwidth]{epsilon_linearscale.eps}
\caption{Plot of approximation errors $\epsilon(\cdot)$ in linear-scale so that the neural network (blue), realized approximations (red), and direct computation (black dashed) provide $\epsilon(\cdot)$-inner approximations.}
\label{fig:2obj-epsilon-linearscale}
\end{subfigure}
~
\begin{subfigure}[t]{0.45\textwidth}
\centering
\includegraphics[width=\textwidth]{epsilon.eps}
\caption{Plot of approximation errors $\epsilon(\cdot)$ in log-scale so that the neural network (blue), realized approximations (red), and direct computation (black dashed) provide $\epsilon(\cdot)$-inner approximations.}
\label{fig:2obj-epsilon}
\end{subfigure}
\caption{Section~\ref{sec:2obj}: Plot of the efficient frontier and approximation errors. Neural networks and direct computation are computed with the same training scalarizations $w \in \{(0,1)^\top,(\frac{1}{3},\frac{2}{3})^\top,(\frac{2}{3},\frac{1}{3})^\top,(1,0)^\top\}$.}
\label{fig:2obj}
\end{figure}
Figure~\ref{fig:2obj-efficient} displays the true weak efficient frontier along with two approaches for approximation.  The true weak efficient frontier is plotted with a solid black curve.  The inner approximation provided by the primal neural network (i.e., the boundary $\operatorname{bd}\operatorname{cl}\operatorname{co}\bigcup_{w \in \mathcal{W}} [f(x(w)) + \mathbb{R}^P_+]$) is plotted as a solid blue line and the outer approximation provided by the dual neural network (i.e., the boundary $\operatorname{bd}\bigcap_{w \in \mathcal{W}} \{y \in \mathbb{R}^P \; | \; w^\top y \geq d(\lambda(w),w)\}$) is plotted as a solid red line.  These machine learning approximations are compared with the direct computation of the primal and dual solutions for the 4 training scalarizations as dashed black lines.  Immediately noticeable, the neural network approximations almost completely overlap with each other except in the south east corner of the figure (i.e., for scalarizations $w \in \mathcal{W}$ with $w_1$ low).  The approximations can be more clearly seen in Figure~\ref{fig:2obj-efficient-zoom} in which the area around the training scalarization $w = (\frac{2}{3},\frac{1}{3})^\top$ is highlighted.
The errors $\epsilon(w)$ can be directly seen in Figures~\ref{fig:2obj-epsilon-linearscale} and~\ref{fig:2obj-epsilon}; these figures display the exact same data but in a linear and logarithmic scale respectively.  For these scalarizations the error for the exact computation drops to 0, whereas the neural network has errors growing to around 0.2. These errors are further improved (see the red line ``NN Realized Error'' in Figures~\ref{fig:2obj-epsilon-linearscale} and~\ref{fig:2obj-epsilon}) when using the realizations of the neural network approximation -- the convex hull in the inner approximation and the intersection in the outer approximation as provided in Lemma~\ref{lemma:approx}.  However, as is clear from all figures, for (nearly) any choice of scalarization outside ``low'' $w_1$, the neural network \emph{outperforms} the direct computation. Further, we wish to remind the reader that the neural networks were trained without hyperparameter tuning and, as such, the errors $\epsilon(w)$ for the machine learning approach can be further improved.

\subsection{Many objective problem}\label{sec:many}
We now wish to consider how the neural networks can consider a problem with many objectives.  In particular we will consider a problem with $P = N = M$ objectives
\begin{align*}
f_i(x) := (x_i - 1)^2 + \sum_{j \neq i} x_j^2
\end{align*}
for every $i = 1,...,P$ and with $M = N$ constraints $g_j(x) = f_j(x) - 1$ for every $j = 1,...,M$.
Note that this problem was considered as ``Test Instance 4'' in~\cite{eichfelder2021approximation}.
Unlike in the prior example of Section~\ref{sec:2obj}, we do not have a closed form solution to the minimizers to this problem.
However, the neural network methodology presented above can still be applied to deduce inner and outer approximations for the weakly efficient frontier.
To consider the outer approximation, we need to consider the Lagrange dual function $d: \mathbb{R}^M \times \mathcal{W} \to \mathbb{R}$ of the weighted-sum scalarizations.  For this problem, the Lagrange dual can be computed as
\[d(\lambda,w) := \sum_{i = 1}^P \left[(w_i + \lambda_i) f_i\left(\frac{w_1+\lambda_1}{\sum_{j = 1}^P [w_j + \lambda_j]} , \cdots , \frac{w_P + \lambda_P}{\sum_{j = 1}^P [w_j + \lambda_j]}\right) - \lambda_i\right]\]
for any $\lambda \in \mathbb{R}^M$ and $w \in \mathcal{W}$.
For test purposes we will vary $P = N = M \in \{2,...,20\}$ within this case study.

With this setup, we can consider specific instantiations of the problem at the different choices of $P = N = M$.
Herein, we train each of these problems with 100 uniformly chosen random scalarizations $w \in \mathcal{W}$ of the correct dimension.
Both primal and dual neural networks $x: \mathcal{W} \to \mathbb{X}$ and $\lambda \to \mathbb{R}^M_+$ are designed with 2 hidden layers, each hidden layer with $500$ hidden nodes.  The terminal (projection) activation function for the primal neural network is considered with a tolerance of $5 \times 10^{-5}$ to guarantee primal feasibility as discussed in Remark~\ref{rem:tolerance} and with strictly feasible point $\bar x := \frac{1}{N} \times {\bf 1}$.
All intermediate activation functions (i.e., all but the terminal activation functions) are chosen to be the hyperbolic tangent function with linear linking between layers.
Finally, we train these networks over 250 epochs.

\begin{figure}[h]
\centering
\begin{subfigure}[t]{0.45\textwidth}
\centering
\includegraphics[width=\textwidth]{boxplot-fixed.eps}
\caption{Box plot of approximation errors $\epsilon(w)$ from machine learning approach. The central line provides the median; the box shows the inter-quartile range; outliers (red plus-signs) are determined to be further from the median than 1.5 times the inter-quartile range.}
\label{fig:boxplot}
\end{subfigure}
~
\begin{subfigure}[t]{0.45\textwidth}
\centering
\includegraphics[width=\textwidth]{NNvsRand.eps}
\caption{Mean (solid) and median (dashed) estimation errors $\epsilon(w)$ from machine learning approach (black) and random primal/dual feasible points (blue).}
\label{fig:NNvsRand}
\end{subfigure}
\caption{Section~\ref{sec:many}: Visualizations of approximation errors $\epsilon(w)$ over 5000 (uniformly random) test scalarizations $w \in \mathcal{W}$ as size of the problem $P = N = M$ changes.}
\label{fig:many}
\end{figure}
In order to compare the performance of the inner and outer approximations as the size of the problem grows, we compute $\epsilon(w)$ for 5000 uniformly sampled scalarizations $w \in \mathcal{W}$ of the correct dimension.  These errors (in logarithmic scale) are displayed in Figure~\ref{fig:boxplot}.  Note that the errors $\epsilon(w)$ become more consistent as the dimension of the problem increases; this is noticeable with the shrinking of the error bars of the box plot.  Additionally, as can be seen in Figure~\ref{fig:NNvsRand}, both the mean and median of these errors decrease as the dimension of the problem grows.
This improvement in the mean and median of the errors is primarily due to the shrinking size of the feasible region as the dimension grows.  For instance, as the problem size $P = N = M$ grows, the feasible region can be proven to exist within a box with rapidly shrinking volume, i.e., $\mathbb{X} \subseteq [-\frac{1}{\sqrt{N-1}},\frac{1}{\sqrt{N-1}}]^N$.  This is made clear in Figure~\ref{fig:NNvsRand} as the inner approximation from selecting a random primal feasible point (uniformly selected in $[-\frac{1}{\sqrt{N-1}},\frac{1}{\sqrt{N-1}}]^N$ and projected onto $\mathbb{X}$ via the projection activation function proposed in Proposition~\ref{prop:primal}) and random dual feasible point (deterministically selecting $\lambda := {\bf 0}$) converges in approximation errors to the neural network approach.  Importantly, the neural networks uniformly outperform the random approximations for all tested problems.
Finally, we wish to highlight that the reported approximation errors $\epsilon(w)$ in Figure~\ref{fig:many} are without considerations of the convex hull for the inner approximation and intersection for the outer approximation as in Lemma~\ref{lemma:approx}.  For this reason, we view the reported errors as an upper bound on those that would be realized by the neural networks constructed in this example.

\subsection{Mean-variance optimization}\label{sec:mean-var}
Consider now the mean-variance optimization problem proposed originally by Markowitz~\cite{markowitz1952} as a bi-objective optimization problem.  Specifically, this problem consists of a portfolio optimizer who is simultaneously attempting to maximize her expected return and minimize her risk (as encoded by the variance of the portfolio returns).  We will present this problem with no-short selling allowed, i.e., the investor can only purchase assets for her portfolio.  As such we will consider the following $P = 2$ objectives
\[f_1(\hat{x}) = -\bar{r}^\top \hat{x} \quad \text{ and } \quad f_2(\hat{x}) = \frac{1}{2}\hat{x}^\top C \hat{x}\]
over $S \geq 2$ assets in which the agent can invest, and where $\bar{r} \in \mathbb{R}^S$ denotes
the mean returns for each asset  and $C \in \mathbb{R}^{S \times S}$ is the covariance structure.  The investor is constrained by $M = S$ inequality constraints $g_j(\hat{x}) = -\hat{x}_j$ for every $j = 1,...,M$ (encoding the no-short selling constraints).
In contrast to the prior examples, herein the optimizer is constrained so that she uses all funds that she has available (i.e., ${\bf 1}^\top \hat{x} = 1$).
As discussed within Remark~\ref{rem:convex}\eqref{rem:convex-equality}, we consider the modification of this problem to guarantee the equality constraint is imposed throughout; specifically, we consider $N = S-1$ variables with $\hat{x} := (x^\top , 1-{\bf 1}^\top x)^\top$ for $x \in \mathbb{R}^N$.
To guarantee that this problem is strictly convex we will restrict the set of scalarizations under consideration to $\mathcal{W}_+ := \{w \in \mathcal{W} \; | \; w_2 > 0\}$ with positive definite covariance matrix $C$ as is discussed in Remark~\ref{rem:convex}\eqref{rem:convex-convex}.

To consider the outer approximation, we need to consider the Lagrange dual function $d: \mathbb{R}^M \times \mathcal{W}_+ \to \mathbb{R}$ of the weighted-sum scalarizations.  For this mean-variance problem, the Lagrange dual can be computed as
\begin{align*}
d(\lambda,w) &:= w^\top f(A_{\perp} \tilde{x}^*(\lambda,w) + e_s) - \lambda^\top(A_{\perp}\tilde{x}^*(\lambda,w) + e_s),\\
\tilde{x}^*(\lambda,w) &:= (w_2 A_{\perp}^\top C A_{\perp})^{-1} \left[A_{\perp}^\top \lambda + w_1(\bar{r}_{-S} - \bar{r}_S{\bf 1}) - w_2 A_{\perp}^\top C e_s\right],\\
A_{\perp} &:= \left(I_{N \times N} \, , \, -{\bf 1}\right)^\top \in \mathbb{R}^{S \times N}, \qquad e_s := ({\bf 0}^\top , 1)^\top \in \mathbb{R}^S
\end{align*}
for any $\lambda \in \mathbb{R}^M$ and $w \in \mathcal{W}_+$.

We will implement this mean-variance optimization problem on $S = 492$ stocks in the S\&P500 index.\footnote{These 492 stocks are currently in the S\&P500 index and had daily data available on Yahoo Finance from January 1, 2018 through December 31, 2021.} Empirical means and covariances for daily log returns from January 1, 2018 through December 31, 2021 are used to construct the mean asset returns $\bar{r}$ and covariance matrix $C$.\footnote{The neural network was considered with adjusted empirical mean returns and covariance matrix so that the objective values have comparable orders of magnitude. Specifically, the mean returns are multiplied by $100^2/6$ and the covariance matrix is multiplied by $100^2$. Without this modification of $1/6$, the training process was found to overly weight the importance of the mean returns and consistently results in the maximum return portfolio. Furthermore, the multiplication of both objectives by $100^2$ was used to eliminate numerical issues with the invertibility of the covariance matrix $C$ for computing the outer approximation.}
Herein, we train this problem with $5$ scalarizations: $w \in \{(0,1)^\top, (\frac{1}{4},\frac{3}{4})^\top, (\frac{1}{2},\frac{1}{2})^\top, (\frac{3}{4},\frac{1}{4})^\top, (1-10^{-5},10^{-5})^\top\}$.  The primal neural network $x: \mathcal{W}_+ \to \mathbb{X}$ is designed with three hidden layers, each with 800 hidden nodes.  Due to the construction of the feasible region, we apply a ReLU activation function prior to the projection activation function, i.e., such that the neural network $z$ in Proposition~\ref{prop:primal} has a terminal ReLU activation function.  The terminal (projection) activation function for this primal neural network is considered with tolerance of $5 \times 10^{-5}$ to guarantee primal feasibility as discussed in Remark~\ref{rem:tolerance} and with strictly feasible point $\bar x := 7.5\times 10^{-5} {\bf 1}$.  The dual neural network $\lambda: \mathcal{W}_+ \to \mathbb{R}^{492}_+$ is designed with three hidden layers, each with 800 hidden nodes.  In contrast to the prior numerical examples, here we use the soft plus activation function (i.e., $\log(1 + \exp(\cdot))$) to guarantee dual feasibility.  All intermediate activation functions (i.e., all but the terminal activation functions) are chosen to be the hyperbolic tangent function with linear linking between layers.  Finally, we train these networks over 5000 epochs.

\begin{figure}[h]
\centering
\begin{subfigure}[t]{0.3\textwidth}
\centering
\includegraphics[width=\textwidth]{mv-efficient.eps}
\caption{Plot of weak efficient frontier (black) and approximations through neural networks (blue and red).}
\label{fig:mv-efficient}
\end{subfigure}
~
\begin{subfigure}[t]{0.3\textwidth}
\centering
\includegraphics[width=\textwidth]{mv-bullet.eps}
\caption{Plot of the Markowitz bullet (black) and approximations through neural networks (blue and red).}
\label{fig:mv-bullet}
\end{subfigure}
~
\begin{subfigure}[t]{0.3\textwidth}
\centering
\includegraphics[width=\textwidth]{mv-epsilon.eps}
\caption{Plot of the realized approximation errors $\epsilon(\cdot)$.}
\label{fig:mv-epsilon}
\end{subfigure}
\caption{Section~\ref{sec:mean-var}: Plot of the efficient frontier, Markowitz bullet, and approximation errors. Neural networks are computed with 5 regularly spaced training scalarizations.}
\label{fig:mv}
\end{figure}

Figure~\ref{fig:mv-efficient} displays the true weak efficient frontier along with the neural network inner and outer approximations.  The true weak efficient frontier is plotted with a solid black curve.  The inner approximation provided by the primal neural network (i.e., the boundary $\operatorname{bd}\operatorname{cl}\operatorname{co}\bigcup_{w \in \mathcal{W}} [f(x(w)) + \mathbb{R}^2_+]$) is plotted as a solid blue line and the outer approximation provided by the dual neural network (i.e., the boundary $\operatorname{bd}\bigcap_{w \in \mathcal{W}_+} \{y \in \mathbb{R}^2 \; | \; w^\top y \geq d(\lambda(w),w)\}$) is plotted as a solid red line.
Using the usual transformation of the the mean-variance problem, the Markowitz bullet approximations providing the relation between the standard deviation and mean of the portfolio return, are displayed in Figure~\ref{fig:mv-bullet}.
As shown in Figure~\ref{fig:mv-epsilon}, the realized approximation errors is on the order of $10^{-4}$ throughout the space of scalarizations $w \in \mathcal{W}_+$.  Note that the approximation errors jump up for $w_2 \approx 0$, i.e., when the problem approaches a linear objective (rather than strictly convex as imposed in Assumption~\ref{ass:slater}).

\section{Extending to non-strictly convex vector optimization problems}\label{sec:convex}
As stated within Remark~\ref{rem:convex}\eqref{rem:convex-convex}, within this section we wish to relax the strict convexity assumption of Assumption~\ref{ass:slater} so that we only require $f_i$ to be convex for every objective function $i = 1,...,P$.
That is, consider the convex vector optimization problem~\eqref{eq:cvop} such that
\begin{itemize}
\item $f_i$ are convex and continuously differentiable;
\item $g_j$ are convex and continuously differentiable;
\item \eqref{eq:cvop} is bounded;
\item Slater's condition holds with strictly feasible point $\bar x \in \mathbb{R}^N$.
\end{itemize}
With this set of assumptions, which differ from Assumption~\ref{ass:slater} only through the non-strict convexity of the objective functions, we again consider the weighted-sum scalarization problem~\eqref{eq:scalarization} and its dual~\eqref{eq:scalarization-dual}.

The impact of the potential non-strict convexity of the objective functions to the machine learning approach can be two-fold.
\begin{itemize}
\item
In contrast to the simple dual feasibility conditions~\eqref{eq:kkt-dual} under Assumption~\ref{ass:slater}, we might need to consider additional implicit dual feasibility conditions to guarantee finiteness of the Lagrange dual function~\eqref{eq:dual}. As this depends on the particular structure of the considered problem, we will discuss this in more detail in the applications where this issue occurs. The case of linear vector optimization problems are, for instance, discussed in Section~\ref{subsec:d} for the case of linear vector optimization problems.
\item
The approach considered in Section~\ref{sec:nn} needs to be modified as neural networks can be unstable for discontinuous functions (e.g., the universal approximation theorem (Theorem~\ref{thm:uat}) only applies for continuous functions).  Specifically, due to the possibility of multiplicity of optimizers for convex problems, there may not exist a continuous selector of optimizers.
\end{itemize}
To address the second point above, the approach taken herein consists of constructing feasible primal and dual neural networks to the original problem~\eqref{eq:cvop} but where the loss function~\eqref{eq:loss} is based on the KKT conditions of the scalarizations of an \emph{augmented} convex vector optimization problem with scalarization
\begin{equation}\label{eq:scalarization-convex}\tag{$\overline{\text{wP}}$}
\min\Big\{(1-\delta) \sum_{i = 1}^P w_i f_i(x) + \delta \bar f(x) \; | \; x \in \mathbb{X}\Big\},
\end{equation}
where $w \in \mathcal{W}$ with $\delta \in (0,1)$ small and $\bar f: \mathbb{R}^N \to \mathbb{R}$ \emph{strictly convex} and continuously differentiable.  In particular, setting $\bar f(x) := \|x\|_{2}^2$ implies that the objective of~\eqref{eq:scalarization-convex} is strongly convex and thus also strictly convex.

As demonstrated by Proposition~\ref{prop:berge} below, this augmentation provides a reliable approximation to the solution of the original weighted-sum scalarization~\eqref{eq:scalarization}.
\begin{proposition}\label{prop:berge}
Assume $\mathbb{X} \subseteq \mathbb{R}^N$ is compact and fix $w \in \mathcal{W}$.
Let $X^*(w) \subseteq \mathbb{X}$ be the set of primal solutions of the scalarization problem~\eqref{eq:scalarization} w.r.t.\ $w$.
Additionally, let $\bar{x}^*(w,\delta)$ be the primal solutions of the scalarization problem~\eqref{eq:scalarization-convex} w.r.t.\ $w$ and $\delta \in (0,1)$.
Then there exists at least one accumulation point of $\{\bar{x}^*(w,\delta) \; | \; \delta \in (0,1)\}$ as $\delta \searrow 0$ and any such accumulation point is an element of $X^*(w)$, i.e., $\lim_{\delta \searrow 0} \bar{x}^*(w,\delta) \subseteq X^*(w)$.
\end{proposition}
\begin{proof}
Consider the extension of~\eqref{eq:scalarization-convex} to include $\delta \in [0,1)$, i.e., such that \eqref{eq:scalarization-convex} reduces to the original weighted-sum scalarization \eqref{eq:scalarization} at $\delta = 0$.  Let $\bar{X}^*(w,\delta) \subseteq \mathbb{X}$ denote the set of optimizers for~\eqref{eq:scalarization-convex} for any $\delta \in [0,1)$, i.e.,
\[\bar{X}^*(w,\delta) := \begin{cases} X^*(w) &\text{if } \delta = 0 \\ \{\bar{x}^*(w,\delta)\} &\text{if } \delta \in (0,1). \end{cases}\]
By the Berge maximum theorem (see, e.g.,~\cite[Theorem 17.31]{AB07}), $\delta \in [0,1) \mapsto \bar{X}^*(w,\delta)$ is upper continuous.  In particular, by the closed graph theorem (see, e.g.,~\cite[Theorem 17.11]{AB07}), the graph of $\delta \in [0,1) \mapsto \bar{X}^*(w,\delta)$ is closed which implies any accumulation point $\bar{x}^*(w)$ of $\{\bar{x}^*(w,\delta) \; | \; \delta \in (0,1)\}$ as $\delta \searrow 0$ is an element of $X^*(w)$.
Finally, since $\{\bar{x}^*(w,\delta) \; | \; \delta \in (0,1)\} \subseteq \mathbb{X}$, there must exist at least one accumulation point as $\delta \searrow 0$.
\end{proof}

\begin{remark}
Under suitable additional assumptions for the uniqueness of the dual solution $\bar{\lambda}^*(w,\delta)$ to the Lagrange dual problem of~\eqref{eq:scalarization-convex}, a comparable statement to Proposition~\ref{prop:berge} can be given for $\bar{\lambda}^*(w,\delta) \to \lambda^*(w)$.
\end{remark}

To demonstrate this approach, we will first illustrate the exact construction for a generic linear vector optimization problems (LVOP) in Section~\ref{sec:convex-lvop} along with a simple numerical example in that setting.  We will then consider in Section~\ref{sec:mr} a version of the mean-risk problem to present a separate, practical problem from financial applications.

\subsection{Linear vector optimization problems}\label{sec:convex-lvop}
Consider the generic inequality-constrained LVOP
\[\min\{Cx \; | \; Ax \leq b\}\]
for $C \in \mathbb{R}^{P \times N},~A \in \mathbb{R}^{M \times N},~b \in \mathbb{R}^M$.  That is, $\mathbb{X} := \{x \in \mathbb{R}^N \; | \; Ax \leq b\}$.\footnote{Equality constraints can be included as discussed within Remark~\ref{rem:convex}\eqref{rem:convex-equality}.}
By construction, this is a CVOP that is \emph{not} strictly convex.  We will use this problem to illustrate the modifications that are necessary to formalize the neural network approach provided within Section~\ref{sec:nn} to guarantee the inner and outer approximation results provided within Section~\ref{sec:approx}.
Briefly, the idea is to construct primal and dual feasible neural networks to the LVOP that approximate the solution to a closely related \emph{strictly convex} CVOP.

Using this LVOP, let us construct our primal and dual neural networks
followed by a discussion of the necessary changes to the loss function for training purposes.  We wish to note that the inner and outer approximation results provided within Section~\ref{sec:approx} can be applied to the constructions provided within this section without modification.

\subsubsection{Primal feasible neural network}\label{subsec:p}
Consider the weighted-sum scalarization problem
\[\min\{w^\top Cx \; | \; Ax \leq b\}\]
for $w \in \mathcal{W}$.
Herein we construct the primal feasible neural network $x: \mathcal{W} \to \mathbb{X}$ exactly as in Section~\ref{sec:nn-primal} and, specifically, as detailed in Proposition~\ref{prop:primal}.
That is, let $z: \mathcal{W} \to \mathbb{R}^N$ define an arbitrary neural network.  Let $x: \mathcal{W} \to \mathbb{R}^N$ such that
\[x(w) := (1-t^*(z(w))) z(w) + t^*(z(w)) \bar x\]
where $t^*: \mathbb{R}^N \to [0,1]$ is defined as in~\eqref{eq:t}.

\subsubsection{Dual feasible neural network}\label{subsec:d}
Consider the Lagrange dual problem to the weighted-sum scalarization problem w.r.t.\ $w \in \mathcal{W}$.  Herein we need to be careful with the domain of the Lagrange dual function $d: \mathbb{R}^M \times \mathcal{W} \to \mathbb{R}$; specifically, we can explicitly compute $d$ as
\begin{align*}
d(\lambda,w) := \inf_{x \in \mathbb{R}^N} \left[w^\top C x + \lambda^\top(Ax - b)\right] = \begin{cases} -b^\top \lambda &\text{if } A^\top\lambda = -C^\top w \\ -\infty &\text{else}\end{cases}
\end{align*}
for any $\lambda \in \mathbb{R}^M,~w \in \mathcal{W}$.  Therefore the dual problem includes the additional equality constraint $A^\top\lambda = -C^\top w$, i.e.,~\eqref{eq:scalarization-dual} can be specified as
\[\max\{-b^\top\lambda \; | \; A^\top\lambda = -C^\top w,~\lambda \in \mathbb{R}^M_+\}.\]
For notational simplicity, define $\Lambda(w) := \{\lambda \in \mathbb{R}^M \; | \; -\lambda \in -\mathbb{R}^M_+,~A^\top \lambda = -C^\top w\}$ to provide the set of feasible dual solutions $\lambda$.

\begin{assumption}\label{ass:dual-feasible}
There exists a strictly feasible point $\bar{\lambda}(w) \in \Lambda(w)$ for every $w \in \mathcal{W}$.  That is, $A^\top\bar{\lambda}(w) = -C^\top w$ with $\bar{\lambda}(w) \in \mathbb{R}^M_{++}$ for every $w \in \mathcal{W}$.
\end{assumption}

\begin{remark}\label{rem:dual_strict_feasible}
A strictly feasible dual point $\bar{\lambda}(w) \in \Lambda(w)$ can be found by solving the linear program
\[\bar{\lambda}(w) \in \argmin\{\sum_{i = 1}^M \lambda_i \; | \; -\lambda + \bar{\epsilon} \in -\mathbb{R}^M_+,~A^\top \lambda = -C^\top w\}\]
for some $\bar{\epsilon} \in \mathbb{R}^M_{++}$.  Due to the existence of a strictly feasible point (Assumption~\ref{ass:dual-feasible}), so long as $\bar{\epsilon}$ is set small enough, this problem will be feasible.
\end{remark}

Herein our goal is to formulate a neural network $\lambda: \mathcal{W} \to \mathbb{R}^M$ that is dual feasible, i.e., $\lambda(w) \in \Lambda(w)$ for every $w \in \mathcal{W}$.
Let $\bar{\lambda}(w) \in \Lambda(w)$ define a \emph{strictly} feasible point for every $w \in \mathcal{W}$ as given in Assumption~\ref{ass:dual-feasible}.
As discussed within Remark~\ref{rem:convex}\eqref{rem:convex-equality}, we want to consider a reformulation of $\Lambda(w)$ to eliminate the need for the equality constraints so that this set takes the general form of $\mathbb{X}$.
Specifically, let $(A^\top)_{\perp} \in \mathbb{R}^{M \times (M-N)}$ provide a basis for the null space of $A^\top$ (assuming $A^\top$ has full rank).  Then $\Lambda(w)$ can be reformulated without the equality constraint as
\[\Lambda(w) = \{(A^\top)_{\perp} z + \bar{\lambda}(w) \; | \; -[(A^\top)_{\perp} z - \bar{\lambda}(w)] \leq {\bf 0},~z \in \mathbb{R}^{M-N}\}.\]

Within the following proposition, we construct an analog of Proposition~\ref{prop:primal} for this dual feasibility constraint.
\begin{proposition}\label{prop:dual}
Let $z: \mathcal{W} \to \mathbb{R}^{M-N}$ define an arbitrary neural network.  Let $\lambda: \mathcal{W} \to \mathbb{R}^M$ such that
\begin{align*}
\lambda(w) &:= (1-\bar{t}^*((A^\top)_\perp z(w) + \bar{\lambda}(w),\bar{\lambda}(w)))(A^\top)_\perp z(w) + \bar{y}(w),\\
\bar{t}^*(\hat{z},\bar{\lambda}) &:= \max_{j = 1,...,M} \left\{\frac{\hat{z}_j^-}{\hat{z}_j^- + \bar{\lambda}_j}\right\},
\end{align*}
where $\bar{t}^*(\hat{z},\bar{\lambda}) = 0$ if $\hat{z}_j \geq 0$ for every $j = 1,...,M$ and $\bar{\lambda}: \mathcal{W} \to \mathbb{R}^M$ provides a strictly feasible point for any $w \in \mathcal{W}$.  Then $\lambda(w) \in \Lambda(w)$ for every $w \in \mathcal{W}$.
\end{proposition}
Due to the similarity of this proposition to Proposition~\ref{prop:primal} we skip the proof of it here.

\subsubsection{Loss function}
Consider now the modified weighted-sum scalarization problem~\eqref{eq:scalarization-convex} with fixed weight $\delta \in (0,1)$.
An application of Theorem~\ref{thm:kkt} implies that $(\bar{x}^*(w,\delta),\bar{\lambda}^*(w,\delta)) \in \mathbb{R}^N \times \mathbb{R}^M$ are a pair of primal ($\bar{x}^*(w,\delta)$) and dual ($\bar{\lambda}^*(w,\delta)$) solutions of the scalarization problem~\eqref{eq:scalarization-convex}, respectively its Lagrange dual problem,  w.r.t.\ $w \in \mathcal{W}$ and $\delta \in (0,1)$ if and only if they jointly satisfy the KKT conditions
\begin{align*}
{\bf 0} &= (1-\delta)C^\top w + \delta \nabla\bar{f}(\bar{x}^*(w,\delta)) + A^\top \bar{\lambda}^*(w,\delta) \\
b &\geq A\bar{x}^*(w,\delta) \\
{\bf 0} &\leq \bar{\lambda}^*(w,\delta) \\
0 &= \bar{\lambda}_i^*(w,\delta) [A\bar{x}^*(w,\delta) - b]_i  \quad i=1,...,m.
\end{align*}
Note that the primal and dual feasibility constraints of the original, unmodified, scalarization problems~\eqref{eq:scalarization} are included in those of the modified scalarization problem~\eqref{eq:scalarization-convex}. Thus, the neural networks ($x(\cdot),\lambda(\cdot)$) constructed in Sections~\ref{subsec:p} and~\ref{subsec:d} from the original, unmodified, scalarization problem~\eqref{eq:scalarization} are automatically also feasible for the modified problem~\eqref{eq:scalarization-convex}. Hence, the only difference in methodology is to use the loss function of the modified problem~\eqref{eq:scalarization-convex} instead of the original loss function~\eqref{eq:loss} of problem~\eqref{eq:scalarization}.
The loss function based on the KKT conditions of the modified problem~\eqref{eq:scalarization-convex} can be constructed in analogy to Section~\ref{sec:nn-loss} as
\begin{equation*}
L(x(\cdot),\lambda(\cdot);\eta) := \sum_{k = 1}^K \left(\begin{array}{l} \|(1-\delta)C^\top w_k + \delta\nabla\bar{f}(x(w_k)) + A^\top \lambda(w_k)\|_2^2 \\ \quad + \eta \|\operatorname{diag}(\lambda(w_k)) [Ax(w_k) - b]\|_2^2\end{array}\right).
\end{equation*}
Recall from Proposition~\ref{prop:berge} that for $\delta \in (0,1)$ small enough, if the loss of the neural networks is small enough, we can confidently conclude that the obtained solutions $(\bar{x}^*(\cdot,\delta),\bar{\lambda}^*(\cdot,\delta))$ are ``close'' to the true solutions of~\eqref{eq:scalarization}.

\subsubsection{Inner and outer approximations}
Following the construction of the primal and dual neural networks for the LVOP, we want to comment on the inner and outer approximations as presented in Section~\ref{sec:approx}.
Formally, the approximation results of Section~\ref{sec:approx} follow directly from feasibility (both primal and dual feasibility) and strong duality.
As we recover feasibility of our neural networks, and strong duality still holds due to Slater's condition, no modifications are necessary to recover these approximations.
In particular, the primal neural network $x: \mathcal{W} \to \mathbb{R}^N$ provides an inner approximation and the dual neural network $\lambda: \mathcal{W} \to \mathbb{R}^M$ provides an outer approximation of the upper set $\mathcal{P}$ via the relation
\[\operatorname{cl}\operatorname{co}\bigcup_{w \in \mathcal{W}}\left[C x(w) + \mathbb{R}^P_+\right] \subseteq \mathcal{P} \subseteq \bigcap_{w \in \mathcal{W}} \left\{y \in \mathbb{R}^P \; | \; w^\top y \geq -b^\top\lambda(w)\right\}.\]

\subsubsection{Discussion and numerical results}

Consider now the bi-objective ($P = 2$) LVOP with $N = 2$ decision variables and $M = 5$ inequality constraints constructed by
\[C := \left(\begin{array}{cc} 1 & 0 \\ 0 & 1 \end{array}\right), \qquad A := \left(\begin{array}{rr} -2 & -1 \\ -1 & -2 \\ 1 & 1 \\ -1 & 0 \\ 0 & -1 \end{array}\right), \quad\text{ and }\quad b := \left(\begin{array}{r} -2 \\ -2 \\ 6 \\ 0 \\ 0 \end{array}\right).\]
To place it within the strictly convex setting, we introduce the augmented form with $\bar f(x) := \|x\|_{2}^2$ included with constant weight $\delta = 10^{-4}$.

As in Section~\ref{sec:cs}, the neural networks for this example were computed with PyTorch on a local machine using the Adam optimizer with learning rate $10^{-4}$.  Here, we weight the complimentary slackness condition~\eqref{eq:kkt-cs} with $\eta = 10^{-4}$ so that the first-order condition error~\eqref{eq:kkt} and complimentary slackness error are, initially, of the same order of magnitude.  We train these neural networks for 500 epochs.  We wish to note that all other hyperparameters -- except the terminal activation functions as presented within this section -- are chosen arbitrarily and were not found via, e.g., a grid search.

Specifically, we have two neural networks to design: the primal ($x: \mathcal{W} \to \mathbb{R}^2$) and dual ($\lambda: \mathcal{W} \to \mathbb{R}^5$) neural networks.  The primal neural network $x$ is designed with four hidden layers, each with 1200 hidden nodes; the terminal (projection) activation function is considered with a tolerance of $5 \times 10^{-5}$ to guarantee primal feasibility as discussed in Remark~\ref{rem:tolerance} and with strictly feasible point $\bar{x} := {\bf 1}$.  The dual neural network $\lambda$ is designed with three hidden layers, each with 1500 hidden nodes.  The terminal (projection) activation function is considered with a tolerance of $5 \times 10^{-5}$ to guarantee dual feasibility as discussed in Remark~\ref{rem:tolerance} and with strictly feasible points $\bar{\lambda}: \mathcal{W} \to \mathbb{R}^M$ constructed via the linear programming approach outlined within Remark~\ref{rem:dual_strict_feasible} with $\bar{\epsilon} := 5 \times 10^{-3} \times {\bf 1}$.
All intermediate activation functions (i.e., all but the terminal activation functions) are chosen to be the hyperbolic tangent function with linear linking between layers.

Figure~\ref{fig:lvop-efficient} displays the true weak efficient frontier along with the neural network inner and outer approximations.  The true weak efficient frontier is plotted with a solid black curve.  The inner approximation provided by the primal neural network (i.e., the boundary $\operatorname{bd}\operatorname{cl}\operatorname{co}\bigcup_{w \in \mathcal{W}} [Cx(w) + \mathbb{R}^2_+]$) is plotted as a solid blue line and the outer approximation provided by the dual neural network (i.e., the boundary $\operatorname{bd}\bigcap_{w \in \mathcal{W}} \{y \in \mathbb{R}^2 \; | \; w^\top y \geq -b^\top\lambda(w)\}$) is plotted as a solid red line.  Immediately noticeable, the neural network approximations almost completely overlap with each other and, in particular, capture the vertices of the true weak efficient frontier.
The errors $\epsilon(w)$ can be directly seen in Figure~\ref{fig:2obj-epsilon} which displays the data in a logarithmic scale.  For these scalarizations, the original neural network approaches provide relatively high errors (see the blue line ``NN Error'').  However, when using the realizations of the neural network approximation -- the convex hull in the inner approximation and the intersection in the outer approximation as provided in Lemma~\ref{lemma:approx} -- these errors drop significantly with most scalarizations having errors around $10^{-2}$.
Further, we wish to remind the reader that the neural networks were trained without hyperparameter tuning and, as such, the errors $\epsilon(w)$ for the machine learning approach can be further improved.

\begin{figure}[h]
\centering
\begin{subfigure}[t]{0.45\textwidth}
\centering
\includegraphics[width=\textwidth]{lvop-efficient.eps}
\caption{Plot of weak efficient frontier (black) and approximations through neural networks (blue and red).}
\label{fig:lvop-efficient}
\end{subfigure}
~
\begin{subfigure}[t]{0.45\textwidth}
\centering
\includegraphics[width=\textwidth]{lvop-epsilon.eps}
\caption{Plot of approximation errors $\epsilon(\cdot)$ in log-scale so that the neural network (blue) and realized approximations (red) provide $\epsilon(\cdot)$-inner approximations.}
\label{fig:lvop-epsilon}
\end{subfigure}
\caption{Section~\ref{sec:convex-lvop}: Plot of the efficient frontier and approximation errors. Neural networks are computed with the 20 regularly spaced training scalarizations.}
\label{fig:lvop}
\end{figure}

\subsection{Mean-risk problem}\label{sec:mr}
Herein we want to revisit a practical example by studying the mean-risk problem in finance with dynamic trading.  Similar to the mean-variance problem presented within Section~\ref{sec:mean-var}, this is a bi-objective problem that allows for easy visualization of the optimal portfolio.  However, in contrast to the mean-variance problem, herein we consider risk as measured by the entropic risk measure \cite[Example 4.34]{FS11}.  We will present this problem with no-short selling allowed, i.e., the investor can only purchase assets for her portfolio.  To simplify this setting we will consider a market with $S = 2$ assets only: asset $0$ is the cash asset and is worth $\$1$ throughout; asset $1$ is a risky stock whose value fluctuates (with equal up or down probability) according to the Cox-Ross-Rubinstein binomial tree~\cite{cox1979option} with $T = 2$ time steps and volatility $\sigma^2 = 0.05$ as depicted in Figure~\ref{fig:mr-tree}.

Within this setting, the investor is optimizing over $6$ variables -- investments in each asset at time $0$ and in each asset in either state at time $1$.
Let $\hat x = (\hat x_{0,0} , \hat x_{0,1} , \hat x_{1,0}^U , \hat x_{1,1}^U , \hat x_{1,0}^D , \hat x_{1,1}^D) \in \mathbb{R}^6$ provide the investments at times $0,1$ and in assets $0,1$.  The investor is seeking to maximize her expected return and minimize her entropic risk at time $T = 2$, i.e.,
\begin{align*}
f_1(\hat x) &= -\frac{1}{2}\left[\left(\hat x_{1,0}^U + \frac{S_2^{UU} + S_2^{UD}}{2}\hat x_{1,1}^U\right) + \left(\hat x_{1,0}^D + \frac{S_2^{UD} + S_2^{DD}}{2}\hat x_{1,1}^D\right)\right] \\
f_2(\hat x) &= \frac{1}{\alpha}\log\left(\frac{1}{4}\left[e^{-\alpha[\hat x_{1,0}^U + S_2^{UU}\hat x_{1,1}^U]} + e^{-\alpha[\hat x_{1,0}^U + S_2^{UD}\hat x_{1,1}^U]} + e^{-\alpha[\hat x_{1,0}^D + S_2^{UD}\hat x_{1,1}^D]} + e^{-\alpha[\hat x_{1,0}^D + S_2^{DD}\hat x_{1,1}^D]}\right]\right)
\end{align*}
for risk-aversion parameter $\alpha > 0$.  Herein we will consider $\alpha = 1/2$ chosen arbitrarily.
As this problem is convex, but not strictly convex, we augment this problem with a 3rd objective function: $f_3(\hat x) = \|\hat x\|_2^2$ which is always included with weight $\delta = 10^{-4}$.

Due to the no-short selling constraint, the investor is constrained by $M = S$ inequality constraints $g(\hat x) = -\hat x \in \mathbb{R}^6$.
In addition to these inequality constraints, this constructed dynamic portfolio is required to satisfy self-financing conditions with initial wealth $W = 10$ such that
\begin{align*}
W &= \hat x_{0,0} + S_0 \hat x_{0,1} \\
0 &= \left(\hat x_{0,0} + S_1^U \hat x_{0,1}\right) - \left(\hat x_{1,0}^U + S_1^U \hat x_{1,1}^U\right) \\
0 &= \left(\hat x_{0,0} + S_1^D \hat x_{0,1}\right) - \left(\hat x_{1,0}^D + S_1^D \hat x_{1,1}^D\right).
\end{align*}
As discussed within Sections~\ref{sec:background-vop} and~\ref{sec:mean-var}, we consider the modification of this problem to guarantee these equality constraints are satisfied.  Denote these self-financing constraints by $A\hat x = (W,0,0)^\top$; then we consider $N = 3$ variables with $\hat{x} := A_{\perp} x + \tilde x$ for $x \in \mathbb{R}^N$ with $A_{\perp} \in \mathbb{R}^{6 \times 3}$ and $\tilde x := \frac{W}{2}{\bf 1} \in \mathbb{R}^6$.

To consider the outer approximation, we need to consider the Lagrange dual function $d: \mathbb{R}^6 \times \mathcal{W} \to \mathbb{R}$ of the weighted-sum scalarizations (prior to augmentation).  For this mean-risk problem, the Lagrange dual function does not have a closed form solution, rather it is computed by solving the unconstrained problem
\begin{align*}
d(\lambda,w) := \inf_{x \in \mathbb{R}^3} \left[w^\top f(A_{\perp} x + \tilde x) - \lambda^\top(A_{\perp} x + \tilde x)\right]
\end{align*}
for any $\lambda \in \mathbb{R}^3$ and $w \in \mathcal{W}$.

With this setting, we have two neural networks to design: the primal ($x: \mathcal{W} \to \mathbb{R}^3$) and dual ($\lambda: \mathcal{W} \to \mathbb{R}^6$) neural entworks.
The primal neural network $x$ is designed with three hidden layers, each with 800 hidden nodes; the terminal (projection) activation function is considered with a tolerance of $5 \times 10^{-5}$ to guarantee primal feasibility as discussed in Remark~\ref{rem:tolerance} and with strictly feasible point $\bar{x} := {\bf 0} \in \mathbb{R}^3$.  The dual neural network $\lambda$ is designed with three hidden layers, each with 800 hidden nodes.  The terminal ReLU activation function is used to guarantee dual feasibility.
All intermediate activation functions (i.e., all but the terminal activation functions) are chosen to be the hyperbolic tangent function with linear linking between layers.
As with the prior examples, the neural networks herein are computed with PyTorch on a local machine using the Adam optimizer with learning rate $10^{-4}$.  Here, we weight the complimentary slackness condition~\eqref{eq:kkt-cs} with $\eta = 10$.  We train these neural networks for $2000$ epochs over 4 training points $\{(0,1)^\top , (\frac{1}{3},\frac{2}{3})^\top , (\frac{2}{3},\frac{1}{3})^\top , (1,0)^\top\}$.  We wish to note that all other hyperparameters -- except the terminal activation functions discussed previously -- are chosen arbitrarily and were not found via, e.g., a grid search.

\begin{figure}[h]
\centering
\begin{subfigure}[t]{0.3\textwidth}
\centering
\begin{tikzpicture}
\node at (0cm,0cm) (n0) {\footnotesize $S_0 = \$1.00$};
\node at (1.5cm,1cm) (nU) {\footnotesize $S_1^U = \$1.25$};
\node at (1.5cm,-1cm) (nD) {\footnotesize $S_1^D = \$0.80$};
\node at (3cm,2cm) (nUU) {\footnotesize $S_2^{UU} = \$1.56$};
\node at (3cm,0cm) (nUD) {\footnotesize $S_2^{UD} = \$1.00$};
\node at (3cm,-2cm) (nDD) {\footnotesize $S_2^{DD} = \$0.64$};

\draw[-,line width=.6mm] (n0) -- (nU);
\draw[-,line width=.6mm] (n0) -- (nD);
\draw[-,line width=.6mm] (nU) -- (nUU);
\draw[-,line width=.6mm] (nU) -- (nUD);
\draw[-,line width=.6mm] (nD) -- (nUD);
\draw[-,line width=.6mm] (nD) -- (nDD);
\end{tikzpicture}
\caption{Depiction of the price process as a tree. All transitions occur with probability $50\%$.}
\label{fig:mr-tree}
\end{subfigure}
~
\begin{subfigure}[t]{0.3\textwidth}
\centering
\includegraphics[width=\textwidth]{mr-efficient.eps}
\caption{Plot of weak efficient frontier (black) and approximations through neural networks (blue and red).}
\label{fig:mr-efficient}
\end{subfigure}
~
\begin{subfigure}[t]{0.3\textwidth}
\centering
\includegraphics[width=\textwidth]{mr-epsilon.eps}
\caption{Plot of the realized approximation errors $\epsilon(\cdot)$.}
\label{fig:mr-epsilon}
\end{subfigure}
\caption{Section~\ref{sec:mr}: Plot of the market model, efficient frontier, and approximation errors.  Neural networks are computed with 4 regularly spaced training scalarizations.}
\label{fig:mr}
\end{figure}

Figure~\ref{fig:mr-efficient} displays the true weak efficient frontier along with the neural network inner and outer approximations.  The true weak efficient frontier is plotted with a solid black curve.  The inner approximation provided by the primal neural network (i.e., the boundary $\operatorname{bd}\operatorname{cl}\operatorname{co}\bigcup_{w \in \mathcal{W}}[f(x(w)) + \mathbb{R}^2_+]$) is plotted as a solid blue line and the outer approximation provided by the dual neural network (i.e., the boundary $\operatorname{bd}\bigcap_{w \in \mathcal{W}}\{y \in \mathbb{R}^2 \; | \; w^\top y \geq d(\lambda(w),w)\}$) is plotted as a solid red line.  In this scenario, the dual neural network uniformly provide the zero output, i.e., $\lambda \equiv {\bf 0}$.  We conjecture this occurs because the true dual value $\lambda^*(w)$ is nearly identically 0 throughout, and only deviates when $w_1 \approx 1$ or $w_2 \approx 0$ with $\lambda^*((1,0)) \approx (0.0257 \, , \, 0 \, , \, 0.0126 \, , \, 0 \, , \, 0.0126 \, , \, 0)^\top$.  This can be observed with the strong fit of the outer approximation except in the upper left corner of Figure~\ref{fig:mr-efficient}.  The strength of these approximations can further be seen in Figure~\ref{fig:mr-epsilon} where the errors $\epsilon(w)$ are plotted directly in a logarithmic scale.  With an improvement in the dual neural network to capture behavior for high values of $w_1$, we would expect these approximation errors to remain on the order of $10^{-3}$ or below. We again wish to remind the reader that the neural network architectures considered herein were chosen arbitrarily and not as the result of hyperparameter tuning.

\section{Conclusion}\label{sec:conclusion}
Within this work we proposed a neural network architecture that provides inner and outer approximations of the weakly efficient frontier for convex vector optimization problems. Through the use of the dual of the scalarization, we provided a functional of the primal and dual neural networks which provides an upper bound to the error of the inner approximation at each point on the weakly efficient frontier. To demonstrate the ability of the proposed method for large-scale problems, we computed the approximations of the efficient frontiers for multiple convex vector optimization problems.

We wish to conclude with an extension to this methodology which we believe would be of great interest.  Within this work, only the final activation function was investigated.  We did not explore the impact of different intermediate layer and activation function structures.  Optimizing this structure can reduce the size of the realized errors and improve the approximations.  We propose investigating such a hyperparameter optimization for specific classes of problems, e.g., for quadratic vector optimization as analytical structures are most likely to be tractable.

\section{Acknowledgement}
We would like to thank Gabriele Eichfelder for her very helpful comments on an earlier version of this paper.

\bibliographystyle{plain}