On Implementing a Two-Step Interior Point Method for Solving Linear Programs

Abstract

A new two-step interior point method for solving linear programs is presented. The technique uses a convex combination of the auxiliary and central points to compute the search direction. To update the central point, we find the best value for step size such that the feasibility condition is held. Since we use the information from the previous iteration to find the search direction, the inverse of the system is evaluated only once every iteration. A detailed empirical evaluation is performed on NETLIB instances, which compares two variants of the approach to the primal-dual log barrier interior point method. Results show that the proposed method is faster. The method reduces the number of iterations and CPU time(s) by $27\%$ and $18\%$, respectively, on NETLIB instances tested compared to the classical interior point algorithm.

1. Introduction

1.1. Background

The following linear optimization (LO) problem is considered:

$$\left(P\right)\hspace{1em}min\{c^{T}x:Ax=b,x\ge 0\},$$

and the corresponding dual

$$\begin{array}{c}\hfill \left(D\right)\hspace{1em}\hspace{1em}max\{b^{T}y:A^{T}y+s=c,s\ge 0\},\end{array}$$

where $A\in {\mathbb{R}}^{m\times n},$$x,c,s\in {\mathbb{R}}^n$ and $y,b\in {\mathbb{R}}^m$.

LO is an essential tool used widely in theoretical and practical science and engineering domains. Its applications extend to various fields such as operations research, engineering, economics and combinatorics. LO finds numerous applications across various fields. For instance, it is used in ${\mathcal{L}}_1$-regularized support vector machines (SVMs) [1], basis pursuit (BP) problems [2], sparse inverse covariance matrix estimation (SICE) [3], non-negative matrix factorization (NMF) [4] and maximum a posteriori (MAP) inference [5].

The Interior-Point Method (IPM) is a powerful optimization technique that has been widely used in the field of LO. It offers an alternative approach to solving LO problems compared to traditional simplex methods. Karmarkar [6] described the first efficient algorithm for the IPM. Nesterov and Nemirovskii [7] introduced a class of self-concordant barrier functions and extended the algorithm to a more general type of optimization problems, including convex programming, non-linear complementarity problem (NCP), semidefinite optimization (SDO) and second-order cone optimization (SOCO) problems.

1.2. Literature Review

IPMs are generally divided into two categories, that is, small-update methods and large-update methods [8]. In practice, large-update methods are significantly more efficient. However, in theory, the best results are associated with small-update methods. From another perspective, IPMs are divided into feasible and infeasible methods [9]. Feasible methods start from a point that satisfies the constraints and seek to find the optimal solution within the feasible region. In contrast, infeasible methods begin with a point that does not meet the problem’s constraints (infeasible). The algorithm first seeks to find a feasible point, and once a feasible point is found, it moves towards the optimal solution [9]. The IPM based on the kernel function is a popular method. The kernel function is used not only to find the search direction for the algorithm but also to estimate a lower bound for the step size [10]. In this type of IPM, the complexity bounds of the interior point algorithm are computed based on the kernel function. Terlaky et al. [10] proposed a new variant of IPMs for solving linear optimization problems in which the self-regular barrier function replaces the logarithmic barrier function and they showed that their method has the best-known complexity bounds. Following this, several different kernel functions with the best-known complexity bounds have been introduced, including the exponential kernel functions [11], trigonometric kernel functions [12,13] and hyperbolic kernel functions [14,15].

Predictor-corrector methods for interior point algorithms are an effective approach for solving both linear and nonlinear optimization problems [16,17]. These methods incorporate two main steps, i.e., the predictor step and the corrector step. In the predictor step, an approximate direction is determined by solving a linearized system of the Karush–Kuhn–Tucker (KKT) conditions [17]. In the subsequent corrector step, this direction is refined to ensure that the iterates remain within the feasible region and maintain the balance between the primal and dual variables. This refinement typically involves solving another system of linear equations that corrects any infeasibilities introduced in the predictor step, thereby enhancing the convergence properties of the algorithm. Mehrotra’s algorithm [17] is a popular predictor-corrector algorithm used for solving optimization problems. In recent years, several new predictor-corrector algorithms have been introduced, further enhancing the effectiveness and efficiency of IPMs for various optimization problems [18,19].

The Newton method is a popular optimization technique that plays a vital role in interior-point methods. In the context of IPMs, the Newton method is used to solve the system of equations that arise from the Karush–Kuhn–Tucker (KKT) optimality conditions, which include both the primal and dual feasibility conditions, as well as the complementary slackness conditions. This approach is known for its fast convergence rate, often exhibiting quadratic convergence. A modification of Newton’s method was proposed by McDougall and Wotherspoon [20] that can significantly reduce the number of iterations when used to the root of functions. Their approach has a convergence of the order of $1+\sqrt{2}$. To compute the search direction, their algorithm uses an auxiliary point and computes the gradient of the objective function only once in each iteration.

McDougall and Wotherspoon [20] applied this method to solve the power flow problem, which is a multivariate problem. They used the improved method to find the roots of a system of multivariable equations and demonstrated that it can reduce the number of iterations compared to Newton’s method. Additionally, Argyros et al. [21] present some theoretical results related to the convergence of the method. In Example 4 of [21], it is shown that the method is effective for solving multivariable problems and can converge to the optimal solution.

1.3. Contribution

Motivated by these ideas, this paper makes the following contributions:

• To solve optimization problems using the IPM, the search direction plays an important role. We use the idea proposed in [20] for determining the search direction. The main goal of this paper is to extend the strategy proposed in [20] to IPMs. This paper introduces a two-step method that calculates the inverse of the matrix only once to determine the search direction in each iteration.
  The next goal of this article is to demonstrate how this method can be integrated with the feasible interior-point algorithm for linear optimization problems to enhance the efficiency of the IPM. We combine the technique described in [20] with the algorithm presented in [8], resulting in the development of two new algorithms.
• To demonstrate the efficiency of the new method, we implement the algorithms and solve a collection of NETLIB test problems. The results show that the proposed method can improve the algorithm’s efficiency by $27\%$ in terms of the number of iterations and $18\%$ in terms of CPU time.
  The implementation on a Graphics Processing Unit (GPU), a specialized processor designed to accelerate graphics rendering and computational tasks, shows a reduction in the number of iterations and GPU time compared to the classical algorithm by $39\%$ and $30\%$, respectively.

The paper is organized as follows: Section 2 briefly recalls the basics of the IPMs and the central path for solving linear programs. We explain the two-step algorithm proposed here next. Section 3 reports the numerical results on a collection of LP instances from NETLIB for CPU and GPU implementations of the algorithm. We end with concluding remarks in Section 4.

In this paper, we use the following notational conventions. The Euclidean norm of a vector is denoted by $\parallel .\parallel$. The non-negative and positive orthants are represented by ${\mathbb{R}}_{+}^n$ and ${\mathbb{R}}_{++}^n$ respectively. The vectors $xs$ and $\frac{x}{s}$ indicate coordinate-wise operations on the vectors x and s, with components $x_i{s}_i$ and $\frac{x_i}{s_i}$ respectively. We also use diagonal matrices X and S with entries x and s respectively. For any variable v belonging to $\{x,y,s\}$, subscripts denote the $n^{th}$ iteration value. For instance, $x_n$ refers to the value of x after n updates. If the current iteration is n, we use v instead of $v_n$ and $v_{+}$ instead of $v_{n+1}$ for $v\in \{x,y,s,\tilde{x},\tilde{y},\tilde{s},\widehat{x},\widehat{y},\widehat{s}\}$. Here, x refers to the current point and $x_{+}$ refers to the updated value of x.

2. Algorithm

In this paper, without loss of generality, we make the following assumptions:

• A is a full rank matrix, that is rank $A=m\le n$.
• Both problems (P) and (D) satisfy the Interior Point Condition (IPC) [8], i.e., there exists $(x^0,s^0)>0$ and $y^0$ so that:

  $$\begin{array}{ccc}\hfill A{x}^0& =& b,\hspace{1em}\hspace{1em}{x}^0>0,\hfill \\ \hfill A^T{y}^0+s^0& =& c,\hspace{1em}\hspace{1em}{s}^0>0.\hfill \end{array}$$

An optimal solution for (P) and (D) can be found by solving the following system [10]:

$$\begin{array}{ccc}\hfill Ax& =& b,\hspace{1em}\hspace{1em}x\ge 0,\hfill \\ \hfill A^{T}y+s& =& c,\hspace{1em}\hspace{1em}s\ge 0,\hfill \\ \hfill xs& =& 0,\hfill \end{array}$$

(1)

where $xs$ stands for the coordinate-wise (Hadamard) product of the vectors x and s. Note that the first and second equations in (1) guarantee the feasibility of x and $(y,s)$ for problems (P) and (D), respectively, and the last equation is satisfied only if the feasible solution is optimal (complementarity conditions).

The algorithm presented here is a two-step primal-dual IPM; therefore, we summarize the key idea behind the generic primal-dual IPMs for solving linear programming problems from the literature, see e.g., [8,10]. The key idea is to replace the complementarity conditions with parameterized equations $xs=\mu \mathbf{e}$, where $\mathbf{e}$ denotes the all-one vector of length n and $\mu$ is a real positive parameter. With these changes, the system of Equation (1) can be expressed as [10]:

$$\begin{array}{ccc}\hfill Ax& =& b,\hspace{1em}\hspace{1em}x\ge 0,\hfill \\ \hfill A^{T}y+s& =& c,\hspace{1em}\hspace{1em}s\ge 0,\hfill \\ \hfill xs& =& \mu \mathbf{e}.\hfill \end{array}$$

(2)

The parametric system (2) has a unique solution for any $\mu >0$ if the IPC holds [10] (Theorem 1.2.4). For a given $\mu >0$, denote the unique solution of (2) by $\left(x\right(\mu ),y(\mu ),s(\mu \left)\right)$, where $x\left(\mu \right)$ is the $\mu$-center of (P) and $\left(y\right(\mu ),s(\mu \left)\right)$ is called the $\mu$-center of (D). The set of all $\mu$-centers, for all $\mu >0$, defines a path the so-called central path for (P) and (D). It has been shown that, as $\mu \to 0$, then the limit of the central path exists and converges to an analytic center of the optimal solutions set of (P) and (D) [22,23].

Given a current point $(x,y,s)$, the goal is to find direction $(\Delta x,\Delta y,\Delta s)$ such that the point $(x+\Delta x,y+\Delta y,s+\Delta s)$ lies on the central path at $\mu$. If it did lie exactly on the central path, then Equation (2) will be satisfied, i.e., [10]:

$$\begin{array}{ccc}\hfill A(x+\Delta x)& =& b\hfill \\ \hfill A^T(y+\Delta y)+(s+\Delta s)& =& c\hfill \\ \hfill (X+\Delta x)(S+\Delta s)& =& \mu \mathbf{e},\hfill \end{array}$$

(3)

where:

$$X:=\mathrm{diag}\left(x\right),\hspace{1em}\hspace{1em}S:=\mathrm{diag}\left(s\right)$$

Assuming that the point $(x,y,s)$ is feasible for the LO and dropping the non-linear term $\Delta x\Delta s$ [10], we get the following system of equations:

$$\begin{array}{ccc}\hfill A\Delta x& =& 0\hfill \\ \hfill A^T\Delta y+\Delta s& =& 0\hfill \\ \hfill S\Delta x+X\Delta s& =& \mu \mathbf{e}-XS\mathbf{e}\hfill \end{array}$$

(4)

System (4) has a unique solution, referred to as the search direction for IPMs. Many interior-point algorithms rely on the Newton method to solve the system (4) and find the search direction. Therefore, improving the speed of convergence of Newton’s method significantly impacts the convergence speed and efficiency of IPMs. To calculate the search direction in Newton’s method, we use the method proposed in [20], which has a convergence speed of $1+\sqrt{2}$. Accordingly, we introduce a two-step method for the IPM based on the approach proposed in [20]. This change can significantly reduce the number of iterations needed to achieve the optimal solution.To apply the mentioned method, we present the system (4) as a problem of finding roots, followed by an explanation of the new approach.

2.1. Relationship with Root Finding

Define:

$$\xi =\left[\begin{array}{c}x\\ y\\ s\end{array}\right]\mathrm{and}F\left(\xi \right)=\left[\begin{array}{c}Ax-b\\ A^{T}y+s-c\\ \mu \mathbf{e}-XS\mathbf{e}\end{array}\right],$$

then a solution to the system (2) is ${\xi }^{*}$ such that $F\left({\xi }^{*}\right)=0$. The root ${\xi }^{*}$ can be computed using Newton’s method (in the neighbourhood of the root), which, given $\xi$, finds $\Delta \xi$ such that $F(\xi +\Delta \xi )=0$ [20]. Using Taylor series approximation (linear approximation) around $\xi$, we have [20]:

$$F\left(\xi \right)+F^{\prime }\left(\xi \right)\Delta \xi \simeq 0$$

where:

$$\mathrm{Jacobian}{F}^{\prime }=\left[\begin{array}{ccc}\frac{\partial F_1}{\partial {\xi }_1}& \frac{\partial F_1}{\partial {\xi }_2}& \frac{\partial F_1}{\partial {\xi }_3}\\ \frac{\partial F_2}{\partial {\xi }_1}& \frac{\partial F_2}{\partial {\xi }_2}& \frac{\partial F_2}{\partial {\xi }_3}\\ \frac{\partial F_3}{\partial {\xi }_1}& \frac{\partial F_3}{\partial {\xi }_2}& \frac{\partial F_3}{\partial {\xi }_3}\end{array}\right]=\left[\begin{array}{ccc}A& 0& 0\\ 0& A^T& I\\ S& 0& X\end{array}\right]\mathrm{and}\Delta \xi =\left[\begin{array}{c}\Delta x\\ \Delta y\\ \Delta s\end{array}\right].$$

This is precisely the system that was derived earlier (4). This approximation yields an iterative algorithm, the Newton–Raphson method, for finding roots. Starting with an initial estimate ${\xi }_0$, the current estimate ${\xi }_n$ of the root is updated using the following rule for some appropriate step size $\alpha$:

$${\xi }_{n+1}={\xi }_n-\alpha {\left(F^{\prime }\left({\xi }_n\right)\right)}^{-1}F\left({\xi }_n\right)$$

The proposed method generalizes the Newton–Raphson method as follows (assume any F). In the first step, an auxiliary point ${\tilde{\xi }}_n$ is computed using a linear approximation to the function around the current estimate of the root ${\xi }_n$. The average ${\widehat{\xi }}_{n+1}$ of the auxiliary point and the current estimate is computed in the second step. In the third step, the current estimate is updated using the linear approximation to F around ${\xi }_n$, but instead of using the Jacobian at point ${\xi }_n$, the Jacobian is evaluated at the average ${\widehat{\xi }}_{n+1}.$ The update rules used in $n^{th}$ iteration can be summarized as [20]:

$$\begin{array}{ccc}\hfill {\tilde{\xi }}_{n+1}& =& {\xi }_n-\alpha {\left(F^{\prime }\left({\widehat{\xi }}_n\right)\right)}^{-1}F\left({\xi }_n\right)\hfill \\ \hfill {\widehat{\xi }}_{n+1}& =& \frac{1}{2}({\tilde{\xi }}_{n+1}+{\xi }_n)\hfill \\ \hfill {\xi }_{n+1}& =& {\xi }_n-\alpha {\left(F^{\prime }\left({\widehat{\xi }}_{n+1}\right)\right)}^{-1}F\left({\xi }_n\right)\hfill \end{array}$$

This sequence of updates requires two evaluations of the Jacobian, once at point ${\widehat{\xi }}_n$ and the second time at point ${\widehat{\xi }}_{n+1}.$ This requirement doubles the computation burden. In practice, we use information regarding $F^{\prime }\left({\widehat{\xi }}_{n+1}\right)$ for the next iteration. Thus only a single evaluation of the Jacobian is needed in each iteration if we save the result from the previous iteration.

2.2. Two-Step IPMs

Here, we extend the idea presented in Section 2.1 to the interior-point algorithm. The extension involves two phases:

• Phase 1: We obtain a feasible starting point $(x,y,s)$ using the method presented in [17]. Then, we create an auxiliary point $(\tilde{x},\tilde{y},\tilde{s})$, which initially has the same values as $(x,y,s)$. The parameter $\mu$ is updated by a factor of $1-\theta$, and $(\tilde{x},\tilde{y},\tilde{s})$ is set to $(x,y,s)$. After that, we set $(\widehat{x},\widehat{y},\widehat{s})=\frac{1}{2}((\tilde{x},\tilde{y},\tilde{s})+(x,y,s))$ and solve the following system to determine the search direction:

  $$\begin{array}{ccc}\hfill A\Delta x& =& 0\hfill \\ \hfill A^T\Delta y-\Delta s& =& 0\hfill \\ \hfill \widehat{S}\Delta x+\widehat{X}\Delta s& =& \mu \mathbf{e}-XS\mathbf{e}\hfill \end{array}$$

  (5)
  To satisfy the feasibility condition, the algorithm performs a line search to find the maximum step size. To obtain the new point, the current point is updated in the following way:

  $$(x,y,s)\leftarrow (x+\alpha \Delta x,y+\alpha \Delta y,s+\alpha \Delta s)$$
• Phase 2: In this phase, starting from $k\ge 1$, we begin by updating the parameter $\mu$ and constructing the auxiliary point. Then, we utilize it to create the principal point. As in Phase 1, we solve the following system to determine the search direction for the auxiliary point:

  $$\begin{array}{ccc}\hfill A\Delta x& =& 0\hfill \\ \hfill A^T\Delta y-\Delta s& =& 0\hfill \\ \hfill \widehat{S}\Delta x+\widehat{X}\Delta s& =& \mu \mathbf{e}-XS\mathbf{e}.\hfill \end{array}$$

  (6)
  Since the left-hand side of (6) is the same as the left-hand side of (5), so there is no need to compute the inverse. Next, a line search is performed to determine the maximum step size, and the auxiliary point is updated as follows:

  $$(\tilde{x},\tilde{y},\tilde{s})\leftarrow (x+\alpha \Delta x,y+\alpha \Delta y,s+\alpha \Delta s)$$
  Next, the middle point is updated as:

  $$(\widehat{x},\widehat{y},\widehat{s})\leftarrow \frac{1}{2}((\tilde{x},\tilde{y},\tilde{s})+(x,y,s))$$
  Finally, the algorithm updates the principal point using the search direction $(\Delta x,\Delta y,\Delta s)$ is obtained by solving the following system:

  $$\begin{array}{ccc}\hfill A\Delta x& =& 0\hfill \\ \hfill A^T\Delta y-\Delta s& =& 0\hfill \\ \hfill \widehat{S}\Delta x+\widehat{X}\Delta s& =& \mu \mathbf{e}-XS\mathbf{e}\hfill \end{array}$$

  (7)
  As earlier, to ensure that the principal point satisfies the feasibility condition, a line search is performed to determine the maximum step size. The principal point is then updated as follows:

  $$(x,y,s)\leftarrow (x+\alpha \Delta x,y+\alpha \Delta y,s+\alpha \Delta s).$$

In order to determine the maximum value for the step size $\alpha$ in each iteration, we utilize the method presented by Cai et al. [24]. This method initially calculates an upper bound based on the positivity constraints of the variables. If $(x+\alpha \Delta x,y+\alpha \Delta y,s+\alpha \Delta s)$ is not feasible, then $\alpha$ is decreased by a certain multiplicative factor until the point is feasible. This approach is more efficient than determining the optimal value for $\alpha$, as it avoids an additional computational step. Algorithm 1 outlines the steps of this new approach.

Algorithm 1 A two-step feasible IPM algorithm for LOs

The method used in Algorithm 1 stands out from existing interior point algorithms. It utilizes a two-step Newton method to determine the search direction, resulting in a higher convergence rate of $1+\sqrt{2}$ compared to traditional methods. This improved rate leads to a significant reduction in the number of iterations needed by the algorithm. Additionally, unlike interior-point predictor-corrector methods that require computing the matrix’s inverse twice per iteration, Algorithm 1 only needs to compute it once per iteration. Specifically, it uses the previous iteration’s matrix inverse to calculate the search direction for the auxiliary point. This computational characteristic aligns Algorithm 1 with one-step methods but requires fewer iterations to converge.

Incorporating the idea proposed by Wang et al. [25], we will now implement Algorithm 2. To begin, they introduced a fresh variable v, defined below:

$$\begin{array}{c}\hfill v:=\sqrt{\frac{xs}{\mu }}\end{array}$$

(8)

Next, they introduced a kernel function as $\psi \left(t\right)=\frac{{(t-1)}^2}{2}$. Based on this kernel function and the fact that $(x,s)>0$, they showed that [25]:

$$\begin{array}{c}\hfill v=\mathbf{e}\iff v^2=\mathbf{e}\iff \frac{xs}{\mu }=\sqrt{\frac{xs}{\mu }}\end{array}$$

Now, we use this idea and rewrite the (6) and (7) as [25]:

$$\begin{array}{ccc}\hfill A\Delta x& =& 0\hfill \\ \hfill A^T\Delta y-\Delta s& =& 0\hfill \\ \hfill \widehat{S}\Delta x+\widehat{X}\Delta s& =& \mu \mathbf{e}-\sqrt{\mu }\sqrt{XS\mathbf{e}}\hfill \end{array}$$

(9)

and

$$\begin{array}{ccc}\hfill A\Delta x& =& 0\hfill \\ \hfill A^T\Delta y-\Delta s& =& 0\hfill \\ \hfill \widehat{S}\Delta x+\widehat{X}\Delta s& =& \mu \mathbf{e}-\sqrt{\mu }\sqrt{XS\mathbf{e}}\hfill \end{array}$$

(10)

The resulting changes lead to Algorithm 2 which is described next:

Algorithm 2 A two-step feasible IPM algorithm for the algorithm proposed in [25]

2.3. Discussion

The idea in the Algorithm is also used in [20] for finding roots of a univariate function f using a modification to the Newton–Raphson method. The method requires two steps a predictor step and a corrector step. The predictor step computes the point ${\tilde{x}}_{+}$ based on a linear approximation of the function around the current point x. In the corrector step, the new point is calculated using the derivative at a point which is the average of the current point and the predictor. This can be expressed as [20]:

$$\begin{array}{c}\hfill {\tilde{x}}_{+}=x-\frac{f\left(x\right)}{f^{\prime }\left(\frac{x+\tilde{x}}{2}\right)}\end{array}$$

If we call $1/2(x+{\tilde{x}}_{+})$ as $\widehat{x}$ and use the derivative at $\widehat{x}$ instead at x in the predictor step, then we have the following scheme:

$$\begin{array}{ccc}\hfill \widehat{x}& =& \frac{x+{\tilde{x}}_{+}}{2}\hfill \\ \hfill x_{+}& =& x-\frac{f\left(x\right)}{f^{\prime }\left(\widehat{x}\right)}\hfill \end{array}$$

In particular applications $f^{\prime }\left(\widehat{x}\right)$ is a good estimate of $f^{\prime }\left(x\right)$. Therefore this reduces the number of derivate evaluations. If the predictor-corrector scheme was not used, then the new point is $\tilde{x}$, and the function value would be $f\left(\tilde{x}\right)=f(x-f\left(x\right)/f^{\prime }\left(\widehat{x}\right))$. Since a corrector is used, the value of f at the new point is $f(x-f\left(x\right)/f^{\prime }\left({\widehat{x}}_{+}\right))$. If there is a favourable relationship between the derivatives of the predictor and the corrector points, then (assuming a minimization problem) we have:

$$\begin{array}{c}\hfill f(x-f\left(x\right)/f^{\prime }\left(\widehat{x}\right))>f(x-f\left(x\right)/f^{\prime }\left({\widehat{x}}_{+}\right))\end{array}$$

(11)

Suppose this inequality is valid in many iterative steps. In that case, there will be a reduction in the number of iterations compared to the Newton–Raphson method. The decrease in iterations required can be credited to the new approach’s ability to attain the optimal point of the Newton method at a rate of $1+\sqrt{2}$. This higher convergence rate means that the new method needs fewer iterations than the traditional approach.

In general, if the function is multivariate denoted F, then the updates can be written as:

$$\begin{array}{c}\hfill {\tilde{x}}_{+}=x-\alpha {(\nabla F\left(\widehat{x}\right))}^{-1}F\left(x\right)\\ \hfill {\widehat{x}}_{+}=1/2(x+{\tilde{x}}_{+})\\ \hfill x_{+}=x-\alpha {(\nabla F\left({\widehat{x}}_{+}\right))}^{-1}F\left(x\right)\end{array}$$

The algorithm presented in this paper is fundamentally distinct from the predictor-corrector method proposed by Mehrotra [17]. Notably, one of the main differences lies in the approach to calculating the search direction. While both algorithms leverage the two-stage search direction concept, Algorithm 1 computes the search direction by utilizing an auxiliary point and taking the average of this auxiliary point and the main point. Another significant contrast arises in the computation of the matrix inverse. Algorithm 1 requires calculating the matrix inverse only once per iteration. Consequently, it falls into the category of one-step methods in terms of time consumption, while achieving efficiency akin to two-step algorithms. These distinctions set Algorithm 1 apart from the predictor-corrector method of [17], highlighting its unique and promising characteristics in solving optimization problems. Moreover, the method in [17] requires tuning many parameters in each iteration, which can be expensive and complex for some problems to find their optimal values. In contrast, our proposed method involves far fewer parameters to tune, simplifying the process and potentially reducing the computational burden.

3. Experimental Evaluation

We use linear programming instances from NETLIB [26]. There are a total of 138 linear programming problems with varying sizes. The smallest one has 22 non-zeros and the largest one has 1,408,073 non-zeros in the coefficient matrix. There are 30 instances that are infeasible. Thirty nine of the instances are not full rank. Further, four test problems have $b=0$. If $b=0$, the initialization scheme [27] not produce an initial feasible solution. We do not use these 73 instances in our experimentation. We define the absolute and relative gap as follows:

$$\begin{array}{cc}\hfill \mathrm{absolute}\mathrm{gap}:& |c^{T}x-b^{T}y|\hfill \\ \hfill \mathrm{relative}\mathrm{gap}:& \frac{|c^{T}x-b^{T}y|}{1 + |c^{T}x| + |b^{T}y|}.\hfill \end{array}$$

For all experiments, we use a relative-gap of at most ${10}^{-6}$, at most 700 iterations as the stopping criteria. Although the absolute gap is generally more stringent than the relative-gap, achieving the absolute gap for problems with large objective function values may require many iterations. In IPMs, significant reductions occur during the initial iterations, but only small reductions happen in the final iterations. Therefore, in this article, we use the relative-gap in our experiments.

3.1. Experiments on CPU

We demonstrate the effectiveness of the new approach by executing Algorithm 1 and comparing the results to the classical interior-point algorithm (the primal-dual logarithmic barrier method Chapter 7) [8]. Note that instead of calculating the search direction using Newton’s method mentioned in the classical algorithm [8], Algorithm 1 finds the search direction using the modified Newton approach. This is the primary difference between Algorithm 1 and the classical algorithm. Specifically, in each iteration, Algorithm 1 first updates the auxiliary point. It then uses the auxiliary point, along with the primary principle, to find the search direction. In contrast, the classical algorithm relies solely on the information from the principal point to determine the search direction. In this regard, we coded the three algorithms (Algorithm 1, Algorithm 2 and Classical approach) in Python 3.10 using BLAS routines for the equation-solving step in all the algorithms. Of the remaining 65 instances the three algorithms solved 50 instances within a total time limit of 24 h on alliance canada cluster CEDAR, 2 × Intel E5-2683 v4 Broadwell CPU with 125 G RAM 2.1 GHz (https://alliancecan.ca/en, accessed on 27 August 2023) [28]. The remaining fifteen instances large-sized instances were not completed due to the time limit for at least one of the algorithms tested.

To demonstrate the efficacy of the new method, we applied Algorithms 1 and 2 to a test problem. For each step of the two-step method, we computed the objective function. Figure 1 illustrates the difference in the objective function between the two points, ${\tilde{x}}_{+}$ (auxiliary point) and $x_{+}$ (principal point), during each iteration. The plot shows that the objective function value in the second stage is consistently lower than in the first stage, indicated by positive values of $c^T{\tilde{x}}_{+}-c^T{x}_{+}$. This validates the proposed method, demonstrating its potential to reduce the number of iterations required. Additionally, it indicates that the auxiliary point effectively guides the algorithm in correcting its search direction, leading to more efficient convergence.

Next, we want to see for the same instance lp_scsd8, (i) the rate at which the absolute gap decreases, (ii) the decrease in the relative gap between primal and dual solution values as a function of the iterations. Figure 2 shows the two plots. Algorithms 1 and 2 achieve a faster reduction in the relative gap and the absolute gap compared to the classical approach, i.e., a lesser number of iterations are needed overall. Also note a sharp reduction in the latter set of iterations for Algorithms 1 and 2 (this phenomenon is not universal though, sometimes there is a sharp reduction in the initial iteration for some instances).

To demonstrate the effectiveness of the two-step method we observe two quantities, (i) the speedup and (ii) the relative reduction in the number of iterations achieved by Algorithms 1 and 2 compared to the classical algorithm. If Algorithm i takes time $t(A_i,j)$ on instance j, classical algorithm C takes time $t(C,j)$ on the same instance j then the speedup of Algorithm i is defined as $t(C,j)/t(A_i,j)$. Similarly, we define the relative reduction in the number of iterations. If Algorithm i takes $n(A_i,j)$ iterations on instance j, classical algorithm C takes $n(C,j)$ iterations on the same instance j then the relative reduction of Algorithm i is defined as $n(C,j)/n(A_i,j)$.

Figure 3 and Figure 4 show the relative reduction and the speedup in the number of iterations as a function of the number of instances. The instances are ordered according to the size in array. The average value of speedup, relative reduction if computed for each suffix of the array. If the $x$-coordinate is n then the y-coordinate is the average speedup, relative reduction for $50-n$ largest instances. The average speedup over all 50 instances is between 1.3 and 1.4 for both Algorithms 1 and 2. The average relative reduction is between 1.20 and 1.25 for both Algorithms 1 and 2. If we focus only on 5 of the largest instances then the speedup and relative reduction are comparably higher. This indicates that the method scales well. Algorithm 1 has a greater relative reduction but it does not translate into a greater speedup. Algorithm 2 has a better speedup. The instances indexed 40–45 in the plots show anomalous behaviour, it appears that these are specific types of instances (staircase instances) and the performance of Algorithms 1 and 2 seems to be different (compared to the average) on these instances. One of the prominent features of the new proposed algorithm is its ability to find an appropriate search direction compared to traditional interior-point algorithms. As demonstrated in [20], the new modified Newton method exhibits a faster convergence rate, approximately $1+\sqrt{2}$. Therefore, the new algorithm leverages this property, allowing it to reach the optimal solution in fewer iterations—a particularly valuable advantage when dealing with high-dimensional test problems. In scenarios where matrix multiplication and inversion can be computationally intensive, the reduced number of iterations provided by the new algorithm significantly decreases the execution time and enhances overall performance.

In

Table 1. The average number of iterations, CPU time (s) and standard deviation.

Methods	Aver. Iter.	Aver. CPU (s)	Std-Iter.
Classical Algorithm	97.82	166.71	70.96
Algorithm 1	71.34	136.29	55.64
Algorithm 2	72.52	135.46	57.55

Note: Bold numbers indicate the lowest values in each column.

, the average number of iterations and CPU time (measured in seconds) for three different methods—Classical Algorithm, Algorithms 1 and 2—are presented. The Classical Algorithm required an average of 97.82 iterations and 166.71 s of CPU time. In contrast, Algorithm 1 demonstrated improved efficiency with an average of 71.34 iterations and 136.29 s of CPU time. Similarly, Algorithm 2 also exhibited promising results with 72.52 iterations and 135.46 s of CPU time. These findings indicate that both Algorithms 1 and 2 are more effective than the Classical Algorithm in reducing the number of iterations and CPU time. Notably, Algorithm 1 achieved the lowest number of iterations, while Algorithm 2 achieved the lowest CPU time. Moreover, the fourth column presents the standard deviation for the number of iterations for the three algorithms, indicating that Algorithm 1 has the lowest variability, followed by Algorithm 2, with the Classical Algorithm having the highest variability.

The details of the problem instances can be found in [29]. The raw data used to build these figures is in the Appendix A and is self-explanatory. The data in Appendix B is divided into blocks of three rows each. The first row in each block is data for the classical algorithm and the other two rows are the data for Algorithms 1 and 2.

Table A2. Relative gap and the norm.

Name	Regap	$\parallel \mathit{Ax}-\mathit{b}\parallel$	$\parallel \mathit{c}-\mathit{s}-{\mathit{A}}^{\mathit{T}}\mathit{y}\parallel$
adlittle	$8.69\times {10}^{-7}$	$5.62\times {10}^{-14}$	$1.29\times {10}^{-10}$
adlittle	$8.82\times {10}^{-7}$	$3.32\times {10}^{-12}$	$1.29\times {10}^{-10}$
adlittle	$1.27\times {10}^{-7}$	$1.31\times {10}^{-13}$	$1.83\times {10}^{-5}$
afiro	$8.37\times {10}^{-7}$	$4.70\times {10}^{-17}$	$1.61\times {10}^{-4}$
afiro	$7.18\times {10}^{-7}$	$2.13\times {10}^{-9}$	$1.61\times {10}^{-4}$
afiro	$9.54\times {10}^{-7}$	$1.64\times {10}^{-9}$	$3.69\times {10}^{-4}$
agg	$1.41\times {10}^{-7}$	$2.38\times {10}^{-12}$	$1.78\times {10}^{-5}$
agg	$7.57\times {10}^{-7}$	$2.02\times {10}^{-11}$	$1.78\times {10}^{-5}$
agg	$7.57\times {10}^{-7}$	$2.02\times {10}^{-11}$	$1.78\times {10}^{-5}$
agg2	$3.23\times {10}^{-8}$	$9.15\times {10}^{-12}$	$3.34\times {10}^{-6}$
agg2	$7.21\times {10}^{-7}$	$5.27\times {10}^{-13}$	$3.34\times {10}^{-6}$
agg2	$2.03\times {10}^{-9}$	$3.29\times {10}^{-14}$	$3.37\times {10}^{-6}$
agg3	$5.27\times {10}^{-7}$	$1.13\times {10}^{-10}$	$2.82\times {10}^{-6}$
agg3	$5.77\times {10}^{-7}$	$4.40\times {10}^{-11}$	$2.82\times {10}^{-6}$
agg3	$2.62\times {10}^{-7}$	$1.01\times {10}^{-11}$	$3.24\times {10}^{-6}$
bandm	$8.64\times {10}^{-7}$	$7.02\times {10}^{-13}$	$1.76\times {10}^{-3}$
bandm	$1.50\times {10}^{-7}$	$4.43\times {10}^{-15}$	$1.76\times {10}^{-7}$
bandm	$7.66\times {10}^{-7}$	$2.86\times {10}^{-13}$	$8.31\times {10}^{-9}$
beaconfd	$6.09\times {10}^{-9}$	$7.13\times {10}^{-5}$	$3.75\times {10}^{-8}$
beaconfd	$4.26\times {10}^{-8}$	$1.25\times {10}^{-8}$	$3.75\times {10}^{-8}$
beaconfd	$2.405\times {10}^{-8}$	$6.97\times {10}^{-9}$	$1.52\times {10}^{-8}$
blend	$5.83\times {10}^{-7}$	$2.87\times {10}^{-16}$	$2.78\times {10}^{-6}$
blend	$5.25\times {10}^{-7}$	$6.43\times {10}^{-10}$	$2.78\times {10}^{-6}$
blend	$9.32\times {10}^{-7}$	$7.54\times {10}^{-14}$	$7.10\times {10}^{-5}$
bnl2	$9.79\times {10}^{-7}$	$1.65\times {10}^{-16}$	$3.37\times {10}^{-3}$
bnl2	$8.55\times {10}^{-7}$	$2.27\times {10}^{-16}$	$3.37\times {10}^{-3}$
bnl2	$6.19\times {10}^{-7}$	$1.97\times {10}^{-16}$	$2.66\times {10}^{-2}$
capri	$2.34\times {10}^{-9}$	$2.57\times {10}^{-16}$	$2.41\times {10}^{-7}$
capri	$1.57\times {10}^{-7}$	$1.45\times {10}^{-10}$	$2.41\times {10}^{-7}$
capri	$4.77\times {10}^{-7}$	$6.37\times {10}^{-10}$	$7.62\times {10}^{-6}$
chemcom	$3.88\times {10}^{-9}$	$1.48\times {10}^{-16}$	$6.08\times {10}^{-9}$
chemcom	$1.52\times {10}^{-7}$	$1.66\times {10}^{-12}$	$6.08\times {10}^{-9}$
chemcom	$1.98\times {10}^{-7}$	$1.72\times {10}^{-12}$	$6.11\times {10}^{-8}$
czprob	$9.42\times {10}^{-7}$	$1.01\times {10}^{-12}$	$1.59\times {10}^{-3}$
czprob	$8.18\times {10}^{-7}$	$8.94\times {10}^{-13}$	$1.59\times {10}^{-3}$
czprob	$7.07\times {10}^{-9}$	$2.98\times {10}^{-13}$	$3.81\times {10}^{-5}$
e226	$2.03\times {10}^{-7}$	$3.85\times {10}^{-15}$	$5.40\times {10}^{-4}$
e226	$2.95\times {10}^{-7}$	$2.60\times {10}^{-14}$	$5.40\times {10}^{-4}$
e226	$7.65\times {10}^{-7}$	$3.31\times {10}^{-16}$	$4.22\times {10}^{-3}$
fffff800	$9.51\times {10}^{-7}$	$1.16\times {10}^{-10}$	$3.45\times {10}^{-7}$
fffff800	$3.96\times {10}^{-7}$	$8.72\times {10}^{-11}$	$3.45\times {10}^{-7}$
fffff800	$1.31\times {10}^{-8}$	$1.65\times {10}^{-10}$	$6.25\times {10}^{-7}$
fit1p	$8.80\times {10}^{-7}$	$1.48\times {10}^{-16}$	$1.02\times {10}^{-5}$
fit1p	$1.25\times {10}^{-7}$	$1.39\times {10}^{-11}$	$1.02\times {10}^{-5}$
fit1p	$1.01\times {10}^{-7}$	$1.26\times {10}^{-11}$	$7.49\times {10}^{-10}$
forest6	$6.54\times {10}^{-8}$	$6.95\times {10}^{-16}$	$1.71\times {10}^{-6}$
forest6	$8.14\times {10}^{-7}$	$5.54\times {10}^{-12}$	$1.71\times {10}^{-6}$
forest6	$8.59\times {10}^{-7}$	$5.57\times {10}^{-12}$	$3.51\times {10}^{-5}$
ganges	$7.87\times {10}^{-7}$	$1.48\times {10}^{-16}$	$2.40\times {10}^{-7}$
ganges	$5.00\times {10}^{-7}$	$6.77\times {10}^{-12}$	$2.40\times {10}^{-7}$
ganges	$4.25\times {10}^{-7}$	$4.97\times {10}^{-12}$	$3.41\times {10}^{-8}$
gfrd_pnc	$1.88\times {10}^{-8}$	$9.77\times {10}^{-17}$	$1.24\times {10}^{-6}$
gfrd_pnc	$2.99\times {10}^{-7}$	$1.69\times {10}^{-16}$	$1.24\times {10}^{-6}$
gfrd_pnc	$2.11\times {10}^{-7}$	$1.57\times {10}^{-16}$	$3.98\times {10}^{-6}$
israel	$8.63\times {10}^{-7}$	$3.42\times {10}^{-16}$	$4.00\times {10}^{-6}$
israel	$2.37\times {10}^{-7}$	$7.63\times {10}^{-14}$	$4.00\times {10}^{-6}$
israel	$2.33\times {10}^{-7}$	$7.32\times {10}^{-14}$	$3.99\times {10}^{-6}$
lotfi	$1.07\times {10}^{-8}$	$2.34\times {10}^{-14}$	$1.89\times {10}^{-8}$
lotfi	$4.72\times {10}^{-7}$	$4.05\times {10}^{-12}$	$1.89\times {10}^{-8}$
lotfi	$2.23\times {10}^{-7}$	$1.84\times {10}^{-12}$	$5.32\times {10}^{-7}$
maros_r7	$2.36\times {10}^{-10}$	$2.94\times {10}^{-16}$	$7.53\times {10}^{-9}$
maros_r7	$6.67\times {10}^{-8}$	$2.92\times {10}^{-12}$	$7.53\times {10}^{-9}$
maros_r7	$5.15\times {10}^{-8}$	$2.059\times {10}^{-12}$	$1.82\times {10}^{-9}$
pilot87	$2.10\times {10}^{-5}$	$7.50\times {10}^{-12}$	$3.58\times {10}^{-7}$
pilot87	$3.86\times {10}^{-5}$	$2.46\times {10}^{-12}$	$3.58\times {10}^{-7}$
pilot87	$1.35\times {10}^{-5}$	$2.76\times {10}^{-12}$	$2.85\times {10}^{-7}$
pilot	$2.26\times {10}^{-7}$	$7.30\times {10}^{-13}$	$5.62\times {10}^{-7}$
pilot	$1.62\times {10}^{-6}$	$1.00\times {10}^{-13}$	$5.62\times {10}^{-7}$
pilot	$1.53\times {10}^{-6}$	$7.56\times {10}^{-14}$	$1.55\times {10}^{-6}$
pilot_we	$9.70\times {10}^{-7}$	$6.77\times {10}^{-15}$	$9.06\times {10}^{-4}$
pilot_we	$2.87\times {10}^{-7}$	$4.01\times {10}^{-15}$	$9.06\times {10}^{-4}$
pilot_we	$8.19\times {10}^{-7}$	$6.18\times {10}^{-15}$	$1.27\times {10}^{-3}$
sc50a	$9.73\times {10}^{-7}$	$7.74\times {10}^{-17}$	$7.90\times {10}^{-3}$
sc50a	$6.08\times {10}^{-7}$	$1.39\times {10}^{-12}$	$7.901\times {10}^{-3}$
sc50a	$7.08\times {10}^{-7}$	$7.57\times {10}^{-16}$	$3.92\times {10}^{-2}$
sc50b	$9.94\times {10}^{-7}$	$1.36\times {10}^{-16}$	$5.45\times {10}^{-3}$
sc50b	$6.46\times {10}^{-7}$	$3.58\times {10}^{-14}$	$5.45\times {10}^{-3}$
sc50b	$4.58\times {10}^{-7}$	$1.90\times {10}^{-16}$	$2.71\times {10}^{-3}$
sc105	$1.02\times {10}^{-7}$	$1.02\times {10}^{-16}$	$1.20\times {10}^{-2}$
sc105	$4.64\times {10}^{-7}$	$8.09\times {10}^{-13}$	$1.20\times {10}^{-2}$
sc105	$6.63\times {10}^{-7}$	$2.12\times {10}^{-16}$	$9.76\times {10}^{-3}$
sc205	$5.80\times {10}^{-9}$	$1.12\times {10}^{-16}$	$3.80\times {10}^{-6}$
sc205	$3.35\times {10}^{-7}$	$9.02\times {10}^{-13}$	$3.80\times {10}^{-6}$
sc205	$5.83\times {10}^{-7}$	$1.69\times {10}^{-16}$	$9.60\times {10}^{-3}$
scagr25	$6.25\times {10}^{-7}$	$1.32\times {10}^{-16}$	$2.79\times {10}^{-6}$
scagr25	$2.11\times {10}^{-7}$	$5.24\times {10}^{-12}$	$2.79\times {10}^{-6}$
scagr25	$2.80\times {10}^{-7}$	$8.46\times {10}^{-11}$	$6.52\times {10}^{-5}$
scfxm1	$9.01\times {10}^{-7}$	$1.82\times {10}^{-10}$	$4.43\times {10}^{-6}$
scfxm1	$1.74\times {10}^{-7}$	$9.71\times {10}^{-10}$	$4.43\times {10}^{-6}$
scfxm1	$9.94\times {10}^{-7}$	$1.45\times {10}^{-9}$	$7.66\times {10}^{-6}$
scfxm2	$7.87\times {10}^{-7}$	$5.19\times {10}^{-10}$	$1.62\times {10}^{-6}$
scfxm2	$4.25\times {10}^{-7}$	$1.41\times {10}^{-10}$	$1.62\times {10}^{-6}$
scfxm2	$9.46\times {10}^{-7}$	$1.25\times {10}^{-10}$	$3.06\times {10}^{-6}$
scfxm3	$9.61\times {10}^{-7}$	$8.06\times {10}^{-9}$	$7.08\times {10}^{-7}$
scfxm3	$1.57\times {10}^{-7}$	$3.51\times {10}^{-10}$	$7.08\times {10}^{-7}$
scfxm3	$1.60\times {10}^{-7}$	$3.86\times {10}^{-10}$	$6.77\times {10}^{-7}$
scrs8	$8.13\times {10}^{-7}$	$5.25\times {10}^{-13}$	$4.79\times {10}^{-8}$
scrs8	$7.91\times {10}^{-7}$	$1.84\times {10}^{-13}$	$4.79\times {10}^{-8}$
scrs8	$2.03\times {10}^{-7}$	$1.82\times {10}^{-14}$	$2.15\times {10}^{-8}$
scagr7	$9.96\times {10}^{-7}$	$1.33\times {10}^{-16}$	$1.01\times {10}^{-6}$
scagr7	$7.30\times {10}^{-7}$	$1.71\times {10}^{-16}$	$1.01\times {10}^{-6}$
scagr7	$3.32\times {10}^{-7}$	$4.26\times {10}^{-14}$	$6.08\times {10}^{-6}$
scsd1	$6.18\times {10}^{-7}$	$1.18\times {10}^{-16}$	$1.18\times {10}^{-7}$
scsd1	$6.90\times {10}^{-7}$	$1.04\times {10}^{-16}$	$1.18\times {10}^{-7}$
scsd1	$7.51\times {10}^{-7}$	$8.19\times {10}^{-17}$	$1.08\times {10}^{-5}$
scsd6	$7.89\times {10}^{-7}$	$1.56\times {10}^{-16}$	$3.14\times {10}^{-5}$
scsd6	$8.16\times {10}^{-7}$	$2.03\times {10}^{-16}$	$3.14\times {10}^{-5}$
scsd6	$6.87\times {10}^{-7}$	$1.83\times {10}^{-16}$	$3.64\times {10}^{-5}$
scsd8	$2.94\times {10}^{-8}$	$2.52\times {10}^{-16}$	$2.64\times {10}^{-5}$
scsd8	$3.12\times {10}^{-7}$	$2.76\times {10}^{-16}$	$2.64\times {10}^{-5}$
scsd8	$5.75\times {10}^{-7}$	$2.41\times {10}^{-16}$	$5.35\times {10}^{-6}$
sctap1	$5.77\times {10}^{-9}$	$2.84\times {10}^{-16}$	$4.66\times {10}^{-6}$
sctap1	$4.76\times {10}^{-7}$	$1.28\times {10}^{-9}$	$4.66\times {10}^{-6}$
sctap1	$2.18\times {10}^{-9}$	$4.64\times {10}^{-12}$	$3.66\times {10}^{-8}$
sctap2	$6.17\times {10}^{-9}$	$5.70\times {10}^{-11}$	$1.33\times {10}^{-7}$
sctap2	$1.49\times {10}^{-7}$	$4.38\times {10}^{-12}$	$1.33\times {10}^{-7}$
sctap2	$1.79\times {10}^{-7}$	$1.94\times {10}^{-11}$	$2.16\times {10}^{-6}$
sctap3	$1.96\times {10}^{-8}$	$8.318\times {10}^{-11}$	$2.61\times {10}^{-6}$
sctap3	$5.08\times {10}^{-7}$	$2.86\times {10}^{-10}$	$2.61\times {10}^{-6}$
sctap3	$2.63\times {10}^{-7}$	$1.41\times {10}^{-10}$	$2.41\times {10}^{-6}$
share1b	$9.68\times {10}^{-7}$	$1.27\times {10}^{-16}$	$1.40\times {10}^{-4}$
share1b	$9.29\times {10}^{-8}$	$3.31\times {10}^{-15}$	$1.40\times {10}^{-4}$
share1b	$8.53\times {10}^{-8}$	$2.59\times {10}^{-15}$	$1.38\times {10}^{-5}$
share2b	$9.60\times {10}^{-7}$	$6.11\times {10}^{-15}$	$4.34\times {10}^{-6}$
share2b	$8.65\times {10}^{-7}$	$1.23\times {10}^{-9}$	$4.34\times {10}^{-6}$
share2b	$1.14\times {10}^{-7}$	$3.10\times {10}^{-10}$	$1.41\times {10}^{-5}$
stair	$9.30\times {10}^{-7}$	$1.90\times {10}^{-11}$	$1.48\times {10}^{-5}$
stair	$6.96\times {10}^{-7}$	$2.04\times {10}^{-10}$	$1.48\times {10}^{-5}$
stair	$2.84\times {10}^{-7}$	$1.22\times {10}^{-10}$	$2.15\times {10}^{-5}$
standata	$7.12\times {10}^{-7}$	$2.66\times {10}^{-12}$	$3.02\times {10}^{-12}$
standata	$3.87\times {10}^{-7}$	$1.22\times {10}^{-13}$	$3.02\times {10}^{-12}$
standata	$9.88\times {10}^{-7}$	$2.35\times {10}^{-15}$	$4.45\times {10}^{-11}$
standmps	$9.68\times {10}^{-7}$	$8.69\times {10}^{-9}$	$5.74\times {10}^{-7}$
standmps	$1.75\times {10}^{-7}$	$1.10\times {10}^{-9}$	$5.74\times {10}^{-7}$
standmps	$3.92\times {10}^{-8}$	$6.92\times {10}^{-11}$	$4.53\times {10}^{-7}$
stocfor1	$1.20\times {10}^{-8}$	$2.06\times {10}^{-16}$	$5.68\times {10}^{-6}$
stocfor1	$1.30\times {10}^{-7}$	$2.69\times {10}^{-14}$	$5.68\times {10}^{-6}$
stocfor1	$1.96\times {10}^{-7}$	$1.98\times {10}^{-15}$	$1.28\times {10}^{-5}$
stocfor2	$1.39\times {10}^{-10}$	$2.15\times {10}^{-16}$	$6.40\times {10}^{-7}$
stocfor2	$5.61\times {10}^{-10}$	$6.68\times {10}^{-16}$	$6.40\times {10}^{-7}$
stocfor2	$2.08\times {10}^{-8}$	$6.16\times {10}^{-15}$	$5.51\times {10}^{-6}$
truss	$5.72\times {10}^{-11}$	$1.40\times {10}^{-15}$	$6.37\times {10}^{-7}$
truss	$2.12\times {10}^{-9}$	$6.87\times {10}^{-16}$	$6.37\times {10}^{-7}$
truss	$1.99\times {10}^{-9}$	$6.67\times {10}^{-16}$	$1.22\times {10}^{-5}$
lp-vtp-base	$9.02\times {10}^{-9}$	$2.12\times {10}^{-9}$	$3.55\times {10}^{-6}$
lp-vtp-base	$3.82\times {10}^{-7}$	$4.49\times {10}^{-9}$	$3.55\times {10}^{-6}$
lp-vtp-base	$5.09\times {10}^{-7}$	$2.42\times {10}^{-8}$	$4.26\times {10}^{-6}$
woodw	$3.85\times {10}^{-5}$	$2.01\times {10}^{-4}$	$4.05\times {10}^{-4}$
woodw	$4.29\times {10}^{-5}$	$1.51\times {10}^{-5}$	$4.05\times {10}^{-4}$
woodw	$1.33\times {10}^{-4}$	$7.42\times {10}^{-6}$	$2.09\times {10}^{-4}$

lists the relative gap and the norms and shows the convergence of the three algorithms to an optimal solution on several instances from NETLIB.

3.2. Experiments on GPU

A study by Świrydowicz et al. [30] investigated the use of GPU-accelerated linear solvers for optimizing large-scale power grid problems. Motivated by their study, we evaluate the performance of the algorithms on massively parallel hardware (GPU with tensor cores) using library routines to perform factoring.

We utilized the Pytorch library to develop a GPU-enabled version of our program. The library assumes that all tensors are located on the same device. During our experiments, we utilized the NVIDIA V100 Volta (on the GRAHAM Compute Canada cluster), which has a maximum size limitation due to its 16GB HBM2 memory. This card features tensor cores that perform certain tensor operations very efficiently in hardware. GPU arithmetic is typically conducted in low precision, ranging from 8 to 32 bits which account for the speed. Unfortunately, the algorithms failed to converge to an optimal solution when we employed low-precision arithmetic on V100. As a result, we opted to use double precision at 64 bits, even though this substantially increased the run time.

To utilize GPU acceleration, we convert all vectors and matrices into tensors, leveraging the PyTorch library instead of Numpy for computations. In contrast, when running on a CPU, we represent these structures as regular vectors and matrices. By performing essential operations like matrix multiplication, inverse calculation and step size computation in tensor mode. This approach enhances the efficiency of calculations and enables faster reaching the optimal point.

The computationally expensive steps in the algorithm proposed here (as in any IPM) are (i) the determination of $\alpha$ and (ii) the computation of the inverse of the Jacobian to determine the step direction. The first expensive step has quadratic complexity ($O\left({(n+m)}^2\right)$ when implemented naturally). The second step is computationally more costly $O\left({(n+m)}^3\right)$. In the approach here, gains are realized asymptotically in the second step. Therefore, we focus on the time the GPU takes to perform matrix inversions. The raw data is shown in Appendix C. The second and third columns are the iterations needed by the classical and the algorithm proposed here. The last two columns in the table are the times (in seconds) needed by the GPU to perform the matrix inversions.

Results for GPU

Table A3. Results for GPU.

	Iteration	Iteration	GPU Time (s)	GPU Time (s)
Name	Classical	Algorithm 1	Classical	Algorithm 1
adlittle	49	31	0.151	0.099
afiro	19	17	0.028	0.027
agg	74	60	1.635	1.339
agg2	84	47	2.108	1.181
agg3	81	54	2.108	1.378
bandm	700	105	10.05	1.529
beaconfd	700	50	5.318	0.392
blend	53	26	0.148	0.073
bnl2	222	236	182.065	192.634
capri	67	51	0.963	0.733
chemcom	76	47	1.719	1.075
d2q06c	306	166	467.293	252.344
fit1p	41	35	3.137	2.670
forest6	42	22	0.125	0.071
ganges	700	41	71.770	4.190
israel	157	128	1.311	1.109
lotfi	68	54	0.642	0.522
pilot87	246	181	462.709	339.440
pilot	128	92	105.27	75.457
pilot_we	399	269	84.727	56.946
sc50a	46	20	0.089	0.042
sc50b	28	20	0.052	0.039
sc105	21	21	0.085	0.082
sc205	45	23	0.397	0.207
scagr25	47	37	1.133	0.890
scfxm1	75	51	1.399	0.960
scfxm2	72	53	3.480	2.551
scfxm3	76	56	7.496	5.498
scagr7	32	26	0.141	0.117
scsd1	26	27	0.526	0.551
scsd6	24	25	1.071	1.110
scsd8	37	22	6.320	3.732
sctap1	60	52	1.229	1.073
3sctap2	90	47	16.300	8.485
sctap3	105	61	36.785	21.257
share1b	700	69	4.096	0.417
share2b	38	27	0.144	0.104
stair	700	82	13.832	1.640
standmps	183	700	8.892	34.161
stocfor1	67	38	0.255	0.151
lp-vtp-base	558	37	0.521	0.354

presents the results for Classical Algorithm and Algorithm 1 on GPU. We report the number of iterations and GPU times for 41 test problems. Note that, Classical algorithm could not reach the optimal solution for 5 test problems and Algorithm 1 could not find the optimal solution for 1 test problem based on the stopping conditions. We remove the results for these 6 test problems and calculate the relative reduction and speed-up. Figure 5 presents the relative reduction and speed-up for 35 test problems. Moreover, we denote the average of the number of iterations and GPU time for these 35 test problems in

Table 2. The average number of iterations, GPU time and standard deviation for the number of iterations.

Methods	Aver. Iter.	Aver. GPU (s)	Std-Iter.
Classical Algorithm	101.77	39.81	113.71
Algorithm 1	61.68	27.83	60.19

Note: Bold numbers indicate the lowest values in each column.

. Based on this Table, we can conclude that the new proposed method reduced the number of iterations and GPU by $39\%$ and $30\%$ respectively. Moreover, the standard deviation for the number of iterations indicates that Algorithm 1 has the lowest variability at 60.19.

3.3. Discussion

Although it has not been theoretically proven that the convergence rate of the modified Newton method is exactly $1+\sqrt{2}$, numerical results show that this rate is achievable. This rate may vary slightly depending on the structure of the problem (see example in [20]). As shown in [20], the modified Newton method can reduce the number of iterations and computational time for solving a system of equations by a factor of 1.22, calculated as $(1+\sqrt{2})/2$, where 2 is the convergence rate of the Newton method and $1+\sqrt{2}$ is the convergence rate of the modified Newton method. However, the average results in Table 1 and Table 2 demonstrate that when we apply the modified Newton method to find the search direction in IPMs, the convergence rate of the proposed method remains close to $1+\sqrt{2}$.

4. Conclusions

This paper presents a two-step method for optimizing linear programming problems using the interior point method. Building on a previous idea to solve the Karush–Kuhn–Tucker (KKT) conditions, we extend this approach to IPMs. Our method re-examines and enhances the classical algorithm proposed by [8] by incorporating a two-step process. In this method, the search direction is calculated using a modified Newton method that employs an auxiliary point. Like the standard Newton method, the modified Newton method computes the matrix inverse only once. We incorporate this approach into IPM and demonstrate its efficiency through implementations on both CPU and GPU, as well as by solving NETLIB test problems. The results show that the proposed method can improve efficiency by $27\%$ in terms of iterations and $18\%$ in terms of CPU time. These results confirm that the convergence rate of the modified Newton method is equal to $1+\sqrt{2}$. This work is pioneering in integrating the modified Newton method into IPM. This approach has the potential to be applied to other interior point algorithms, enhancing their performance as well.

One of the key areas for future work is proving the complexity bounds of the algorithm, which could enhance the overall complexity bounds of IPMs. Additionally, demonstrating the effectiveness of the algorithm on real-world problems modeled as linear optimization problems is another important avenue for exploration. Another focus will be evaluating the algorithm’s efficiency on large datasets, which is particularly relevant given the increasing dimensions of real-world problems. Further, we plan to optimize the implementation on GPUs, which includes developing parallel implementations of the algorithm and leveraging parallel processing to calculate the inverse of the Jacobian matrix. This work can yield significant performance improvements and provide valuable insights into the practical applications of the algorithm.