Optimal Piecewise Polynomial Approximation for Minimum Computing Cost by Using Constrained Least Squares

1. Introduction

In many applications of scientific and engineering fields, approximation for complex function or data set is important such as compression of ECG signal [1], various voice processing applications speech recognition, synthesis, and conversion [2,3], and correction of sensor data [4,5]. Numerous methods have been studied for approximation, with the least squares method being the preferred choice. Because the least squares method offers high accuracy by minimizing the residuals and is also simple to implement in a computer program. Formulating the optimal coefficients as least squares problems dates back to Gendre (1805) and Gauss (1809). The first application can be found by Gergonne in 1815 [6]. Therefore, in this paper, the least squares method is employed for approximation.

The approximation error at the sample points generally tends to decrease as high-order polynomials are used, while this causes the overfitting which gives large oscillation between sample points. To resolve the overfitting, it is recommended to use lower-order polynomials by splitting the sections. In the case of using piecewise polynomials, it takes less computation time than using a single high-order polynomial with the same approximation error. Since the maximum number of intervals and orders depend on the number of samples and complexity, several researches have been proposed to find the appropriate order and number of intervals. In [7,8,9,10], piecewise polynomial fitting was proposed. [7] proposed Least Squares Piece fitting using a cubic polynomial. [9] used a method of adjusting the boundary of the segments to increase the operation speed, and [10] used Piecewise Polynomial Modeling to lower the error rate. However, in [8,9,10], there is a limitation that the order of the polynomial must be determined by the user. In [4], the whole domain is evenly divided into several intervals, and each interval was approximated by a cubic polynomial using the least squares method with the constraint that it had continuity in all boundaries. [11,12,13] proposed an approximation method using Piecewise Linear Approximation (PLA). [14] proposed a method using several linear functions and then moving the endpoints of the interval appropriately to reduce the approximation error of the interval.

There are two considerations when approximating with piece-wised polynomials: the order of polynomials and the number of intervals. The higher order of polynomials needs more computational cost. Further, this may cause overfitting. Therefore, it is necessary to determine the appropriate order of the polynomials. When using approximated function, i.e. piece-wised function, it requires additional steps to determine the subdomain corresponding to input value. If there are many intervals, the approximation precision increases. But more time also needed to determine subdomain. So, the order and the number of polynomials should be balanced and optimized. In order to take this, optimization scheme is proposed in contrast with [4,5,6,7,8,9,10,11,12,13,14] which utilize predetermined polynomial order and number of polynomials.

Several researches are adopted optimization methods for curve fitting. Among the optimization methods, hp-adaptive pseudo-spectral method, proposed in [15], is most similar to the algorithm proposed in this paper. The method determines the locations of segment breaks and the order of the polynomial in each segment. To be more specific, when the error within a segment displays a consistent pattern, the number of collocation points is increased. On the other hand, if there are isolated points with significantly higher errors than the rest of the segments, then the segment is split at those points. It produced solutions with better accuracy than global pseudo-spectral collocation while utilizing less computational time and resources. The hp-adaptive pseudo-spectral method and the optimal piecewise polynomial (OPP) function-approximation algorithm, proposed in this paper, are similar in that they increase accuracy by increasing the number of orders and intervals. While the OPP function-approximation algorithm compares computational costs to obtain an approximation of the minimum cost while satisfying the error norm construct for approximation.

As seen in Figure 1, approximation is processed by finding the optimized order and number of intervals that satisfy the error norm constraint with minimum computational cost. To do this, the algorithm adopts two cost functions. One is for the approximation error that is used by the least squares. The calculated approximation error is compared with the error norm construct for approximation. The other is computational cost function. It means program run time to calculate value of function for input value, x, and is composed of two runtime costs: cost for polynomial function call and cost for binary search tree to determine the relevant interval. In order words, the computational cost is calculated given the order of polynomial and the number of intervals. The costs for all possible combinations are calculated and sorted in ascending order offline. This is used to search for the order and number of intervals with the minimum cost satisfying the error constraint. Therefore, the OPP function-approximation algorithm is useful in systems that must be efficient to compute, such as real-time systems and embedded systems. When approximating, if the approximation functions of each interval are placed independently, discontinuity occurs over the entire interval. To address this, a constraint has been introduced using the Lagrange multiplier method to smoothly connect the approximation functions of each interval.

This paper is organized as follows. In Section 2, previous knowledge, related to the algorithms proposed Section 3, is described. In Section 4, the algorithm is examined with representative test functions. Subsequently, the OPP function-approximation algorithm is discussed and further work, to improve the performance, is proposed in Section 5.

2. Preliminaries

2.1. Nomenclature

The algorithm and all formulas in this study use the notation in Table 1. The symbols necessary to understand the algorithm are as follows: $a_{k,i}$means the coefficient of the $i^{th}$ order term of the $k^{th}$ polynomial. ${\mathit{q}}_k$ represents the coefficients of the approximation polynomial in the $k^{th}$ interval as a vector, and $\mathit{q}$ is a $(m+1)n$ vector concatenated all ${\mathit{q}}_k$. $c_a$ is the approximation cost and $c_c$ is the computational cost at runtime. $m^{*}$ and $n^{*}$ mean the optimal polynomial order and the number of intervals, respectively. Where $*$ describes optimized value for $c_a$ within ${ϵ}_a$ and minimum $c_c$. Scalars, vectors and matrices are written in lowercase, lowercase boldface and uppercase boldface, respectively.

2.2. Formulation of Constrained Least Squares

The least squares method is a representative approach in regression analysis to approximate functions or discrete samples by minimizing the sum of the squares of the residuals made in the results of each individual equation. The residual, is the difference between the given sample and the approximate value, represented by weighted squares as shown in (1). The matrix form of the least squares algorithm is represented as follows:

$$e=\frac{1}{2}{\sum }_{k=1}^n\left({\left({\mathit{y}}_k-{\mathit{F}}_k{\mathit{q}}_k\right)}^T{\mathit{W}}_k\right({\mathit{y}}_k-{\mathit{F}}_k{\mathit{q}}_k)$$

(1)

where ${\mathit{F}}_k=\left[\begin{array}{ccccc}1& x_{k,1}& x_{k,1}^2& \cdots & x_{k,1}^m\\ 1& x_{k,2}& x_{k,2}^2& \cdots & x_{k,2}^m\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& x_{k,\eta \left(k\right)}& x_{k,\eta \left(k\right)}^2& \cdots & x_{k,\eta \left(k\right)}^m\end{array}\right]$, ${\mathbf{y}}_k={\left[\begin{array}{ccc}{y}_{k,1}& \cdots & y_{k,\eta \left(k\right)}\end{array}\right]}^T$, and${\mathit{q}}_k={\left[a_{k,0}{a}_{k,1}\cdots a_{k,m}\right]}^T$. The least squares algorithm is accomplished by finding $\mathit{q}$ which makes the partial differential of error, i.e. the gradient, to be 0. Consequently, this leads to optimal polynomials with the minimum average of the residual sum of squares at sample points. For the sake of simplicity, Equation (1) can be shortened into a single monomial expression using sparse diagonal matrix $\mathit{F}$ as follows:

$$e=\frac{1}{2}{\left(\mathit{y}-\mathit{F}\mathit{q}\right)}^T\mathit{W}\left(\mathit{y}-\mathit{F}\mathit{q}\right)$$

(2)

where $\mathit{F}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left({\mathit{F}}_k\right),\mathit{y}={\left[\begin{array}{ccc}{\mathit{y}}_1^T& \cdots & {\mathit{y}}_n^T\end{array}\right]}^T$, and $\mathit{q}={\left[\begin{array}{cc}{\mathit{q}}_1^T& \begin{array}{cc}\cdots & {\mathit{q}}_n^T\end{array}\end{array}\right]}^T$. The polynomial, $p_k$, should be smooth at the boundaries, ${\mathit{v}}_{k-1}$ and ${\mathit{v}}_k$, with adjacent polynomials, $p_{k-1}$ and $p_{k+1}$, respectively. Without loss of generality, in the proposed algorithm the 1st order differentiability was appended as constraint, $\mathit{G}\mathit{q}=0$ for the smoothness. Note that if necessary, it is possible to increase the order for differentiability at the intersection of the adjacent intervals. Consequently, Equation (2) was modified with Lagrange Multiplier, $\mathit{\lambda }$ as follows:

$$e=\frac{1}{2}{\left(\mathit{y}-\mathit{F}\mathit{q}\right)}^T\mathit{W}\left(\mathit{y}-\mathit{F}\mathit{q}\right)+{\mathit{\lambda }}^T\mathit{G}$$

(3)

where $\mathit{G}=\mathrm{diag}\left({\mathit{G}}_k\right)$ and ${\mathit{G}}_l=\left[\begin{array}{cccccccc}1& {\mathit{v}}_k& \cdots & {\mathit{v}}_k^m& -1& -{\mathit{v}}_k& \cdots & -{\mathit{v}}_k^m\\ 0& 1& \cdots & m{\mathit{v}}_k^{m-1}& 0& -1& \cdots & -m{\mathit{v}}_k^{m-1}\end{array}\right]$. By partial differentiation of (3) with respect to $\mathit{q}$ and $\mathit{\lambda }$, two equations are obtained, that is $\frac{\partial }{\partial q}e=0$ and $\mathit{G}\mathit{q}=0$. By solving the above equation, the optimal coefficients, which minimize $c_a$, are obtained.

$$\mathit{q}=\mathit{H}(-{\mathit{G}}^T{\left(\mathit{G}\mathit{H}{\mathit{G}}^T\right)}^{-1}\mathit{G}\mathit{H}{\mathit{F}}^T{\mathit{W}}^T\mathit{y}+{\mathit{F}}^T{\mathit{W}}^T\mathit{y})$$

(4)

where $\mathit{H}={\left({\mathit{F}}^T{\mathit{W}}^T\mathit{F}\right)}^{-1}$.

3. OPP Approximation

3.1. Overall Algorithm Flow

The proposed algorithm is aimed at use in systems where computing time is critical, such as real-time systems and/or embedded systems. At algorithm runtime, it is particularly important to reduce function call times. Since the approximate function is composed of several polynomials, it needs to know how much time is required to execute the function at runtime. To do this, the computational cost function was defined from the four arithmetic operations and binary search tree. This is explained in more detail in the next subsection. Offline, based on the computational cost function, an approximation polynomial with the smallest calculation time is obtained, then used to the system.

The overall flow for the OPP approximation is described in Figure 2. First, the maximum values of polynomial order, $m_{max}$, and the number of intervals, $n_{max}$, to explore was set. Afterward, computational cost for all possible combinations of the number of intervals and the order is calculated and tabulated in ascending order ($m\in \left[2,m_{max}\right],n\in [1,n_{max}]$). Subsequently, $\mathit{q}$ and $c_a$ are calculated using the constrained least squares for the $k^{th}$ pair ${\left(m,n\right)}_k$, the sample points $(x,y)$, and the polynomial intervals $\mathit{v}$. This is repeated with randomly selected $\mathit{v}$ for $N$ times and the case with the smallest $c_a$ was selected. If $c_a$ is greater than ${ϵ}_a$, the next pair ${\left(m,n\right)}_{k+1}$ is examined and, $\mathit{q}$ and $c_a$ are calculated again. The loop is repeated until $c_a$ is less than ${ϵ}_a$. When $c_a$ becomes smaller than ${ϵ}_a$, the piecewise polynomials with the order and the number of intervals are determined as optimized approximate. If $c_a$ of the pair ${\left(m,n\right)}_{max}$ is greater than ${ϵ}_a$ about the pair, it is considered that finding the piecewise polynomial approximation function is failed in the specified range. In this case, it is possible to continue searching by larger $m_{max}$ and/or $n_{max}$ or by relaxing ${ϵ}_a$.

3.2. Computational Cost

In this section, the computational cost is described. Since the approximated function is piece-wised with several polynomials, at runtime it takes two kinds of computation times: the time for a single polynomial computation and the time for selecting specific polynomial. Therefore, the computational cost function is defined as follows:

$$c_c\left(m,n\right)=c_p+c_b$$

(5)

Computing time for a single polynomial depends on the order of polynomial. This is defined as $c_p\left(m\right)$ as follows:

$$c_p\left(m\right)=n_a{r}_a+n_r{r}_{r}m+n_f{r}_{f}m$$

(6)

where $n_a$, $n_r$, and $n_f$ are the number of assignments, arithmetic, and ‘for’ instructions, respectively. In addition, $r_a$, $r_r$, and $r_f$ are cycle counts of assignment, arithmetic, and ‘for’ instructions. $m$, the order of polynomial, means repetition count.

Since the approximate function is implemented with several polynomials, the clock cycles are required for the binary search which is to find the suitable intervals for a certain input $x$. This is defined as $c_b$ and is determined as the average of the short path case and long path case in a binary search tree. Two examples are illustrated in Figure 3.

Since the average clock cycles to determine the corresponding polynomial depends on the probability for which the input is in a certain interval, $c_b$ is a probability distribution function for input value. For the sake of simplicity, the probability that an input is in a particular interval is the same for all intervals. Consequently, the costs associated with the number of intervals are as follows:

$$c_b=c_{bs}{p}_s+c_{bl}{p}_l$$

(7)

where $c_{bs}=(d_m-1)(n_a{r}_a+n_d{r}_d+n_r{r}_r+n_w{r}_w+n_i{r}_i)$, $p_s=n_s/n$, $c_{bl}=d_m(n_a{r}_a+n_d{r}_d+n_r{r}_r+n_w{r}_w+n_i{r}_i)$, and $p_l=n_l/n$. $c_{bs}$ and $c_{bl}$ represent binary search tree cost for the short path and the long, respectively. Similarly, $n_s$ and $n_l$ are the number of the short path and the long path, respectively. Note that $n_s+n_l=n$. $d_m$ means the depth of the tree and is an integer value satisfying $2^{d_m-1}\le n<2^{d_m}$. In addition, $n_d$, $n_h$, and $n_i$ are the number of divisions, ‘while’, and ‘if’ instructions, respectively. Finally, $r_d$, $r_h$, and $r_i$ are cycle counts of division, ‘while’ and ‘if’ instruction for ARM Cortex-M7 core, respectively. The calculated $c_c(m,n)$ is organized into a table by $m$ and $n$ (see Figure 4) and it is sorted in ascending order. In this paper, it is a priority to increase $m$ when the costs are the same.

4. Experimental Results

In order to verify the performance of the proposed OPP algorithm, it was implemented using MATLAB with several nonlinear functions. In this experiment, we set the number of samples to 100, $m_{max}=10$, $n_{max}=20$, and ${ϵ}_a=0.0001$. And the boundary values are determined as the value when the approximate error is the smallest after changing 100 times randomly. Figure 5 is the result of approximation for the logarithmic function. The logarithmic function was chosen as the test function because there are many sensors whose output results are logarithmic, such as temperature sensors and distance sensors. Figure 6 is approximation results for the sine function that is frequently used in many applications. Figure 7 and Figure 8 show approximate results for triangular and square functions, respectively. The triangular function is adopted to examine an approximation to the continuous nonlinear function with undifferentiable points. In Figure 7, approximate with a smooth curve at a sharp point. And the square function is used to curve fitting the function which has discontinuous points. Figure 8 shows that overfit occurs in the section where the value changes rapidly, resulting in oscillation. Figure 9 is the result of approximation of the sigmoid function. The sigmoid function, also known as the logistic function, is commonly used in various fields such as neural networks, machine learning, and statistics. Figure 10 is the result of approximation for the ReLU function which is one of the most commonly used activation functions in deep learning models. Figure 11 is the result of the user specified function for testing a nonlinear function that is more complex than the previous differential functions across all intervals. The function is represented as follows:

$$y=\frac{1}{1+e^{-x}}+\sin \left(x\right)-\ln \left(x\right)$$

(8)

Overall, the approximation for continuous and differentiable functions over all intervals is a satisfactory result. However, in functions with continuous but undifferentiable points and functions with discontinuous points, they have less accurate and higher $c_c(m,n)$ than continuous functions. This can be solved by increasing the sample point. Another way is to vibrate the boundary value, which will be studied later.

5. Conclusions

To approximate complex function or discrete sample points, Optimal Piecewise Polynomial (OPP) function-approximation algorithm was proposed. The maximum values of polynomial order and the number of intervals to explore are preset, and the computational cost is calculated to determine the order and number of intervals for the approximation function. Then the combination of the order and the number of intervals is sorted in ascending order based on computational cost offline. Subsequently, the coefficient of the polynomial and the approximation error at the sample points are determined using constrained least squares. If the error is greater than the given error bound, the next combination is examined until the error is less than the bound. Ultimately, the OPP function-approximation algorithm determines the fastest approximation function at runtime within permissible error. The performance of the proposed algorithm has been demonstrated through several nonlinear functions.

There are some further works to obtain a more accurate approximation function. In this paper, the optimal intervals were obtained by random sampling, i.e. Mote Carlo method. To improve the performance, gradient-based method will be also applied to determine the optimal boundaries which minimizes the approximation error. Moreover, the performance will be verified and analyzed by experiments with actual embedded systems.