This article is an expanded and revised version of the conference paper [1]. The following motivation encouraged us to delve into this subject. The rapid development of big data storage and processing technologies [2,3] has led to the emergence of optimization mathematical models in the form of multidimensional linear programming (LP) problems [4]. Such LP problems arise in industry, economics, logistics, statistics, quantum physics and other fields [5,6,7]. An important class of LP applications are non-stationary problems related to optimization in dynamic environments [8]. For a non-stationary LP problem, the objective function and/or constraints change over the computational process. Examples of non-stationary problems are the following: decision support in high-frequency trading [9,10], hydro-gas-dynamics problems [11], optimal control of technological processes [12,13,14], transportation [15,16,17], scheduling [18,19].

One of the standard approaches to solving non-stationary optimization problems is to consider each change as the appearance of a new optimization problem that needs to be solved from scratch [8]. However, this approach is often impractical because solving a problem from scratch without reusing information from the past can take too much time. Thus, it is desirable to have an optimization algorithm capable of continuously adapting the solution to the changing environment, reusing information obtained in the past. This approach is applicable for real-time processes if the algorithm tracks the trajectory of the moving optimal point fast enough. In the case of large-scale LP problems, the last requirement necessitates the development of scalable methods and parallel algorithms of LP.

Until now, one of the most common ways to solve the LP problem was a class of algorithms proposed and developed by Danzig based on the simplex method [20]. The simplex method was found to be efficient for solving a large class of LP problems. In particular, the simplex method easily takes advantage of any hyper-sparsity in a LP problem [21]. However, the simplex method has some fundamental features that limit its use for the solution of large LP problems. First, in certain cases, the simplex method has to iterate through all the vertices of the simplex, which corresponds to exponential time complexity [22,23,24]. Second, in most cases, the simplex method successfully solves LP problems containing up to 50,000 variables. However, a loss of precision is observed when the simplex method is used for solving large LP problems [25]. Such a loss of precision cannot be compensated even by applying such computational intensive procedures as “affine scaling” or “iterative refinement” [26]. Third, in the general case, the sequential nature of simplex method makes it difficult to parallelize on multiprocessor systems with distributed memory [27]. Numerous attempts have been made to make a scalable parallel implementation of the simplex method, but all of them were unsuccessful [28]. In all cases, the scalability boundary was from 16 to 32 processor nodes (see, for example, [29]). Khachiyan proved [30], using a variant of the ellipsoid method (proposed in the 1970’s by Shor [31] and Yudin and Nemirovskii [32]), that LP problems can be solved in polynomial time. However, attempts to apply this approach in practice were unsuccessful, since in the vast majority of cases, the ellipsoid method demonstrated much worse efficiency compared to the simplex method. Later, Karmarkar [33] proposed a polynomial-time interior-point algorithm, which can be used in practice. This algorithm generated the whole field of modern interior-point methods [34], which are able to solve large-scale LP problems with millions of variables and millions of equations [35,36,37,38,39]. Moreover, these methods are self-correcting, and therefore provide high precision of calculations. A general lack of interior-point methods is the need to find some feasible point satisfying all the constraints of the LP problem before starting calculations. Finding such an interior point can be reduced to solving an additional LP problem [40]. Another method for finding an interior point is the pseudo-projection method [41], which uses Fejér mappings [42]. One more significant disadvantage of the interior-point method is its poor scalability on multiprocessor systems with distributed memory. There are several successful parallel implementations of the interior-point method for particular cases (see, for example, [43]), but, in the general case, an efficient parallel implementation on multiprocessor systems could not be built. In accordance with this, the development and research of new approaches to solving multidimensional non-stationary LP problems in real time is an urgent direction.

One of the most promising approaches to solving complex problems in real time is the use of neural network models [44]. Artificial neural networks (ANN) are a powerful universal tool that is applicable to solving problems in almost all areas. The most popular neural network model is the feedforward neural network. Training and operation of such networks can be implemented very efficiently on GPUs [45]. An important property of a feedforward neural network is that the time to solve a problem is a known constant that does not depend on the problem parameters. This feature is necessary for real-time mode. Pioneering work on the use of neural networks to solve LP problems is the article by Tank and Hopfield [46]. The article describes a two-layer recurrent neural network. The number of neurons in the first layer is the number of variables of the LP problem. The number of neurons in the second layer coincides with the number of constraints of the LP problem. The first and second layers are fully connected. The weights and biases are uniquely determined by the coefficients and the right-hand sides of the linear inequalities defining the constraints, and the coefficients of the linear objective function. Thus, this network does not require training. The state of the neural network is described by the differential equation $\dot{x}\left(t\right)=\nabla E\left(x\left(t\right)\right)$, where $E\left(x\right(t\left)\right)$ is an energy function of a special type. Initially, an arbitrary point of the feasible region is fed to the input of the neural network. Then, the signal of the second layer is recursively fed to the first layer. Such process leads to convergence to a stable state in which the output remains constant. This state corresponds to the minimum of the energy function, and the output signal is a solution to the LP problem. The Tank and Hopfield approach has been expanded and improved in numerous works (see, for example, [47,48,49,50,51]). The main disadvantage of this approach is the unpredictable number of work cycles of the neural network. Therefore, a recurrent network based on an energy function cannot be used to solve large-scale LP problems in real time.

In the recent paper [52], a n-dimensional mathematical model for visualizing LP problems was proposed. This model makes it possible to use feedforward neural networks, including convolutional networks [53], to solve multidimensional LP problems, the feasible region of which is a bounded non-empty set. However, there are practically no works in scientific periodicals devoted to the use of convolutional neural networks for solving LP problems [54]. The reason is that convolutional neural networks focus on image processing, but there are no methods for constructing training datasets based on a visual representation of the multidimensional LP problems.

This article describes a new scalable iterative method for solving multidimensional LP problems. This method called “apex method”. The apex method allows you to generate training datasets for the development of feedforward neural networks capable of finding a solution to a multidimensional LP problem based on its visual representation. The apex method is based on the predictor-corrector framework. At the prediction step, a point belonging to the feasible region of the LP problem is calculated. The corrector step calculates a sequence of points converging to the exact solution of the LP problem. The rest of the paper is organized as follows. Section 2 provides a review of iterative projection-type methods and algorithms for solving linear feasibility problems and LP problems. Section 3 includes a theoretical basis of the apex method. Section 4 presents a formal description of apex method. Section 4.1 considers the implementation of the pseudoprojection operation in the form of sequential and parallel algorithms and provides an analytical estimation of the scalability of parallel implementation. Section 4.2 describes the quest stage. Section 4.3 describes the target stage. Section 5 presents an information about the software implementation of the apex method and describes the results of large-scale computational experiments on a cluster computing system. Section 6 discusses the issues related to the main contribution of this article, the advantages and disadvantages of the proposed approach, possible applications and some other aspects of using the apex method. Finally, in Section 7, we present our conclusions and comment on possible further studies of the apex method.

This section provides an overview of works devoted to iterative projection-type methods used to solve linear feasibility and LP problems. The linear feasibility problem can be stated as follows. Consider the system of linear inequalities in matrix form where $A\in {\mathbb{R}}^{m\times n}$, and $b\in {\mathbb{R}}^n$. To avoid triviality, we assume that $m>1$. The linear feasibility problem consists of finding a point $\tilde{x}\in {\mathbb{R}}^n$ satisfying matrix inequality system (1). We assume from now on that such a point exists.

Projection-type methods rely on the following geometric interpretation of the linear feasibility problem. Let $a_i\in {\mathbb{R}}^n$ be a vector formed by the elements of the ith row of the matrix A. Then, the matrix inequality $Ax⩽b$ is represented as a system of inequalities

Here and further on, $\langle ·,·\rangle$ stands for the dot product of vectors. We assume from now on that for all $i=1,\dots ,m$. For each inequality $〈a_i,x〉⩽b_i$, define the closed half space and its bounding hyperplane

For any point $x\in {\mathbb{R}}^n$, the orthogonal projection$\pi \left(x\right)$ of point x onto the hyperplane $H_i$ can be calculated by the equation

Here and below, $∥·∥$ denotes the Euclidean norm. Let us define

We assume that $M\ne \varnothing$, i.e. the solution of system (1) exists. In this case, M is a convex closed polytope in the space ${\mathbb{R}}^n$. In geometric interpretation, the linear feasibility problem consists of finding a point $\tilde{x}\in M$.

The forefathers of the iterative projection-type methods for solving linear feasibility problems are Kaczmarz and Cimmino. In [55] (English translation in [56]), Kaczmarz presented a sequential projections method for solving a consistent system of linear equations

His method, starting from an arbitrary point $x^{(0,m)}\in \mathbb{R}$, calculates the following sequence of point groups: for $k=1,2,3,\dots$ Here, ${\pi }_i$$(i=1,\dots ,m)$ is the orthogonal projection onto hyperplane $H_i$ defined by equation (6). This sequence converges to the solution of system (8). Geometrically, the method can be interpreted as follows. The initial point $x^{(0,m)}$ is projected orthogonally onto hyperplane $H_1$. The projection is the point $x^{(1,1)}$, which now is thrown onto $H_2$. The resulting point $x^{(1,2)}$ is then thrown onto $H_3$ and gives the point $x^{(1,3)}$, etc. As a result, we obtain the last point $x^{(1,m)}$ from the first point group. The second point group is constructed in the same way, starting from the point $x^{(1,m)}$. The process is repeated for $k=2,3,\dots$

Cimmino proposed in [57] (English description in [58]) a simultaneous projection method for the same problem. This method uses the following orthogonal reflection operation which calculates the point ${\rho }_i\left(x\right)$ symmetric to the point x with respect to the hyperplane $H_i$. For the current approximation $x^{\left(k\right)}$, the Cimmino method simultaneously calculates reflections with respect to all hyperplanes $H_i$$(i=1,\dots ,m)$, and then a convex combination of these reflections is used to form the next approximation: where $w_i>0$$(i=1,\dots ,m)$, ${\displaystyle \sum _{i=1}^m}{w}_i=1$. When $w_i=\frac{1}{m}$$(i=1,\dots ,m)$, equation (11) is transformed into the following equation:

Agmon [59] and Motzkin and Schoenberg [60] generalized the projection method from equations to inequalities. To solve problem (1), they introduce the relaxed projection where $0<\lambda <2$. It is obvious that ${\pi }_i^1\left(x\right)={\pi }_i\left(x\right)$. To calculate the next approximation, relaxed projection method uses the following equation: where

Informally, the next approximation $x^{(k+1)}$ is a relaxed projection of the previous approximation $x^{\left(k\right)}$ with respect to the furthest hyperplane $H_l$ bounding the half-space ${\widehat{H}}_l$ not containing $x^{\left(k\right)}$. Agmon in [59] showed that sequence $x^{\left(k\right)}$ converges, as $k\to \infty$, to a point on the boundary of M.

Censor and Elfving, in [61], generalized the Cimmino method to the case of linear inequalities. They consider the relaxed projection onto the half-space ${\widehat{H}}_i$ defined as follows: that gives the equation Here, $0<\lambda <2$, and $w_i>0$$(i=1,\dots ,m)$, ${\displaystyle \sum _{i=1}^m}{w}_i=1$. In [62], De Pierro proposed an approach to convergence proof for this method, which differs from the approach of Censor and Elfving. De Piero’s approach is also acceptable for the case when the underlying system of linear inequalities is infeasible. In this case, for $\lambda =1$, sequence (17) converges to the point that is the minimum of the function $f\left(x\right)={\sum }_{i=1}^m{w}_i{∥{\widehat{\pi }}_i\left(x\right)-x∥}^2$, i.e., it is a weighted (with the weights $w_i$) least squares solution of system (1).

The Cimmino-like methods allow efficient parallelization, since orthogonal projections (reflections) can be calculated simultaneously and independently. The article [63] investigates the efficiency of parallelization of the Cimmino-like method on Xeon Phi manycore processors. In [64], the scalability of the Cimmino method for multiprocessor systems with distributed memory is evaluated. The applicability of the Cimmino-like method for solving non-stationary systems of linear inequalities on computing clusters is considered in [65].

As a recent work, we can mention article [66], which extends the Agmon-Motzkin-Schoenberg relaxation method for the case of semi-infinite inequality systems. The authors consider the system with an infinite number of inequalities in the finite-dimensional Euclidean space ${\mathbb{R}}^n$: where I is an arbitrary infinite index set. The main idea of the method is as follows. Let the hyperplane $H_x^{(\infty )}=sup\left\{\left.max\left(〈a_i,x〉-b_i,0\right)\right|i\in I\right\}$ be the biggest violation with respect to x. Let $x^{\left(0\right)}$ be an arbitrary initial point. If the current iteration $x^{\left(k\right)}$ is not a solution of system (18), then let $x^{(k+1)}$ be the orthogonal projection of $x^{\left(k\right)}$ onto a hyperplane $H_i$ ($i\in I$) near the biggest violation $H_{x^{\left(k\right)}}^{(\infty )}$. If system (18) is consistent, then the sequence $\left\{x^{\left(k\right)}\left|k=1,2,\dots \right.\right\}$ generated by the described method converges to the solution of this system.

Solving systems of linear inequalities is closely related to LP problems, so projection-type methods can be effectively used to solve this class of problems. The equivalence of the linear feasibility problem and the LP problem is based on the primal-dual LP problem. Consider the primal LP problem in the matrix form: where $c,x\in {\mathbb{R}}^n$, $b\in {\mathbb{R}}^m$, $A\in {\mathbb{R}}^{m\times n}$, and $c\ne \mathbf{0}$. Here and below, $〈·,·〉$ stands for the dot product of vectors. Let us construct the dual problem with respect to problem (19): where $u\in {\mathbb{R}}^m$. The following primal-dual equality holds:

In [67,68], Eremin proposed the following method based on the primal-dual approach. Let the inequality system define the feasible region of primal problem (19). This system is obtained by adding to the system $Ax⩽b$ the vector inequality $-x⩽0$. In this case, $A^{\prime }\in {\mathbb{R}}^{(m+n)\times n}$, and $b^{\prime }\in {\mathbb{R}}^{m+n}$. Let $a_i^{\prime }$ stand the ith row of the matrix $A^{\prime }$. For each inequality $〈a_i^{\prime },x〉⩽b_i^{\prime }$, define the closed half space and its bounding hyperplane

Let ${\pi }_i^{\prime }\left(x\right)$ stand the orthogonal projection of point x onto the hyperplane $H_i^{\prime }$:

Let us define the projection onto the half-space ${\widehat{H}}_i^{\prime }$:

This projection has the following two properties:

Define ${\phi }_1:{\mathbb{R}}^n\to {\mathbb{R}}^n$ as follows:

In the same way, define the feasible region of dual problem (20) as follows: where $D=A^T\in {\mathbb{R}}^{n\times m}$, $D^{\prime }\in {\mathbb{R}}^{(m+n)\times m}$, and $c^{\prime }\in {\mathbb{R}}^{n+m}$. Denote and define ${\phi }_2:{\mathbb{R}}^m\to {\mathbb{R}}^m$ as follows:

Now, define ${\phi }_3:{\mathbb{R}}^{n+m}\to {\mathbb{R}}^{n+m}$ as follows: which is corresponding to equation (21). Here, $[·,·]$ stands for the concatenation of vectors.

Finally, define $\phi :{\mathbb{R}}^{n+m}\to {\mathbb{R}}^{n+m}$ as follows:

If the feasible region of the primal problem is a bounded and nonempty set, then the sequence converges to the point $[\overline{x},\overline{u}]$, where $\overline{x}$ is the solution of primal problem (19), and $\overline{u}$ is the solution of dual problem (20).

Article [69] proposes a method for solving non-degenerate LP problems based on calculating the orthogonal projection of some special point, independent of the main part of data describing the LP problem, onto a problem-dependent cone generated by the constraint inequalities. Actually, this method solves a symmetric positive definite system of linear equations of a special kind. The author demonstrates a finite algorithm of active-set family that is capable of calculating orthogonal projections for problems with up to thousands of rows and columns. The main drawback of this method is a significant increasing the dimension of the primary problem.

In article [70], Censor proposes the linear superiorization (LinSup) method as a tool for solving LP problems. The LinSup method does not guarantee finding the minimum point of the LP problem, but it directs the linear feasibility-seeking algorithm that it uses toward a point with a decreasing value of the objective function. This process is not identical with that employed by LP solvers but it is a possible alternative to the Simplex method for problems of huge size. The basic idea of LinSup is to add an extra term, called perturbation term, to the iterative equation of the projection method. The perturbation term steers the feasibility-seeking algorithm toward reduced the objective function values. In the case of LP problem (19), the objective function is $f\left(x\right)=〈c,x〉$, and LinSup adds $\left(-\eta \frac{c}{∥c∥}\right)$ as a perturbation term to iterative equation (17):

Here, $0<\eta <1$ is a perturbation parameter.

Article [52] proposes a mathematical model for the visual representation of multidimensional LP problems. To visualize a feasible LP problem, an objective hyperplane $H_c$ is introduced, the normal to which is the gradient of the objective function $f\left(x\right)=〈c,x〉$. In the case of seeking the maximum, the objective hyperplane is positioned in such a way that the value of the objective function at all its points is greater then the value of the objective function at all points of the convex polytope M, which is the feasible region of the LP problem. For any point $g\in H_c$, the objective projection ${\gamma }_M\left(g\right)$ onto M is defined as follows:

Here, ${\pi }_{H_c}\left(x\right)$ denotes the orthogonal projection onto $H_c$. On the objective hyperplane $H_c$, a rectangular lattice of points $\mathfrak{G}\in {\mathbb{R}}^n\times {\mathbb{R}}^{K^{(n-1)}}$ is constructed, where K is the number of lattice points in one dimension. Each point $g\in \mathfrak{G}$ is mapped to the real number $∥{\gamma }_M\left(g\right)-g∥$. This mapping generates a matrix of dimension $(n-1)$, which is an image of the LP problem. This approach opens up the possibility of using feed-forward artificial neural networks, including convolutional neural networks, to solve multidimensional LP problems. One of the main obstacles to the implementation of this approach is the problem of generating a training set. The literature review shows that there is no suitable method capable of constructing such a training set compatible with the described approach. In the next sections, we present such a method.

In this section we present a theoretical basis used to construct the apex method. Consider the LP problem in the following form: where $c\in {\mathbb{R}}^n$, $b\in {\mathbb{R}}^m$, $A\in {\mathbb{R}}^{m\times n}$, $m>1$, and $c\ne \mathbf{0}$. We assume that the constraint $x⩾\mathbf{0}$ is also included in the system $Ax⩽b$ in the form of the following inequalities:

Let $\mathcal{P}$ stand for the set of row indices in matrix A:

Let $a_i\in {\mathbb{R}}^n$ be a vector formed by the elements of the ith row of the matrix A, and $a_i\ne \mathbf{0}$ for all $i\in \mathcal{P}$. We denote by ${\widehat{H}}_i$ the closed half space defined by the inequality $〈a_i,x〉⩽b_i$, and by $H_i$ the hyperplane bounding ${\widehat{H}}_i$:

The geometric meaning of this definition is that a ray outgoing from a point belonging to a half-space in the direction of vector c belongs to this half-space.

The following proposition provides the necessary and sufficient condition for the c-recessivity of the half-space.

Denote i.e., $e_c$ stands for the unit vector parallel to vector c.

Define i.e., ${\mathcal{I}}_c$ is the set of indices for which the half-space ${\widehat{H}}_i$ is c-recessive. We assume from now on that

The following proposition specifies the region containing a solution of LP problem (38).

Let ${\pi }_i\left(x\right)$ stand for the orthogonal projection of point x onto hyperplane $H_i$:

The next proposition provides a continuous M-Fejér mapping, which will be used in the apex method.

The correctness of this definition is ensured by propositions 4 and 5.

In this section, we describe a new scalable iterative method for solving LP problem (38), called the “apex method”. The apex method is based on the predictor-corrector framework and proceeds in two stages: quest (predictor) and target (corrector). The quest stage calculates a rough initial approximation of LP problem (38). The target stage refines the initial approximation with a given precision. The main operation used in the apex method is an operation that calculates a pseudoprojection according to Definition 4. This operation is used both in the quest stage and in the target stage. In the next section, we describe and investigate a parallel algorithm for calculating a pseudoprojection.