# Pruning Multilayered ELM Using Cholesky Factorization and Givens Rotation Transformation

Authors: Jingyi Liu; Xinxin Liu; Chongmin Liu; Ba Tuan Le; Dong Xiao
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

Extreme learning machine is originally proposed for the learning of the single hidden layer feedforward neural network to overcome the challenges faced by the backpropagation (BP) learning algorithm and its variants. Recent studies show that ELM can be extended to the multilayered feedforward neural network in which the hidden node could be a subnetwork of nodes or a combination of other hidden nodes. Although the ELM algorithm with multiple hidden layers shows stronger nonlinear expression ability and stability in both theoretical and experimental results than the ELM algorithm with the single hidden layer, with the deepening of the network structure, the problem of parameter optimization is also highlighted, which usually requires more time for model selection and increases the computational complexity. This paper uses Cholesky factorization strategy and Givens rotation transformation to choose the hidden nodes of MELM and obtains the number of nodes more suitable for the network. First, the initial network has a large number of hidden nodes and then uses the idea of ridge regression to prune the nodes. Finally, a complete neural network can be obtained. Therefore, the ELM algorithm eliminates the need to manually set nodes and achieves complete automation. By using information from the previous generation’s connection weight matrix, it can be evitable to re-calculate the weight matrix in the network simplification process. As in the matrix factorization methods, the Cholesky factorization factor is calculated by Givens rotation transform to achieve the fast decreasing update of the current connection weight matrix, thus ensuring the numerical stability and high efficiency of the pruning process. Empirical studies on several commonly used classification benchmark problems and the real datasets collected from coal industry show that compared with the traditional ELM algorithm, the pruning multilayered ELM algorithm proposed in this paper can find the optimal number of hidden nodes automatically and has better generalization performance.


---

## Body

## 1. Introduction

Extreme learning machine (ELM) was first proposed by Huang and has attracted extensive attentions for its extremely fast learning speed, least human intervention, and easy implementation [1–3]. In the ELM algorithm, the input weights and hidden biases are randomly generated from any continuous probability distribution, and then, the output weights can be solved using the generalized Moore–Penrose inverse. Compared with the BP neural network, this algorithm has a good performance network in regression [4–6], classification [7–9], feature learning [10–12], and cluster tasks [13–15]. Different from conventional gradient-based neural network learning algorithms, which are sensitive to the combination of parameters and easy to trap in local optimum, ELM not only has a faster training speed but also has a smaller training error. The learning theories of ELM show that it can maintain the universal approximation and classification capability of SLFNs even if it works with randomly generated hidden nodes. However, because the essence of the extreme learning machine is the empirical risk minimization model and the randomness of the input weights and hidden biases, the extreme learning machine requires more hidden nodes than the BP neural network. Inspired by the structural risk minimization (SRM) principle of statistical learning theory and the weighted least square method, Deng et al. [16] proposed a regularized extreme learning machine (RELM). On the basis of regularization, Zhou et al. [17] used Cholesky factorization to solve the regularized extreme learning machine and recursive least squares method for online learning. Chou’s method implements online training by adding new training samples one by one and introduces the forgetting factor to reduce the influence of old training samples, so as to highlight the role of new training samples.Although RELM has better robustness, it does not solve how to obtain the most suitable network structure like ELM. When there are too few hidden nodes in the neural network, the phenomenon of under fitting will appear, and when there are too many hidden nodes, there will be over fitting. Therefore, suitable hidden nodes are particularly important for neural networks. Rong et al. [18] proposed a pruning extreme learning machine. This method first defines a large amount of hidden nodes for the network, then prunes the label according to the correlation degree of each node, and finally obtains a suitable network. Martinez-Martinez et al. [19] used regularization to prune the hidden nodes. The initial stage of this method is the same as the pruning extreme learning machine, which is to first define a large number of hidden nodes for the network. The difference is that the method proposed by Martinez uses ridge regression to prune the network structure, thereby pruning redundant nodes.When the input of the extreme learning machine contains a lot of noise, the fitting effect is often not ideal. Therefore, Miche et al. [20] proposed an optimal pruning extreme learning machine (OP-ELM) based on the correlation of input variables. This method uses the leave-one-out criterion to determine whether the node is left or not, which can effectively reduce the sensitivity to noise and obtain a more concise network model. Miche et al. [21] proposed a new and improved ELM (TROP-ELM) on the basis of the optimal pruning extreme learning machine. TROP-ELM adds L1 and L2 penalty terms to solve the OP-ELM instability in the process of solving the residual sum of squares. Indeed, the first step in the improvement process performs a L1 penalty on the output layer to rank the hidden nodes by implementing a LARS and then followed sequentially by a L2 penalty applied on the L1 penalized production from the LARS, thus avoiding numerical instabilities and efficient pruning of the network accordingly.In recent years, due to the successful application of deep learning in many practical problems, academic research on ELM has shifted from single-layer networks to multilayer networks. The multilayered extreme learning machine (MELM) [22] is a new framework that attempts to integrate the characteristics of deep neural networks with a multilayer structure into ELM, but it also brings huge hidden layers of MELM that need to be predefined manually. Therefore, the exploration of the neural network structure is still an important issue. Moreover, the Moore–Penrose generalized inverse matrix is involved in the process of using the least squares method to solve the output matrix. If the matrix is not full rank, it will often cause numerical calculation errors.Aiming at the above open problem of MELM, this paper proposes a pruning multilayer ELM algorithm framework (P-MELM). The algorithm adds the Tikhonov regularization factor in the process of calculating the output weights to achieve efficient pruning of nodes and then uses fast matrix decomposition strategy to automatically learn the optimal number of nodes. Since a lot of operations are involved in the pruning process of hidden nodes, this paper proposes a more stable numerical operation method, which uses Cholesky factorization of the covariance matrix and Givens rotation transformation to calculate pseudo inverse. This method can avoid large-scale matrix inversion, so compared with the direct calculation of the least square solution, this method needs less computing time and memory.The chapters of this article are arranged as follows. Section1 reviews the improvements made by previous generations to ELM. Section 2 mainly introduces MELM. Section 3 proposes and introduces our algorithm P-MELM. Section 4 establishes a model and compares it with other algorithms. Section 5 summarizes the full text.


## 2. Multilayered ELM (MELM)

Suppose there areN arbitrarily different samples X,T=xi,tii=1,2,…,N, where X=x1,x2,…,xNT∈RN×n is the input sample, T=t1,t2,…,tNT∈RN×m is the corresponding sign sample, and Gx=1/1+e−x is the selected activation function. First, using the principle of random sampling, a value is randomly selected within the range of [0-1] to become the input weight W1,L=w1,1,w1,2,…,w1,LT∈RL×n of the first hidden layer and the threshold B1,L=b11,b12,…,b1LT∈RL×N of the first hidden layer. The standard mathematical formula for MELM is as follows:(1)∑j=1Lβ1,LjGw1,jxi+b1j=ti,i=1,2,…,N.By using the ELM algorithm to write formula (1) in the matrix form, we can get(2)H1,Lβ1,L=T,where(3)H1,L=GW1,LX+B1,L,=Gw1,1x1+b11⋯Gw1,Lx1+b1L⋮⋮⋮Gw1,1xN+b11⋯Gw1,LxN+b1L,=h1,1⋯h1,L,(4)β1,L=β1,L1,β1,L2,…,β1,LLT.As described in RMELM [23], the output matrix of the first hidden layer in MELM is defined as H1,L∈RN×L, whose elements h1,ji=Gw1,jxi+b1ji=1,2,…,Nandj=1,2,…,L may be explained as the output of the jth node in the first hidden layer in reference to xi, where h1,j=gw1,jx1+b1j⋯gw1,jxN+b1jTj=1,2,…,L.Unlike the conventional learning algorithm for SLFNs, the MELM approach is thus to generate randomly the weightsw1,j and the bias b1j in the first hidden layer and determine analytically the connection weights β1,L∈RL×m between the first hidden layer and the output layer using Moore–Penrose generalized inverse, which is identical to finding the minimum norm least square solution of the linear system, and therefore, representation (2) can be written as(5)β1,L=H1,L+T,where H1,L+ is the Moore–Penrose (MP) pseudoinverse of the matrix H1,L. If H1,LTH1,L is nonsingular, then H1,L+=H1,LTH1,L−1H1,LT, or alternatively H1,L+=H1,LTH1,LH1,LT−1 if H1,LH1,LT is nonsingular.According to the solving process of the MELM algorithm, we can get(6)H2,L=gW2,LH1,L+B2,L,where W2,L represents the weights’ matrix that links the first hidden layer to the second hidden layer. The matrices B2,L and H2,L, respectively, denote the bias and the expected output of the second hidden layer. Moreover, the output matrix H2,L is made on the basis of the multilayer theory and can be rewritten as(7)H2,L=Tβ1,L+,where β1,L+ is the MP pseudoinverse of the matrix β1,L, obtained using the approach discussed before, namely, β1,L+=β1,LTβ1,L−1β1,L if β1,LTβ1,L is nonsingular; otherwise, β1,L+=β1,LTβ1,Lβ1,LT−1 if β1,Lβ1,LT is nonsingular.The augmented matrix is then defined asW2HE,L=B2,LW2,L and computed as follows:(8)W2HE,L=g−1H2∗,LH2E,L+,where H2E,L+ represents the pseudoinverse of the matrix H2E,L=1H1,LT, 1=1,1,…,1T∈RN, and g−1x represents the inverse function of gx. The calculation of H2E,L+ is the same as the previous description of H1,L+. Subsequently, we can calculate the actual output matrix of the second hidden layer as(9)H2,L=gW2,LH1,L+B2,L=gW2HE,LH2E,L.Therefore,β2,L can be updated as follows:(10)β2,L=H2,L+T.The calculating method ofH2,L proceeds in the form discussed before. If MELM has only two hidden layers, the final output of MELM can be expressed as(11)tx=H2,Lβ2,L.When the number of hidden layers of MELM is greater than 2, as indicated in (5)–(10), an iterative solution scheme is adopted to implement the calculation process. Then, all the related parameters of these hidden layers can be obtained after the (M-2)-iteration operation. To make the final actual hidden layer output close to the expected hidden layer output, in the training stage, the MELM designs a novel parameter-setting step for the newly added hidden layer, as described in (8), which is the key step to guarantee the stability and feasibility of the MELM model.


## 3. Pruning Multilayered ELM (P-MELM)

This newly modified approach is in the form of a L2 regularization penalty applied within the MELM. First, the regularization method can constrain the norm of the output weights to avoid the possible numerical instabilities and prevent the model from overfitting. Second, the regularization term makes the coefficient matrix of linear equations positive definite and thus uses a fast matrix factorization strategy to recursively update the output weights in the calculation of the pruning criterion, which provides more accurate and reliable computation results. The pruning multilayered ELM algorithm based on the Cholesky factorization and Givens rotation transformation are deduced in this section, which avoids the complexities of matrix inversion.


### 3.1. The Solutions of P-MELM by Cholesky Factorization and Givens Rotation Transformation

Previous research has shown that due to the randomly generated parameters of the hidden nodes, MELM typically requires more hidden nodes than conventional gradient descent methods to provide desired performance levels. This essence of MELM clearly brings faster learning speed, least human intervention, and easy implementation. However, it suffers from performance problems connected with too many hidden nodes chosen at random, that is, the lack of sparsity and adjustment ability in the hidden layers. Furthermore, another common risk is that the network will produce some redundant nodes that add little or nothing to the existing efforts to respond to the model training. As a result of all these trends, the initial output matrix may not be a full column rank or even ill-conditioned; then, the exact analytical solution of the output weights can hardly be found, which directly affects the generalization capability and stability of the network. Aiming at this problem, it is proposed here to modify the pruning strategy for the selection of the optimal number of nodes by taking heuristic techniques to tentatively prune out the redundant nodes of the model that will decrease the computational complexity of providing an accurate ranking of the hidden layer nodes and raise the efficiency of the network learning and the adaptability to real problem resolving. In the following, the most well-known algorithms used to perform regularization are reviewed, using aL2 penalty in the MELM, to regularize the network.It is assumed that the initial value of hidden nodes is set asL=Lmax, the output matrix of the first hidden layer is H1,L, and the connection weights’ matrix between the first hidden layer and the output layer is β1,L. Through studying the training dataset, the optimal number of hidden nodes and hidden layers are determined, which use the P-MELM prediction model for a framework and employ the Cholesky factorization method and Givens rotation transformation to recursively calculate the factorization factor UL. As the number of hidden nodes decreases, the new connection weights’ matrix β1,L+1 can be obtained rapidly according to the parameter information of the existing network. More specifically, if we reduce the total number of hidden nodes from L to L − 1, the first hidden layer output matrix of P-MELM will change from H1,L∈RN×L to H1,L−1∈RN×L−1, which can be intuitively represented as the following:(12)H1,L=h1,1⋮h1,2⋯h1,L=h1,1⋮H1,L−1,where h1,j=Gw1,jx1+b1j⋯Gw1,jxN+b1jTj=1,2,…,L. According to (12), we can get(13)AL=CIL+H1,LTH1,L,=C100IL−1+h1,1TH1,L−1Th1,1H1,L−1,=C+h1,1Th1,1h1,1TH1,L−1H1,L−1Th1,1CIL−1+H1,L−1TH1,L−1,=PL−1QL−1QL−1TAL−1,where IL∈RL×L and IL−1∈RL−1×L−1 denote the identity matrices of order L and (L-1), respectively, QL−1=h1,1Th1,2⋯h1,1Th1,L, and PL−1=C+h1,1Th1,1. Due to the tight coupling that exists between AL−1 and AL, as shown in (13) above, the Cholesky factorization UL−1 of matrix AL−1 can be analytically and easily determined based on the Cholesky factorization theorem [24]. The recursive calculation of matrix factorizations is implemented as follows. First, the matrices AL and ULT are partitioned into four submatrices, and their corresponding block matrices have the same size. For convenience, the matrix ULT can thus be rewritten as(14)ULT=vVT0TWT,where VT is an (L − 1)-dimensional row vector, WT is an (L − 1) × (L − 1) square matrix, and v is a scalar. Substituting formula (14) to (13), AL−1 can be solved as(15)AL−1=VVT+WWT=UL−1UL−1T.Next, the rank one update algorithm [25] is used to calculate the Cholesky factorization factor, so UL−1T can be obtained from the equation of the form:(16)JL−1JL−2⋯Jk⋯J1VWT=0UL−1TL×L−1.Here,Jk is the Givens rotation transformation matrix, which can be expressed in the matrix form as(17)Jk=I+ck−1δ0δ0T+δkδkT+skδ0δkT−δkδ0T,where ck=cosθk=WkkT/WkkT2+VkT2, sk=sinθk=VkT/WkkT2+VkT2, and k=1,2,…,L−1.It is worth noting that, in (16), the matrix V and W change with the rank one update algorithm proceeding, so the calculation process of (17) does not need to explicitly record the matrix Jk. Meanwhile, as shown in (12), BL can be represented as(18)BL=H1,LTT=h1,1⋮H1,L−1TT,=h1,1TTH1,L−1TT=h1,1TTBL−1.If the number of hidden nodes decreases, the connection weights’ matrixβ1,L−1 between the first hidden layer and the output layer can be directly and analytically determined by solving the following equations:(19)AL−1β1,L−1=BL−1.Then, substitute formula (15) into (19) and multiply UL−1−1 of both sides, and the problem seeking for β1,L−1 can thus be formulated as the following:(20)UL−1Tβ1,L−1=FL−1,where FL−1=UL−1−1BL−1. For convenience, FL−1=UL−1−1BL−1 is intuitively transformed as(21)UL−1FL−1=BL−1,and it can be shown that FL−1 is associated with UL−1 and BL−1, which are indeed obtained by solving the linear equation system. Accordingly, the optimal output weight matrix β1,L−1 can be calculated simultaneously using the corresponding entries of matrices UL−1 and FL−1.


### 3.2. Decremental Learning Procedure of P-MELM

The process builds an optimum network structure, known as decremental learning of P-MELM, through the calculation of the corresponding Cholesky factorization factorUL. It can be seen from formulas (14)–(20) that the P-MELM algorithm adopts a fast calculation format and can achieve a decreasing update of the connection weight matrix. The solution method of β1,L−1 based on Cholesky factorization and Givens rotation transformation makes full use of the information stored in the calculation of β1,L, where UL−1 can be obtained through Givens rotation transformation on the basis of UL, and BL−1 can be obtained directly on the basis of BL. Therefore, if the number of hidden nodes decreases successively, β1,L−1 can be rapidly calculated by simple matrix arithmetic operation on the basis of computing β1,L. In this case, if β1,L−1 is calculated by using the method shown in (5), it must be recalculated in the way of finding the inverse of the higher-order matrix, and it cannot be directly solved on the basis of computing β1,L. Therefore, the P-MELM algorithm can guarantee the learning accuracy and further improve the training speed. We have summarized the steps of the P-MELM algorithm as follows:Step 1: first, we limit the number of nodes. The maximum number of nodes isLmax, the minimum number of nodes is Lmin, and the iteration stopping criterion is ξL. Then, we set the initial number of hidden nodes as L=Lmax and calculate AL and BL.Step 2: calculate the Cholesky factorizationUL of AL according to (13), and calculate FL using the matrices UL and BL according to (21).Step 3: calculateβ1,L with the matrices UL and FL according to (20). Then, build a P-MELM model according to β1,L. The network structure of the model has M hidden layers, and each hidden layer has L hidden layer nodes.Step 4: the empirical risk and structural risk of the P-MELM algorithm are solved, and the sum of these two parts is carried out:(22)R1,L=H1,Lβ1,L−T2+Cβ1,L2,=H1,Lβ1,L−TTH1,Lβ1,L−T+Cβ1,LTβ1,L.Step 5: letL=L−1, calculate UL and BL on the basis of UL+1 and BL+1 according to the principle shown in (14)–(18), and then, go to Step 3. Start from L=Lmax−1 to determine whether (23) is satisfied:(23)R1,L+1−R1,LR1,max>ξL,whereR1,max is the resulting maximum value of R1,L,⋯,R1,Lmax. If the conditions of (23) are met, then modeling iteration is completed, the optimal number value of hidden nodes is confirmed as L, and the corresponding prediction model of P-MELM is built. Otherwise, if L>Lmin, continue to decrease L until the condition L=Lmin is satisfied.In the P-MELM algorithm flow, first, the neural network has a larger number of hidden nodesLmax, and then, the number of nodes is gradually reduced until R1,L changes significantly. If the number of hidden nodes continues to decrease at this time, the generalization ability and fitting effect of P- P-MELM ELM will decrease along with the progress, and a large number of effective nodes will also be lost in the network. Therefore, when R1,L changes significantly, P-P-MELM ELM has the optimal number of hidden nodes.


### 3.3. Proposed P-MELM

P-MELM shows an obvious superiority in faster computations and numerical stability compared with the conventional MELM by introducing Cholesky factorization and Givens rotation transformation. This technique uses the basic arithmetic operation on matrix elements to construct weights matrix and avoids the complicated matrix inverse. In addition, when the hidden node is reduced one by one, the dynamically changing network structure needs to retrain the model by the existing training data, and this will lead to significantly longer model training times. Faced with this situation, a better alternative is to fully utilize information gained by the process of P-MELM training and implement the decremental updating of network parameters recursively.First, train P-MELM with the dataset {X, T} and then, use trained P-MELM to provide the prediction of the verifying data Z. The chosen neural network architecture consisted of M hidden layers, one input layer, and output layer; meanwhile, each of the hidden layers contains L hidden nodes. Then, the output of the standard MELM classifier with respect to the input x can be modeled as(24)tM,Lx=∑j=1LβM,LjGwM,j,bMj,x,x∈Rn.If the number of training samples isN >> L, we set the initial value of hidden nodes as Lmax, the minimum number of hidden nodes as Lmin, the initial number of hidden layers as 1, and the maximum number of hidden layers as Mmax. ξL is expressed as the best accuracy of neural network learning.Phase 1: initialization of the P-MELM network.(1)
LetL=Lmax and M=1. Randomly select the input weight w1,j and bias b1j of the first hidden layer, j=1,⋯,L.(2)
Solve the output matrixH1,L of the first hidden layer according to (3), H1,L=Gw1,1,…,w1,L,b11,…,b1L,x1,…,xN.(3)
Calculate the matrixAL and BL according to (13) and (18).(4)
Cholesky factorization ofAL according to (14) to get UL.(5)
According to (21), calculate FL by using the matrices UL and BL.(6)
According to (20), calculate β1,L by using the matrices UL and FL. Then, build a P-MELM model according to β1,L. The network structure of the model has M hidden layers, and each hidden layer has L hidden layer nodes.(7)
According to (22), the empirical risk and structural risk of the P-MELM algorithm are obtained, and the two parts are summed.Phase 2: Update P-MELM neural network sturctural parameters (the number of hidden nodes and hidden layers). WhileL≥Lmin.(1)
LetL=L−1. Calculate UL and BL on the basis of UL+1 and BL+1 according to the principle shown in (14)–(18), and then, go to Step 1 (5). If (23) is satisfied, then L is the optimal number of hidden nodes at this time. If (23) is not satisfied, then L=Lmax−1 continues to solve.(2)
LetM=M+1. Calculate the expected output matrix HM,L of the Mth layer according to (7), HM,L=TβM−1,L+.(3)
Calculate the learning parameterWMHE,L in the Mth hidden layer according to (8), WMHE,L=g−1HM,LHME,L+.(4)
Update the output matrix of theMth hidden layer, HM,L=gWMHE,LHME,L.(5)
According to (19)–(21), the connection weights between the Mth hidden layer and the output layer βM,L are updated.(6)
IfM=Mmax, then stop training; otherwise, return to 2 (1), and continue to increase the number of hidden layer M until the condition is met.(7)
Update the output weight matrixβM,L and establish a P-MELM model. And, calculate the final output tM,Lx of the P-MELM algorithm according to (24).Therefore, the predicted output of test instance Z can be represented as tM,LZ.If we setMmax=1 in the P-MELM approach, P-MELM becomes the pruning extreme learning machine (P-ELM). The flow chart of the P-MELM algorithm is shown in Figure 1. It can be seen from Figure 1 that when the number of hidden layer nodes in the P-MELM algorithm gradually decreases, P-MELM does not need to retrain the network and recalculate the weights, but updates the weights on the basis of the original network. MELM needs to recalculate the output weight every time once the number of hidden nodes is reduced. Therefore, when the number of hidden layer nodes decreases gradually, P-MELM is simpler and takes less time than MELM.Figure 1
Flow chart of the P-MELM algorithm.


## 3.1. The Solutions of P-MELM by Cholesky Factorization and Givens Rotation Transformation

Previous research has shown that due to the randomly generated parameters of the hidden nodes, MELM typically requires more hidden nodes than conventional gradient descent methods to provide desired performance levels. This essence of MELM clearly brings faster learning speed, least human intervention, and easy implementation. However, it suffers from performance problems connected with too many hidden nodes chosen at random, that is, the lack of sparsity and adjustment ability in the hidden layers. Furthermore, another common risk is that the network will produce some redundant nodes that add little or nothing to the existing efforts to respond to the model training. As a result of all these trends, the initial output matrix may not be a full column rank or even ill-conditioned; then, the exact analytical solution of the output weights can hardly be found, which directly affects the generalization capability and stability of the network. Aiming at this problem, it is proposed here to modify the pruning strategy for the selection of the optimal number of nodes by taking heuristic techniques to tentatively prune out the redundant nodes of the model that will decrease the computational complexity of providing an accurate ranking of the hidden layer nodes and raise the efficiency of the network learning and the adaptability to real problem resolving. In the following, the most well-known algorithms used to perform regularization are reviewed, using aL2 penalty in the MELM, to regularize the network.It is assumed that the initial value of hidden nodes is set asL=Lmax, the output matrix of the first hidden layer is H1,L, and the connection weights’ matrix between the first hidden layer and the output layer is β1,L. Through studying the training dataset, the optimal number of hidden nodes and hidden layers are determined, which use the P-MELM prediction model for a framework and employ the Cholesky factorization method and Givens rotation transformation to recursively calculate the factorization factor UL. As the number of hidden nodes decreases, the new connection weights’ matrix β1,L+1 can be obtained rapidly according to the parameter information of the existing network. More specifically, if we reduce the total number of hidden nodes from L to L − 1, the first hidden layer output matrix of P-MELM will change from H1,L∈RN×L to H1,L−1∈RN×L−1, which can be intuitively represented as the following:(12)H1,L=h1,1⋮h1,2⋯h1,L=h1,1⋮H1,L−1,where h1,j=Gw1,jx1+b1j⋯Gw1,jxN+b1jTj=1,2,…,L. According to (12), we can get(13)AL=CIL+H1,LTH1,L,=C100IL−1+h1,1TH1,L−1Th1,1H1,L−1,=C+h1,1Th1,1h1,1TH1,L−1H1,L−1Th1,1CIL−1+H1,L−1TH1,L−1,=PL−1QL−1QL−1TAL−1,where IL∈RL×L and IL−1∈RL−1×L−1 denote the identity matrices of order L and (L-1), respectively, QL−1=h1,1Th1,2⋯h1,1Th1,L, and PL−1=C+h1,1Th1,1. Due to the tight coupling that exists between AL−1 and AL, as shown in (13) above, the Cholesky factorization UL−1 of matrix AL−1 can be analytically and easily determined based on the Cholesky factorization theorem [24]. The recursive calculation of matrix factorizations is implemented as follows. First, the matrices AL and ULT are partitioned into four submatrices, and their corresponding block matrices have the same size. For convenience, the matrix ULT can thus be rewritten as(14)ULT=vVT0TWT,where VT is an (L − 1)-dimensional row vector, WT is an (L − 1) × (L − 1) square matrix, and v is a scalar. Substituting formula (14) to (13), AL−1 can be solved as(15)AL−1=VVT+WWT=UL−1UL−1T.Next, the rank one update algorithm [25] is used to calculate the Cholesky factorization factor, so UL−1T can be obtained from the equation of the form:(16)JL−1JL−2⋯Jk⋯J1VWT=0UL−1TL×L−1.Here,Jk is the Givens rotation transformation matrix, which can be expressed in the matrix form as(17)Jk=I+ck−1δ0δ0T+δkδkT+skδ0δkT−δkδ0T,where ck=cosθk=WkkT/WkkT2+VkT2, sk=sinθk=VkT/WkkT2+VkT2, and k=1,2,…,L−1.It is worth noting that, in (16), the matrix V and W change with the rank one update algorithm proceeding, so the calculation process of (17) does not need to explicitly record the matrix Jk. Meanwhile, as shown in (12), BL can be represented as(18)BL=H1,LTT=h1,1⋮H1,L−1TT,=h1,1TTH1,L−1TT=h1,1TTBL−1.If the number of hidden nodes decreases, the connection weights’ matrixβ1,L−1 between the first hidden layer and the output layer can be directly and analytically determined by solving the following equations:(19)AL−1β1,L−1=BL−1.Then, substitute formula (15) into (19) and multiply UL−1−1 of both sides, and the problem seeking for β1,L−1 can thus be formulated as the following:(20)UL−1Tβ1,L−1=FL−1,where FL−1=UL−1−1BL−1. For convenience, FL−1=UL−1−1BL−1 is intuitively transformed as(21)UL−1FL−1=BL−1,and it can be shown that FL−1 is associated with UL−1 and BL−1, which are indeed obtained by solving the linear equation system. Accordingly, the optimal output weight matrix β1,L−1 can be calculated simultaneously using the corresponding entries of matrices UL−1 and FL−1.


## 3.2. Decremental Learning Procedure of P-MELM

The process builds an optimum network structure, known as decremental learning of P-MELM, through the calculation of the corresponding Cholesky factorization factorUL. It can be seen from formulas (14)–(20) that the P-MELM algorithm adopts a fast calculation format and can achieve a decreasing update of the connection weight matrix. The solution method of β1,L−1 based on Cholesky factorization and Givens rotation transformation makes full use of the information stored in the calculation of β1,L, where UL−1 can be obtained through Givens rotation transformation on the basis of UL, and BL−1 can be obtained directly on the basis of BL. Therefore, if the number of hidden nodes decreases successively, β1,L−1 can be rapidly calculated by simple matrix arithmetic operation on the basis of computing β1,L. In this case, if β1,L−1 is calculated by using the method shown in (5), it must be recalculated in the way of finding the inverse of the higher-order matrix, and it cannot be directly solved on the basis of computing β1,L. Therefore, the P-MELM algorithm can guarantee the learning accuracy and further improve the training speed. We have summarized the steps of the P-MELM algorithm as follows:Step 1: first, we limit the number of nodes. The maximum number of nodes isLmax, the minimum number of nodes is Lmin, and the iteration stopping criterion is ξL. Then, we set the initial number of hidden nodes as L=Lmax and calculate AL and BL.Step 2: calculate the Cholesky factorizationUL of AL according to (13), and calculate FL using the matrices UL and BL according to (21).Step 3: calculateβ1,L with the matrices UL and FL according to (20). Then, build a P-MELM model according to β1,L. The network structure of the model has M hidden layers, and each hidden layer has L hidden layer nodes.Step 4: the empirical risk and structural risk of the P-MELM algorithm are solved, and the sum of these two parts is carried out:(22)R1,L=H1,Lβ1,L−T2+Cβ1,L2,=H1,Lβ1,L−TTH1,Lβ1,L−T+Cβ1,LTβ1,L.Step 5: letL=L−1, calculate UL and BL on the basis of UL+1 and BL+1 according to the principle shown in (14)–(18), and then, go to Step 3. Start from L=Lmax−1 to determine whether (23) is satisfied:(23)R1,L+1−R1,LR1,max>ξL,whereR1,max is the resulting maximum value of R1,L,⋯,R1,Lmax. If the conditions of (23) are met, then modeling iteration is completed, the optimal number value of hidden nodes is confirmed as L, and the corresponding prediction model of P-MELM is built. Otherwise, if L>Lmin, continue to decrease L until the condition L=Lmin is satisfied.In the P-MELM algorithm flow, first, the neural network has a larger number of hidden nodesLmax, and then, the number of nodes is gradually reduced until R1,L changes significantly. If the number of hidden nodes continues to decrease at this time, the generalization ability and fitting effect of P- P-MELM ELM will decrease along with the progress, and a large number of effective nodes will also be lost in the network. Therefore, when R1,L changes significantly, P-P-MELM ELM has the optimal number of hidden nodes.


## 3.3. Proposed P-MELM

P-MELM shows an obvious superiority in faster computations and numerical stability compared with the conventional MELM by introducing Cholesky factorization and Givens rotation transformation. This technique uses the basic arithmetic operation on matrix elements to construct weights matrix and avoids the complicated matrix inverse. In addition, when the hidden node is reduced one by one, the dynamically changing network structure needs to retrain the model by the existing training data, and this will lead to significantly longer model training times. Faced with this situation, a better alternative is to fully utilize information gained by the process of P-MELM training and implement the decremental updating of network parameters recursively.First, train P-MELM with the dataset {X, T} and then, use trained P-MELM to provide the prediction of the verifying data Z. The chosen neural network architecture consisted of M hidden layers, one input layer, and output layer; meanwhile, each of the hidden layers contains L hidden nodes. Then, the output of the standard MELM classifier with respect to the input x can be modeled as(24)tM,Lx=∑j=1LβM,LjGwM,j,bMj,x,x∈Rn.If the number of training samples isN >> L, we set the initial value of hidden nodes as Lmax, the minimum number of hidden nodes as Lmin, the initial number of hidden layers as 1, and the maximum number of hidden layers as Mmax. ξL is expressed as the best accuracy of neural network learning.Phase 1: initialization of the P-MELM network.(1)
LetL=Lmax and M=1. Randomly select the input weight w1,j and bias b1j of the first hidden layer, j=1,⋯,L.(2)
Solve the output matrixH1,L of the first hidden layer according to (3), H1,L=Gw1,1,…,w1,L,b11,…,b1L,x1,…,xN.(3)
Calculate the matrixAL and BL according to (13) and (18).(4)
Cholesky factorization ofAL according to (14) to get UL.(5)
According to (21), calculate FL by using the matrices UL and BL.(6)
According to (20), calculate β1,L by using the matrices UL and FL. Then, build a P-MELM model according to β1,L. The network structure of the model has M hidden layers, and each hidden layer has L hidden layer nodes.(7)
According to (22), the empirical risk and structural risk of the P-MELM algorithm are obtained, and the two parts are summed.Phase 2: Update P-MELM neural network sturctural parameters (the number of hidden nodes and hidden layers). WhileL≥Lmin.(1)
LetL=L−1. Calculate UL and BL on the basis of UL+1 and BL+1 according to the principle shown in (14)–(18), and then, go to Step 1 (5). If (23) is satisfied, then L is the optimal number of hidden nodes at this time. If (23) is not satisfied, then L=Lmax−1 continues to solve.(2)
LetM=M+1. Calculate the expected output matrix HM,L of the Mth layer according to (7), HM,L=TβM−1,L+.(3)
Calculate the learning parameterWMHE,L in the Mth hidden layer according to (8), WMHE,L=g−1HM,LHME,L+.(4)
Update the output matrix of theMth hidden layer, HM,L=gWMHE,LHME,L.(5)
According to (19)–(21), the connection weights between the Mth hidden layer and the output layer βM,L are updated.(6)
IfM=Mmax, then stop training; otherwise, return to 2 (1), and continue to increase the number of hidden layer M until the condition is met.(7)
Update the output weight matrixβM,L and establish a P-MELM model. And, calculate the final output tM,Lx of the P-MELM algorithm according to (24).Therefore, the predicted output of test instance Z can be represented as tM,LZ.If we setMmax=1 in the P-MELM approach, P-MELM becomes the pruning extreme learning machine (P-ELM). The flow chart of the P-MELM algorithm is shown in Figure 1. It can be seen from Figure 1 that when the number of hidden layer nodes in the P-MELM algorithm gradually decreases, P-MELM does not need to retrain the network and recalculate the weights, but updates the weights on the basis of the original network. MELM needs to recalculate the output weight every time once the number of hidden nodes is reduced. Therefore, when the number of hidden layer nodes decreases gradually, P-MELM is simpler and takes less time than MELM.Figure 1
Flow chart of the P-MELM algorithm.


## 4. Algorithm Verification

In order to verify the stability of the P-MELM algorithm and P-ELM algorithm proposed in this paper, we have carried out experiments in MATLAB 2016a environment and carried out 100 experiments on each algorithm. Our experiment is divided into two parts. The first part is a simple benchmark classification [26], and the second part is to classify coal by using the spectral information of coal [27]. In addition, in order to comprehensively evaluate the performance of the P-ELM algorithm and P-MELM algorithm, the initial value of the number of nodes in a single hidden layer is set as 100 in the experiment, and the step size is set as 1 to delete redundant nodes. The number of layers of the multilayer extreme learning machine is set from 2 to 10.


### 4.1. Data Acquisition and Processing

We download two commonly used datasets of Diabetic Retinopathy Debrecen and image segmentation from the UCI Machine Learning Library to verify the reliability of the algorithm proposed in simple benchmark classification. In addition, when faced with extremely time-consuming and complex classification problems, researchers now tend to use machine learning methods to solve complex problems. Therefore, this paper also uses the spectral information of coal to classify coal to verify the reliability of P-ELM and P-MELM. The collected coal mine information is shown in Table1.Table 1
Type and quantity of coal collected.
Mineral speciesAttributeTraining dataTesting dataCoalAnthracite4427Bituminous coal5030Lignite3622Coal gangue6640Total196119We use the SVC-HR-1024 spectrometer produced by American Spectra Vista Company as the experimental instrument in the spectrum experiment. The spectral bandwidth of this instrument is 1024, and the spectral range is 350–2500 nm. In the experiment, firstly, the sample is cleaned. The spectrometer probe is 380 mm away from the surface of the ore sample, and as far as possible, it is set perpendicular to the sample surface. Secondly, the experimenters should reduce the walking process and wear dark clothes during the experiment to reduce the interference of environmental factors on the spectrum experiment. Then, the spectrometer was used to conduct 6 spectroscopic experiments on each sample of coal. Finally, summarize all the experimental data, and calculate the average of the six spectral experimental data of the same coal mine as the final spectral data. Figure2 shows the process of data collection and processing, during which whiteboard measurement and calibration are carried out every 10 minutes. Because the original spectral data included 1024 bands, some of which are overlapped bands, it is necessary to remove the overlapped bands from the data, and finally, 973 useful bands are obtained.Figure 2
Spectrum experiment and spectrum curve of coal.


### 4.2. Evaluation of Testing Accuracy

We use P-ELM and P-MELM to model the classification problem. Table2 shows the parameters of P-ELM and P-MELM. Figures 3 and 4, respectively, show the accuracy of data classification by P-ELM and P-MELM under different hidden nodes. From Figures 3 and 4, we can see that P-ELM and P-MELM can obtain good learning accuracy and have a relatively accurate classification effect. If the redundant hidden nodes in the neural network are pruned, the testing accuracy of P-ELM and P-MELM enhance or remain unchanged, and this phenomenon is especially evident in Diabetic Retinopathy Debrecen, image segmentation, and coal spectral datasets. More importantly, the classification accuracy surfaces of P-MELM keep relatively smooth, but the curves of P-ELM have some oscillations in Diabetic Retinopathy Debrecen and image segmentation cases when the redundant hidden nodes are pruned. By comparing the algorithm process of P-ELM and P-MELM, it can be seen that P-MELM changes from a single hidden layer to multiple hidden layers on the basis of P-ELM, and the calculation method of the parameters in the P-MELM algorithm is different. So, P-MELM can eliminate the fluctuations caused by P-ELM to some extent. As a consequence, the P-MELM algorithm is more accurate in classification than the P-ELM algorithm, and the fitting effect is better.Table 2
Parameters of P-ELM and P-MELM under different classification problems.
DatasetsAlgorithmInitial nodesEnd nodesHidden layersDiabetic Retinopathy DebrecenP-ELM100501P-MELM100802–10Image segmentationP-ELM100801P-MELM80552–10Coal spectral dataP-ELM100751P-MELM80552–10Figure 3
The average testing accuracy for the algorithm P-ELM using (a) Diabetic Retinopathy Debrecen, (b) image segmentation, and (c) coal spectral data, while the hidden node is decreased one by one.
(a)(b)(c)Figure 4
The average testing accuracy for the algorithm P-MELM with different hidden layers using (a) Diabetic Retinopathy Debrecen, (b) image segmentation, and (c) coal spectral data, while the hidden node is decreased one by one. Here, M and L are the number of hidden layers and the number of redundant nodes removed, respectively.
(a)(b)(c)


### 4.3. Evaluation of Average Training Time

In Section 4.3, we compare the training time of MELM and P-MELM under the same parameters. The time required for MELM and P-MELM algorithms is shown in Figure5. From Figure 5, we can see that whether it is coal mine classification, diabetic retinopathy classification, or image segmentation classification, and the training time required for P-MELM is much less than the training time required for MELM. This is because when the number of hidden nodes changes, MELM needs to retrain the network, and P-MELM does not need to recalculate the network parameters, but only updates the weights in the model. Therefore, under the same network structure, P-MELM has a faster convergence speed than MELM. Furthermore, the training time curves for several cases shown in Figure 4 also indicate that the spent learning time is still monotonically decreasing with the learning steps, which is consistent with MELM’s conclusions. Compared with P-MELM, MELM needs more time to retrain the network and recalculate the output weights. However, with the increase of hidden layer nodes, the time difference of training time will gradually decrease.Figure 5
Comparison of training time between the algorithms MELM and P-MELM using (a) Diabetic Retinopathy Debrecen, (b) image segmentation, and (c) coal spectral data, while the hidden node is decreased one by one.
(a)(b)(c)


### 4.4. Comparison of Testing Accuracy by Different Methods

The classification accuracy of the testing data is selected as the performance evaluation criterion for the algorithms investigated, namely, ELM, MELM, P-ELM, P-MELM, ELM-AE, and SELM. The best learning accuracy that six algorithms can achieve can be observed from Table3. On the basis of these experimental results, we can draw the following conclusions. If the initial network structure is fixed, P-ELM and P-MELM always get higher testing accuracy than ELM and MELM, respectively. This means that the P-ELM and P-MELM techniques can get more stable networks in most cases. What is more, a higher classification accuracy indicates a better classification algorithm performance. As observed from Table 3, it is obvious that the proposed P-MELM technique achieves the highest classification accuracy among all classification datasets, relative to the original ELM, MELM, and P-ELM algorithms. Finally, these results further prove that the P-ELM and P-MELM methods have a better classification effect than the ELM method and the training speed can meet the requirements.Table 3
The best testing accuracy that the algorithms ELM, MELM, P-ELM, P-MELM, ELM-AE, and SELM can achieve.
DatasetsBest testing accuracy (%)ELMMELMP-ELMP-MELMELM-AESELMDiabetic Retinopathy Debrecen63.6665.9767.5968.7568.0666.20Image segmentation75.0875.8981.3183.0482.493.0Coal spectral data89.9189.9191.6092.2491.6067.23


## 4.1. Data Acquisition and Processing

We download two commonly used datasets of Diabetic Retinopathy Debrecen and image segmentation from the UCI Machine Learning Library to verify the reliability of the algorithm proposed in simple benchmark classification. In addition, when faced with extremely time-consuming and complex classification problems, researchers now tend to use machine learning methods to solve complex problems. Therefore, this paper also uses the spectral information of coal to classify coal to verify the reliability of P-ELM and P-MELM. The collected coal mine information is shown in Table1.Table 1
Type and quantity of coal collected.
Mineral speciesAttributeTraining dataTesting dataCoalAnthracite4427Bituminous coal5030Lignite3622Coal gangue6640Total196119We use the SVC-HR-1024 spectrometer produced by American Spectra Vista Company as the experimental instrument in the spectrum experiment. The spectral bandwidth of this instrument is 1024, and the spectral range is 350–2500 nm. In the experiment, firstly, the sample is cleaned. The spectrometer probe is 380 mm away from the surface of the ore sample, and as far as possible, it is set perpendicular to the sample surface. Secondly, the experimenters should reduce the walking process and wear dark clothes during the experiment to reduce the interference of environmental factors on the spectrum experiment. Then, the spectrometer was used to conduct 6 spectroscopic experiments on each sample of coal. Finally, summarize all the experimental data, and calculate the average of the six spectral experimental data of the same coal mine as the final spectral data. Figure2 shows the process of data collection and processing, during which whiteboard measurement and calibration are carried out every 10 minutes. Because the original spectral data included 1024 bands, some of which are overlapped bands, it is necessary to remove the overlapped bands from the data, and finally, 973 useful bands are obtained.Figure 2
Spectrum experiment and spectrum curve of coal.


## 4.2. Evaluation of Testing Accuracy

We use P-ELM and P-MELM to model the classification problem. Table2 shows the parameters of P-ELM and P-MELM. Figures 3 and 4, respectively, show the accuracy of data classification by P-ELM and P-MELM under different hidden nodes. From Figures 3 and 4, we can see that P-ELM and P-MELM can obtain good learning accuracy and have a relatively accurate classification effect. If the redundant hidden nodes in the neural network are pruned, the testing accuracy of P-ELM and P-MELM enhance or remain unchanged, and this phenomenon is especially evident in Diabetic Retinopathy Debrecen, image segmentation, and coal spectral datasets. More importantly, the classification accuracy surfaces of P-MELM keep relatively smooth, but the curves of P-ELM have some oscillations in Diabetic Retinopathy Debrecen and image segmentation cases when the redundant hidden nodes are pruned. By comparing the algorithm process of P-ELM and P-MELM, it can be seen that P-MELM changes from a single hidden layer to multiple hidden layers on the basis of P-ELM, and the calculation method of the parameters in the P-MELM algorithm is different. So, P-MELM can eliminate the fluctuations caused by P-ELM to some extent. As a consequence, the P-MELM algorithm is more accurate in classification than the P-ELM algorithm, and the fitting effect is better.Table 2
Parameters of P-ELM and P-MELM under different classification problems.
DatasetsAlgorithmInitial nodesEnd nodesHidden layersDiabetic Retinopathy DebrecenP-ELM100501P-MELM100802–10Image segmentationP-ELM100801P-MELM80552–10Coal spectral dataP-ELM100751P-MELM80552–10Figure 3
The average testing accuracy for the algorithm P-ELM using (a) Diabetic Retinopathy Debrecen, (b) image segmentation, and (c) coal spectral data, while the hidden node is decreased one by one.
(a)(b)(c)Figure 4
The average testing accuracy for the algorithm P-MELM with different hidden layers using (a) Diabetic Retinopathy Debrecen, (b) image segmentation, and (c) coal spectral data, while the hidden node is decreased one by one. Here, M and L are the number of hidden layers and the number of redundant nodes removed, respectively.
(a)(b)(c)


## 4.3. Evaluation of Average Training Time

In Section 4.3, we compare the training time of MELM and P-MELM under the same parameters. The time required for MELM and P-MELM algorithms is shown in Figure5. From Figure 5, we can see that whether it is coal mine classification, diabetic retinopathy classification, or image segmentation classification, and the training time required for P-MELM is much less than the training time required for MELM. This is because when the number of hidden nodes changes, MELM needs to retrain the network, and P-MELM does not need to recalculate the network parameters, but only updates the weights in the model. Therefore, under the same network structure, P-MELM has a faster convergence speed than MELM. Furthermore, the training time curves for several cases shown in Figure 4 also indicate that the spent learning time is still monotonically decreasing with the learning steps, which is consistent with MELM’s conclusions. Compared with P-MELM, MELM needs more time to retrain the network and recalculate the output weights. However, with the increase of hidden layer nodes, the time difference of training time will gradually decrease.Figure 5
Comparison of training time between the algorithms MELM and P-MELM using (a) Diabetic Retinopathy Debrecen, (b) image segmentation, and (c) coal spectral data, while the hidden node is decreased one by one.
(a)(b)(c)


## 4.4. Comparison of Testing Accuracy by Different Methods

The classification accuracy of the testing data is selected as the performance evaluation criterion for the algorithms investigated, namely, ELM, MELM, P-ELM, P-MELM, ELM-AE, and SELM. The best learning accuracy that six algorithms can achieve can be observed from Table3. On the basis of these experimental results, we can draw the following conclusions. If the initial network structure is fixed, P-ELM and P-MELM always get higher testing accuracy than ELM and MELM, respectively. This means that the P-ELM and P-MELM techniques can get more stable networks in most cases. What is more, a higher classification accuracy indicates a better classification algorithm performance. As observed from Table 3, it is obvious that the proposed P-MELM technique achieves the highest classification accuracy among all classification datasets, relative to the original ELM, MELM, and P-ELM algorithms. Finally, these results further prove that the P-ELM and P-MELM methods have a better classification effect than the ELM method and the training speed can meet the requirements.Table 3
The best testing accuracy that the algorithms ELM, MELM, P-ELM, P-MELM, ELM-AE, and SELM can achieve.
DatasetsBest testing accuracy (%)ELMMELMP-ELMP-MELMELM-AESELMDiabetic Retinopathy Debrecen63.6665.9767.5968.7568.0666.20Image segmentation75.0875.8981.3183.0482.493.0Coal spectral data89.9189.9191.6092.2491.6067.23


## 5. Conclusions and Discussion

Based on MELM, this paper introduces the L2 penalty factor between the hidden layer and the output layer, uses the regularization parameter C to measure the structural risk and empirical risk of the MELM model, and proposes the P-MELM algorithm. P-MELM uses the Cholesky factorization method to effectively reduce the amount of calculation for a single solution of the output weights’ matrix, and when pruning hidden nodes, P-MELM can also update the original connection weights without the need for network recalculation, so P-MELM greatly reduces the amount of calculation and improves the calculation efficiency. P-MELM introduces the L2 regularization term into MELM to reduce the sensitivity of nodes. Starting from the initial large number of hidden nodes, by minimizing its prediction error and regularization control item, redundant nodes are rationally pruned to obtain a more concise network. Finally, data simulation indicates that P-MELM features are fast, precise, and reliable.



---

*Source: 5588426-2021-04-27.xml*