RAYEN: Imposition of Hard Convex Constraints on Neural Networks
This paper presents RAYEN, a framework to impose hard convex constraints on
the output or latent variable of a neural network. RAYEN guarantees that, for
any input or any weights of the network, the constraints are satisfied at all
times. Compared to other approaches, RAYEN does not perform a
computationally-expensive orthogonal projection step onto the feasible set,
does not rely on soft constraints (which do not guarantee the satisfaction of
the constraints at test time), does not use conservative approximations of the
feasible set, and does not perform a potentially slow inner gradient descent
correction to enforce the constraints. RAYEN supports any combination of
linear, convex quadratic, second-order cone (SOC), and linear matrix inequality
(LMI) constraints, achieving a very small computational overhead compared to
unconstrained networks. For example, it is able to impose 1K quadratic
constraints on a 1K-dimensional variable with an overhead of less than 8 ms,
and an LMI constraint with 300x300 dense matrices on a 10K-dimensional variable
in less than 12 ms. When used in neural networks that approximate the solution
of constrained optimization problems, RAYEN achieves computation times between
20 and 7468 times faster than state-of-the-art algorithms, while guaranteeing
the satisfaction of the constraints at all times and obtaining a cost very
close to the optimal one.
Jesus Tordesillas, Jonathan P. How, Marco Hutter
# RAYEN: Imposition of Hard Convex Constraints on Neural Networks

###### Abstract

This paper presents RAYEN, a framework to impose hard convex constraints on the output or latent variable of a neural network. RAYEN guarantees that, for any input or any weights of the network, the constraints are satisfied at all times. Compared to other approaches, RAYEN does not perform a computationally-expensive orthogonal projection step onto the feasible set, does not rely on soft constraints (which do not guarantee the satisfaction of the constraints at test time), does not use conservative approximations of the feasible set, and does not perform a potentially slow inner gradient descent correction to enforce the constraints. RAYEN supports any combination of linear, convex quadratic, second-order cone (SOC), and linear matrix inequality (LMI) constraints, achieving a very small computational overhead compared to unconstrained networks. For example, it is able to impose 1K quadratic constraints on a 1K-dimensional variable with an overhead of less than 8 ms, and an LMI constraint with \(300\times 300\) dense matrices on a 10K-dimensional variable in less than 12 ms. When used in neural networks that approximate the solution of constrained optimization problems, RAYEN achieves computation times between 20 and 7468 times faster than state-of-the-art algorithms, while guaranteeing the satisfaction of the constraints at all times and obtaining a cost very close to the optimal one.

 Constraints, convex, neural network, optimization, differentiable.

## I Introduction and Related Work

Convex constraints are widely used in areas such as computer vision [1], control [2], trajectory planning [3], supply chain [4], signal processing [5], and economics [6]. In recent years, there has been an extensive use of neural networks in all these applications due to their expressive power. However, their lack of constraint satisfaction guarantees greatly limits their applicability, especially for safety-critical applications.

There are some constraints that can be easily imposed on the output (or latent variable) of a neural network. For instance, box constraints can be imposed using sigmod functions, nonnegative constraints can be enforced using relu functions, simplex constraints can be guaranteed using softmax functions, and some spherical or ellipsoidal constraints can be imposed using normalization. However, these methods are not applicable to more general types of constraints.

One common way to bias the network towards the satisfaction of the constraints is via soft constraints, which consist of the addition of terms in the training loss to penalize the violation of the constraints [7]. Similar penalties are also used in physics-informed neural networks [8]. The main disadvantage of this approach is that there are no constraint satisfaction guarantees at test time.

Figure 1: RAYEN applied to a batch of 500 samples with a feasible set (���⃝1) defined by linear, convex quadratic, SOC, and LMI constraints. For each sample in the batch, RAYEN lets the corresponding latent variable of the network be the vector that defines the step to take from an interior point of the feasible set (�⃝2). The length of this vector is then adjusted to ensure that the endpoint lies within the set (�⃝3). For visualization purposes, a section of the set has been removed in the right plots.
When the constraints are homogeneous linear inequality constraints (i.e., \(\mathbf{Ay}\leq\mathbf{0}\)), Frerix et al. [9] leverage the Minkowski-Weyl theorem to guarantee the satisfaction of these linear constraints at all times.1 The main disadvantage of this approach is that it is only applicable to linear inequality constraints. Moreover, this method requires running offline the double description method [10, 11] to obtain the V-representation (vertices and rays) of the polyhedron defined by its H-representation (intersection of half-spaces). Even if done offline, the double description method quickly becomes intractable for high dimensions. Other related methods that focus on linear constraints include [12, 13, 14]. In contrast to these works, which focus only on linear constraints, RAYEN supports any combination of linear, convex quadratic, SOC, and LMI constraints.

Footnote 1: The Minkowski-Weyl theorem can be leveraged for any (convex) polyhedron defined by \(\mathbf{Ay}\leq\mathbf{b}\), but [9] focuses on the case \(\mathbf{Ay}\leq\mathbf{0}\)

There have also been many recent advances in implicit layers that solve different types of optimization problems [15, 16, 17, 18], with some works focusing on quadratic programming [19] or, more generally, convex optimization problems [20]. These implicit layers can be leveraged to enforce constraints on the network during training and/or testing by obtaining the orthogonal projection onto the feasible set. This orthogonal projection could also be obtained using the Dykstra's Projection Algorithm [21, 22, 23]. While the use of these orthogonal projections guarantees the satisfaction of the constraints, it is typically at the expense of very high computation times. By leveraging analytic expressions of the distance to the boundaries of the convex set, RAYEN is able to avoid this slow orthogonal projection step and guarantee the constraints in a much more computationally-efficient manner.

With the goal of reducing the computation time, Donti et al. proposed DC3 [24], an algorithm that first uses completion to guarantee the equality constraints, and then enforces the inequality constraints by using an inner gradient decent procedure that takes steps along the manifold defined by the equality constraints. This inner gradient descent procedure is performed both at training and testing time. While DC3 is typically less computationally expensive that projection-based methods, this inequality correction may still require many steps, as we will show in Section V-A. Moreover, this inner gradient correction may also suffer from convergence issues for general convex constraints. Compared to DC3, RAYEN does not need to rely on an inner gradient descent correction, avoiding therefore any convergence issues and substantially reducing the computation time.

The contributions of this work are therefore summarized as follows:

* Framework to impose by construction hard convex constraints on the output or latent variable of a neural network. The constraints are guaranteed to be satisfied at all times, for any input and/or weights of the network. No conservative approximations of the set are used.
* Any combination of linear, convex quadratic, SOC, and LMI constraints is supported. For example, RAYEN can impose 1K quadratic constraints on a 1K-dimensional variable with a computation overhead of only 8 ms. Similarly, a \(300\times 300\) dense LMI constraint can be imposed on a 10K-dimensional variable with an overhead of less than 11 ms.
* When used in neural networks that approximate the solution of optimization problems, RAYEN showcases computation times between 20 and 7468 times faster than other state-of-the-art algorithms, while generating feasible solutions whose costs are very close to the optimal value.

The notation used throughout the paper is available in Table I.

## II Problem Setup

This work addresses the problem of how to ensure that the output (or a latent variable) \(\mathbf{y}:=\left[\begin{smallmatrix}\mathbf{y}_{[0]}&\cdots&\mathbf{y}_{[k-1]}\end{smallmatrix} \right]^{T}\in\mathbb{R}^{k}\) of a neural network lies in a convex set \(\mathcal{Y}\) defined by the following constraints:

\[\mathbf{A}_{1}\mathbf{y}\leq\mathbf{b}_{1} \tag{1}\] \[\mathbf{A}_{2}\mathbf{y}=\mathbf{b}_{2}\] (2) \[g_{i}\left(\mathbf{y}\right):=\frac{1}{2}\mathbf{y}^{T}\mathbf{P}_{i}\mathbf{y }+\mathbf{q}_{i}^{T}\mathbf{y}+r_{i}\leq 0\;\;i=0,\ldots,\eta-1\] (3) \[h_{j}\left(\mathbf{y}\right):=\|\mathbf{M}_{j}\mathbf{y}+\mathbf{s}_{j}\|-\mathbf{c} _{j}^{T}\mathbf{y}-d_{j}\leq 0\;\;j=0,\ldots,\mu-1\] (4) \[\mathbf{W}(\mathbf{y}):=\mathbf{y}_{[0]}\mathbf{F}_{0}+\ldots+\mathbf{y}_{[k-1]}\mathbf{F }_{k-1}+\mathbf{F}_{k}\succeq\mathbf{0} \tag{5}\]

where \(\mathbf{P}_{i}\succeq\mathbf{0}\;\forall i=0,...,\eta-1\), and where \(\mathbf{F}_{0},...,\mathbf{F}_{k}\) are symmetric matrices. This set \(\mathcal{Y}\), which can be bounded or unbounded, is defined by linear constraints (Eqs. 1 and 2), convex quadratic constraints (Eq. 3), SOC constraints (Eq. 4), and LMI constraints (Eq. 5, also known as semidefinite constraints). Constraints 1, 2, 3, and 4 could also be written as an LMI constraint, but for notational convenience

\begin{table}
\begin{tabular}{|c|l|} \hline
**Symbol** & **Mearing** \\ \hline \hline \(c,c,\mathbf{C},\mathbf{C},\mathbf{C}\) & Scalar, column vector, matrix, and set \\ \hline \(\|c\|\) & Euclidean norm of the vector \(\mathbf{c}\) \\ \hline \(\mathbf{C}\succeq\mathbf{0}\) & \(\mathbf{C}\) is a (symmetric) positive semidefinite matrix \\ \hline \(\mathbf{C}\succeq\mathbf{0}\) & \(\mathbf{C}\) is a (symmetric) positive definite matrix \\ \hline \(\mathbf{I}\) & Identity matrix \\ \hline \(\mathbf{c}\) & \(\mathbf{c}^{-}\equiv\frac{1}{2\pi}\mathbf{c}\) \\ \hline \(\mathbf{C}^{\dagger}\) & Pseudoinverse of a matrix \(\mathbf{C}\) \\ \hline max \((\mathbf{c})\) & Maximum of the elements of the column vector \(\mathbf{c}\) \\ \hline max \((\mathbf{c})\) & Maximum of the elements of the set \(\mathbf{C}\) \\ \hline inf \((\mathbf{C})\) & Infimum of the set \(\mathbf{C}\) \\ \hline \(\mathbf{c}\) & Frontier of the set \(\mathbf{C}\) \\ \hline \(\mathbf{a}\leq\mathbf{b}\), \(\mathbf{a}=\mathbf{b}\) & Element-wise inequality, element-wise equality \\ \hline \(\mathbf{A}\odot\mathbf{B}\) & Element-wise division \\ \hline \multirow{3}{*}{\(\mathbf{0},\mathbf{1}\)} & Matrix or vector of zeros/(means dimensions given by the \\  & context). If the dimensions need to be specified, the subscript \(a\)\(\times\) rows \(\times\) columns will be empty. \\ \hline \(\mathbf{0}\) & Empty set \\ \hline inf \((\mathbf{5})\) & Affine hull of a set \(\mathbf{S}\) (i.e., smallest affine set containing \(\mathbf{S}\)) \\ \hline \(\mathbf{c}_{[i]}\) & \(\mathbf{i}\)-th element of the column vector \(\mathbf{c}\in\mathbb{R}^{n}\). \(i\in\{0,\ldots,a-1\}\). If \(\mathbf{C}\in\mathbb{R}^{n\times n}\), then \(\mathbf{c}_{[i]}\in\mathbb{R}^{n\times n}\) (i.e., \(\mathbf{x}\in\mathbf{0}\) vector). \\ \hline \(\mathbf{C}_{[S,i]}\cdot\mathbf{c}_{[S]}\) & \(\mathbf{C}_{[S,i]}\) is a matrix whose rows are the rows of \(\mathbf{C}\) whose \\  & indices are in the set \(\mathbf{S}\leq\mathbf{1}\). An analogous definition for \(\mathbf{c}_{[S]}\) \\ \hline abs \((\mathbf{a})\) & Element-wise absolute value \\ \hline softmax \((\mathbf{a})\) & Softmax function, in \(\mathbf{x}\) \(\mathbf{b}\) - softmax \((\mathbf{a})\), then \(\mathbf{b}_{[i]}=\mathbf{c}^{-}\equiv\mathbf{c}^{-}\equiv\mathbf{c}^{-}\equiv\mathbf{c}^{-}\equiv \mathbf{c}^{-}\equiv\mathbf{c}^{-}\equiv\mathbf{c}^{-}\equiv\mathbf{c}^{-}\equiv\mathbf{c}^{-}\) \\ \hline cig \((\mathbf{C})\) & Column vector containing all the eigenvalues of the matrix \(\mathbf{C}\) \\ \hline SOC & Second-order Cone \\ \hline LMI & Linear Matrix Inequality \\ \hline \end{tabular}
\end{table} TABLE I: Notation used in this paper
throughout the paper (and without any loss of generality), we explicitly distinguish them. Note also that several LMIs can be converted into one LMI by simply stacking the matrices appropriately (see, e.g., [25, Section 4.6.2]).

Let us also introduce the following definitions:

\[\mathcal{Y}_{L} :=\big{\{}\mathbf{y}\in\mathbb{R}^{k}|\mathbf{A}_{1}\mathbf{y}\leq\mathbf{b}_{1} \big{\}}\cap\big{\{}\mathbf{y}\in\mathbb{R}^{k}|\mathbf{A}_{2}\mathbf{y}=\mathbf{b}_{2}\big{\}}\] \[\mathcal{Y}_{Q} :=\big{\{}\mathbf{y}\in\mathbb{R}^{k}|g_{i}\left(\mathbf{y}\right)\leq 0 \ \forall i=0,\ldots,\eta-1\big{\}}\] \[\mathcal{Y}_{S} :=\big{\{}\mathbf{y}\in\mathbb{R}^{k}|h_{j}\left(\mathbf{y}\right)\leq 0 \ \forall j=0,\ldots,\mu-1\big{\}}\] \[\mathcal{Y}_{M} :=\big{\{}\mathbf{y}\in\mathbb{R}^{k}|\mathbf{W}(\mathbf{y})\succeq\mathbf{0} \big{\}}\]

Hence, we have that \(\mathcal{Y}:=\mathcal{Y}_{L}\cap\mathcal{Y}_{Q}\cap\mathcal{Y}_{S}\cap \mathcal{Y}_{M}\). For simplicity, we will assume throughout the paper that \(\text{aff}\left(\mathcal{Y}_{L}\right)=\text{aff}\left(\mathcal{Y}\right)\), where \(\text{aff}\left(\cdot\right)\) denotes the affine hull. Some degenerate cases do not satisfy this assumption, but in those cases the same set \(\mathcal{Y}\) can be parametrized with a different set of constraints for which this assumption holds.2

Footnote 2: An example of a degenerate case in \(\mathcal{Y}\) (i.e., \(k=3\)) with only linear and quadratic constraints would be when \(\mathcal{Y}_{L}\) is the unit cube, and \(\mathcal{Y}_{Q}\) is a cylinder tangent to one of its faces such that \(\mathcal{Y}=\mathcal{Y}_{L}\cap\mathcal{Y}_{Q}\) is the segment they have in common. By simply reparametrizing \(\mathcal{Y}\) with only linear constraints (the segment itself), then the assumption \(\text{aff}\left(\mathcal{Y}_{L}\right)=\text{aff}\left(\mathcal{Y}\right)\) is satisfied.

Overall, RAYEN works as follows (see also Fig. 2): In an offline phase (Section III), the affine hull of \(\mathcal{Y}\) is obtained, and the constraints are expressed in that linear subspace. Letting \(\mathcal{Z}\) denote the feasible set in this subspace, an interior point \(\mathbf{z}_{0}\) of this set \(\mathcal{Z}\) is also found. In the online phase (Section IV), a linear map is applied to an upstream latent variable of the network to obtain a vector \(\mathbf{v}\) of the same dimension of the subspace, that will determine the step direction from \(\mathbf{z}_{0}\). The length \(\lambda\geq 0\) of this step is adjusted to ensure that \(\mathbf{z}_{0}+\lambda\frac{\mathbf{v}}{\|\mathbf{v}\|}\) lies in \(\mathcal{Z}\). After lifting to the original dimension, the output is therefore guaranteed to lie in \(\mathcal{Y}\).

## III RAYEN: Offline

### _Affine Hull of \(\mathcal{Y}\)_

To find the affine hull of the set \(\mathcal{Y}\), we first express all the linear constraints (Eqs. 1 and 2) as linear inequality constraints by simply stacking the matrices:

\[\tilde{\mathbf{A}} :=\big{[}\,\mathbf{A}_{1}^{T}\,\mathbf{A}_{2}^{T}-\mathbf{A}_{2}^{T}\,\big{]} ^{T}\] \[\tilde{\mathbf{b}} :=\big{[}\,b_{1}^{T}\,\mathbf{b}_{2}^{T}-\mathbf{b}_{2}^{T}\,\big{]}^{T}\]

and therefore \(\mathcal{Y}_{L}\) is now \(\Big{\{}\mathbf{y}\in\mathbb{R}^{k}|\tilde{\mathbf{A}}\mathbf{y}\leq\tilde{\mathbf{b}}\Big{\}}\). To reduce the computation time during the online phase, the redundant constraints of the inequality system \(\tilde{\mathbf{A}}\mathbf{y}\leq\tilde{\mathbf{b}}\) (if there are any) are then deleted. This is done by solving a sequence of Linear Programs (see, e.g. [26, Eq. 1.5]), which obtains the reduced system of inequalities \(\mathbf{A}\mathbf{y}\leq\mathbf{b}\) such that
\(\mathcal{Y}_{L}=\left\{\mathbf{y}\in\mathbb{R}^{k}|\mathbf{A}\mathbf{y}\leq\mathbf{b}\right\}\). Then, another sequence of Linear Programs is solved to find the set

\[\mathcal{E}:=\left\{i|\mathbf{A}_{[i,:]}\mathbf{y}=\mathbf{b}_{[i]}\ \forall\mathbf{y}\in \mathcal{Y}_{L}\right\}\,\]

which contains the indexes of the constraints that are always active for all the points in \(\mathcal{Y}_{L}\) (see, e.g. [27, Alg. 1.3.1] or [28, Section 5.2]). Defining \(\mathcal{I}:=\{0,...,a-1\}\backslash\mathcal{E}\) (where \(a\) is the number of rows of \(\mathbf{A}\)), let us now introduce the following notation:3

Footnote 3: Note that if \(\mathcal{I}=\emptyset\), then we need \(\mathbf{b}_{\mathcal{I}}=1\) (instead of \(\mathbf{b}_{\mathcal{I}}=0\)) to make sure that there exists a point in the interior of \(\mathcal{Z}_{L}:=\left\{\mathbf{z}\in\mathbb{R}^{n}|\mathbf{0}_{1\times n}\mathbf{z}\leq 1 \right\}\equiv\mathbb{R}^{n}\) (see Section III-C).

\[\mathbf{A}_{\mathcal{I}}:=\begin{cases}\mathbf{0}_{1\times k}&\mathcal{I}=\emptyset\\ \mathbf{A}_{[\mathcal{I},:]}&\text{otherwise}\end{cases},\ \mathbf{b}_{ \mathcal{I}}:=\begin{cases}1&\mathcal{I}=\emptyset\\ \mathbf{b}_{[\mathcal{I}]}&\text{otherwise}\end{cases}\]

Then, the affine hull of \(\mathcal{Y}\) is given by

\[\text{aff}\left(\mathcal{Y}\right)=\left\{\mathbf{y}\in\mathbb{R}^{k}|\mathbf{A}_{ \mathcal{E}}\mathbf{y}=\mathbf{b}_{\mathcal{E}}\right\}\,\]

where we have used the assumption that \(\text{aff}\left(\mathcal{Y}\right)=\text{aff}\left(\mathcal{Y}_{L}\right)\) (see Section II). The dimension of this affine hull is then

\[n:=k-\text{rank}\left(\mathbf{A}_{\mathcal{E}}\right). \tag{6}\]

Note also that

\[\mathcal{Y}_{L}=\underbrace{\left\{\mathbf{y}\in\mathbb{R}^{k}|\mathbf{A}_{\mathcal{E }}\mathbf{y}=\mathbf{b}_{\mathcal{E}}\right\}}_{=\text{aff}\left(\mathcal{Y}_{L} \right)=\text{aff}\left(\mathcal{Y}\right)}\cap\left\{\mathbf{y}\in\mathbb{R}^{k}| \mathbf{A}_{\mathcal{I}}\mathbf{y}\leq\mathbf{b}_{\mathcal{I}}\right\}\.\]

In general, \((\mathbf{A}_{\mathcal{E}},\mathbf{b}_{\mathcal{E}})\) may be different than \((\mathbf{A}_{\mathcal{I}},\mathbf{b}_{2})\), and \((\mathbf{A}_{\mathcal{I}},\mathbf{b}_{\mathcal{I}})\) may be different than \((\mathbf{A}_{\mathcal{I}},\mathbf{b}_{1})\). This can happen, for example, when \(\text{aff}\left(\left\{\mathbf{y}\in\mathbb{R}^{k}|\mathbf{A}_{1}\mathbf{y}\leq\mathbf{b}_{1} \right\}\right)\neq\mathbb{R}^{k}\) (intuitively this means that the inequality constraint defined by Eq. 1_hides_ equality constraints), or when there are redundant constraints in Eqs. 1 and/or 2.

### _Set \(\mathcal{Z}\)_

Any point in \(\mathbf{y}\in\text{aff}\left(\mathcal{Y}\right)\) can be expressed as

\[\mathbf{y}=\mathbf{f}(\mathbf{z}):=\mathbf{N}\mathbf{z}+\underbrace{\mathbf{A}_{\mathcal{E}}^{\dagger }\mathbf{b}_{\mathcal{E}}}_{:=\mathbf{y}_{p}}\, \tag{7}\]

where \(\mathbf{z}\in\mathbb{R}^{n}\), \(\mathbf{N}\) is a \(k\times n\) matrix whose columns form an orthonormal basis for the null space of \(\mathbf{A}_{\mathcal{E}}\), and \(\mathbf{A}_{\mathcal{E}}^{\dagger}\) is the pseudoinverse of \(\mathbf{A}_{\mathcal{E}}\). Using Eq. 7, we have that \(\mathbf{A}_{\mathcal{I}}\mathbf{y}\leq\mathbf{b}_{\mathcal{I}}\) is equivalent to:

\[\mathbf{A}_{\mathcal{I}}\left(\mathbf{N}\mathbf{z}+\mathbf{y}_{p}\right)\leq\mathbf{b}_{ \mathcal{I}}\] \[\underbrace{\mathbf{A}_{\mathcal{I}}\mathbf{N}\mathbf{z}}_{:=\mathbf{A}_{p}}\leq \underbrace{\mathbf{b}_{\mathcal{I}}-\mathbf{A}_{\mathcal{I}}\mathbf{y}_{p}}_{:=\mathbf{b}_{p}} \tag{8}\]

Defining now these sets (see also Fig. 2):

\[\mathcal{Z}_{L}:= \left\{\mathbf{z}\in\mathbb{R}^{n}|\mathbf{A}_{p}\mathbf{z}\leq\mathbf{b}_{p}\right\}\] \[\mathcal{Z}_{Q_{i}}:= \left\{\mathbf{z}\in\mathbb{R}^{n}|g_{i}\left(\mathbf{f}(\mathbf{z})\right) \leq 0\right\}\qquad\mathcal{Z}_{Q}:=\bigcap_{i=0}^{\eta-1}\mathcal{Z}_{Q_{i}}\] \[\mathcal{Z}_{S_{j}}:= \left\{\mathbf{z}\in\mathbb{R}^{n}|h_{j}\left(\mathbf{f}(\mathbf{z})\right) \leq 0\right\}\qquad\mathcal{Z}_{S}:=\bigcap_{j=0}^{\mu-1}\mathcal{Z}_{S_{j}}\] \[\mathcal{Z}_{M}:= \left\{\mathbf{z}\in\mathbb{R}^{n}|\mathbf{W}\left(\mathbf{f}(\mathbf{z})\right) \succeq\mathbf{0}\right\}\] \[\mathcal{Z}:= \mathcal{Z}_{L}\cap\mathcal{Z}_{Q}\cap\mathcal{Z}_{S}\cap \mathcal{Z}_{M}\]

Any point \(\mathbf{z}\in\mathcal{Z}\) can therefore be mapped to the corresponding point \(\mathbf{y}\in\mathcal{Y}\) using Eq. 7.

### _Interior point of \(\mathcal{Z}\)_

Finally, we obtain a point \(\mathbf{z}_{0}\) in the interior of \(\mathcal{Z}\), by solving the following convex program:

\[(\mathbf{z}_{0},\epsilon^{*})=\underset{\mathbf{z},\epsilon}{\text{argmax}} \epsilon\] s.t.: \[\mathbf{A}_{p}\mathbf{z}-\mathbf{b}_{p}\leq-\epsilon\mathbf{1}\] \[g_{i}\left(\mathbf{f}\left(\mathbf{z}\right)\right)\leq-\epsilon\ \ i=0,...,\eta-1\] \[h_{j}\left(\mathbf{f}\left(\mathbf{z}\right)\right)\leq-\epsilon\ \ j=0,...,\mu-1\] \[\mathbf{W}\left(\mathbf{f}(\mathbf{z})\right)\succeq\epsilon\mathbf{I}\] \[\epsilon>0\]

## IV RAYEN: Online

As shown in Fig. 3, to obtain an \(n\)-dimensional latent variable \(\mathbf{v}\) (where \(n\) is the dimension of \(\text{aff}\left(\mathcal{Y}\right)\), see Eq. 6), we first apply a linear layer to the upstream latent variable of the network \(\mathbf{x}\in\mathbb{R}^{m}\). This mapping can be omitted if \(m=n\). Assuming now \(\mathbf{v}\neq\mathbf{0}\),4 let us define the inverse distance to the frontier of \(\mathcal{Z}\) along \(\bar{\mathbf{v}}:=\frac{\mathbf{v}}{\|\mathbf{v}\|}\) as

Footnote 4: If \(\mathbf{v}=\mathbf{0}\), we simply take \(\mathbf{z}_{1}:=\mathbf{z}_{0}\), and therefore \(\kappa\) does not need to be computed.

\[\kappa:=\text{inf}\left\{\lambda^{-1}|\mathbf{z}_{0}+\lambda\bar{\mathbf{v}}\in\mathcal{Z },\ \lambda>0\right\}\.\]

Then, it is clear that

\[\mathbf{z}_{1}:=\left(\mathbf{z}_{0}+\text{min}\left(\frac{1}{\kappa},\|\mathbf{v}\|\right) \bar{\mathbf{v}}\right)\in\mathcal{Z}\, \tag{9}\]

Figure 3: Neural network equipped with RAYEN. \(\kappa\) is computed as shown in Table II, and the map \(\mathbf{f}\left(\mathbf{z}_{0}+\text{min}\left(\frac{1}{\kappa},\|\mathbf{v}\|\right)\bar{ \mathbf{v}}\right)\) is used to guarantee \(\mathbf{y}\in\mathcal{Y}\). Here, L(\(m,n\)) denotes a linear layer with input size \(m\) and output size \(n\). This linear layer is not needed if \(m=n\). After the RAYEN module, there may be more downstream layers. During the training procedure, gradients can then be backpropagated through RAYEN.
and that therefore \(\mathbf{f}\left(\mathbf{z}_{1}\right)\in\mathcal{Y}\). Note that no conservatism is introduced here. In other words, and given that \(\mathcal{Y}\) is a convex set, for any point \(\mathbf{y}\in\mathcal{Y}\), there exists at least one \(\mathbf{v}\) such that \(\mathbf{f}\left(\mathbf{z}_{1}\right)\in\mathcal{Y}\) holds. To obtain \(\kappa\), first note that

\[\kappa=\max\left[\,\kappa_{L}\;\kappa_{Q}\;\kappa_{S}\;\kappa_{M}\,\right]^{T}\,\]

where \(\kappa_{L}\), \(\kappa_{Q}\), \(\kappa_{S}\), and \(\kappa_{M}\) are defined in Table II and are, respectively, the inverses of the distances to \(\partial\mathcal{Z}_{L}\), \(\partial\mathcal{Z}_{Q}\), \(\partial\mathcal{Z}_{S}\), and \(\partial\mathcal{Z}_{M}\) along \(\bar{\mathbf{v}}\) (see also Figs. 2 and 4). In the following subsections we explain how \(\kappa_{L}\), \(\kappa_{Q}\), \(\kappa_{S}\), and \(\kappa_{M}\) can be obtained.

### _Computation of \(\kappa_{L}\)_

Given a face \(\mathcal{F}_{j}:=\left\{\mathbf{z}\in\mathbb{R}^{n}|\mathbf{A}_{p_{[j,:]}}\mathbf{z}=\mathbf{b }_{p_{[j]}}\right\}\) of the polyhedron \(\mathcal{Z}_{L}\), the inverse distance from \(\mathbf{z}_{0}\) to \(\mathcal{F}_{j}\) along the direction \(\bar{\mathbf{v}}\) can be easily obtained as

\[\kappa_{L,j}=\text{relu}\left(\frac{\mathbf{A}_{p_{[j,:]}}\bar{\mathbf{v}}}{\mathbf{b}_{p _{[j]}}-\mathbf{A}_{p_{[j,:]}}\mathbf{z}_{0}}\right)\,\]

where the \(\text{relu}(\cdot)\) operator is needed to enforce moving along \(\bar{\mathbf{v}}\). Taking into account now all the faces of \(\mathcal{Z}_{L}\), we have that \(\kappa_{L}\) is given by

\[\kappa_{L}:=\text{relu}\left(\max\left(\underbrace{\mathbf{A}_{p}\oslash((\mathbf{b}_ {p}-\mathbf{A}_{p}\mathbf{z}_{0})\,\mathbf{1}_{1\times n})}_{:=\mathcal{D}}\bar{\mathbf{v}} \right)\right)\,\]

where \(\oslash\) denotes the element-wise division between two matrices, and where \(\mathbf{D}\) can be computed offline because it does not depend on \(\bar{\mathbf{v}}\). Note that if \(\kappa_{L}=0\), then \(\mathcal{Z}_{L}\) is unbounded along the ray that starts at \(\mathbf{z}_{0}\) and follows \(\bar{\mathbf{v}}\).

### _Computation of \(\kappa_{Q}\)_

The inverse of the distance from \(\mathbf{z}_{0}\) to \(\partial\mathcal{Z}_{Q_{i}}\) along the direction \(\bar{\mathbf{v}}\) can be obtained by computing the nonnegative \(\kappa_{Q,i}\) that satisfies:

\[g_{i}\left(\mathbf{f}\left(\mathbf{z}_{0}+\frac{1}{\kappa_{Q,i}}\bar{\mathbf{v}}\right) \right)=0\,\]

where \(\mathbf{f}(\cdot)\) is defined in Eq. 7. Multiplying both sides of this equation by \(\kappa_{Q,i}^{2}\) yields a quadratic equation on \(\kappa_{Q,i}\), from which its nonnegative root \(\kappa_{Q,i}^{+}\) can be easily obtained. Taking now into account all the quadratic constraints, \(\kappa_{Q}\) is therefore given by:

\[\kappa_{Q}:=\max\left(\left[\,\kappa_{Q,0}^{+}\cdots\kappa_{Q,\eta-1}^{+} \right]^{T}\right)\]
### _Computation of \(\kappa_{S}\)_

Similar to the previous case, the inverse of the distance from \(\mathbf{z}_{0}\) to \(\partial\mathcal{Z}_{S_{j}}\) along the direction \(\bar{\mathbf{v}}\) can be obtained by computing the nonnegative \(\kappa_{S,j}\) that satisfies:

\[h_{j}\left(\mathbf{f}\left(\mathbf{z}_{0}+\frac{1}{\kappa_{S,j}}\bar{\mathbf{v}}\right) \right)=0\]

\[\left\|\mathbf{M}_{j}\mathbf{f}\left(\mathbf{z}_{0}+\frac{1}{\kappa_{S,j}}\bar{\mathbf{v}} \right)+\mathbf{s}_{j}\right\|=\mathbf{c}_{j}^{T}\mathbf{f}\left(\mathbf{z}_{0}+\frac{1}{ \kappa_{S,j}}\bar{\mathbf{v}}\right)+d_{j} \tag{10}\]

Squaring both sides of Eq. 10, and then multiplying them by \(\kappa_{S,j}^{2}\) yields a quadratic equation on \(\kappa_{S,j}\), from which its two roots \(\kappa_{S,j}^{(0)}\) and \(\kappa_{S,j}^{(1)}\) can be easily obtained. The inverse of the distance from \(z_{0}\) to \(\partial\mathcal{Z}_{S_{j}}\) along \(\bar{\mathbf{v}}\) will therefore be given by5

Footnote 5: The \(\text{relu}(\cdot)\) operator is needed because both roots can be negative due to the fact that we have squared both sides of Eq. 10. Both roots being negative means that in that case \(\mathcal{Z}_{S_{j}}\) is unbounded in the direction \(\bar{\mathbf{v}}\) and therefore we have \(\kappa_{S,j}=0\). When both roots are positive, we need to select the largest one, which is why the \(\max\left(\cdot\right)\) operator is needed.

\[\text{relu}\left(\max\left(\left[\begin{array}{cc}\kappa_{S,j}^{(0)}& \kappa_{S,j}^{(1)}\end{array}\right]^{T}\right)\right)\.\]

Taking into account all the SOC constraints, \(\kappa_{S}\) is therefore given by

\[\kappa_{S}:=\text{relu}\left(\max\left(\left[\begin{array}{cc}\kappa_{S,0}^ {(0)}&\kappa_{S,0}^{(1)}\cdots\kappa_{S,\mu-1}^{(0)}&\kappa_{S,\mu-1}^{(1)} \end{array}\right]^{T}\right)\right)\.\]

### _Computation of \(\kappa_{M}\)_

The inverse of the distance from \(\mathbf{z}_{0}\) to \(\partial\mathcal{Z}_{M}\) along the direction \(\bar{\mathbf{v}}\) is defined as:

\[\kappa_{M}:=\text{inf}\left\{\lambda^{-1}|\mathbf{z}_{0}+\lambda\bar{\mathbf{v}}\in \mathcal{Z}_{M},\ \lambda>0\right\}\]

Noting that

\[\mathbf{f}\left(\mathbf{z}_{0}+\lambda\bar{\mathbf{v}}\right)=\lambda\underbrace{\mathbf{N} \bar{\mathbf{v}}}_{:=\mathbf{p}}+\underbrace{\mathbf{N}\mathbf{z}_{0}+\mathbf{y}_{p}}_{:=\mathbf{y}_{0 }}=\left(\lambda\mathbf{\rho}+\mathbf{y}_{0}\right)\,\]

then the condition \(\mathbf{z}_{0}+\lambda\bar{\mathbf{v}}\in\mathcal{Z}_{M}\) is equivalent to

\[\mathbf{W}\left(\mathbf{f}\left(\mathbf{z}_{0}+\lambda\bar{\mathbf{v}}\right)\right)=\mathbf{F}_{ k}+\sum_{\alpha=0}^{k-1}\left(\lambda\mathbf{\rho}_{[\alpha]}+\mathbf{y}_{0_{[\alpha]}} \right)\mathbf{F}_{\alpha}\succeq\mathbf{0}\.\]

Dividing everything by \(\lambda\) (recall that \(\lambda>0\)) and rearranging the terms we have

\[\lambda^{-1}\underbrace{\left(\mathbf{F}_{k}+\sum_{\alpha=0}^{k-1}\left(\mathbf{y}_{ 0_{[\alpha]}}\mathbf{F}_{\alpha}\right)\right)}_{:=\mathbf{H}\succ\mathbf{0}}+\underbrace {\sum_{\alpha=0}^{k-1}\left(\mathbf{\rho}_{[\alpha]}\mathbf{F}_{\alpha}\right)}_{:=\mathbf{ S}}\succeq\mathbf{0}\,\]

where \(\mathbf{H}\succ\mathbf{0}\) derives from the fact that \(\mathbf{z}_{0}\) is an interior point of \(\mathcal{Z}_{M}\). Hence:

\[\kappa_{M}:=\text{inf}\left\{\delta|\delta\mathbf{H}+\mathbf{S}\succeq\mathbf{0},\ \delta>0\right\} \tag{11}\]

Let us now distinguish two cases:

* If \(\mathbf{S}\succeq\mathbf{0}\), then it is clear that \(\kappa_{M}=0\).
* Otherwise, we need \(\kappa_{M}\mathbf{H}+\mathbf{S}\in\partial\mathcal{Z}_{M}\), and hence \[\kappa_{M}=\text{max}\left(\left\{\tau_{i}|\left(\tau_{i}\mathbf{H}+\mathbf{S}\right) \mathbf{v}_{i}=0\mathbf{v}_{i},\mathbf{v}_{i}\neq\mathbf{0}\right\}\right)\,\]

where we have used the fact that the matrices that belong to \(\partial\mathcal{Z}_{M}\) have at least one zero eigenvalue (see also Fig. 5). Note also that: \[\left(\tau_{i}\mathbf{H}+\mathbf{S}\right)\mathbf{v}_{i}=0\mathbf{v}_{i}\iff-\mathbf{H}^{-1}\mathbf{S }\mathbf{v}_{i}=\tau_{i}\mathbf{v}_{i}\] (12) Letting \(\mathbf{L}\) denote the lower triangular matrix with positive diagonal entries such that \(\mathbf{H}^{-1}=\mathbf{L}\mathbf{L}^{T}\) (i.e., the Cholesky decomposition of \(\mathbf{H}^{-1}\)), and using the fact that \(\text{eig}\left(-\mathbf{H}^{-1}\mathbf{S}\right)=\text{eig}\left(\mathbf{L}^{T}\left(-\mathbf{S }\right)\mathbf{L}\right)\) (see proof in Appendix A), we have that6 Footnote 6: Eq. 12 can also be written as \(-\mathbf{S}\mathbf{v}_{i}=\tau_{i}\mathbf{H}\mathbf{v}_{i}\), and hence \(\kappa_{M}\) could also be found solving the generalized eigenvalue problem [25, 29] of the pair \((-\mathbf{S},\mathbf{H})\).

\[\kappa_{M}=\text{relu}\left(\max\left(\text{eig}\left(\mathbf{L}^{T}\left(-\mathbf{S} \right)\mathbf{L}\right)\right)\right). \tag{13}\]

As only the maximum eigenvalue is needed, we can use methods such as power iteration [31, 32, 33], LOBPCG [34], or the Lanczos algorithm [35, 36] to avoid computing the whole spectrum of the matrix. Note also that the matrix \(\mathbf{L}\) can be computed offline, since it does not depend on \(\mathbf{v}\).

### _Remarks_

We can have these two cases in Eq. 9:

* If \(\mathbf{z}_{0}+\mathbf{v}\in\mathcal{Z}\), then \(\left\|\mathbf{v}\right\|\leq\frac{1}{\kappa}\), and therefore Eq. 9 becomes \(\mathbf{z}_{1}=\mathbf{z}_{0}+\mathbf{v}\).
* If \(\mathbf{z}_{0}+\mathbf{v}\notin\mathcal{Z}\), then \(\left\|\mathbf{v}\right\|>\frac{1}{\kappa}\) and therefore Eq. 9 becomes \(\mathbf{z}_{1}=\mathbf{z}_{0}+\frac{1}{\kappa}\bar{\mathbf{v}}\). This operation performs a (nonorthogonal) projection of \(\mathbf{z}_{0}+\mathbf{v}\) onto \(\mathcal{Z}\) along the ray that starts at \(\mathbf{z}_{0}\) and follows \(\bar{\mathbf{v}}\).

Note also that, if \(\mathcal{Z}\) (\(\mathcal{Z}_{L}\), \(\mathcal{Z}_{Q}\), \(\mathcal{Z}_{S}\), \(\mathcal{Z}_{M}\) respectively) is unbounded along the ray that starts at \(\mathbf{z}_{0}\) and follows \(\bar{\mathbf{v}}\), then
\(\kappa\) (\(\kappa_{L}\), \(\kappa_{Q}\), \(\kappa_{S}\), \(\kappa_{M}\) respectively) will be zero. If \(\kappa=0\), then we have that \(\frac{1}{\kappa}=\infty\) and therefore \(\mathbf{z}_{1}=\mathbf{z}_{0}+\mathbf{v}\).

Fig. 6 shows the results of RAYEN applied to 12 different sets \(\mathcal{Y}\). For each set, we generate 12K samples \(\mathbf{v}\in\mathbb{R}^{n}\) sampled uniformly in the box \([-2.5,2.5]^{n}\) (see Fig. 3), compute the corresponding \(\kappa\), and plot the resulting point \(\mathbf{y}=\mathbf{f}\left(\mathbf{z}_{1}\right)\in\mathbb{R}^{k}\) for each of those samples. All these points \(\mathbf{y}\) are guaranteed to lie in the set \(\mathcal{Y}\). Note also that the density of the produced points is higher in the frontiers of the set due to the \(\text{min}\left(\frac{1}{\kappa},\left\|\mathbf{v}\right\|\right)\) operation performed.

## V Results

In this Section, we first apply RAYEN to approximate the optimal solution of a constrained optimization problem (Section V-A), and then analyze its computation time when applied to different constraints of varying dimensions (Section V-B). All the results were obtained using a desktop equipped with an Intel Core i9-12900 \(\times\) 24, 64GB of RAM, and an NVIDIA GeForce RTX 4080. The double-precision floating-point format float64 was used for all the results.

### _Optimization problems_

One of the (many) applications of being able to impose constraints on neural networks is when they are used to approximately solve optimization problems much faster than standard solvers. In this section, we analyze how a neural network equipped with RAYEN is able to approximate the optimal solution of an optimization problem while satisfying all the constraints. We use the optimization problems shown in Table IV: one with only linear constraints, and another one with both linear and convex quadratic constraints. Both are convex optimization problems in which the objective function depends on a parameter \(\mathbf{\gamma}\). Nonconvex objective functions could also be used, but we use convex objective functions to enable a direct comparison with the globally optimal solution obtained using a state-of-the-art convex optimization solver,

Figure 6: RAYEN applied to 12 different constraints. For each of the constraints, we generate 12K samples \(\mathbf{v}\in\mathbb{R}^{n}\) sampled uniformly in the box \([-2.5,2.5]^{n}\), compute the corresponding \(\kappa\), and plot the resulting \(\mathbf{y}=\mathbf{f}\left(\mathbf{z}_{0}+\text{min}\left(\frac{1}{\kappa},\left\|\mathbf{v} \right\|\right)\bar{\mathbf{v}}\right)\in\mathbb{R}^{k}\) for each of those samples (see also Fig. 3). For each of the subfigures, the left plot shows the feasible set \(\mathcal{Y}\), while the right plot shows both \(\mathcal{Y}\) and the output samples \(\mathbf{y}\). Note how all these output samples satisfy the constraints. The constraints of the first row have \(k=n=2\). The first two constraints of the second row have \(k=3\) and \(n=2\). The rest of the constraints have \(k=n=3\). RAYEN works in any dimension, but, for visualization purposes, only 2D and 3D examples are shown here.
such as Gurobi [37]. A detailed description of these problems is available in Appendix B.

We compare RAYEN with the following approaches (see also Table III):

* **UU:** Both training and testing are Unconstrained. This method produces directly \(\mathbf{y}=\mathbf{v}\). A soft cost, described below, is added to the training loss to penalize the violation of the constraints.
* **UP:** When training, this method simply outputs \(\mathbf{f}(\mathbf{v})\) (see Eq. 7), meaning that \(\mathbf{y}\) is Unconstrained with respect to the inequality constraints. When testing, this method outputs \(\mathbf{f}(\mathbf{v}^{\prime})\), where \(\mathbf{v}^{\prime}\) is the orthogonal Projection of \(\mathbf{v}\) onto \(\mathcal{Z}\). The projection is computed using [20]. This method also has a soft cost in the training loss to penalize the violation of the constraints.
* **PP:** The orthogonal Projection onto \(\mathcal{Z}\) is performed during both training and testing. These projections are implemented using [20].
* **DC3:** Algorithm proposed by [24]. Completion is used to enforce the equality constraints, and then an inner gradient decent method (which is also performed in the testing phase) enforces the inequality constraints. This method also includes a soft cost in the training loss to penalize the violation of the constraints.
* **Bar:** Generalization of the barycentric coordinates method [9] (proposed for constraints \(\mathbf{A}\mathbf{y}\leq\mathbf{0}\)) for any constraint of the type \(\mathbf{A}\mathbf{y}\leq\mathbf{b}\). This method only supports linear constraints. The matrices of vertices \(\mathbf{V}\) and rays \(\mathbf{R}\) (see Table III) are computed offline using the double description method [10, 11] on the H-representation (half-space representation) of \(\mathcal{Y}_{L}\). As detailed in Table III, the \(\text{softmax}(\cdot)\) operator produces the weights for the convex combination of the vertices, while the \(\text{abs}(\cdot)\) operator produces the weights for the conical combination of the rays.

Defining the following soft costs:

\[p_{\text{soft},L} := \left\|\text{relu}\left(\mathbf{A}_{1}\mathbf{y}-\mathbf{b}_{1}\right) \right\|^{2}+\left\|\mathbf{A}_{2}\mathbf{y}-\mathbf{b}_{2}\right\|^{2} \tag{14}\] \[p_{\text{soft},Q} := \sum_{i=0}^{\eta-1}\left\|\text{relu}\left(g_{i}\left(\mathbf{y} \right)\right)\right\|^{2} \tag{15}\]

the training and test losses in each algorithm are defined in Table IV. In these losses, the soft costs are weighted with
\(\omega\geq 0\). The methods **Bar**, **PP**, and **RAYEN** do not need a soft cost since all the constraints are guaranteed to be satisfied both during training and testing. Note also that in **DC3**, the term \(\left\|\mathbf{A}_{2}\mathbf{y}-\mathbf{b}_{2}\right\|^{2}\) of \(p_{\text{soft},L}\) will always be zero since this algorithm satisfies the equality constraints by construction.

All the methods are implemented in Pytorch[38], and the Adam optimizer [39] with a learning rate of \(10^{-4}\) is used for training. A total of 2000 epochs are performed for training, and the policy with the best validation loss is selected for testing. For each optimization problem, a total of 1216 samples of \(\mathbf{\gamma}\in\mathbb{R}^{n_{\mathbf{\gamma}}}\) are drawn from uniform distributions detailed in Appendix C, and \(\approx 70\%\) of them are used in the training set and \(\approx 30\%\) in the validation set. The batch size for training is 256.

For each optimization problem \(j\in\{1,2\}\), we use two testing sets (each one with 512 samples \(\mathbf{\gamma}\in\mathbb{R}^{n_{\mathbf{\gamma}}}\)): \(\text{TE}_{j,\text{in}}\) and \(\text{TE}_{j,\text{out}}\). \(\text{TE}_{j,\text{in}}\) uses the same distribution as the one used for training, while \(\text{TE}_{j,\text{out}}\) does not. Details on these distributions are available in Appendix C. The goal of using the test set \(\text{TE}_{j,\text{out}}\) is to study how well the trained policy generalizes

\begin{table}
\begin{tabular}{|c|c|c|c|c|c|c|c|} \hline \multirow{2}{*}{**Method**} & \multirow{2}{*}{\begin{tabular}{c} **Num. of** \\ **params** \\ \end{tabular} } & \multicolumn{6}{c|}{**Optimization 1 (Linear Constraints)**} \\ \cline{3-8}  & & \multicolumn{2}{c|}{**Same distribution as training: \(\mathbf{\gamma}\in\text{TE}_{1,\text{in}}\)**} & \multicolumn{2}{c|}{**Different distribution as training: \(\mathbf{\gamma}\in\text{TE}_{1,\text{out}}\)**} \\ \cline{3-8}  & & **Normalized Loss** & **Violation** & **Inference Time (\(\mu s\))** & **Normalized Loss** & **Violation** & **Inference Time (\(\mu s\))** \\ \hline \multirow{2}{*}{\begin{tabular}{c} **UU**, \(\omega=0\) \\ **UU**, \(\omega=100\) \\ \end{tabular} } & 0.0006 & 1945.58 & 0.34 & 0.0013 & 1937.43 & 0.43 \\  & \(\omega=100\) & 2.4431 & 2.68 & 0.36 & 1.8362 & 14.02 & 0.43 \\  & \(\omega=100\) & 9872 & 1.0278 & 0.00418 & 0.36 & 1.0220 & 0.08535 & 0.43 \\  & \(\omega=1000\) & 1.4391 & 0.00255 & 0.36 & 1.2263 & 0.19527 & 0.44 \\  & \(\omega=5000\) & 2.1784 & 0.00325 & 0.36 & 1.4767 & 0.02354 & 0.43 \\ \hline \multirow{2}{*}{
\begin{tabular}{c} **UP**, \(\omega=0\) \\ **UP**, \(\omega=100\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=10000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(\omega=1000\) \\ **UP**, \(
to unseen inputs. The results for these test sets are shown in Table V and Fig. 7. Defining the normalized loss as the obtained loss divided by the globally-optimal loss found by Gurobi [37], the following conclusions can be extracted from these results:

**Same distribution as in training:**

* Optimization 1: **RAYEN** is able to obtain the normalized loss closest to the globally-optimal one, while guaranteeing a zero violation of the constraints. The methods **UP**, **PP**, **DC3** with \(\omega=5000\), **Bar**, and **RAYEN** generate feasible solutions. In terms of computation time, **RAYEN** is between 20 and 7468 times faster than these methods.
* Optimization 2: Again, **RAYEN** is able to obtain the normalized loss closest to the globally-optimal one while guaranteeing a zero violation of the constraints. Compared to the feasible approaches (**UP** and **PP**), **RAYEN**'s computation time is between 2094 and 2417 times faster.

**Generalization: Different distribution as in training:**

* Optimization 1: Both **RAYEN** and **Bar** are able to obtain normalized losses very close to the globally-optimal one, while guaranteeing a zero violation of the constraints. Although **Bar**'s loss is slightly better than **RAYEN**'s (1.0070 vs. 1.0076), **Bar**'s method is 16 times slower, and it requires 197 times as many parameters as **RAYEN**. This high number of parameters required by **Bar** is due to the high number of vertices of the polyhedron \(\mathcal{Y}\). Note also that **Bar** only supports linear constraints.
* Optimization 2: **RAYEN** is the algorithm that has the best generalization (smallest normalized loss) while achieving a zero violation of the constraints.

It is also important to note that, when tested in the same distribution as the training set, the violations of both **UU** and **DC3** tend to decrease with higher values of \(\omega\). These violations however tend to become quite large when using a different distribution than the one used in training.

### _Computation time_

To study the computation time of the proposed approach, we generate random dense constraints of varying dimensions, and measure the computation time required by RAYEN to impose these constraints on \(\mathbf{y}\in\mathbb{R}^{k}\) for a random \(\mathbf{v}\in\mathbb{R}^{n}\) (see Fig. 3). All the details of how these random constraints are

Figure 8: Computation times of RAYEN for different dense random constraints of varying size. In these plots, \(\mathbf{y}\in\mathbb{R}^{k}\), \(\eta\) is the number of quadratic constraints (see Eq. 3), and \(\mu\) is the number of SOC constraints (see Eq. 4). The plots show the inference time for each sample in a batch of 2000 samples. The double-precision floating-point format float64 was used.

Figure 7: Comparison of the computation time vs. the normalized loss for the methods that generate feasible solutions. The dashed line corresponds to a normalized cost of 1.
generated are available in the attached code. The computation times per sample (in a batch of 2000 samples) achieved by RAYEN are available in Fig. 8. Note that they are short, even for high-dimensional variables, as shown in the following statistics:

* 2000 linear constraints on a 10K-dimensional variable: \(<0.9\) ms.
* 1000 dense convex quadratic constraints on a 1K-dimensional variable: \(<8\) ms.
* 500 dense SOC constraints on a 1K-dimensional variable with matrices \(\mathbf{M}_{j}\) of 300 rows: \(<8\) ms.
* A dense \(300\times 300\) LMI constraint on a 10K-dimensional variable: \(<11\) ms.

## VI Discussion and Limitations

### _Nonconvex Constraints_

Compared to methods such as DC3 [24], one limitation of the proposed framework is that right now it does not support generic nonconvex constraints. However, and as we have seen from the results in Section V-A, when the constraints are convex the performance of our method is much better than [24]. As future work (see also Section VII), we plan to extend RAYEN to nonconvex constraints, leveraging techniques such as sequential convexification.

### _Non-fixed Constraints_

This paper has focused on the case where the parameters that define the constraints (i.e., \(\mathbf{A}_{1},\mathbf{b}_{1},\mathbf{A}_{2},\mathbf{b}_{2},\mathbf{P}_{i},...\)) are fixed and known beforehand. A natural question to ask is whether or not RAYEN could be used when these parameters are not fixed, and depend on the previous layers (or input) of the network. The answer turns out to be yes, as long as some extra care is taken:

* The constraints must define a nonempty set at all times, and \(\mathcal{Z}\) must have interior points. There are many different ways to ensure this. For instance, if there are no equality constraints, one way would be to ensure these conditions: \[\mathbf{b}_{1}>\mathbf{0}\] (16) \[r_{i}<0\;\;i=0,...,\eta-1\] (17) \[d_{j}>\|\mathbf{s}_{j}\|\;\;\;j=0,...,\mu-1\] (18) \[\mathbf{F}_{k}\succ\mathbf{0}\] (19) Then, we can take \(\mathbf{z}_{0}=\mathbf{0}\), as it is guaranteed to be in the interior of \(\mathcal{Z}\). Eqs. 16, 17, and 18 can be easily enforced by, e.g., using sigmoid functions. Eq. 19 can be enforced by setting \[\mathbf{F}_{k}:=\mathbf{\Gamma}^{T}\mathbf{\Gamma}+\epsilon\mathbf{I}\;,\] where \(\epsilon>0\) is a predefined scalar, and where the matrix \(\mathbf{\Gamma}\) depends on the previous layers (or input) of the network. When there are only linear constraints, other ways to guarantee a nonempty feasible set are explained in [24, Section 4.1] and [19, Appendix A].
* \(\mathbf{P}_{i}\succeq\mathbf{0}\) (\(i=0,...,\eta-1\)) is required. Without any loss of generality, this can be enforced by setting \(\mathbf{P}_{i}:=\mathbf{V}_{i}^{T}\mathbf{V}_{i}\), where \(\mathbf{V}_{i}\) is a matrix that depends on the previous layers (or input) of the network.
* The matrices \(\mathbf{F}_{\alpha}\) (\(\alpha=0,...,k\)) must be symmetric. One way to guarantee this is by setting \(\mathbf{F}_{\alpha}:=(\mathbf{\Delta}_{\alpha}+\mathbf{\Delta}_{\alpha}^{T})/2\), where \(\mathbf{\Delta}_{\alpha}\) is a matrix that depends on the previous layers (or input) of the network.
* \(\text{aff}\left(\left\{\mathbf{y}\in\mathbb{R}^{k}|\mathbf{A}_{1}\mathbf{y}\leq\mathbf{b}_{1} \right\}\cap\mathcal{Y}_{\mathcal{Q}}\cap\mathcal{Y}_{\mathcal{S}}\cap\mathcal{ Y}_{M}\right)=\mathbb{R}^{k}\). Intuitively, this means that the inequality constraints must not _hide_ equality constraints. We can then take \(\mathbf{A}_{\mathcal{I}}\equiv\mathbf{A}_{1}\), \(\mathbf{b}_{\mathcal{I}}\equiv\mathbf{b}_{1}\), \(\mathbf{A}_{\mathcal{E}}\equiv\mathbf{A}_{2}\), and \(\mathbf{b}_{\mathcal{E}}\equiv\mathbf{b}_{2}\).

If these conditions are satisfied, then the offline phase (Section III) is not needed, since \(\mathbf{z}_{0}\) is already known, and the matrices \(\mathbf{A}_{p}\) and \(\mathbf{b}_{p}\) can be easily found online by simply computing the nullspace of \(\mathbf{A}_{\mathcal{E}}\) (Eq. 8).

## VII Conclusion and Future work

This work presented RAYEN, a framework to impose hard convex constraints on the output or latent variable of a neural network. RAYEN is able to guarantee by construction the satisfaction of any combination of linear, convex quadratic, SOC, and LMI constraints. When applied to approximate the solution of constrained optimization problems, RAYEN showcases computation times between 20 and 7468 times faster than state-of-the-art algorithms, while guaranteeing the satisfaction of the constraints at all times and obtaining a loss very close to the optimal one.

Future work includes the incorporation of more types of convex constraints, like for example exponential cone constraints. RAYEN could still be used for any convex constraint as long as the corresponding \(\kappa_{\square}\) for that constraint can be found. Even if an analytic solution for \(\kappa_{\square}\) does not exist, one could always leverage implicit layers that find the roots of nonlinear equations [40] to numerically find it. We also plan to study how RAYEN could be used to impose generic nonconvex constraints [24, 41] (via, e.g., sequential convexification), to impose Lipschitz constraints [42], to impose generic matrix manifolds constraints [43], or leveraged for safe robot control [44].

_Proof of \(\text{eig}\left(-\mathbf{H}^{-1}\mathbf{S}\right)=\text{eig}\left(\mathbf{L}^{T}\left(- \mathbf{S}\right)\mathbf{L}\right)\)_

Given that \(\mathbf{H}^{-1}\succ\mathbf{0}\) (because \(\mathbf{H}\succ\mathbf{0}\)), we can write its Cholesky decomposition as \(\mathbf{H}^{-1}=\mathbf{L}\mathbf{L}^{T}\), where \(\mathbf{L}\) is a lower triangular matrix with positive diagonal entries. Hence:

\[\text{eig}\left(-\mathbf{H}^{-1}\mathbf{S}\right)=\text{eig}\left(\mathbf{L}\mathbf{L}^{T} \left(-\mathbf{S}\right)\right)=\text{eig}\left(\mathbf{L}^{T}\left(-\mathbf{S}\right) \mathbf{L}\right)\;,\]

where in the last step we have use the fact that \(\text{eig}\left(\mathbf{A}\mathbf{B}\right)=\text{eig}\left(\mathbf{B}\mathbf{A}\right)\) for square matrices \(\mathbf{A}\) and \(\mathbf{B}\) (see, e.g., [45]). As \(\mathbf{L}^{T}\left(-\mathbf{S}\right)\mathbf{L}\) is clearly a real symmetric matrix, this also proves that all the eigenvalues of \(-\mathbf{H}^{-1}\mathbf{S}\) are real, even though \(-\mathbf{H}^{-1}\mathbf{S}\) may not be a symmetric matrix [46].
### _Planning problem_

The optimization problems used in Section V-A are trajectory optimization problems that aim at obtaining the clamped uniform B-Spline that minimizes a weighted sum of the velocity smoothness, the acceleration smoothness, and the jerk smoothness, while ensuring that the whole trajectory remains within a sequence of overlapping polyhedra that represent the free space. This problem (or small variations of it) appears extensively in many trajectory planning problems in Robotics, such as [3, 47, 52, 53]. Using the notation shown in in Table VI, the optimization problems are available in Table VII. The constraints are these:

* **Initial and final constraints:** The initial and final states are stop conditions. The initial position is fixed, while the final position of the trajectory is included in the cost as a soft constraint. This final position is taken inside \(\mathcal{P}_{n_{\text{polyh}}-1}\) (the last polyhedron of the corridor).
* **Corridor constraints:** Each interval of the trajectory is assigned to one polyhedron, and the convex hull property of the MINVO basis [54] is leveraged to ensure that each interval remains inside that assigned polyhedron. Compared to approaches that only impose constraints on discretization points along the trajectory, this approach guarantees safety for the whole trajectory \(t\in[0,t_{f}]\).
* **Dynamic limit constraints:** Optimization 1 uses linear constraints to impose a maximum velocity \(v_{\text{max}}\) and acceleration \(a_{\text{max}}\). Optimization 2 uses quadratic constraints instead, and also imposes a maximum jerk \(j_{\text{max}}\). Similar to the corridor constraints, we also leverage here the convex hull property of the MINVO basis.

Using RAYEN, some examples of the trajectories generated by the trained network for 100 different \(\boldsymbol{\gamma}\) taken from the testing set \(\text{TE}_{1,\text{in}}\) (see Appendix C) are shown in Fig. 9. All the constraints are guaranteed to be satisfied at all times.8

Footnote 8: As the degree of the spline of the Optimization 1 is 2, then the term \(\int_{0}^{t_{f}}\left\|\mathbf{j}(t)\right\|^{2}dt\) is clearly 0, and therefore \(\alpha_{\text{j}}\) does not affect the optimal solution. We keep this term simply for consistency purposes with Optimization 2, but it could also be removed as well.

### _Training, validation, and test sets_

The \(\boldsymbol{\gamma}\) used in Table III (the input of the network) is \(\boldsymbol{\gamma}:=\left[\begin{array}{cc}\alpha_{\text{v}}&\alpha_{\text {a}}&\alpha_{\text{j}}\end{array}\mathbf{p}_{f}^{T}\right]^{T}\), which contains the nonnegative weights \(\{\alpha_{\text{v}},\alpha_{\text{a}},\alpha_{\text{j}}\}\) and the final position \(\mathbf{p}_{f}\). Letting \(j\in\{1,2\}\) denote the optimization problem, we use the following training, validation, and test sets:

* **Training and Validation sets:**\(\left[\begin{array}{cc}\alpha_{\text{v}}&\alpha_{\text{a}}&\alpha_{\text{j}} \end{array}\right]^{T}\) is taken uniformly within the cube \([0.0,1.0]^{3}\), while \(\mathbf{p}_{f}\) is taken uniformly within \(\mathcal{P}_{n_{\text{polyh}}-1}\).
* **Testing sets:*
* \(\text{TE}_{j,\text{in}}\): Same distributions as above.
* \(\text{TE}_{j,\text{out}}\): \(\left[\begin{array}{cc}\alpha_{\text{v}}&\alpha_{\text{a}}&\alpha_{\text{j}} \end{array}\right]^{T}\) is taken uniformly within the cube \([1.0,2.0]^{3}\), while \(\mathbf{p}_{f}\) is taken uniformly within \(\mathcal{P}_{n_{\text{polyh}}-1}\).

\begin{table}
\begin{tabular}{|c|c|} \hline
**Symbol** & **Meaning** \\ \hline \hline \(\mathbf{p},\mathbf{v},\mathbf{a},\mathbf{j}\) & Position, velocity, acceleration, and jerk \\ \hline \(\mathbf{p}_{\text{o}},\mathbf{p}_{f}\) & Initial and final positions \\ \hline \(n_{\text{inter}}\) & Number of intervals of the B-Spline \\ \hline \(n_{\text{rep}}\) & Number of position control points of the B-Spline \\ \hline \(a_{\text{j}}\) & \(l\)-th acceleration control point of the B-Spline \\ \hline \(j_{\text{j}}\) & \(l\)-th jerk control point of the B-Spline \\ \hline \(j,j\) & \(j\) is the index of the interval of the trajectory, \\ \hline \(\mathcal{Q}_{j}^{\text{MV}},j_{j}^{\text{MV}}\) & \(\frac{\mathcal{Q}_{j}^{\text{MV}}}{\text{s}}\) is the set of position control points of the interval \\ \(j\) using the MINVO basis. Analogous definition for the velocity control points \(\mathcal{V}_{\text{MV}}^{\text{MV}}\) \\ \hline \(\mathcal{P}_{0},\ldots,\mathcal{P}_{n_{\text{polyh}}-1}\) & Sequence of overlapping polyhedra. They can be \\ \cline{1-1}  & obtained using methods such as [47, 48, 49, 50, 51] \\ \hline \([a]\) & Floor function (i.e., rounding to the nearest integer \(\leq a\)) \\ \hline \end{tabular}
\end{table}
Table VI: Notation used in Appendix B.

Figure 9: Trajectories generated by the trained network equipped with RAYEN. For each optimization \(j\in\{1,2\}\), each plot shows the trajectories obtained using 100 different \(\boldsymbol{\gamma}\) from the testing set \(\text{TE}_{j,\text{in}}\). \(\mathbf{p}_{0}\) is the initial condition, and \(\mathcal{P}_{i}\) denotes each polyhedron. The final position \(\mathbf{p}_{f}\) is taken randomly inside the last polyhedron.
## Acknowledgment

The authors would like to thank Victor Klemm and Dr. Kasra Khosoussi for helpful insights and discussions. Research supported in part by the Air Force Office of Scientific Research MURI FA9550-19-1-0386.
