A novel fixed-point based two-step inertial algorithm for convex minimization in deep learning data classification

1.   Introduction

Convex minimization problems are widespread in optimization. They are applied in various fields such as machine learning, signal processing, economics, and engineering. These problems are important because they are often assured to have a global minimum and thus give robust and reliable solutions to complex systems ^[1]. Convex minimization problems are fundamental in machine learning since they play an important role in training models, in particular for the support vector machines and logistic regression, where the goal is to minimize a loss function over additional regularizing terms ^[2]. These problems occur in signal processing tasks such as image reconstruction and denoising, where we seek to recover signals from noisy observations ^[3]. In finance as well, portfolio optimization is also based upon convex optimization, where one must find the optimal asset allocation that optimizes returns while minimizing risk ^[4].

The general form of a convex minimization problem can be expressed as:

where $f: \mathbb{R}^n \rightarrow \mathbb{R}$ is a smooth, convex function representing data fidelity, and $g: \mathbb{R}^n \rightarrow \mathbb{R}$ is a proper, lower semi-continuous convex function representing regularization or constraints. This formulation is widely used because it encompasses a broad range of practical problems where balancing accuracy and simplicity (or sparsity) of the solution is necessary ^[5]. Problem (1.1) has been studied by many researchers, and a number of algorithms have been proposed to solve this problem; see, for example, ^[6,7,8,9].

A powerful class of methods for solving convex minimization problems are fixed-point algorithms. These algorithms deal with finding $x^{*}$ such that $x^{*} = Tx^{*}$, where $T = prox_{\lambda g}(I – \lambda \nabla f)$ is a proximal operator, $f$ and $g$ are as in Eq (1.1), and $\lambda > 0$ is a step size parameter. In general, the proximal operator is defined for a function $g$ as:

This operator is used to deal efficiently with functions g that are nonsmooth, a common characteristic of regularization terms in convex optimization. It is well known that if $x$ is a fixed-point of the forward-backward operator

which involves the proximal operator, then $x$ is a minimizer of $f + g$. This relationship between fixed-points and minimizers is a key insight in convex minimization. It serves as the theoretical foundation for using fixed-point iterations to solve convex minimization problems ^[10,11].

Many researchers have proposed fixed-point algorithms that leverage the properties of the proximal operator ^{[12,13,14,15]}. Algorithms such as the proximal gradient method or forward-backward splitting (FBS) algorithm ^[7] update iterates using:

This method combines gradient descent for the smooth function $f$ with the proximal step for the nonsmooth function $g$, making it particularly useful for large-scale optimization problems. It is noted that when $L$ is a Lipschitz constant of $\nabla f$, then $T = \text{prox}_{\lambda g}(I - \lambda \nabla f)$ is nonexpansive if $\lambda \in (0, 2 / L)$. Solving variational inequalities and equilibrium problems is especially well suited to such an approach.

One of the most prominent forward-backward type algorithms with inertial step and design is the fast iterative shrinkage-thresholding algorithm (FISTA) introduced by Beck et al. ^[16]. The authors proved the algorithm's convergence and showed its practical use by applying it to image restoration problems, demonstrating its efficacy in real-world tasks. The FISTA algorithm is defined as follows:

where $n \geq 1$, $T : = prox_{\frac{1}{L}g}(I - \frac{1}{L}\nabla f)$, $x_1 = y_0 \in \mathbb{R}^n$, $t_1 = 1$, and $\theta_n$ represents the inertial step size introduced by Nesterov ^[17].

Bussaban et al. ^[12] introduced the parallel inertial S-iteration forward-backward algorithm (PISFBA), a new method designed to solve convex minimization problems efficiently. The algorithm is defined as:

where $n \geq 1$, $T_n = prox_{c_{n}g}(I-c_{n}\nabla f)$, $x_0 = x_1 \in H$, $0 < q < \alpha_n \leq 1$, $0 < s < \beta_n < r < 1$, and $\sum_{n = 1}^{\infty} \theta_{n} \|x_n - x_{n-1}\| < \infty$. The authors proved a weak convergence theorem for PISFBA under the assumption of Lipschit continuity of $\nabla f$. In addition, they showed their method's practical applicability by applying it to regression and data classification problems, showing that PISFBA can effectively solve real-world machine learning problems. After that, various fixed-point algorithms with a one-step inertial technique were introduced for solving the convex minimization problem (1.1) and employed to solving the data classification problem and the image restoration problem; see ^{[18,19,20,21]}.

To relax the continuity assumption of $\nabla f,$ D. Reem et al. ^[22] introduced a new variant of the proximal gradient method that does not impose the above-mentioned global Lipschitz continuity assumption.

The FISTA and PISFBA algorithms both employ a single inertial parameter, $\theta$. Such single-step algorithms are most effective when the problems are high-dimensional, and under more general conditions can be unstable and slow to converge. Two-step inertial algorithms incorporate two inertial parameters, resulting in improved convergence and stability. The basic idea is to update iterates using a combination of current and previous step information, formalized as:

where $\theta > 0$ and $\delta < 0$ are inertial parameters that incorporate momentum from previous iterates. This approach enhances motion modeling, improves stability, and increases redundancy, making it suitable for a broader range of applications. The incorporation of two inertial parameters also provides flexibility and adaptability in algorithm design, allowing fine-tuning for specific optimization problems. Recent research showcases the improved stability of two-step inertial algorithms over their one-step counterparts. Izuchukwu et al. (2023) ^[23] introduced a two-step inertial forward-reflected-anchored-backward splitting algorithm for monotone inclusion problems, demonstrating strong convergence with fewer computational steps compared to one-step methods. Iyiola and Shehu (2022) ^[24] proposed a two-step inertial proximal point algorithm for convex minimization, establishing a non-asymptotic $O(1/n)$ convergence rate. Recently, Thong et al. (2025) ^[25] introduced a double-inertial-step algorithm for split common fixed-point problems, demonstrating strong convergence with an application to signal processing.

Motivated by the literature discussed above, we propose a new two-step inertial algorithm that incorporates two inertial parameters and exhibits enhanced convergence. The proposed method is versatile and can be applied to various complex hierarchical optimization tasks. The remainder of the paper is organized as follows: In Section 2, we recall some basic definitions and results that are crucial for understanding the proposed method. In Section 3, we present our algorithm and the relevant convergence analysis, and we discuss the algorithm's application to solving convex minimization problems. Section 4 delves into the underlying principles of the machine learning models used for data classification and demonstrates the application of our algorithm by reformulating these models as convex minimization problems. Section 5 presents numerical experiments that demonstrate the performance of the proposed method. Finally, Section 6 concludes the paper with a summary of findings and potential future research directions.

2.   Preliminaries

In this section, we introduce the notation and give some fundamental definitions and lemmas that are used in the following sections. Let $H$ be a real Hilbert space with an inner product $\langle \cdot, \cdot \rangle$ and the induced norm $\|\cdot\|$. The set of positive integers, the set of real numbers, the set of non-negative real numbers, and the set of positive real numbers are denoted by $\mathbb{N}$, $\mathbb{R}$, $\mathbb{R}_{+}$, and $\mathbb{R}_{ > 0}$, respectively. For a mapping $T: H \to H$, the set of fixed-points of $T$ is denoted by $F(T)$ i.e., $F(T) = \{x \in H : Tx = x\}.$ Let $\psi$ be a family of mappings from $H$ into itself; $F(\psi)$ denotes the set of all common fixed-points of $\psi$, that is, $F(\psi) = \{x \in H: T(x) = x \, \, \, \text{for all}\, \, \, T \in \psi\}$. A mapping $T:H \to H$ is said to be nonexpansive if $\|Tx - Ty\| \leq \|x-y\|$ for all $x, y \in H.$

Definition 2.1. ^[26] Let $\{T_{n}\}$ and $\psi$ be countable families of nonexpansive mappings of $H$ into itself such that $\varnothing\neq F(\psi)\subset\cap_{n = 1}^{\infty}F(T_{n})$. We say that the sequence $\{T_{n}\}$ satisfies the NST-condition (I) with respect to $\psi$ if for any bounded sequence $\{u_{n}\}$ in $C$ and all $T\in \psi$, the following holds:

In the case that $\psi = \{T\},$ $\{T_{n}\}$ is said to satisfy the NST-condition (I) with respect to $T$.

Definition 2.2. Let $f: \mathbb{R}^n \rightarrow \mathbb{R}$ be a smooth convex function and let $g: \mathbb{R}^n \rightarrow \mathbb{R}$ be a proper, lower semi-continuous, convex function. Following Moreau ^[10], define, for $\lambda > 0$, the proximity operator with respect to $\lambda$ and $g$ by

see ^[11,27]. Define the forward-backward operator $T$ as

with $\lambda > 0$, where $\nabla f$ represents the gradient of the function $f$. If $\lambda \in (0, 2/L)$, where $L$ is a Lipschitz constant of $\nabla f$, then $T$ is nonexpansive. More related results on nonexpansive projections and resolvents of accretive operators in Banach spaces are discussed in ^[28].

The lemmas stated below play a significant role in proving our main result.

Lemma 2.1. ^[29] Suppose that $\{a_{n}\}$ and $\{b_{n}\}$ are two sequences of nonnegative numbers such that $a_{n+1}\leq a_{n}+b_{n}$ for all $n\geq 1$. If $\sum_{n = 1}^{\infty}b_{n}$ converges, then $\lim_{n\to\infty}a_{n}$ exists.

Lemma 2.2. ^[30] For a real Hilbert space $H$, the following results hold:

(i) For any vectors $x, y \in H$ and any scalar $\gamma\in [0, 1],$

(ii) For any vectors $x, y \in H,$

Lemma 2.3. ^[31] Let $H$ be a Hilbert space and $\{u_{n}\}$ be a sequence in $H$ such that for some nonempty subset $\Upsilon \subset H$, the following conditions are satisfied:

(i) For every $p\in\Upsilon$, $\lim_{n\to\infty}\|u_{n}-p\|$ exists.

(ii) Any weak-cluster point of the sequence $\{u_{n}\}$ belongs in $\Upsilon$.

Then, there exists $v\in \Upsilon$ such that $\{u_{n}\}$ weakly converges to $v$.

Lemma 2.4. ^[12] Let $H$ be a real Hilbert space. Consider a proper, lower semi-continuous convex function $g: H \rightarrow \mathbb{R} \cup \{\infty\}$ and a convex differentiable function $f: H \rightarrow \mathbb{R}$ whose gradient $\nabla f$ is Lipschitz continuous with constant $L > 0$. Let $T$ denote the forward-backward operator associated with $f$ and $g$. Suppose that $\{T_n\}$ is a sequence of forward-backward operators corresponding to step sizes $\{c_n\}$ such that $c_n \rightarrow c$ with $c_n, c \in (0, 2/L)$. Then, the sequence $\{T_n\}$ satisfies the NST-condition (I) with respect to $T$.

3.   Main results

In what follows, $T: H \rightarrow H$ is a nonexpansive mapping with $F(T) \neq \emptyset$, and $\{T_n : H \rightarrow H\}$ is family of nonexpansive mappings such that $F(T) \subset \Gamma$, where $\Gamma : = \cap_{n = 1}^{\infty}F(T_{n})$.

We now present Algorithm 1 and demonstrate its weak convergence.

Remark 3.1. Allowing $\delta_n$ to take nonpositive values provides a damping or "pull-back" effect that helps to stabilize the forward inertial term. Empirically, combining a positive $\theta_n$ and a negative $\delta_n$ strikes an effective balance between acceleration and stability in a two-step inertial algorithm; see ^[23,24,25] for more details.

Assume that the control sequences $\{\tau_{n}\}$, $\{\theta_{n}\}$, $\{\delta_{n}\}$, and $\{\beta_{n}\}$ in Algorithm 1 satisfy the following conditions:

(C1) $\sum_{n = 1}^{\infty}\tau_{n} < \infty$;

(C2) $0 < a_{1}\leq\beta_{n}\leq a_{2} < 1$ for some $a_{1}, a_{2}\in \mathbb{R}$.

The following lemmas provide a crucial estimate in the proof of Theorem 3.1.

Lemma 3.1. Let $\{x_{n}\}$ be a sequence generated by Algorithm 1. Then, for any $x^{*} \in \Gamma,$ we have

(i) $\|x_{n+1}-x^{*}\|\leq\|x_{n}-x^{*}\|+2\tau_n,$

(ii) $\lim_{n\to\infty}\|x_{n}-x^{*}\|$ exists.

Proof. Let $x^{*}\in\Gamma$. By Algorithm 1, we obtain

and

From the choice of $\theta_n$ and $\delta_n$, we note that

By (3.1)–(3.3), we obtain

Hence, by Lemma 2.1, it follows that $\lim_{n\to\infty}\|x_n - x^{*}\|$ exists, as required. □

Lemma 3.2. Let $\{x_{n}\}$ be a sequence generated by Algorithm 1. Then $\lim_{n\to\infty}\|T_{n}x_{n}-x_{n}\| = 0$.

Proof. By the definition of $y_{n}$ and Lemma 2.2, we obtain

This implies that, for $n\geq 1,$

By Lemma 3.1, we let

From (3.1), we obtain

It follows that

Using (3.2) and (3.4), we obtain

Since (C1) and the definitions of $\theta_{n}$ and $|\delta_{n}|$, we obtain

From (3.6) and (3.8), we obtain

It follows from (3.7) and (3.9) that

So, we get from (C3), (3.5), (3.6) and (3.10) that

as required. □

We now establish the main convergence property of Algorithm 1.

Theorem 3.1. Let $\{x_{n}\}$ be a sequence generated by Algorithm 1. Suppose that $\{T_{n}\}$ satisfies the NST-condition (I) with respect to $T$. Then, the sequence $\{x_{n}\}$ converges weakly to an element $x^{*}\in\Gamma,$ where $\Gamma: = \cap_{n = 1}^{\infty}F(T_{n})\neq\emptyset.$

Proof. Let $x^{*}\in\Gamma$. Using Lemma 3.1, we obtain $\lim_{n\to\infty}\|x_{n}-x^{*}\|$ exists. From Lemma 3.2 and $\{T_{n}\}$ satisfying the NST-condition (I) with respect to $T$, we obtain

Since $I-T$ is demiclosed at $0$, we obtain $\omega_{w}(x_{n})\subset F(T)$, where $\omega_{w}(x_{n})$ is the set of all weak cluster points of $\{x_{n}\}$. So, based on Lemma 2.3, we can conclude that $\{x_{n}\}$ converges weakly to $x^{*}\in F(T) \subset \Gamma$. □

Remark 3.2. From Lemma 3.1 and Theorem 3.1, we know that the sequence $\{x_n\}$ generated by Algorithm 1 converges weakly to $x^* \in \bigcap_{n = 1}^\infty F(T_n)$ and

where $\nu_n = 2\tau_n$.For each $n \in \mathbb{N}$, put $e_n = \|x_n -x^*\|.$ From the condition (C2), we can choose $\nu_n$ in such a way that $\lim_{n\to\infty} \frac{\nu_n}{e_n} = 0.$ One can show that $\limsup_{n\to\infty} \frac{e_{n+1}}{e_n} \leq 1.$ In this case, the sequence $\| x_{n} - x^* \|$ seems to be linear or sublinear convergence.

Next, we present Algorithm 2, which adapts Algorithm 1 by setting $T : = \operatorname{prox}_{c g}(I - c \nabla f)$ and $T_n : = \operatorname{prox}_{c_n g}(I - c_n \nabla f)$. This modification allows us to solve the convex minimization problem and establish the weak convergence of the generated sequence to a minimizer of $f + g$.

Theorem 3.2. Let $f, g: \mathbb{R}^n \rightarrow \mathbb{R} \cup \{\infty\}$ be functions such that $f$ is convex and differentiable, having Lipschitz continuous $\nabla f$ with constant $L > 0$, and $g$ is a proper, lower semicontinuous, convex function. Consider the sequence $\{x_n\}$ generated by Algorithm 2, with parameters $\{\tau_n\}$, $\{\theta_n\}$, $\{\delta_n\}$, and $\{\beta_n\}$ as defined in Algorithm 1. Then, the sequence $\{x_n\}$ converges to an element of $\operatorname{Argmin}(f + g)$.

Proof. We begin by noting that both $T$ and each $T_n$ are nonexpansive operators. For all $n$, the fixed-points of $T_n = \operatorname{prox}_{c_n g}(I - c_n \nabla f)$ are exactly the minimizers of $f + g$. Therefore,

By Lemma 2.4, this means that the sequence $\{T_n\}$ satisfies the NST-condition (I) with respect to $T$. Applying Theorem 3.1, we conclude that the sequence $\{x_n\}$ generated by Algorithm 2 converges to a point in $F(T) = \operatorname{Argmin}(f + g)$. □

4.   Applications

Artificial intelligence (AI) and its subfields, like machine learning (ML) and deep learning (DL), have been applied to many domains, such as image processing and data classification. The adoption of these technologies has made possible huge advances in solving complex problems that were previously intractable. AI, for example, has played a key role in the analysis of images for making medical diagnoses, helping diagnose diseases earlier. In image processing, DL programs have remarkably advanced object detection and increased the accuracy of recognition tasks. Numerous additional examples could be given from a variety of industries, where ML and DL models have improved predictive capabilities and helped solve difficult classification problems ^[32,33].

One of the major building blocks for the development of AI is a type of ML algorithm called a neural network, which is a model inspired by the structure of the human brain, consisting of interconnected layers of nodes (or 'neurons') that process data. In the very early days of neural networks, it was found that multiple-layered neural networks, called deep neural networks, could be used to model complex nonlinear relationships in data. However, traditional training algorithms, such as backpropagation, were computationally intensive and prone to overfitting and local minima ^[34].

To address these challenges, Huang et al. (2006) ^[35] introduced the extreme learning machine (ELM), simplifying the training of single hidden-layer, feed-forward neural networks (SLFNs) significantly. ELM randomly initializes the weights and biases that exist between the input and hidden layers and then, using a least squares solution, figures out the output weights. Training time is reduced, and generalization performance is improved using this approach. After their success with ELM, the researchers studied the deeper architectures in order to improve the learning capabilities further. To allow ELM to capture more complex data patterns, Qu et al. (2016) ^[36] proposed the two-hidden-layer ELM (TELM). The TELM algorithm exhibited improved performance over traditional ELM when working with nonlinear and high-dimensional data.

Xiao et al. (2017) ^[37] take this concept and expand it by developing the multihidden-layer extreme learning machine (MELM), which uses multiple hidden layers. For the first hidden (random) layer, MELM follows the characteristic random initialization of ELM, whereas the parameters of subsequent hidden layers are iteratively adjusted. This multilayer structure helps the network to better model important data relationships. Thus, it is an excellent network algorithm for regression and classification problems.

Let $\{ (\mathbf{X}, \mathbf{T}): \mathbf{X} \in \mathbb{R}^{N \times n}, \mathbf{T} \in \mathbb{R}^{N \times m} \}$ be a training set of $N$ distinct samples, where $\mathbf{X}$ is the training data matrix, $n$ is the number of input nodes (features), $\mathbf{T}$ is the target matrix, and $m$ is the number of output nodes (classes).

4.1.   Extreme Learning Machine (ELM)

In the context of ELM (see Figure 1 for the ELM structure), the training set consists of $N$ samples. The training data matrix $\mathbf{X}$, the target matrix $\mathbf{T}$, the input weight matrix $\mathbf{W}$, and the bias matrix $\mathbf{B}$ are defined as follows:

First, the input data matrix $\mathbf{X}$ passes through the model, where it is multiplied by the weight matrix $\mathbf{W}$ that connects the input layer to the hidden layer. After multiplying by $\mathbf{W}$, the bias matrix $\mathbf{B}$ is added, where $L$ is the number of hidden nodes. The core principle of ELM is that with a sufficiently large number of hidden neurons, the random initialization of the weight matrix $\mathbf{W}$ and bias vector $\mathbf{B}$ has minimal impact on the model's accuracy. As such, in ELM, $\mathbf{W}$ and $\mathbf{B}$ are typically initialized randomly without significant consequences on performance ^[35].

The output of this operation is passed through an activation function $g(\cdot)$ at the hidden layer. The purpose of the activation function is to introduce nonlinearity into the model, enabling it to learn complex patterns in the data. Common activation functions include the sigmoid, ReLU, and hyperbolic tangent functions. The choice of an activation function depends on the specific application and the nature of the data.

The hidden layer output matrix $\mathbf{H} \in \mathbb{R}^{N \times L}$ is computed as:

The ELM model uses the equation $\mathbf{H} \mathbf{\beta} = \mathbf{T}$, where the weight matrix between the hidden and output layers (output weights) $\mathbf{\beta} \in \mathbb{R}^{L \times m}$ is defined by

To determine the output weights $\mathbf{\beta}$, the ELM model minimizes the least-squares error between the predicted output $\mathbf{Y}$ and the target output $\mathbf{T}$. This is achieved using the Moore-Penrose pseudoinverse:

where $\mathbf{H}^\dagger$ is the Moore-Penrose pseudoinverse of $\mathbf{H}$. The pseudoinverse $\mathbf{H}^\dagger$ is computed as:

assuming $\mathbf{H}^T \mathbf{H}$ is invertible. The output of the network is then given by:

where $\mathbf{Y} \in \mathbb{R}^{N \times m}$ is the predicted output matrix.

4.2.   Two-Hidden-Layer Extreme Learning Machine (TELM)

The TELM algorithm tries to calculate and update the weights and biases between the first and second hidden layers, as well as between the second hidden layer and the output layer (see Figure 2 for the TELM structure). Initially, the computation is simplified by merging the two hidden layers into one equivalent hidden layer. Output weights are calculated for the combined hidden layer, and then the combined layer is separated, and the model is refined further by updating the weights and biases via additional steps.

The TELM algorithm starts by treating the two hidden layers as one hidden layer to simplify the initialization. The output of this combined hidden layer, denoted by $\mathbf{H} \in \mathbb{R}^{N \times L}$, is expressed as:

where $\mathbf{X} \in \mathbb{R}^{N \times n}$ is the input matrix, $\mathbf{W} \in \mathbb{R}^{n \times L}$ is the weight matrix of the first hidden layer, $\mathbf{B} \in \mathbb{R}^{N \times L}$ is the bias matrix of the first hidden layer (broadcasted to match the dimensions), $L$ is the number of hidden nodes, and $g(\cdot)$ is the activation function applied element-wise. The weight matrix $\mathbf{W}$ and bias matrix $\mathbf{B}$ are initialized randomly.

Once $\mathbf{H}$ is obtained, the output weight matrix $\mathbf{\beta}$ between the second hidden layer and the output layer is determined using:

where $\mathbf{H}^{\dagger}$ is the Moore-Penrose pseudoinverse of $\mathbf{H}$, and $\mathbf{T} \in \mathbb{R}^{N \times m}$ is the target matrix.

After the initial output weight matrix $\mathbf{\beta}$ is determined, the TELM algorithm separates the two hidden layers, which were previously merged, so that the network now contains two distinct hidden layers. To refine the model, the expected output, $\mathbf{H}_1$, of the second hidden layer is computed in a way that minimizes error, with the previously calculated output weight matrix beta acting as a constraint. This is done using $\mathbf{\beta}$ and the target matrix $\mathbf{T}$:

where $\mathbf{\beta}^{\dagger}$ is the Moore-Penrose pseudoinverse of the matrix $\mathbf{\beta}$.

Next, the algorithm updates the weights between the first and second hidden layers in a way that minimizes error, with $\mathbf{H}_1$ and $\mathbf{H}$ acting as constraints. Observe that, according to the typical algorithm, the expected output of the second hidden layer would be calculated using

where $\mathbf{W}_1 \in \mathbb{R}^{L \times L}$ is the weight matrix between the first and second hidden layers, and $\mathbf{B}_1 \in \mathbb{R}^{N \times L}$ is the bias matrix for the second hidden layer. Using the expected output $\mathbf{H}_1$ and the inverse function of the activation function g(·), the updated weight matrix $\mathbf{W}_{HE}$ is calculated as:

where $\mathbf{W}_{HE} = [\mathbf{B}_1 \quad \mathbf{W}_1]^{T} \in \mathbb{R}^{(L+1) \times L},$ $\mathbf{H}_E^{\dagger}$ is the Moore-Penrose pseudoinverse of $\mathbf{H}_E = [1 \quad \mathbf{H}] \in \mathbb{R}^{N \times (L+1)}$, and $g^{-1}(\cdot)$ is the inverse of the activation function $g(\cdot)$.

With the updated weight matrix $\mathbf{W}_{HE}$, the actual output of the second hidden layer $\mathbf{H}_2 \in \mathbb{R}^{N \times L}$ is updated as follows:

where $\mathbf{H}_2$ represents the final output of the second hidden layer after weight and bias adjustments.

Finally, the output weight matrix $\mathbf{\beta}_{new} \in \mathbb{R}^{L \times m}$ between the second hidden layer and the output layer is updated as follows:

where $\mathbf{H}_2^{\dagger}$ is the Moore-Penrose pseudoinverse of $\mathbf{H}_2$. The final output (predicted output matrix) $\mathbf{Y} \in \mathbb{R}^{N \times m}$ of the TELM network is then computed by:

Briefly summarized, the TELM first simplifies the network by merging the hidden layers, calculates the output weight matrix as in the one-hidden layer case, and then separates the layers, iteratively refining the weight and bias by a series of updates. This process guarantees that the network can replicate complex relationships in the data by using multiple hidden layers and fine-tuning their parameters.

4.3.   Multihidden-Layer Extreme Learning Machine (MELM)

The multihidden-layer extreme learning machine (MELM) architecture extends ELM and TELM by adding additional hidden layers, improving the learning capabilities of the model. The incorporation of the additional layers enables the network to leverage the data in finer detail, making the model especially effective for complicated tasks. An examination of a model with three hidden layers, as shown in Figure 3, serves to illustrate how MELM works. The strategy is to reduce the three-layer network by replacing the second and third hidden layers with a single combined equivalent hidden layer, thus reducing it to a TELM framework. In this way, the TELM framework is extended for use on more complex network structures.

A 3-layer MELM begins by combining the second and third hidden layers into one. The TELM methodology is then used to derive the output weight matrix, $\mathbf{\beta}_{new}$, between the combined second hidden layer and the output layer. After obtaining the initial output weight matrix, $\mathbf{\beta}_{new}$, the MELM algorithm then separates the merged hidden layers, resulting in three distinct hidden layers. To further refine the model, the expected output of the third hidden layer, $\mathbf{H}_3 \in \mathbb{R}^{N \times L}$, is computed using the updated weight matrix $\mathbf{\beta}_{new}$ as a constraint:

where $\mathbf{\beta}_{new}^{\dagger} \in \mathbb{R}^{L \times m}$ is the Moore-Penrose pseudoinverse of the weight matrix $\mathbf{\beta}_{new}$.

Observe that, according to the typical algorithm, the expected output of the third hidden layer, $\mathbf{H}_3$, would be calculated using

where $\mathbf{W}_2 \in \mathbb{R}^{L \times L}$ is the weight matrix between the second and third hidden layers, $\mathbf{H}_2$ is the output from the second hidden layer, and $\mathbf{B}_2 \in \mathbb{R}^{N \times L}$ is the bias term for the third hidden layer. The algorithm updates the weights between the second and third hidden layers using the inverse function of the activation function g(·) and the expected output $\mathbf{H}_3$. The weight matrix $\mathbf{W}_{HE1}$ is calculated as:

where $\mathbf{W}_{HE1} = [\mathbf{B}_2 \quad \mathbf{W}_2]^T \in \mathbb{R}^{(L+1) \times L}$, $\mathbf{H}_{E1}^{\dagger}$ is the Moore-Penrose pseudoinverse of $\mathbf{H}_{E1} = [1 \quad \mathbf{H}_2] \in \mathbb{R}^{N \times (L + 1)}$, and $g^{-1}(\cdot)$ is the inverse of the activation function $g(\cdot)$.

With the updated weight matrix $\mathbf{W}_{HE1}$, the actual output of the third hidden layer, $\mathbf{H}_4 \in \mathbb{R}^{N \times L}$, is calculated as follows:

where $\mathbf{H}_4$ represents the final output of the third hidden layer before reaching the output layer.

Finally, the output weight matrix $\mathbf{\beta}_{final} \in \mathbb{R}^{L \times m}$ between the second hidden layer and the output layer is updated as follows:

The final output (predicted output matrix) $\mathbf{Y} \in \mathbb{R}^{N \times m}$ of the three-hidden-layer ELM network is given by:

This process can be generalized to networks with more than three hidden layers in a natural inductive way. For a network with $n$ hidden layers, the network structure is simplified by merging the last two layers, applying the algorithm for $n-1$ layers to obtain an output weight matrix, then separating the merged layers, updating the weights between layers $n-1$ and $n$, and completing the calculation for the last layer, as shown above in the discussions of TELM and MELM. This approach guarantees that the network can handle still more complicated data patterns with little additional computational cost. Specifically, the computational complexity of this process is $\mathcal{O}(n)$, as each additional layer contributes a fixed computational cost, ensuring scalability for deeper networks.

4.4.   Reformulation as a convex minimization problem and application of TIFBSPA

When the Moore-Penrose pseudoinverse matrices in expressions (4.2)–(4.5) and (4.7) exist, the computation of such matrices is usually straightforward. Unfortunately, in cases in which the Moore-Penrose pseudoinverse does not exist, computing the actual inverse directly is impractical or numerically unstable, making the usual approach nontrivial to apply. To overcome this problem, we employ regularization techniques using the least absolute shrinkage and selection operator (lasso) ^[38], which adds a regularization term to perform stable inverse calculations. This reformulates the calculation as a convex minimization problem. The resulting computation is more robust, and the solution exists.

We can reformulate Eqs (4.2)–(4.5) and (4.7) as follows:

where $\lambda$ is the regularization parameter controlling the trade-off between fitting the data and penalizing the complexity of the patterns that the model is able to exploit (i.e., the L1 norm of the coefficients). This regularization ensures that even when the Moore-Penrose pseudoinverse does not exist or the problem is ill-posed, a solution can still be obtained by minimizing the above convex objective functions.

We note that Eqs (4.15)–(4.18) are each in the form of a sum of two convex functions.

We employ our algorithm (TIFBSPA) to solve these convex minimization problems.

5.   Numerical experiments

We present numerical experiments to illustrate the performance of TIFBSPA when applied to convex minimization problems arising from standard classification tasks. In these experiments, the main goal is to demonstrate the performance of our algorithm as compared to the FISTA and PISFBA methods from the literature. The experimental setup, datasets utilized, parameter settings for the algorithms, and the results are discussed in detail.

The specific convex minimization problems we deal with in our experiments are presented in Table 1. Each row of the table shows the functions $f(\cdot)$ and $g(\cdot)$ that arise in the MELM algorithm (see Section 4.4). In managing the trade-off between data fitting and model complexity, the value $\lambda = 10^{-5}$ is used. In solving these problems with TIFBSPA, we show how TIFBSPA can address both the accuracy and complexity of the problem we are solving.

We test the performance of our algorithm by applying it to several datasets that vary in their dimensions (i.e., number of features), number of samples, and class distributions. Applying the algorithm to data sets with diverse characteristics allows for a comprehensive study of the algorithm's capability to handle different kinds and amounts of data. Each dataset is divided into training and testing sets. Table 2 provides an overview.

The datasets were chosen for their relevance to common classification tasks. Each dataset is briefly described below.

● Iris ^[39]: A well-known dataset in the ML community, this dataset is used for basic classification algorithm testing. It has 4 features and 3 classes.

● Ionosphere ^[40]: This is a binary classification using 34 features on a radar signal classification dataset.

● Hypertension: The healthcare related classification task obtained with 9 features and binary outcomes, which is particularly relevant for real-world applications, is presented in this dataset. The data was collected by Sripat Medical Center, Faculty of Medicine, Chiang Mai University, so it was real-world relevant and authentic.

● Weather Type ^[41]: Test the performance of the algorithm on multi-class classification problems with a more complex dataset containing 10 features and 4 classes. Because this dataset contains outliers, it can be used to test a machine learning model's accuracy and resiliency and provides a good platform to test different types of models.

In the first experiment, we measured the accuracy of each model in solving data classification, defining accuracy as follows:

To assess the effectiveness of different machine learning models, we conduct a series of experiments comparing the performance of the MELM (3-hidden-layer ELM), TELM, and ELM. For each dataset, the models are trained and tested across a range of hidden nodes, varying from 50 to 300, to examine how the number of hidden nodes influences mean accuracy. Our proposed fixed-point algorithm (TIFBSPA) is applied to each model, running 500 iterations across 10 executions. The mean accuracy from these executions is calculated and presented to highlight the algorithm's performance consistency. The variation in accuracy of classification results from the randomness of the initial weights and biases in each execution. Hence, it is more accurate to average the above results to give a better view of the model performance.

The results of the first experiment, summarized in Figure 4, demonstrate that in most of the tested datasets, the MELM performs better than both the traditional ELM and the TELM models. The most transparently observed advantage is in the Weather Type dataset (Figure 4d). Here, the MELM's better performance reflects that model's increased ability to capture complicated data trends. High variability and the presence of outliers in the weather data pose problems for models with fewer hidden layers. The three hidden layers of MELM enable the flowability of the model among the data and the acceptance of a greater number of features. This allows modeling of more complicated relationships, increasing the algorithm's accuracy. However, for the Hypertension dataset (Figure 4c), the performance of MELM, TELM, and ELM is more similar. This is likely because the Hypertension dataset is well-cleaned and preprocessed, reducing the need for a deeper model like MELM to handle the data's complexity. As a result, the simpler models perform comparably well. Because it exhibited the best performance in most datasets and offered the possibility to deal with more complex data, MELM was chosen as the model for the second experiment of the study: a comparison of TIFBSPA with FISTA and PISFBA.

In this experiment, we implemented the MELM model using TIFBSPA, FISTA, and PISFBA in order to compare the algorithms' performance across four datasets. All algorithms are run for a maximum of 1000 iterations, and the results report the best performance achieved during this process. The datasets used, Iris (L = 50), Ionosphere (L = 80), Hypertension (L = 100), and Weather Type (L = 60), are configured with optimal hidden nodes based on Figure 4. In particular, the number of hidden nodes for each dataset was selected to balance accuracy, computational efficiency, and the risk of overfitting, as larger numbers of nodes provide diminishing returns in accuracy while increasing complexity.

For the numerical experiments, it is crucial to carefully select the algorithm parameters, as they directly influence convergence behavior and the quality of the solution. Table 3 summarizes the parameters used for TIFBSPA, FISTA, and PISFBA. Careful selection of the parameters ensures that each algorithm performs optimally under the same conditions, allowing for fair comparison. The parameters for TIFBSPA are chosen to ensure a balance between convergence speed and solution accuracy. The choices for FISTA and PISFBA are based on recommendations from the literature.

The measures being used are accuracy, precision, recall, and F1 score. Since accuracy was already defined in the first experiment, the focus here is on the additional metrics.

● Recall: The proportion of actual positive cases that are correctly predicted.

Recall is all about identifying all the positive cases and therefore is important in circumstances where false negatives are costly.

● Precision: The proportion of predicted positive cases that are actually positive.

Precision focuses on minimizing false positives, which is essential when false alarms are undesirable.

● F1 Score: The harmonic mean of Precision and Recall, balancing the trade-off between them.

F1-Score takes both Precision and Recall and is useful where there is an unequal distribution of data and both false positives and false negatives are costly.

TP, TN, FP, and FN stands for true positives, true negatives, false positives, and false negatives, respectively. The performance of these algorithms, measured using these metrics, is summarized in Tables 4 and 5.

The results show that TIFBSPA performs better than the other methods in terms of test accuracy, with fewer iterations required to achieve an optimal solution. For example, with the Iris dataset, TIFBSPA achieved an optimal solution in 69 iterations, as compared to 97 iterations for FISTA and 556 iterations for PISFBA. Although TIFBSPA does not outperform the other algorithms in every single metric across all cases, it offers a reasonable compromise between speed and accuracy. In some cases, measures like precision and F1 score are a little lower, but they are still reasonable. These results demonstrate that TIFBSPA converges faster and performs comparably or better than popular, existing methods.

6.   Conclusions

In this paper, we proposed a novel two-step inertial algorithm (TIFBSPA) to tackle convex minimization problems when they arise in hierarchical optimization tasks. We proved a weak convergence theorem that the fixed-points of the algorithm are indeed minimizers of the associated convex minimization problems. This algorithm was then applied to a series of data classification tasks using the ELM, TELM and MELM models. All models were reformulated in the form of a series of convex minimization problems, which enabled effective use of the TIFBSPA algorithm. Numerical experiments performed on multiple datasets showed that the MELM model plus our proposed algorithm achieved the best classification accuracy and required fewer iterations to converge compared to FISTA and PISFBA. These results highlight the superior efficiency and performance of TIFBSPA.

We demonstrated that the proposed TIFBSPA algorithm is both theoretically robust and practically effective, making it a valuable tool for addressing machine learning and optimization challenges. This method can be extended to other optimization problems, and future research may investigate how this method can be integrated with other advanced machine learning approaches.

Author contributions

Conceptualization, R. W.; Formal analysis, K. J. and R. W.; Investigation, K. J. and R. W.; Methodology, R. W.; Software, K. J.; Supervision, S. S.; Validation, S. S.; Visualization, R. W.; Writing-original draft, K. J.; Writing-review and editing, R. W. and S. S.

Use of Generative-AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The authors thank the referees for their helpful comments and valuable suggestions. K. Janngam was supported by the CMU Proactive Researcher Program, Chiang Mai University, under Contract No. 568/2567. S. Suantai would like to thank Chiang Mai University, Fundamental Fund 2025, Chiang Mai University. R. Wattanataweekul would like to thank Ubon Ratchathani University.

Conflict of interest

The authors declare no conflict of interest.