FRACTIONAL ORDER LOG BARRIER INTERIOR POINT ALGORITHM FOR POLYNOMIAL REGRESSION IN THE ℓ _p -NORM

ABSTRACT

Fractional calculus is the branch of mathematics that studies the several possibilities of generalizing the derivative and integral of a function to noninteger order. Recent studies found in literature have confirmed the importance of fractional calculus for minimization problems. However, the study of fractional calculus in interior point methods for solving optimization problems is still new. In this study, inspired in applications of fractional calculus in many fields, was developed the so-called fractional order log barrier interior point algorithm by replacing some integer derivatives for the corresponding fractional ones on the first order optimality conditions of Karush-Kuhn-Tucker to solve polynomial regression models in the ℓ _p −norm for 1 < p < 2. Finally, numerical experiments are performed to illustrate the proposed algorithm.

Keywords:
nonlinear programming; polynomial regression; ℓp−norm; interior point method; fractional derivative

1 INTRODUCTION

It is fundamentally important to make predictions based upon scientific data. The problem of fitting curves to data points has many practical applications Bard (1974); Sevaux & Mineur (2007); Chatterjee & Hadi (2012).

Given a set of m data points in ${\mathrm{ℝ}}^2:{\left\{\left(a_i,b_i\right)\right\}}_{i=1}^m$, where a _i is an argument value and b _i a corresponding dependent value, with $a_i\ne a_j$ for all $i\ne j$. The curve fitting procedure tries to build a linear or nonlinear function y = f (x), defined for all possible choices of x, that approximately fits the data set. The fitted curves to the data by f are most often chosen to be polynomials Süli & Mayers (2003).

Let y = f (x) be a polynomial function of degree n− 1 of the form $f\left(x\right)=x_0+x_{1}x+\dots +x_{n-1}{x}^{n-1}$, the candidate function to fit data. The procedure to fit the polynomial function to data, in the ℓ _p −norm, determines the vector $x={\left(x_0,x_1,x_2,\dots ,x_{n-1}\right)}^{⊺}\in {\mathrm{ℝ}}^n$, where the superscript T represents transpose, that minimizes the ℓ _p −norm of the residual error as follows

$\min \Phi \left(x\right)={||b-Ax||}_p,$ (1)

where ${||·||}_p$ denotes ℓ _p −norm, $b={\left(b_1,b_2,\dots ,b_m\right)}^{⊺}\in {\mathrm{ℝ}}^m$ and $A\in {\mathrm{ℝ}}^{m\times n}$ is a Vandermonde matrix, with rank (A) = m, which can be written as

$A=\left(\begin{array}{ccccc}1& a_1& {a_1}^2& \cdots & {a_1}^{n-1}\\ 1& a_2& {a_2}^2& \cdots & {a_2}^{n-1}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& a_m& {a_m}^2& \cdots & {a_m}^{n-1}\end{array}\right).$ (2)

This is a nonlinear (or linear) regression model if n > 2 (or n = 2). If m ≤ n, there is an (n − 1)th degree polynomial satisfying $f\left(a_i\right)=b_i\left(i=1,2,\dots ,m\right)$. If m > n, the problem cannot be exactly solved and one needs to find the vector x in (1).

In many industrial applications, missing data or anomalous values can arise as errors of measurement instruments during the data generation. In this cases, polynomial regression models in the ℓ _p −norm, with $p\ne 2$, are more robust Forsythe (1972). It is important to choose an appropriate value for p and several criteria for the choice of p have been studied Rice & White (1964).

The calculus of integrals and derivatives of arbitrary order, named as fractional calculus, was conceptualized in connection with the infinitesimal calculus since 1695 Oldham & Spanier (1974). Some of the application areas include: viscoelasticity, signal processing, probability, statistics, electrochemistry, diffusion in porous media, fluid flow, backpropagation training of neural networks, fuzzy control method, and so on Kilbas et al. (2006); Dalir & Bashour (2010); Mohammadzadeh & Kayacan (2020); Wang et al. (2017a); Chen et al. (2017a); Grigoletto & Oliveira (2020). Also, its importance in minimization problems has been confirmed by recent works found in literature. For example, Chen et al. Chen et al. (2017b) presented the fractional order gradient methods (FOGMs) by writing the Riemann-Liouville and Caputo fractional derivatives as Taylor series, and Wang et al. Wang et al. (2017b) proposed the fractional gradient descent method employing the Caputo derivative for the backpropagation training of neural networks. However, the study of fractional calculus in the field of interior point methods for solving optimization problems is still new.

The goal of this study is to investigate the fractional order log barrier interior point algorithm, involving the Caputo fractional derivative. It is based on replacing some integer derivatives for the corresponding fractional ones on the first order optimality conditions for solving polynomial regression models in the ℓ _p −norm for 1 < p < 2. Functions of the form $g\left(x\right)={\left(x+\kappa \right)}^p$, with κ > 0 and 1 < p < 2, arise when solving problem 1), according to the approach which will be discussed throughout this study. The Caputo fractional derivative, in addition to generalize the integer order derivative, is a useful tool to obtain the derivative for functions as g(x), as p is on the interval 1 < p < 2, then it will be a non integer number.

This paper is organized as follows. Preliminary concepts of fractional calculus are presented in Section 2. In Section 3, the fractional order log barrier interior point algorithm for solving polynomial regression models in the ℓ _p −norm is discussed. Numerical experiments are performed to illustrate the proposed algorithm in Section 4, and the Section 5 contains the conclusions.

2 PRELIMINARY CONCEPTS OF FRACTIONAL CALCULUS

Some basic concepts and definitions involving special functions and the Caputo fractional derivatives will be presented in this section.

Definition 1. The gamma function Γ(z) is originally defined by Erdélyi et al. (1981); Andrews et al. (1999):

$\Gamma \left(z\right)={\int }_0^{\infty }{t}^{z-1}{e}^{-t}\text{d}t\left(z>0\right).$ (3)

For gamma function, the reduction formula

$\Gamma \left(z+1\right)=z\Gamma \left(z\right),$ (4)

holds. In particular, when $z=n\in \mathrm{\mathbb{N}}$, then $\Gamma \left(n+1\right)=n!$.

Definition 2. The beta function B(z, w) is defined by Erdélyi et al. (1981); Andrews et al. (1999):

$\text{B}\left(z,w\right)={\int }_0^{\infty }{\xi }^{z-1}{\left(1+\xi \right)}^{-z-w}\text{d}\xi \left(z>0;w>0\right).$ (5)

It is connected with the gamma function by the following relation

$\text{B}\left(z,w\right)=\frac{\Gamma \left(z\right)\Gamma \left(w\right)}{\Gamma \left(z+w\right)}\left(z>0;w>0\right).$ (6)

The fractional calculus basically defines several integral and derivative operators of arbitrary order. The Caputo fractional derivative is one of the fractional derivative operators Kilbas et al. (2006).

Definition 3. The Caputo left-sided fractional derivative with respect to x, of order 0 < α < 1, of functions f on a subset Ω = [a, b] of real axis $\mathrm{ℝ}=\left(-\infty ,\infty \right)$, where $f\in {\text{C}}^1\left[a,b\right]$, denoted by $\left({}_x{}^{\mathrm{C}}\mathrm{D}_{a+}^{\alpha }f\right)\left(x\right)$, is given by

$\left({}_x{}^{\mathrm{C}}\mathrm{D}_{a+}^{\alpha }f\right)\left(x\right):=\frac{1}{\Gamma (1-\alpha )}{\int }_a^x{\left(x-\tau \right)}^{-\alpha }f\text{'}\left(\tau \right)\text{d}\tau \left(x\in \Omega \right).$ (7)

Note that the Caputo fractional derivative is a nonlocal fractional derivative, because it depends on the choice of the order α, function f (x), x, and also relies on the total effects on the interval [a, x]. Usually, this is called memory effect.

In particular, $\left({}_x{}^{\mathrm{C}}\mathrm{D}_{a+}^{1}f\right)\left(x\right):=f\text{'}\left(x\right)$.

Property 1. Let $g\left(x\right)={\left(x+\kappa \right)}^p$ be defined for $-\kappa \le x<\infty$, with κ ≥ 0, and 1 < p < 2, and let 0 < α < 1. Then

$\left[{}_x{}^{\mathrm{C}}\mathrm{D}_{-\kappa +}^{\alpha }{\left(x+\kappa \right)}^p\right]\left(x\right)=\frac{\Gamma \left(p+1\right)}{\Gamma \left(p-\alpha +1\right)}{\left(x+\kappa \right)}^{p-\alpha }.$ (8)

Proof. Applying the Caputo left-sided fractional derivative (7) to the function g(x) with a = −κ,

$\left[{}_x{}^{\mathrm{C}}\mathrm{D}_{-\kappa +}^{\alpha }{\left(x+\kappa \right)}^p\right]\left(x\right)\frac{1}{\Gamma (1-\alpha )}{\int }_{-\kappa }^x{\left(x-\tau \right)}^{-\alpha }p{\left(\tau +\kappa \right)}^{p-1}\text{d}\tau .$ (9)

Using the change of variables $\xi =\frac{\tau +\kappa }{x-\tau }$, then:

$\tau =\frac{\xi x-\kappa }{1+\xi },d\tau =\frac{\left(x+\kappa \right)}{{\left(1+\xi \right)}^2}d\xi ,\tau \to -\kappa \Rightarrow \xi \to 0\mathrm{and}\tau \to x\Rightarrow \xi \to \infty .$

The integral (9) under this change of variables and by means of (4)-(6) takes the form

$$\begin{gathered} \left[{}_x{}^{C}D_{-\kappa +}^{\alpha }{\left(x+\kappa \right)}^p\right]\left(x\right)=\frac{p}{\Gamma (1-\alpha )}{\int }_{-\kappa }^x{\left(x-\tau \right)}^{-\alpha }{\left(\tau +\kappa \right)}^{p-1}\text{d}\tau \\ =\frac{p}{\Gamma (1-\alpha )}{\int }_0^{\infty }{\left(\frac{x+\kappa }{1+\xi }\right)}^{-\alpha }{\left(\frac{\xi (x+\kappa )}{1+\xi }\right)}^{p-1}\frac{\left(x+\kappa \right)}{{\left(1+\xi \right)}^2}\text{d}\xi \\ =\frac{p{\left(x+\kappa \right)}^{p-\alpha }}{\Gamma (1-\alpha )}{\int }_0^{\infty }{\xi }^{p-1}{\left(1+\xi \right)}^{\alpha -p-1}\text{d}\xi \\ =\frac{p{\left(x+\kappa \right)}^{p-\alpha }}{\Gamma (1-\alpha )}\text{B}\left(p,1-\alpha \right) \\ =\frac{p{\left(x+\kappa \right)}^{p-\alpha }}{\Gamma (1-\alpha )}\frac{\Gamma \left(p\right)\Gamma \left(1-\alpha \right)}{\Gamma \left(p-\alpha +1\right)} \\ =\frac{\Gamma \left(p+1\right)}{\Gamma \left(p-\alpha +1\right)}{\left(x+\kappa \right)}^{p-\alpha }. \end{gathered}$$

□

3 FRACTIONAL ORDER LOG BARRIER INTERIOR POINT ALGORITHM

Raising $\Phi \left(x\right)={||b-Ax||}_p$ to the power p, the problem (1) can be rewritten in the alternative form as follows

$\min {||b-Ax||}_p^p,$ (10)

where x ∈ ℝⁿ . The problem (10) is similar to the following nonlinear optimization problem

$$\begin{gathered} \min {||r||}_p^p \\ \text{s.t.}Ax+r=b, \\ x\in {\mathrm{ℝ}}^n, \end{gathered}$$ (11)

where $r={\left(r_1,r_2,\dots ,r_m\right)}^{⊺}\in {\mathrm{ℝ}}^m$ is the residual vector of regression and ℓ _p −norm is defined of the form ${||r||}_p={\left(\sum _{i=1}^m{\left|r_i\right|}^p\right)}^{\frac{1}{p}}\left(1\le p<\infty \right)$ and ${||r||}_{\infty }=\max \left|r_i\right|$.

The parameter values p most commonly used to the problem (11) are p = 1, p = 2 and p → ∞. Linear programming procedures Charnes et al. (1955); Oliveira et al. (2000); Oliveira & Lyra (2004) can be used for p = 1 and p → ∞ cases. For p = 2, the direct solution is given by x = (A ^⊺ A)⁻¹ A ^⊺ b. For other values of p, unconstrained minimization procedures can be used Dennis Jr. & Schnabel (1983); Li (1993); Cantante et al. (2012).

Considering the unrestricted residual term r _i as the difference between two nonnegative variables u _i and v _i: r _i = u _i −v _i , with u _i and v _i being defined by

$u_i=\left\{\begin{array}{l}{r}_i,\text{if}{r}_i\ge 0,\\ 0,\text{if}{r}_i<0,\end{array}\right.v_i=\left\{\begin{array}{l}0,\text{if}{r}_i\ge 0,\\ -r_i,\text{if}{r}_i<0,\end{array}\right.$ (12)

for all i = 1, 2, . . . , m, then |r _i | = u _i + v _i and the problem (11) can be converted into a convex programming problem Charnes et al. (1955); Cantante et al. (2012):

$$\begin{gathered} \min \sum _{i=1}^m{\left(u_i+v_i\right)}^p \\ \text{s.t.}Ax+u-v=b, \\ x\in {\mathrm{ℝ}}^n,u\in {\mathrm{ℝ}}_{+}^m,v\in {\mathrm{ℝ}}_{+}^m, \end{gathered}$$ (13)

where ${\mathrm{ℝ}}_{+}^m$ denote the set of m-dimensional nonnegative vectors.

Interior point methods are widely used to solve convex optimization problems because of their good performance in practice Biegler (2010); Gondzio (2012); Lilian et al. (2016). By adding a logarithmic barrier function to the objective function in (13), the barrier problem is given by

$$\begin{gathered} \min \sum _{i=1}^m\left[{\left(u_i+v_i\right)}^p-\mu \ln u_i-\mu \ln v_i\right] \\ \text{s.t.}Ax+u-v=b, \\ x\in {\mathrm{ℝ}}^n, \end{gathered}$$ (14)

where µ > 0 is the barrier parameter. An optimal solution of (13) can be found by solving a series of barrier problems of the form (14) while µ is decreasing and going to zero. The Lagrangian function associated with the problem (14) is

$\text{L}\left(x,\lambda ,u,v\right)=\sum _{i=1}^m\left[{\left(u_i+v_i\right)}^p-\mu \ln u_i-\mu \ln v_i\right]+{\lambda }^{⊺}\left(Ax+u-v-b\right),$ (15)

where $\lambda ={\left({\lambda }_1,{\lambda }_2,\dots ,{\lambda }_m\right)}^{⊺}\in {\mathrm{ℝ}}^m$ is a Lagrange multiplier vector.

Let ϕ and φ be the functions given by

$\varphi \left(x,\lambda ,u,v\right)=\sum _{i=1}^m\left[{\left(u_i+v_i\right)}^p\right],$ (16)

and

$\phi \left(x,\lambda ,u,v\right)=-\mu \sum _{i=1}^m\left(\ln u_i+\ln v_i\right)+{\lambda }^{⊺}\left(Ax+u-v-b\right),$ (17)

then the Lagrangian function can be written as

$\text{L}\left(x,\lambda ,u,v\right)=\varphi \left(x,\lambda ,u,v\right)+\phi \left(x,\lambda ,u,v\right).$ (18)

The fractional order log barrier interior point algorithm, involving the Caputo fractional derivative, is based on replacing some integer derivatives for the corresponding fractional ones on the first order optimality conditions (𝛻L = 𝛻ϕ + 𝛻φ = 0), of the following form

${\nabla }^{\alpha }\varphi +\nabla \phi =\left[\begin{array}{c}{\nabla }_x\varphi +{\nabla }_x\phi \\ {\nabla }_{\lambda }\varphi +{\nabla }_{\lambda }\phi \\ {\nabla }_u^{\alpha }\varphi +{\nabla }_u\phi \\ {\nabla }_v^{\alpha }\varphi +{\nabla }_v\phi \end{array}\right]=\left[\begin{array}{c}{A}^{⊺}\lambda \\ Ax+u-v-b\\ g_{\alpha }-\mu U^{-1}{1}_m+\lambda \\ g_{\alpha }-\mu V^{-1}{1}_m-\lambda \end{array}\right]=\left[\begin{array}{r}0\\ 0\\ 0\\ 0\end{array}\right],$ (19)

where 𝛻_x , 𝛻_λ , 𝛻_u and 𝛻_v represent the gradient with respect to x, λ , u and v, respectively, ${\nabla }_u^{\alpha }$ and ${\nabla }_v^{\alpha }$ represent the fractional order gradient with respect to u and v, respectively, given by Caputo fractional derivative (7).

Furthermore, $1_m={\left(1,1,\dots ,1\right)}^{⊺}\in {\mathrm{ℝ}}^m,U=\text{diag}\left(u\right),V=\text{diag}\left(v\right)$, where diag(w) denote the diagonal matrix from a vector w, and $g_{\alpha }\in {\mathrm{ℝ}}^m$. The ith component of the vector g _α is ${\left(g_{\alpha }\right)}_i=\left[{}_{u_i}{}^C{\text{D}}_{-v_i+}^{\alpha }{\left(u_i+v_i\right)}^p\right]\left(u_i\right)=\left[{}_{v_i}{}^C{\text{D}}_{-u_i+}^{\alpha }{\left(u_i+v_i\right)}^p\right]\left(v_i\right)$. By Property 1, for i = 1, 2, . . . , m, (g _α )_i is given by

${\left(g_{\alpha }\right)}_i=\frac{\Gamma \left(p+1\right)}{\Gamma \left(p+1-\alpha \right)}{\left(u_i+v_i\right)}^{p-\alpha }\left(1<p<2;0<\alpha <1\right).$ (20)

The nonlinear system of equations (19) can be rewritten in the alternative form as follows

$\left[\begin{array}{c}{A}^{⊺}\lambda \\ Ax+u-v-b\\ U\left(g_{\alpha }+\lambda \right)\\ V\left(g_{\alpha }-\lambda \right)\end{array}\right]=\left[\begin{array}{r}0\\ 0\\ {\mu }_m\\ {\mu }_m\end{array}\right],$ (21)

where ${\mu }_m={\left(\mu ,\mu ,\dots ,\mu \right)}^{⊺}\in {\mathrm{ℝ}}^m$.

For given µ > 0, if α = 1 in the system (21), then it recovers the first order optimality conditions, since $\underset{\alpha \to 1}{\lim }\left({\nabla }^{\alpha }\varphi +\nabla \phi \right)=\nabla \varphi +\nabla \phi =\nabla \text{L}$, and has a unique solution $\left(x_{\mu }^{*},{\lambda }_{\mu }^{*},u_{\mu }^{*},v_{\mu }^{*}\right)$. When the classical gradient is replaced by the fractional one (0 < α < 1), then if the system (21) has a solution $\left({x_{\mu }^{\alpha }}^{*},{{\lambda }_{\mu }^{\alpha }}^{*},{u_{\mu }^{\alpha }}^{*},{v_{\mu }^{\alpha }}^{*}\right)$, it can be called fractional solution and $\underset{\alpha \to 1}{\lim }\left({x_{\mu }^{\alpha }}^{*},{{\lambda }_{\mu }^{\alpha }}^{*},{u_{\mu }^{\alpha }}^{*},{v_{\mu }^{\alpha }}^{*}\right)=\left(x_{\mu }^{*},{\lambda }_{\mu }^{*},u_{\mu }^{*},v_{\mu }^{*}\right)$.

For given µ > 0 and 0 < α ≤ 1, suppose that the system (21) has a solution. Given an initial point (x ⁰, λ ⁰, u ⁰, v ⁰) such that $u^0\in {\mathrm{ℝ}}_{++}^m,v^0\in {\mathrm{ℝ}}_{++}^m$, where ${\mathrm{ℝ}}_{++}^m$ denote the set of m-dimensional positive vectors, and Ax ⁰ + u ⁰ −v ⁰ = b, by applying Newton method to the system (21), then the search direction (Δx, Δλ, Δu, Δv) can be found by solving the following Newton system

$\left[\begin{array}{cccc}0& A^{⊺}& 0& 0\\ A& 0_m& I_m& -I_m\\ 0& U& D_u& U{H}_{\alpha }\\ 0& -V& V{H}_{\alpha }& D_v\end{array}\right]\left[\begin{array}{r}\Delta x\\ \Delta \lambda \\ \Delta u\\ \Delta v\end{array}\right]=\left[\begin{array}{r}{r}_1\\ r_2\\ r_3\\ r_4\end{array}\right],$ (22)

where

$\begin{array}{rcl}{r}_1& =& -A^{⊺}\lambda ,\\ r_2& =& -Ax-u+v+b,\\ r_3& =& -U\left(g_{\alpha }+\lambda \right)+{\mu }_m,\\ r_4& =& -V\left(g_{\alpha }-\lambda \right)+{\mu }_m,\end{array}$ (23)

I_m is the identity matrix of order m, D _u = G _α + λI _m +UH _α , D _v = G _α −λI _m +VH _α , where G _α = diag(g _α ), H _α = diag(h _α ), where h _α ∈ ℝ_m and for i = 1, 2, . . . , m, the ith component of the vector h _α is given by

${\left(h_{\alpha }\right)}_i=\frac{\Gamma \left(p+1\right)}{\Gamma \left(p-\alpha \right)}{\left(u_i+v_i\right)}^{p-\alpha -1}.$ (24)

The (h _α )_i component is obtained by evaluating the derivative

${\left(h_{\alpha }\right)}_i=\frac{\partial }{\partial u_i}\left[{\left(g_{\alpha }\right)}_i\right]=\frac{\left(p-\alpha \right)\Gamma \left(p+1\right)}{\Gamma \left(p+1-\alpha \right)}{\left(u_i+v_i\right)}^{p-\alpha -1}=\frac{\partial }{\partial v_i}\left[{\left(g_{\alpha }\right)}_i\right],$ (25)

and taking (4) into account, Γ(p + 1 −α) can be rewritten in the form (p−α)Γ(p−α). So, (h _α )_i in equation (24) can be obtained.

For given µ > 0, the new iterate barrier parameter $\hat{\mu }$ is updated in the following form Wright (1996):

$\hat{\mu }=\frac{\mu }{\beta },$ (26)

where β > 1 is added to control your decay and to improve the convergence process.

To keep off the next û and $\hat{v}$ in the border region is required to shorten the step size α _uv by introducing a parameter σ, with 0 < σ < 1, which is often a value close to 1, as follows Biegler (2010); Vanderbei (2020):

${\alpha }_{uv}=\min \left\{\sigma \left(\underset{{\left(\Delta u\right)}_i<0}{\min }\frac{u_i}{\left|{\left(\Delta u\right)}_i\right|}\right),\sigma \left(\underset{{\left(\Delta v\right)}_i<0}{\min }\frac{v_i}{\left|{\left(\Delta v\right)}_i\right|}\right),\sigma \right\}.$ (27)

Given a point (x, λ, u, v), the new iterate $\left(\hat{x},\hat{\lambda },\hat{u},\hat{v}\right)$ is given by

$\hat{x}=x+{\alpha }_{uv}\Delta x,\hat{\lambda }=\lambda +{\alpha }_{uv}\Delta \lambda ,\hat{u}=u+{\alpha }_{uv}\Delta u,\hat{v}=v+{\alpha }_{uv}\Delta v.$

3.1 Search direction

The search direction (∆x, ∆λ, ∆u, ∆v) is unique and can be obtained by solving the Newton system (22):

$A^{⊺}\Delta \lambda =r_1,$ (28)

$A\Delta x+\Delta u-\Delta v=r_2,$ (29)

$U\Delta \lambda +D_u\Delta u+U{H}_{\alpha }\Delta v=r_3,$ (30)

$-V\Delta \lambda +V{H}_{\alpha }\Delta u+D_v\Delta v=r_4.$ (31)

Isolating ∆u in (30),

$\Delta u=D_u^{-1}\left(r_3-U\Delta \lambda -U{H}_{\alpha }\Delta v\right).$ (32)

Substituting (32) into (31) and isolating ∆v,

$\Delta v={D_1}^{-1}\left[r_4+V\Delta \lambda +V{H}_{\alpha }{D}_u^{-1}\left(U\Delta \lambda -r_3\right)\right],$ (33)

where

$D_1=D_v-UV{\left(H_{\alpha }\right)}^2{D}_u^{-1}.$ (34)

Substituting (32) into (29),

$A\Delta x+D_u^{-1}\left(r_3-U\Delta \lambda \right)-\left(U{H}_{\alpha }{D}_u^{-1}+I_m\right)\Delta v=r_2.$ (35)

Now, substituting (33) into (35),

$\Delta \lambda ={D_2}^{-1}\left[r_2-A\Delta x-D_u^{-1}{r}_3+{D_1}^{-1}\left(U{H}_{\alpha }{D}_u^{-1}+I_m\right)\left(r_4-V{H}_{\alpha }{D}_u^{-1}{r}_3\right)\right],$ (36)

where

$D_2=-U{D}_u^{-1}-V{D_1}^{-1}{\left(U{H}_{\alpha }{D}_u^{-1}+I_m\right)}^2.$ (37)

Then, ∆x can be obtained by replacing ∆λ given by (36) into (28). In this case, ∆x can be written as follows

$\Delta x={\left(A^{⊺}{D_2}^{-1}A\right)}^{-1}\overline{r},$ (38)

where

$\overline{r}=-r_1+A^{⊺}{D_2}^{-1}\left[r_2-D_u^{-1}{r}_3+{D_1}^{-1}\left(U{H}_{\alpha }{D}_u^{-1}+I_m\right)\left(r_4-V{H}_{\alpha }{D}_u^{-1}{r}_3\right)\right].$ (39)

3.2 Initialization

The initial point (x ⁰, λ ⁰, u ⁰, v ⁰) is chosen so that $u^0\in {\mathrm{ℝ}}_{++}^m,v^0\in {\mathrm{ℝ}}_{++}^m$ and $A{x}^0+u^0-v^0=b$, and is given by Coleman & Li (1992); Oliveira & Cantante (2004):

$\begin{array}{rcl}{x}^0& =& {\left(A^{⊺}A\right)}^{-1}{A}^{⊺}b,\\ r^0& =& b-A{x}^0,\\ {\lambda }^0& =& \frac{\delta r^0}{{||r^0||}_{\infty }}\text{,with}\delta =0.975\text{,}\\ u_i^0& =& \left\{\begin{array}{ll}{r}_i^0+{ϵ}_1,& \text{if}{r}_i^0\ge 0,\\ {ϵ}_1,& \text{if}{r}_i^0<0,\end{array}\right.\left(\text{for}i=1,2,\dots ,m\right),\\ v_i^0& =& \left\{\begin{array}{ll}{ϵ}_1,& \text{if}{r}_i^0\ge 0,\\ {ϵ}_1-r_i^0,& \text{if}{r}_i^0<0,\end{array}\right.\left(\text{for}i=1,2,\dots ,m\right),\end{array}$ (40)

where ε ₁ is a positive value close to zero.

The x ⁰ choice is the direct solution to the problem (11) for p = 2, and has proved to be a good choice in numerical experiments performed for polynomial regression models in the ℓ _p −norm for values of p reasonably closer to 2 Cantante et al. (2012); Oliveira & Cantante (2004); Cantane (2004); Grigoletto (2011).

3.3 Termination criteria

The fractional order log barrier interior point algorithm terminates if at least one of the two following convergence criteria is satisfied:

$N=\frac{||\left({\nabla }^{\alpha }\varphi +\nabla \phi \right)||}{\left(1+||x||+||u||+||v||+||\lambda ||\right)\left(2m\right)}\le {ϵ}_2,$ (41)

where m is the number of rows of the matrix A, and

$\left|\hat{N}-N\right|\le {ϵ}_3,$ (42)

where ε ₂ and ε ₃ are positive values close to zero, and $\hat{N}$ is the new iterate, given by

$\hat{N}=\frac{||\left({\nabla }^{\alpha }\hat{\varphi }+\nabla \hat{\phi }\right)||}{\left(1+||\hat{x}||+||\hat{u}||+||\hat{v}||+||\hat{\lambda }||\right)\left(2m\right)},$

where ${\nabla }^{\alpha }\hat{\varphi }+\nabla \hat{\phi }$ is obtained by considering the new iterates $x:=\hat{x},\lambda :=\hat{\lambda },u:=\hat{u},v:=\hat{v}$ and $\mu :=\hat{\mu }$ in equation (19).

3.4 Algorithm

The fractional order log barrier interior point algorithm works as follows. For given p, α, µ, β, σ, ε ₁, ε ₂, ε ₃, δ, and initial point (x, λ, u, v) from (40), the Newton system (22) is solved, and its solution (Δx, Δλ, Δu, Δv) is obtained. The step size α _uv is determined by (27) and the new iterate $\left(\hat{x},\hat{\lambda },\hat{u},\hat{v}\right)$ is obtained. The new barrier parameter $\hat{\mu }$ is updated by (26). This procedure repeat until the termination criterion is satisfied.

In particular, when α = 1, the fractional order log barrier interior point algorithm recover the classical log barrier interior point algorithm.

Algorithm 1
Fractional order log barrier interior point algorithm for polynomial regression models in the ℓ _p −norm.

4 NUMERICAL EXPERIMENTS

In order to illustrate the proposed fractional order log barrier interior point algorithm to solve polynomial regression models in the ℓ _p −norm, numerical experiments were performed to compare it with the classical log barrier interior point algorithm. The implementation of the fractional order log barrier interior point algorithm was performed under Windows 10 and Matlab (R2016a) running on a desktop with 2.20 GHz Intel Core i5-5200 central processing unit (CPU) and 4G random-access memory (RAM).

A data set containing daily interest rates observed over 40 years was used for the analysis of polynomial regressions. The data set contains 10958 observations: ${\left\{\left(a_i,b_i\right)\right\}}_{i=1}^{10958}$, where a _i represents the day of the week for a given date and b _i the interest rate (in percentage) for the specific day a _i . The ai values were normalized to the [0, 1] interval of the real line ℝ to avoid numerical stability problems of the algorithm Oliveira & Cantante (2004); Cantane (2004). Figure 1 shows the data set.

Figure 1
Daily changes in the interest rate.

For the numerical results below, a comparison of the fractional order log barrier interior point algorithm for different values of the order α considers: “It.”, the number of iterations; “Res. Err.”, the residual error, given by ∥r∥_p = ∥b−Ax∥_p ; and “Time (s)”, the CPU time in seconds.

The number of iterations and the value for the residual error for different values of α and p, obtained from the fractional order log barrier interior point algorithm for linear regression (n = 2) in the ℓ _p −norm, are shown in Tables 1-3. If the algorithm does not converge (divergence of Newton’s method for solving the nonlinear system of equations (21): ${||b-A{x}^{k+1}||}_p\ge {||b-A{x}^k||}_p)$, or if the algorithm fails (ill conditioned matrix), then the results will be shown in Tables 1-3 as “-” or “*”, respectively. When α = 1, the results of the classical log barrier interior point algorithm are recovered.

Table 1
Numerical results for linear regression in the ℓ _p −norm (p = 1.1, 1.2, 1.3) for different values of α.

	p = 1.1			p = 1.2			p = 1.3
α	It.	Res. Err.		It.	Res. Err.		It.	Res. Err.
0.1	−	−		−	−		−	−
0.2	5	10250		−	−		−	−
0.3	19	10118		5	5293		−	−
0.4	11	10297		20	5251		5	3044
0.5	∗	∗		14	5313		20	3024
0.6	3	10357		17	5324		11	3055
0.7	3	10358		3	5339		6	3064
0.8	5	10358		3	5339		3	3065
0.9	3	10359		5	5339		3	3065
1	2	10360		3	5340		5	3066

Table 2
Numerical results for linear regression in the ℓ _p −norm (p = 1.4, 1.5, 1.6) for different values of α.

	p = 1.4			p = 1.5			p = 1.6
α	It.	Res. Err.		It.	Res. Err.		It.	Res. Err.
0.1, 0.2, 0.3, 0.4	−	−		−	−		−	−
0.5	5	1905		−	−		−	−
0.6	18	1899		5	1274		−	−
0.7	11	1910		18	1272		5	900.8
0.8	123	1909		9	1278		19	900.3
0.9	3	1915		3	1280		11	902
1	3	1915		3	1280		3	903.5

Table 3
Numerical results for linear regression in the ℓ _p −norm (p = 1.7, 1.8, 1.9) for different values of α.

	p = 1.7			p = 1.8			p = 1.9
α	It.	Res. Err.		It.	Res. Err.		It.	Res. Err.
0.1, 0.2, . . . , 0.7	−	−		−	−		−	−
0.8	5	665.5		−	−		−	−
0.9	18	665.3		6	510.2		−	−
1	11	666.2		19	510.1		26	403.4

The algorithm failed for α = 0.5 and p = 1.1, as can be seen in Table 1, but for example, for α = 0.55 and p = 1.1, the algorithm converges.

According to the Tables 1-3, the fractional order log barrier interior point algorithm for linear regression models in the ℓ _p −norm (p = 1.1, 1.2, . . . , 1.7) yield smaller residual error when α = p − 0.8. It is a relationship between p and α, which can also be written as p−α = 0.8. In this cases, the algorithm takes more iterations until convergence.

For linear regression in the ℓ _1.2 −norm, a smaller residual error, given by 5242, was obtained with 16 iterations when α = 0.39 (p−α = 0.81). For α = 0.36, the residual error is 5256 and the number of iterations is 22. For α = 0.37, the residual error is 5252 and the number of iterations is 28. For α = 0.38, the residual error is 5247 and the number of iterations is 14.

For linear regression in the ℓ _1.8 −norm, a smaller residual error, given by 509.967, was obtained with 26 iterations when α = 0.96 (p − α = 0.84). For α = 0.94, the residual error is 510.011 and the number of iterations is 28. For α = 0.95, the residual error is 509.987 and the number of iterations is 16. For α = 0.97, the residual error is 510.015 and the number of iterations is 17.

For linear regression in the ℓ _1.9 −norm, the algorithm does not converge for α ≤ 0.92. For α = 0.93, the residual error is 403.48 and the number of iterations is 4. For α = 0.95, the residual error is 403.44 and the number of iterations is 5. For α = 0.97, the residual error is 403.42 and the number of iterations is 6. For α = 0.99, the residual error is 403.41 and the number of iterations is 6. When α = 0.99, the algorithm takes fewer iterations than when α = 1 (See Table 3), and the residual error is very close to the residual error obtained when α = 1.

The number of iterations, the residual error and the CPU time for (n − 1)th degree polynomial regressions (n = 2, 5, 10, 15), in the ℓ _1.3 −norm, for different values of α, obtained from the fractional order log barrier interior point algorithm, are given in the following tables.

The smallest residual errors for n = 2 and n = 5 were obtained with a greater number of iterations when α = 0.5 (see Table 4), while for n = 10 and n = 15 (see Table 5), the residual errors were smallest for α = 0.4.

Table 4
Numerical results for polynomial regressions (n = 2, 5) in the ℓ _1.3 −norm.

	n = 2				n = 5
α	It.	Time (s)	Res. Err.		It.	Time (s)	Res. Err.
0.1, 0.2, 0.3	−	−	−		−	−	−
0.4	5	0.0690	3044.3		16	0.1839	2433.6
0.5	20	0.1561	3024.8		20	0.2306	2423.5
0.6	11	0.1033	3055		17	0.2090	2431.8
0.7	6	0.0727	3064.1		5	0.0800	2441.4
0.8	3	0.0379	3065.7		5	0.0849	2443.3
0.9	3	0.0492	3065.8		3	0.0679	2443.7
1	5	0.0725	3065.9		3	0.0718	2443.8

Table 5
Numerical results for polynomial regressions (n = 10, 15) in the ℓ _1.3 −norm.

	n = 10				n = 15
α	It.	Time (s)	Res. Err.		It.	Time (s)	Res. Err.
0.1, 0.2, 0.3	−	−	−		−	−	−
0.4	2	0.0622	2213.4169		1	0.0555	1953.1475
0.5	1	0.0317	2220.5557		1	0.0583	1953.4556
0.6	1	0.0337	2221.8459		1	0.055	1953.9432
0.7	1	0.0320	2222.0088		1	0.0698	1954.1043
0.8	1	0.0312	2222.2097		1	0.0526	1953.8845
0.9	1	0.0318	2222.2880		1	0.0584	1953.8840
1	1	0.0314	2222.2941		1	0.0697	1953.8973

The performance of the fractional order log barrier interior point algorithm appears to be computationally consistent. For example, the values of the residual error at each iteration k for linear regression in the ℓ _1.3 −norm are shown in Figure 2.

Figure 2
Linear regression in the ℓ1.3−norm.

In the Figure 2, one can observe that the use of any one of the fractional order derivatives (α = 0, 4, 0.5, . . . , 0.9) produces smaller residual errors than with the use of the integer order derivative (α = 1). The use of the fractional derivatives does not interfere in the iterations running times, that compare similarly.

The approximate solution for (n − 1)th degree polynomial regression in the ℓ _p −norm, obtained at the kth iteration of the fractional order log barrier interior point algorithm and given by $x^k=\left(x_0^k,x_1^k,x_2^k,\dots ,x_{n-1}^k\right)$, provides the coefficients of the (n − 1)th degree polynomial that approximately fits the data set, given by $f\left(x\right)=x_0+x_{1}x+\dots +x_{n-1}{x}^{n-1}$. In particular, the approximate solutions x ^k , obtained at the kth iteration of the fractional order log barrier interior point algorithm, for (n − 1)th degree polynomial regressions (n = 2, 5, 10, 15) in the ℓ _1.3 −norm, with smaller residual errors are:

$\begin{array}{ll}\begin{array}{rl}{x}^{20}=& \left[\begin{array}{r}{x}_0^{20}\\ x_1^{20}\end{array}\right]=\left[\begin{array}{r}4.5E+00\\ 5.9E+00\end{array}\right]\end{array}& \left(n=2;\alpha =0.5\right),\\ & \end{array}$

i.e., at the 20th iteration of the algorithm with α = 0.5, the linear regression (n = 2) provides the polynomial coefficients from the above vector x ²⁰. So, the linear regression is given by f (x) = 4.5 + 5.9 x.

$$\begin{gathered} x^{20}=\left[\begin{array}{r}{x}_0^{20}\\ x_1^{20}\\ x_2^{20}\\ x_3^{20}\\ x_4^{20}\end{array}\right]=\left[\begin{array}{r}5.3E+00\\ -3.4E+01\\ 1.9E+02\\ -2.8E+02\\ 1.2E+02\end{array}\right]\left(n=5;\alpha =0.5\right), \\ {x}^2=\left[\begin{array}{r}{x}_0^2\\ x_1^2\\ x_2^2\\ x_3^2\\ x_4^2\\ x_5^2\\ x_6^2\\ x_7^2\\ x_8^2\\ x_9^2\end{array}\right]=\left[\begin{array}{r}1.8E+00\\ 1.7E+02\\ -3.7E+03\\ 3.5E+04\\ -1.6E+05\\ 4.6E+05\\ -7.6E+05\\ 7.2E+00\\ -3.7E+05\\ 7.9E+04\end{array}\right](n=10;\alpha =0.4), \\ {x}^1=\left[\begin{array}{r}{x}_0^1\\ x_1^1\\ x_2^1\\ x_3^1\\ x_4^1\\ x_5^1\\ x_6^1\\ x_7^1\\ x_8^1\\ x_9^1\\ x_{10}^1\\ x_{11}^1\\ x_{12}^1\\ x_{13}^1\\ x_{14}^1\\ \end{array}\right]=\left[\begin{array}{r}3.6E+00\\ 2.6E+00\\ -1.9E+02\\ 1.3E+04\\ -2.1E+05\\ 1.4E+06\\ -5.6E+06\\ 1.1E+07\\ -1.3E+07\\ 5.9E+06\\ -2.6E+06\\ 1.2E+07\\ -1.9E+07\\ 1.2E+07\\ -2.8E+06\end{array}\right]\left(n=15;\alpha =0.4\right). \end{gathered}$$

The polynomials from the solutions described are shown in Figure 3.

Figure 3
Polynomial regressions in the ℓ _1.3 −norm.

5 CONCLUSIONS

The fractional order log barrier interior point algorithm for polynomial regression models in the ℓ _p −norm (1 < p < 2), obtained by replacing some integer derivatives for the corresponding fractional ones on the first order optimality conditions, was investigated in this paper. This algorithm appears to be computationally consistent but depending on the fractional order α, the algorithm does not converge or fails. The numerical experiments showed that the use of the fractional derivatives can be beneficial for solving optimization problems. In future, further studies on this subject could be undertaken.

Acknowledgments

The author thanks the reviewers for their numerous helpful suggestions.

References

• ANDREWS GE, ASKEY R & ROY R. 1999. Special Functions. Cambridge: Cambridge University Press.
• BARD Y. 1974. Nonlinear Parameter Estimation. New York: Academic Press.
• BIEGLER LT. 2010. Nonlinear Programming. Philadelphia: SIAM.
• CANTANE DR. 2004. Métodos de Pontos Interiores Aplicados ao Problema de Regressão pela Norma Lp. São Carlos: Dissertação de Mestrado, Instituto de Ciências Matemáticas e de Computação, Universidade de São Paulo.
• CANTANTE DR, GRIGOLETTO EC & OLIVEIRA ARL. 2012. Método de pontos interiores barreira logarítmica preditor-corretor especializado para o problema de regressão pela norma Lp. Tendências em Matemática Aplicada e Computacional - TEMA, 13: 219-231.
• CHARNES A, COOPER WW & FERGUSON R. 1955. Optimal estimation of executive compensation by linear programming. Management Sci., 1: 138-151.
• CHATTERJEE S & HADI AS. 2012. Regression Analysis by Example. Hoboken: Wiley.
• CHEN Y, GAO Q, WEI Y & WANG Y. 2017a. Study on fractional order gradient methods. Appl. Math. Comput., 314: 310-321.
• CHEN Y, GAO Q, WEI Y & WANG Y. 2017b. Study on fractional order gradient methods. Appl. Math. Comput. , 314: 310-321.
• COLEMAN TF & LI Y. 1992. A globally and quadratically convergent affine scaling method for linear ℓ₁ problems. Math. Programming, 56: 189-222.
• DALIR M & BASHOUR M. 2010. Applications of fractional calculus. Appl. Math. Sci., 4: 1021-1032.
• DENNIS JR JE & SCHNABEL RB. 1983. Numerical Methods for Unconstrained Optimization. NJ: Prentice-Hall.
• ERDÉLYI A, MAGNUS W, OBERHETTINGER F & TRICOMI F. 1981. Higher Transcendental Functions. vol. I-III. Melbourne, Florida: Krieger Pub.
• FORSYTHE AB. 1972. Robust estimation of straight line regression coefficients by minimizing pth power deviations. Technometrics, 14: 159-166.
• GONDZIO J. 2012. Interior point methods 25 years later. Eur. J. Oper. Res., 218: 587-601.
• GRIGOLETTO EC. 2011. Implementação Eficiente dos Métodos de Pontos Interiores Especializados para o Problema de Regressão pela Norma Lp. Campinas: Dissertação de Mestrado, Instituto de Matemática, Estatística e Computação Científica, Universidade Estadual de Campinas.
• GRIGOLETTO EC & OLIVEIRA ARL. 2020. Fractional order gradient descent algorithm. In: Proceeding Series of the Brazilian Society of Computational and Applied Mathematics - XXXIX CNMAC. p. 010387. São Carlos: SBMAC.
• KILBAS AA, SRIVASTAVA HM & TRUJILLO JJ. 2006. Theory and Applications of Fractional Differential Equations. Amsterdam: Elsevier.
• LI Y. 1993. A globally convergent method for ℓ_p problems. SIAM J. Optimiz., 3: 609-629.
• LILIAN FB, OLIVEIRA ARL & GHIDINI CTLS. 2016. Use of continued iteration on the reduction of iterations of the interior point method. Pesqui. Oper., 36: 487-501.
• MOHAMMADZADEH A & KAYACAN E. 2020. A novel fractional-order type-2 fuzzy control method for online frequency regulation in ac microgrid. Eng. Appl. Artif. Intell., 90: 103483.
• OLDHAM KB & SPANIER J. 1974. The Fractional Calculus. New York: Academic Press .
• OLIVEIRA ARL & CANTANTE DR. 2004. Método de pontos interiores aplicados ao problema de regressão pela norma Lp. Tendências em Matemática Aplicada e Computacional - TEMA , 5: 281-291.
• OLIVEIRA ARL & LYRA C. 2004. Interior point methods for the polynomial ℓ_∞ fitting problems. Int. Trans. Oper. Res., 11: 309322.
• OLIVEIRA ARL, NASCIMENTO MA & LYRA C. 2000. Efficient implementation and benchmark of interior point methods for the polynomial ℓ₁ fitting problem. Comput. Statist. Data Anal., 35: 119-135.
• RICE JR & WHITE JS. 1964. Norms for Smoothing and Estimation. SIAM Review, 6: 243-256.
• SEVAUX M & MINEUR Y. 2007. A curve-fitting genetic algorithm for a styling application. Eur. J. Oper. Res. , 179: 895-905.
• SÜLI E & MAYERS DF. 2003. An Introduction to Numerical Analysis. New York: Cambridge University Press.
• VANDERBEI RJ. 2020. Linear Programming: Foundations and Extensions. Springer Nature.
• WANG J, WEN Y, GOU Y, YE Z & CHEN H. 2017a. Fractional-order gradient descent learning of BP neural networks with Caputo derivative. Neural Netw., 89: 19-30.
• WANG J, WEN Y, GOU Y, YE Z & CHEN H. 2017b. Fractional-order gradient descent learning of BP neural networks with Caputo derivative. Neural Netw. , 89: 19-30.
• WRIGHT SJ. 1996. Primal-Dual Interior-Point Methods. Philadelphia, USA: SIAM.

Publication Dates

• Publication in this collection
  25 Nov 2022
• Date of issue
  2022

History

• Received
  18 June 2021
• Accepted
  23 June 2022