# Incremental Multiple Hidden Layers Regularized Extreme Learning Machine Based on Forced Positive-Definite Cholesky Factorization

Authors: Jingyi Liu; Ba Tuan Le
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

The theory and implementation of extreme learning machine (ELM) prove that it is a simple, efficient, and accurate machine learning method. Compared with other single hidden layer feedforward neural network algorithms, ELM is characterized by simpler parameter selection rules, faster convergence speed, and less human intervention. The multiple hidden layer regularized extreme learning machine (MRELM) inherits these advantages of ELM and has higher prediction accuracy. In the MRELM model, the number of hidden layers is randomly initiated and fixed, and there is no iterative tuning process. However, the optimal number of hidden layers is the key factor to determine the generalization ability of MRELM. Given this situation, it is obviously unreasonable to determine this number by trial and random initialization. In this paper, an incremental MRELM training algorithm (FC-IMRELM) based on forced positive-definite Cholesky factorization is put forward to solve the network structure design problem of MRELM. First, an MRELM-based prediction model with one hidden layer is constructed, and then a new hidden layer is added to the prediction model in each training step until the generalization performance of the prediction model reaches its peak value. Thus, the optimal network structure of the prediction model is determined. In the training procedure, forced positive-definite Cholesky factorization is used to calculate the output weights of MRELM, which avoids the calculation of the inverse matrix and Moore-Penrose generalized inverse of matrix involved in the training process of hidden layer parameters. Therefore, FC-IMRELM prediction model can effectively reduce the computational cost brought by the process of increasing the number of hidden layers. Experiments on classification and regression problems indicate that the algorithm can be effectively used to determine the optimal network structure of MRELM, and the prediction model training by the algorithm has excellent performance in prediction accuracy and computational cost.


---

## Body

## 1. Introduction

The neural network is a complex nonlinear system interconnected by a large number of neurons, and it is based on the research of human brain information processing ability by modern neurobiology and cognitive science. The neural network is also a mathematical simulation form of human brain physiological structure, with strong adaptability, self-learning ability, and nonlinear mapping, and it has been widely used by researchers in many scientific fields [1–3]. However, the above prediction models are all based on the traditional neural networks, and the network training process needs to modify the network weights repeatedly according to the training objectives and gradient information. The entire network training process usually takes hundreds or even thousands of iterations before it can be finally completed, which requires a large amount of calculation.Extreme learning machine (ELM) is a novel single hidden layer feedforward neural network. It transforms the iterative adjustment process of traditional neural network parameter training into solving linear equations. According to Moore-Penrose generalized inverse matrix theory, the least squares solution with the minimum norm is obtained analytically as the network weights. The whole training process can be completed in one time without iteration. Compared with the traditional neural network training algorithm, which requires several iterations to determine the network weights, the training speed of ELM is significantly improved [4, 5]. This advantage enables ELM to be successfully applied in pattern recognition [6, 7] and regression estimation [8–10]. In order to improve the generalization ability of ELM, the literature [11] draws on the principle of structural risk minimization in statistical learning theory and proposes a regularized extreme learning machine (RELM). RELM has a better generalization ability by introducing parameters λ to weigh structural risks and empirical risks [12–15]. For the single hidden layer RELM model with multiple input and single output, the literature [16] designed the Cholesky factorization method for regularized output weight matrix. In the learning and forgetting process of the sample sequence, the Cholesky factorization factor is calculated recursively by adding and deleting samples one by one, and then the output weights are adjusted, and the network structure is fixed. However, if dealing with input data with complex noise signals and high-dimensional information, or with more classification categories, RELM also shows its own shortcomings, and the accuracy of the established model is greatly reduced.In order to improve the embarrassing situation of RELM, the literature [17, 18] starts from improving its network structure. On the basis of the traditional RELM three-layer structure, the number of hidden layers is increased to form a neural network with one input layer, multiple hidden layers, and one output layer, that is, the multiple hidden layers RELM network model (MRELM), in which the neuron nodes of each hidden layer are fully connected. MRELM inherits the idea that RELM randomly initializes the weights matrix between the input layer and the hidden layer as well as the bias vector of the hidden layer. By forcing the actual output of the hidden layer to be as close as possible to the expected output, the weights matrix and the bias vector of the added hidden layers are calculated; thereby a neural network model with multiple hidden layers is established. The parameter training process needs to calculate the inverse matrix and the MP generalized inverse matrix, in which the first hidden layer parameters are randomly initialized, and the remaining hidden layers parameters are obtained by minimizing the error between the actual output and the expected output of the corresponding hidden layer. Compared with the traditional RELM model, MRELM can effectively improve the prediction accuracy through the layer-by-layer optimization of network parameters between different hidden layers. Moreover, it has the advantages of strong generalization ability and fast computing speed and is not easy to fall into local optimum [19–22]. However, since the initial parameter values of MRELM are randomly initialized, although this avoids the situation that the algorithm falls into local optimum and overfitting, it also leads to the failure of some hidden layers or the reduction of the effect on the neural network during the modeling process. As a result, there are some redundant hidden layers in the MRELM network, which often require more reasonable selection methods and theories for the number of hidden layers. Meanwhile, the network structure of MRELM is determined by the users based on their own practical experience, but this empirical choice is not reasonable, and it is difficult to guarantee the optimality. In practical applications, users often need to carry out repeated experiments for many times and choose the network structure with the least time consuming and the highest accuracy from the complex results comparison as the optimal network model for the training and prediction of the actual data.In order to realize the effective design of MRELM network structure, select the number of hidden layers reasonably, and achieve the desired accuracy requirements, an incremental MRELM training algorithm based on forced positive definite Cholesky factorization (FC-IMRELM) [23, 24] is put forward in this paper. The algorithm can adjust the number of hidden layers in the network adaptively according to the predicted data, so as to determine the optimal network structure of FC-IMRELM. At the same time, a novel method is adopted to calculate the parameters of the newly added hidden layers, that is, the connection weight matrix and the bias vector of the hidden layers. Based on previous research, MRELM typically requires fewer hidden neurons than ELM to achieve a desirable performance level. This is a basic requirement for considering the multiple-hidden-layers structure presented. The foundational ideas for the FC-IMRELM algorithm are far simpler to produce and more stable by comparing and contrasting its characteristics with other ELM variants. Experimental results for classification and regression problems show that the proposed FC-IMRELM algorithm has more advantages in terms of average accuracy compared to the traditional RELM model and other improved models of MRELM.The rest of this paper is organized as follows: Section2 presents a brief review of the basic concepts and related work of multiple hidden layers RELM, Section 3 describes the proposed incremental FC-IMRELM technique, Section 4 reports and analyzes the experimental results, and, finally, Section 5 summarizes key conclusions of the present study.


## 2. Brief Review of Multiple Hidden Layers Regularized Extreme Learning Machine

The MRELM algorithm tries to find a mapping relationship that makes the output predicted by the ELM neural network with multiple hidden layers infinitely close to the actual given result. This mapping relationship will be embodied in the solution process of the weight and bias parameters of the hidden layers. The number of hidden layers in the MRELM neural network needs to be selected according to the change of the predicted data. Therefore, in the training process of network parameters, in order to ensure that the final hidden layer output is closer to the expected hidden layer output, in addition to the random initialization of the parameters of the first hidden layer, the parameter training process starts from the second hidden layer to optimize the network parameters until all the network parameters are completed. Furthermore, during the establishment of the neural network, the weight matrix and bias vector of each hidden layer are acquired and recorded, so as to obtain the final predicted output result of the MRELM neural network. The solving process of the network parameters will be explained in detail in the following algorithm flow.Suppose that a set of training sample dataset given in MRELM neural network isX,T=xi,ti(i=1,2,⋯,N), where X=[x1,x2,⋯,xN]T is the input samples, xi=xi1,xi2,⋯,xinT∈Rn is the n×1 input vector, T=[t1,t2,⋯,tN]T is the corresponding labeled samples, ti=ti1,ti2,⋯,timT∈Rm is the m×1 observation vector, and N is the total number of training samples. Meanwhile, it is assumed that all hidden layers in the MRELM model contain the same number of hidden nodes L, and each hidden node chooses the same activation function g(x). In the modeling process of MRELM algorithm, multiple hidden layers in the neural network are first treated as a single hidden layer, and then the hidden layer parameters in the MRELM network containing only a single hidden layer are randomly initialized, namely, the input weights matrix W1=[(W1)1,(W1)2,⋯,(W1)L]T∈RL×n connecting the input layer and the first hidden layer, and the bias vector B1=[b11,b12,⋯,b1L]T∈RL of the first hidden nodes. Thus, the output matrix of the first hidden layer H1 can be calculated as follows:(1)H1=gW1X+B1whose scalar entries (h1)ij=g(W1jxi+b1j)(i=1,2,⋯,N,j=1,2,⋯,L) are interpreted as the output of the j-th hidden node in the first hidden layer with respect to xi and (W1)j=[(W1)j1,(W1)j2,⋯,(W1)jn]T is the vector of connection weights between n input nodes and the j-th hidden node in the first hidden layer. To better balance the empirical risk and structural risk, the MRELM adjusts the proportion of the two risks by introducing parameter λ, which can be expressed as the following constrained optimization problem.(2)min12β12+λ2ε2s.t.T=H1β1-εwhere β1=[(β1)1,(β1)2,⋯,(β1)L]T∈RL×m is the connection weights matrix between the first hidden layer and the output layer, with vector components (β1)j=[(β1)j1,(β1)j2,⋯,(β1)jm]T(j=1,2,⋯,L) that denote the connection weights between the j-th hidden node in the first hidden layer and m output nodes, ε=[ε1,ε2,⋯,εN]T denotes the training error, and λ>0 is the regularization parameter.According on the KKT theorem, the constrained optimization of (2) can be transformed into the following dual optimization problem:(3)Lα,β1,ε=12β12+λ2ε2-αH1β1-T-εwhere α=[α1,α2,⋯,αN] is the Lagrange multipliers vector. Utilizing KKT optimality conditions, the following equations can be obtained:(4a)∂L∂β1=0⇒β1T=αH1(4b)∂L∂ε=0⇒λεT+α=0(4c)∂L∂α=0⇒H1β1-T-ε=0Finally, β1 can be gotten as follows:(5a)β1=λ-1I+H1TH1-1H1TTor(5b)β1=H1Tλ-1I+H1H1T-1TIn order to reduce the computational costs, if N>L, one may prefer to apply the solution (5a), and if N<L, one may prefer to apply the solution (5b).Now the second hidden layer is added to the MRELM neural network, the network structure with two hidden layers is restored, and the two hidden layers are fully connected, so the prediction output of the second hidden layerH2 can be obtained as follows:(6)H2=gW2H1+B2where W2 denotes the weights matrix between the first hidden layer and the second hidden layer. We suppose that the first and second hidden layers have the same number of nodes, and thus W2 is a square matrix. The matrix B2 represents the bias of the second hidden layer. The expected output of the second hidden layer H2∗ can be calculated as(7)H2∗=Tβ1+where β1+ is the MP generalized inverse of the matrix β1, which can be calculated using the orthogonal projection method. Namely, if (β1)Tβ1 is nonsingular, then (β1)+=((β1)Tβ1)-1(β1)T; otherwise (β1)+=(β1)T(β1(β1)T)-1 if β1(β1)T is nonsingular. To make the predicted output of the hidden layer in the MRELM neural network infinitely close to the expected output, we may set H2=H2∗. Subsequently, we define the augmented matrix (W2)HE=[B2W2], and it can be gotten as(8)W2HE=g-1H2∗H2E+where (H2E)+ is the MP generalized inverse of the matrix H2E=[1H1]T, and 1 represents a one-column vector of size N whose elements are the scalar unit 1. The solving method of (H2E)+ is the same as previously discussed for (β1)+. The notation g-1(x) indicates the inverse of the activation function g(x). For classification and regression problems, we all invoke the widely used logistic sigmoid function g(x)=1/(1+e-x). The predicted output of the second hidden layer H2 is obtained as (9)H2=gW2H1+B2=gW2HEH2ETherefore, the connection weights matrix β2 between the second hidden layer and the output layer is calculated as(10a)β2=λ-1I+H2TH2-1H2TTor(10b)β2=H2Tλ-1I+H2H2T-1TThe solving method of β2 is chosen according to what is previously discussed for β1.According to the MRELM algorithm flow, the third hidden layer is added to the MRELM network, and restore the network structure with three hidden layers. Since the nodes between each hidden layer are all connected together, the prediction output of the third hidden layer can be obtained as(11)H3=gW3H2+B3where W3 represents the weights matrix between the second hidden layer and the third hidden layer, and the vector B3 denotes the bias of the third hidden layer. Thus, the expected output of the third hidden layer can be gotten as (12)H3∗=Tβ2+where β2+ is the MP generalized inverse of the weights matrix β2, obtained using the approach described before. To meet the requirement that the predicted output of the third hidden layer is infinitely close to the expected output, let H3=H3∗. Accordingly, the augmented matrix can be defined as (W3)HE=[B3W3], and we can solve it as follows.(13)W3HE=g-1H3∗H3E+where (H3E)+ is the MP generalized inverse of the matrix H3E=[1H2]T, the specific meaning of the symbol 1 is described above, and the calculation of (H3E)+ also proceeds in the manner discussed before. Therefore, we can update the predicted output of the third hidden layer as(14)H3=gW3H2+B3=gW3HEH3EFinally, the connection weight matrix β3 between the third hidden layer and the output layer can be calculated as(15a)β3=λ-1I+H3TH3-1H3TTor(15b)β3=H3Tλ-1I+H3H3T-1TThe calculation approach of β3 is still selected according to the principle of β1 discussed previously. The final output of the MRELM network with three hidden layers after training can be expressed as(16)fx=H3β3If the number of hidden layersl in the MRELM network is more than 3, an iterative format can be adopted to realize the calculation process. In other words, the iterative calculation of formula (6) to formulas (15a) and (15b) is performed for l-3 times until all hidden layer parameters are solved. Emphasized finally, this algorithm does not add all hidden layers to the network at one time, nor does it calculate all hidden layer parameters at one time, but one hidden layer after another is added to the network. Every time a new hidden layer is added, the weights matrix and the bias vector of the hidden layers are calculated immediately to prepare for the parameter calculation of the hidden layer to be added next time.


## 3. Solutions of IMRELM by the Forced Positive-Definite Cholesky Factorization

For the single hidden layer feed-forward neural network, the literature [15] puts forward the regularized extreme learning machine algorithm based on Cholesky factorization (CF-FORELM), introduces the Cholesky factorization of positive definite matrix into the solving process of RELM, and designs a recursive solution method for the calculation of the regularized output matrix Cholesky factorization factor. The advantages of CF-FORELM algorithm prompted us to introduce the forced positive-definite Cholesky factorization method into the framework of MRELM algorithm with multiple hidden layers, and we then proposed an MRELM neural network training algorithm based on forced positive definite Cholesky factorization (FC-IMRELM). Compared with the inverse matrix calculation of invertible matrix λ-1I+(Hl)THl in traditional RELM algorithm and the calculation of MP generalized inverse of matrix (Hl)+ in MELM algorithm (l denotes the number of hidden layers), the algorithm effectively reduces the computational cost and complexity brought by the matrix inverse process. Meanwhile, the numerical stability of the forced positive definite Cholesky factorization method also greatly weakens the randomness effect of the ELM algorithm on the prediction results.


### 3.1. Forced Positive-Definite Cholesky Factorization (FC)

The main difficulty of the MRELM algorithm is the calculation of the inverse matrix and the MP generalized inverse matrix involved in the training process, including the inverse calculation of symmetric positive semidefinite matrix. In this case, the improved MRELM mode based on the traditional Cholesky factorization could not be realized, because the Cholesky factorization of the symmetric positive semidefinite matrix might not exist. Even if such a factorization exists, the calculation process is generally numerically unstable for the elements of the matrix factorization factor may be unbounded. In order to overcome these difficulties, we put forward a modified approach based on the forced positive definite Cholesky factorization for the MRELM algorithm with multiple hidden layers, which is a numerical stability approach.When the forced positive definite strategy is adopted to improve the MRELM algorithm, the key problem is how to form the positive definite matrix from the modified Cholesky decomposition of the undetermined matrix. If the matrixG is not a positive definite matrix, the Cholesky factorization method, which forces the matrix to have positive definite property, is to find a unit lower triangular matrix L and a positive definite diagonal matrix D for the general symmetric matrix G, so that the matrix G-=LDLT is positively definite, and it is only one diagonal matrix E away from the matrix G.(17)G-=LDLT=G+EIn fact, the Cholesky factorization of symmetric positive definite matrix can be described as follows:(18a)djj=gjj-∑s=1j-1dssljs2(18b)lij=1djjgij-∑s=1j-1dssljslis,i≥j-1where gij represents the element of matrix G and djj denotes the main diagonal element of matrix D. Here, the Cholesky factorization factors L and D are required to satisfy two requirements: one is that all elements of D are strictly positive, and the other is that the elements of the factorization factor L are uniformly bounded. That is, for k=1,⋯,n and a positive number ρ, the formula (19) is required:(19)dkk>δ,rik≤ρ,i>kwhere the auxiliary quantity rik=likdkk, δ is a given small positive number. The matrix satisfying the above conditions is said to be sufficiently positive definite, where E is a zero matrix.Next, we describe thej-th step of this factorization. Suppose the j-1 column of the forced positive-definite Cholesky factorization has been calculated. For k=1,⋯,j-1, equation (19) holds. First calculate(20)γj=ξj-∑s=1j-1dssljs2,where ξj is taken as gjj, and the test value d- is defined as(21)d-=max⁡γj,δwhere δ is a small positive number. In order to determine whether d- can accept as the j-th element of D, we check whether rij=lijd- satisfies the formula (19). If so, let djj=d-, and get the j-th column of L from lij=rij/djj. Otherwise,(22)djj=ξj-∑s=1j-1dssljs2let ξj=gjj+ejj, select positive number ejj to make max⁡rij=ρ, and produce the j-th column of L.If the above process is completed, we obtain the Cholesky factorization formula (17) of the positive definite matrix G-, where E is a nonnegative diagonal matrix and the diagonal element is ejj. For the given matrix G, this nonnegative diagonal matrix E depends on ρ. If n>1, (23)Eρ∞≤ξρ+n-1ρ2+2γ+n-1ρ2δwhere ξ is the maximum norm of the nondiagonal elements of G and γ is the maximum norm of the diagonal elements of G. If ρ2=ξ/n2-1, the upper bound is minimized. So, let ρ satisfy formula (24):(24)ρ2=max⁡γ,ξn2-1,εMwhere εM represents the machine precision. We increase εM to prevent G from being small.Finally, we present the forced positive-definite Cholesky factorization algorithm, where the auxiliary quantitycis=lisdss, s=1,⋯,j, i=j,⋯,n. These values need not be stored separately; they can be stored in the matrix G.Algorithm 1.
Forced Positive-Definite Cholesky Factorization (FC)
S tep 1. Calculate the bounds of the elements of the factorization factor. Let ρ2=max⁡γ,ξ/v,εM, where v=max⁡1,n2-1 and γ and ξ are the maximum norm of the diagonal and nondiagonal elements of G, respectively.
S tep 2. Initialization. Let j=1, cii=gii, i=1,⋯,n.
S tep 3. Determine the minimum index q, so that cqq=maxj≤i≤n⁡cii, and exchange the information of q rows and j rows, q columns and j columns of G.
S tep4. Calculate the j-th row of L, and solve the maximum norm of lijdjj. Let ljs=cjs/dss, s=1,⋯,j-1, calculate cij=gij-∑s=1j-1ljscis, i=j+1,⋯,n, and let θj=maxj+1≤i≤n⁡cij. If j=n, let θj=0.
S tep 5. Calculate the j-th diagonal element djj=max⁡δ,cjj,θj2/ρ2 of D. The diagonal element of E is modified to Ej=djj-cjj. If j=n, stop.
S tep 6. Correct the diagonal elements and column index, and let cii=cii-cij2/djj, i=j+1,⋯,n and j=j+1; jump to Step 3.


### 3.2. Process of Matrices Decomposition for FC-IMRELM

According to the MRELM training process shown in equations (1) to (16), its essence is to solve the connection weight matrix βl between the hidden layer and the output layer. However, it can be seen from equations (5a), (5b), (10a), (10b), (15a), and (15b) that the solution method given in literature [18] involves matrix inversion, and the solution process of each hidden layer's parameters (Wl)HE (l is the number of hidden layers) involves the calculation of the MP generalized inverse of matrix. The problem of large amount of calculation reduces the modeling efficiency of the MRELM prediction model. In order to solve the above problems effectively, we propose a solution method of the weights matrix βl and hidden layer parameters (Wl)HE based on the forced positive definite Cholesky factorization.First, on the basis of equations (4a), (4b), and (4c), β1=(H1)TαT can be obtained from equation (4a), and ε=-αT/λ can be obtained from equation (4b), and then substituting β1 and ε into equation (4c), we can obtain(25)Aβ1k=bAk,k=1,2,⋯,m.where A=λ-1I+(H1)TH1, ((bA)(1),(bA)(2),⋯,(bA)(m))=(H1)TT, and β1=((β1)(1),(β1)(2),⋯,(β1)(m)), and(26)BUBk=Tk,k=1,2,⋯,mwhere B=λ-1I+H1(H1)T, T=(T(1),T(2),⋯,T(m)), and β1=(H1)T((UB)(1),(UB)(2),⋯,(UB)(m)).Therefore, the process of solvingβ1 according to equation (5a) can be transformed into solving m linear equations in the form of equation (25), and the process of solving β1 according to equation (5b) can be transformed into solving m linear equations in the form of equation (26); m is the dimension of observation vector.If the number of training samples is greater than the number of hidden nodes, that is,N>L, the solution process of β1 based on the Cholesky factorization is as follows. First, calculate the Cholesky factorization of matrix A:(27)A=SASATwhere SA is a lower triangular matrix with positive diagonal elements. The nonzero element (sA)ij in SA can be calculated by the element aij of A according to equation (28):(28)sAij=aii-∑n=1i-1sAin2i=j,aij-∑n=1j-1sAinsAjnsAjji>jwhere i=1,2,⋯,L, j=1,2,⋯,L. Substitute equation (27) into equation (25) and multiply both sides of the equation by (SA)-1: (29)SATβ1k=VAk,k=1,2,⋯,m,where (VA)(k)=(SA)-1(bA)(k). Solving β1=((β1)(1),(β1)(2),⋯,(β1)(m)) is equivalent to solving equation (29). Because (VA)(k)=(SA)-1(bA)(k) is equivalent to SA(VA)(k)=(bA)(k), the calculation formula for the elements (vA)i(k) of (VA)(k) can be obtained by comparing the elements on both sides of the equation (30)vAik=bAiksAiii=1,bAik-∑n=1i-1sAnivAnksAiii>1where i=1,2,⋯,L, (bA)i(k) is the element of (bA)(k). Finally, on the basis of obtaining SA and (VA)(k), the element (β1)i(k) of (β1)(k) can be calculated by using the elements of SA and (VA)(k): (31)β1ik=vAiksAiii=L,vAik-∑n=1L-isAi+n,iβ1i+nksAiii<Lwhere i=1,2,⋯,L, k=1,2,⋯,m, β1=((β1)(1),(β1)(2),⋯,(β1)(m)).If the number of hidden nodes is greater than the number of training samples, that is,N<L, the solution process of β1 based on the Cholesky factorization is as follows. First, calculate the Cholesky factorization of B:(32)B=SBSBTwhere SB is a lower triangular matrix with positive diagonal elements. The nonzero element (sB)ij in SB can be calculated by the element bij of B according to equation (33): (33)sBij=bii-∑n=1i-1sBin2i=j,bij-∑n=1j-1sBinsBjnsBjji>jwhere i=1,2,⋯,N, j=1,2,⋯,N. Substitute equation (32) into equation (26) and multiply both sides of the equation by (SB)-1: (34)SBTUBk=VBk,k=1,2,⋯,m,where (VB)(k)=(SB)-1T(k). Solving β1=(H1)T((UB)(1),(UB)(2),⋯,(UB)(m)) is equivalent to solving equation (34). Because (VB)(k)=(SB)-1T(k) is equivalent to SB(VB)(k)=T(k), the calculation formula for the element (vB)i(k) of (VB)(k) can be obtained by comparing the elements on both sides of the equation (35)vBik=TiksBiii=1,Tik-∑n=1i-1sBnivBnksBiii>1where i=1,2,⋯,N, Ti(k) is the element of T(k). Finally, on the basis of obtaining SB and (VB)(k), the element (uB)i(k) of (UB)(k) can be calculated by using the elements of SB and (VB)(k). (36)uBik=vBiksBiii=L,vBik-∑n=1L-isBi+n,iuBi+nksBiii<Lwhere i=1,2,⋯,N, k=1,2,⋯,m, β1=(H1)T((UB)(1),(UB)(2),⋯,(UB)(m)).At this point, we get the connection weight matrixβ1 between the first hidden layer and the output layer. The connection weight matrix βl between the other hidden layer and output layer can also be calculated by the above method. Compared with the solution method of the connecting weight matrix βl as shown in equations (5a), (5b), (10a), (10b), (15a), and (15b), the solution of βl based on Cholesky factorization does not involve the inverse operation of the matrix, and it can be achieved by using simple algebraic operations.The MRELM model contains multiple hidden layers, and the solving process of each hidden layer parameter(Wl)HE needs to calculate the MP generalized inverse of the corresponding matrices. However, it can be concluded from equations (7), (8), (11), and (13) that the solution method of (βl)+ and (HlE)+ using orthogonal projection method [5] involves the inverse calculation of symmetric semipositive definite matrixes (βl)Tβl, βl(βl)T, (HlE)THlE, and HlE(HlE)T, which has the problem of large computational cost and numerical instability, and if the condition number of the above matrices is too large, the calculation results of MP generalized inverse of matrices (βl)+ and (HlE)+ are usually unable to be obtained. This not only affects the modeling efficiency and prediction effect of MRELM model, but also may make the modeling process impossible to complete. However, the traditional Cholesky factorization method can only be used to solve the calculation of symmetric positive definite matrix. To effectively overcome the above difficulties, we use the forced positive definite Cholesky factorization to solve the MP generalized inverse of matrices (βl)+ and (HlE)+, and we then get the hidden layer parameters (Wl)HE.If(βl)Tβl is nonsingular, then take (βl)+=((βl)Tβl)-1(βl)T and substitute it into (7) and (11): (37)Cβl+k=βlTk,k=1,2,⋯,L,where C=(βl)Tβl, H(l+1)∗=T(((βl)+)(1),((βl)+)(2),⋯,((βl)+)(L)).Ifβl(βl)T is nonsingular, then take (βl)+=(βl)T(βl(βl)T)-1 and substitute it into (7): (38)DTβl+Tk=βlk,k=1,2,⋯,m,where D=βl(βl)T, H(l+1)∗=T((((βl)+)T)(1),(((βl)+)T)(2),⋯,(((βl)+)T)(m))T.Therefore, the process of solvingH(l+1)∗ can be transformed into solving L linear equations in the form of equation (37), or solving m linear equations in the form of equation (38), L is the number of hidden nodes, and m is the dimension of the observation vector.If(βl)Tβl is nonsingular, the solution process of H(l+1)∗ based on forced positive definite Cholesky factorization is as follows. First, calculate the modified Cholesky factorization result of matrix C:(39a)C=PCMCPCTwhere PC is the unit lower triangular matrix and MC is the positive definite diagonal matrix. (pC)ij is the nonzero element below the main diagonal in PC, and (mC)ii is the main diagonal element in MC; it can be calculated according to Algorithm 1 by using element cij of C, where i=1,2,⋯,m, j=1,2,⋯,m. If we set MC=(MC)1/2(MC)1/2, then (39a) can be written as(39b)C=QCQCTwhere QC=PC(MC)1/2. Substitute equation (39b) into equation (37) and multiply both sides of the equation by (QC)-1. (40)QCTβl+k=VCk,k=1,2,⋯,L,where (VC)(k)=(QC)-1((βl)T)(k). Solving H(l+1)∗=T(((βl)+)(1),((βl)+)(2),⋯,((βl)+)(L)) is equivalent to solving equation (40). Since (VC)(k)=(QC)-1((βl)T)(k) is equivalent to QC(VC)(k)=((βl)T)(k), by comparing the elements on both sides of the equation, the calculation formula for the element (vC)i(k) of (VC)(k) can be obtained: (41)vCik=βlTikqCiii=1,βlTik-∑n=1i-1qCnivCnkqCiii>1where ((βl)T)i(k) is the element of ((βl)T)(k), and (qC)ij is the element of QC, i=1,2,⋯,m, j=1,2,⋯,m. Finally, on the basis of obtaining QC and (VC)(k), the element ((βl)+)i(k) of ((βl)+)(k) can be calculated by using the elements of QC and (VC)(k). (42)βl+ik=vCikqCiii=L,vCik-∑n=1L-iqCi+n,iβl+i+nkqCiii<Lwhere i=1,2,⋯,m, k=1,2,⋯,L, (βl)+=(((βl)+)(1),((βl)+)(2),⋯,((βl)+)(L)); consequently, we can get H(l+1)∗=T(((βl)+)(1),((βl)+)(2),⋯,((βl)+)(L)).Ifβl(βl)T is nonsingular, the solution process of H(l+1)∗ based on forced positive definite Cholesky factorization is as follows. First, calculate the modified Cholesky factorization results of matrix D: (43a)D=PDMDPDTwhere PD is the unit lower triangular matrix and MD is the positive definite diagonal matrix. (pD)ij is the nonzero element below the main diagonal in PD, and (mD)ii is the main diagonal element in MD; it can be calculated according to Algorithm 1 by using the element dij of D, i=1,2,⋯,L, j=1,2,⋯,L. If we set MD=(MD)1/2(MD)1/2, then (43a) can be written as(43b)D=QDQDTwhere QD=PD(MD)1/2. Substitute equation (43b) into equation (38) and multiply both sides of the equation by (QD)-1. (44)QDTβl+Tk=VDk,k=1,2,⋯,m,where (VD)(k)=(QD)-1(βl)(k). Solving H(l+1)∗ = T((((βl)+)T)(1),(((βl)+)T)(2),⋯,(((βl)+)T)(m))T is equivalent to solving equation (44). Since (VD)(k)=(QD)-1(βl)(k) is equivalent to QD(VD)(k)=(βl)(k), by comparing the elements on both sides of the equation, the calculation formula for the element (vD)i(k) of (VD)(k) can be gotten as (45)
v
D
i
k
=
β
l
i
k
q
D
i
i
i
=
1
,
β
l
i
k
-
∑
n
=
1
i
-
1
q
D
n
i
v
D
n
k
q
D
i
i
i
>
1where (βl)i(k) is the element of (βl)(k) and (qD)ij is the element of QD, i=1,2,⋯,L, j=1,2,⋯,L. Finally, on the basis of obtaining QD and (VD)(k), the element (((βl)+)T)i(k) of (((βl)+)T)(k) can be obtained by using elements of QD and (VD)(k). (46)βl+Tik=vDikqDiii=L,vDik-∑n=1L-iqDi+n,iβl+Ti+nkqDiii<Lwhere i=1,2,⋯,L, k=1,2,⋯,m, (βl)+=((((βl)+)T)(1),(((βl)+)T)(2),⋯,(((βl)+)T)(m))T; therefore, we can get H(l+1)∗=T((((βl)+)T)(1),(((βl)+)T)(2),⋯,(((βl)+)T)(m))T.If(HlE)THlE is nonsingular, then let (HlE)+=((HlE)THlE)-1(HlE)T and substitute it into (8) and (13): (47)EUEk=g-1Hl∗Tk,k=1,2,⋯,L,where E=(HlE)THlE, (Wl)HE=((UE)(1),(UE)(2),⋯,(UE)(L))T(HlE)T.IfHlE(HlE)T is nonsingular, then let (HlE)+=(HlE)T(HlE(HlE)T)-1 and substitute it into (8) and (13): (48)FWlHETk=bFk,k=1,2,⋯,L,where (Wl)HE = ((((Wl)HE)T)(1),(((Wl)HE)T)(2), ⋯, (((Wl)HE)T)(L))T, ((bF)(1),(bF)(2),⋯(bF)(L))T = g-1(Hl∗)(HlE)T, F=HlE(HlE)T.Therefore, the process of solving(Wl)HE can be transformed into solving L linear equations in the form of equation (47), or solving L linear equations in the form of equation (48), where L is the number of hidden nodes.If(HlE)THlE is nonsingular, the solution process of (Wl)HE based on the forced positive definite Cholesky factorization is as follows. First, calculate modified Cholesky factorization results of matrix E: (49a)E=PEMEPETwhere PE is the unit lower triangular matrix and ME is the positive definite diagonal matrix. (pE)ij is the nonzero element below the main diagonal in PE, and (mE)ii is the main diagonal element in ME; it can be calculated according to Algorithm 1 by using element Eij of E, i=1,2,⋯,N, j=1,2,⋯,N. If we set ME=(ME)1/2(ME)1/2, then (49a) can be written as(49b)E=QEQETwhere QE=PE(ME)1/2. Substitute equation (49b) into equation (47) and multiply both sides of the equation by (QE)-1. (50)QETUEk=VEk,k=1,2,⋯,L,where (VE)(k)=(QE)-1((g-1(Hl∗))T)(k). Solving (Wl)HE=((UE)(1),(UE)(2),⋯,(UE)(L))T(HlE)T is equivalent to solving equation (46). For (VE)(k)=(QE)-1((g-1(Hl∗))T)(k) is equivalent to QE(VE)(k)=((g-1(Hl∗))T)(k), by comparing the elements on both sides of the equation, the calculation formula for the element (vE)i(k) of (VE)(k) can be calculated as (51)vEik=g-1Hl∗TikqEiii=1,g-1Hl∗Tik-∑n=1i-1qEnivEnkqEiii>1where ((g-1(Hl∗))T)i(k) is the element of ((g-1(Hl∗))T)(k) and (qE)ij is the element of QE, i=1,2,⋯,N, j=1,2,⋯,N. Finally, on the basis of obtaining QE and (VE)(k), the element (UE)i(k) of (UE)(k) can be calculated by using the elements of QE and (VE)(k). (52)UEik=vEikqEiii=L,vEik-∑n=1L-iqEi+n,iUEi+nkqEiii<Lwhere i=1,2,⋯,N, k=1,2,⋯,L, UE=((UE)(1),(UE)(2),⋯,(UE)(L)); accordingly, we can obtain the matrix (Wl)HE=((UE)(1),(UE)(2),⋯,(UE)(L))T(HlE)T.IfHlE(HlE)T is nonsingular, the solution process of (Wl)HE based on the forced positive definite Cholesky factorization is as follows. First, calculate modified Cholesky factorization results of matrix F:(53a)F=PFMFPFTwhere PF is the unit lower triangular matrix and MF is the positive definite diagonal matrix. (pF)ij is the nonzero element below the main diagonal in PF, and (mF)ii is the main diagonal element in MF; it can be calculated according to Algorithm 1 by using element Fij of F, i=1,2,⋯,L+1, j=1,2,⋯,L+1. If we set MF=(MF)1/2(MF)1/2, then (53a) can be written as(53b)F=QFQFTwhere QF=PF(MF)1/2. Substitute equation (53b) into equation (48) and multiply both sides of the equation by (QF)-1. (54)QFTWlHETk=VFk,k=1,2,⋯,L,where (VF)(k)=(QF)-1(bF)(k). Because solving (Wl)HE=((((Wl)HE)T)(1),(((Wl)HE)T)(2),⋯,(((Wl)HE)T)(L))T is equivalent to solving equation (54) and (VF)(k)=(QF)-1(bF)(k) is equivalent to QF(VF)(k)=(bF)(k), by comparing the elements on both sides of the equation, the calculation formula for the element (vF)i(k) of (VF)(k) can be calculated as (55)vFik=bFikqFiii=1,bFik-∑n=1i-1qFnivFnkqFiii>1where (bF)i(k) is the element of (bF)(k) and (qF)ij is the element of QF, i=1,2,⋯,L+1, j=1,2,⋯,L+1. Finally, on the basis of obtaining QF and (VF)(k), the element (((Wl)HE)T)i(k) of (((Wl)HE)T)(k) can be calculated by using the elements of QF and (VF)(k). (56)WlHETik=vFikqFiii=L,vFik-∑n=1L-iqFi+n,iWlHETi+nkqFiii<Lwhere i=1,2,⋯,L+1, k=1,2,⋯,L. So we can get (Wl)HE=((((Wl)HE)T)(1),(((Wl)HE)T)(2),⋯,(((Wl)HE)T)(L))T.So far we have obtained the parameters(Wl)HE of each hidden layer. Compared with the solving method of hidden layer parameters (Wl)HE shown in equations (7), (8), (11), and (13), the calculation approach of (Wl)HE based on the forced positive definite Cholesky factorization does not involve the operation of inverse matrix, which can be realized only by simple algebraic calculations. Meanwhile, the method of forcing the matrix to be positive definite matrix can guarantee the numerical stability of MRELM neural network training process.


### 3.3. FC-IMRELM Training Algorithm

Studies have shown that the number of hidden layers determines the learning accuracy and generalization ability of the MRELM model [25], and it is also a key factor that must be determined in advance when designing the MRELM network structure. Due to the complexity of various training samples applied in MRELM prediction model, it is difficult to accurately determine the optimal number of hidden layers by human experience, so that the MRELM prediction model has enough hidden layers to ensure its learning accuracy, while, at the same time, it has as few hidden layers as possible to maintain its contracted network structure. To avoid the disadvantages and difficulties of artificially selecting the number of hidden layers, we propose an incremental MRELM training algorithm based on forced positive Cholesky factorization, which can automatically determine the optimal number of hidden layers in MRELM, and the training process is as follows.Algorithm 2.
Incremental MRELM based on Forced Positive-Definite Cholesky Factorization (FC-IMRELM)
S tep 1. Suppose the training sample dataset is X,T=xi,ti(i=1,2,⋯,N), where X=[x1,x2,⋯,xN]T is the input samples, T=[t1,t2,⋯,tN]T is the labeled samples, each hidden layer contains the number of hidden nodes L, and the activation function is g(x).
S tep 2. Let the number of hidden layers in MRELM l=1, randomly initialize the input weights matrix W1 between the input layer and the first hidden layer and the bias vector B1 of the first hidden layer, and let (W1)HE=[B1W1], H1E=[1X]T.
S tep 3. Calculate the output matrix of the first hidden layer H1=g((W1)HEH1E), and calculate the matrices A, B, and (bA)(k)(k=1,2,⋯,m).
S tep 4. Calculate the connection weights matrix between the first hidden layer and the output layer β1=(λ-1I+(H1)TH1)-1(H1)TT or β1=(H1)T(λ-1I+H1(H1)T)-1T.
(1) IfN>L, the Cholesky factorization factor SA of matrix A is calculated according to the formula (28), and (VA)(k)(k=1,2,⋯,m) is calculated according to the formula (30) using matrices SA and (bA)(k)(k=1,2,⋯,m). Then matrices SA and (VA)(k)(k=1,2,⋯,m) are used to calculate (β1)(k) according to equation (31), and we can get β1=((β1)(1),(β1)(2),⋯,(β1)(m)).
(2) IfN<L, the Cholesky factorization factor SB of matrix B is calculated according to the formula (33), and (VB)(k)(k=1,2,⋯,m) is calculated according to the formula (35) using SB and T(k)(k=1,2,⋯,m). Then the matrices SB and (VB)(k)(k=1,2,⋯,m) are used to calculate (UB)(k)(k=1,2,⋯,m) according to the equation (36); subsequently, we can obtain β1=(H1)T((UB)(1),(UB)(2),⋯,(UB)(m)).
S tep 5. Calculate the expected output H(l+1)∗=T(βl)+ of the (l+1)-th hidden layer.
(1) If(βl)Tβl is nonsingular, the forced positive-definite Cholesky factorization factor QC of matrix C is calculated according to equations (39a) and (39b), and (VC)(k)(k=1,2,⋯,L) is calculated according to equation (41) using QC and ((βl)T)(k)(k=1,2,⋯,L). Then the matrices QC and (VC)(k)(k=1,2,⋯,L) are used to calculate ((βl)+)(k)(k=1,2,⋯,L) according to formula (42), so we can obtain H(l+1)∗=T(((βl)+)(1),((βl)+)(2),⋯,((βl)+)(L)).
(2) Ifβl(βl)T is nonsingular, the forced positive-definite Cholesky factorization factor QD of D is calculated according to equations (43a) and (43b), and (VD)(k)(k=1,2,⋯,m) is calculated according to equation (45) using QD and (βl)(k)(k=1,2,⋯,m). Then the matrices QD and (VD)(k)(k=1,2,⋯,m) are used to calculate (((βl)+)T)(k)(k=1,2,⋯,m) according to formula (46); thus we can get H(l+1)∗=T((((βl)+)T)(1),(((βl)+)T)(2),⋯,(((βl)+)T)(m))T.
S tep 6. Let l=l+1, calculate the parameter (Wl)HE=g-1(Hl∗)(HlE)+ of the l-th hidden layer, set (Wl)HE=[BlWl], Wl is the weights matrix between the (l-1)-th hidden layer and the l-th hidden layer, and Bl is the bias vector of the l-th hidden layer.
(1) If(HlE)THlE is nonsingular, the forced positive-definite Cholesky factorization factor QE of matrix E is calculated according to equations (49a) and (49b), and (VE)(k)(k=1,2,⋯,L) is calculated according to equation (51) using QE and ((g-1(Hl∗))T)(k)(k=1,2,⋯,L). Then the matrices QE and (VE)(k)(k=1,2,⋯,L)are used to calculate (UE)(k)(k=1,2,⋯,L) according to formula (52); therefore we can obtain (Wl)HE=((UE)(1),(UE)(2),⋯,(UE)(L))T(HlE)T.
(2) IfHlE(HlE)T is nonsingular, the forced positive-definite Cholesky factorization factor QF of matrix F is calculated according to equations (53a) and (53b), and (VF)(k)(k=1,2,⋯,L) is calculated according to equation (55) using QF and (bF)(k)(k=1,2,⋯,L). Then the matrices QF and (VF)(k)(k=1,2,⋯,L) are used to calculate (((Wl)HE)T)(k) according to formula (56); thereby we can get (Wl)HE=((((Wl)HE)T)(1),(((Wl)HE)T)(2),⋯,(((Wl)HE)T)(L))T.
S tep 7. Calculate the predicted output matrix of the l-th hidden layer according to the equation Hl=g(WlHl-1+Bl)=g((Wl)HEHlE).
S tep 8. Recalculate the connection weight matrix βl=(λ-1I+(Hl)THl)-1(Hl)TT or βl=(Hl)T(λ-1I+Hl(Hl)T)-1T between the l-th hidden layer and the output layer according to the principle shown in (28)-(36).
S tep 9. On the basis of βl, the MRELM prediction model with l hidden layers is established, and the final output result f(x)=Hlβl is calculated.
S tep 10. Calculate the sum of the empirical risk and the structural risk of the MRELM prediction model. (57)Rl=12βlTβl+λ2εTε.Then jump to Step 5. Judging the formula (58) from l=5,(58)Rl-i-Rl-i-1Rmax≤ξ,where ξ is the learning precision and Rmax is the maximum value of R1,⋯Rl, i=0,⋯,3. If the formula (58) is satisfied, the training process terminates, determine l as the optimal number of hidden layers, and establish the corresponding MRELM prediction model; otherwise continue to increase l until the condition is met.The number of hidden layers in MRELM increases successively from the initial value, and the expansion stops whenRl is no longer significantly reduced. At this time, even if the hidden layer is continued to be added, the Rl representing the learning accuracy and generalization ability of MRELM model will not be significantly improved but will lead to a large number of redundant hidden layers in MRELM. Therefore, the MRELM model at this time has the optimal number of hidden layers.


### 3.4. Proof of Positive Definiteness for FC-IMRELM

In the implementation of FC-IMRELM algorithm, the process of solving the connection weights matrixβl and the learning parameters (Wl)HE in the l-th hidden layer can be transformed into solving linear equations in the form of equations (25), (26), (37), (38), (47), and (48). The premise of applying standard Cholesky factorization to solving linear equations is that its coefficient matrix must be a symmetric positive definite matrix, and thus we need to prove that the matrices A=λ-1I+(H1)TH1 and B=λ-1I+H1(H1)T are symmetric positive definite matrices. The precondition of applying the forced positive definite Cholesky factorization to solving linear equations is that the coefficient matrix must be a symmetric matrix; therefore it is necessary to prove that the matrices C=(βl)Tβl, D=βl(βl)T, E=(HlE)THlE, and F=HlE(HlE)T are symmetric semipositive definite matrices. The following theorem shows that the matrices (Hl)THl and Hl(Hl)T are symmetric semipositive definite matrices, and λ-1I is symmetric positive definite matrix; hence the matrices λ-1I+(Hl)THl and λ-1I+Hl(Hl)T are symmetric positive definite matrices, too. Consequently, Eq. (25) and Eq. (26) can be solved by the standard Cholesky factorization, and the connection weights matrix βl in the l-th hidden layer can be gotten. In addition, the matrices (βl)Tβl, βl(βl)T, (HlE)THlE, and HlE(HlE)T are symmetric semipositive definite matrices, so Eq. (37), Eq. (38), Eq. (47), and Eq. (48) can be calculated by force positive-definite Cholesky factorization, and the learning parameters (Wl)HE in the l-th hidden layer can be obtained.Theorem 3.
LetHl be a matrix. Then λ-1I+(Hl)THl and λ-1I+Hl(Hl)T are all symmetric positive definite matrices.Proof.
Obviously, for the matrixA=λ-1I+(Hl)THl, then by properties of transpose we have (59)AT=λ-1I+HlTHlT=λ-1I+HlTHlTT=λ-1I+HlTHl=A.So, A is symmetric.
Similarly, for the matrixB=λ-1I+Hl(Hl)T, applying properties of transpose again, we get (60)BT=λ-1I+HlHlTT=λ-1I+HlTTHlT=λ-1I+HlHlT=B.Hence, B is symmetric.
On the other hand, letη=(η1,η2,⋯,ηL)T be a nonzero vector in RL. We can easily verify that (61)ηTλ-1Iη=λ-1∑i=1Lηi2>0.In addition,(62)ηTHlTHlη=HlηTHlη=Hlη2≥0.Thus, the matrix A=λ-1I+(Hl)THl is positive definite
Also, letμ=(μ1,μ2,⋯,μN)T be a nonzero vector in RN. Then we deduce that (63)μTλ-1Iμ=λ-1∑i=1Lμi2>0.Additionally,(64)μTHlHlTμ=HlTμTHlTμ=HlTη2≥0.Consequently, the matrix B=λ-1I+Hl(Hl)T is positive definite.Theorem 4.
Letβl be a matrix. Then (βl)Tβl and βl(βl)T are all symmetric positive semidefinite matrices.Proof.
First, for the matrixC=(βl)Tβl, it is not hard to see from properties of transpose that (65)CT=βlTβlT=βlTβlTT=βlTβl=C.Thus, C is symmetric.
Second, by the same method we have forD=βl(βl)T(66)DT=βlβlTT=βlTTβlT=βlβlT=D.Therefore, D is symmetric.
Finally, it is clear that(βl)Tβl is a square symmetric matrix. For any nonzero vector σ=(σ1,σ2,⋯σm)T in Rm, we have(67)σTβlTβlσ=βlσTβlσ=βlσ2≥0.Thus, the matrix C=(βl)Tβl is positive semidefinite.
It is also clear thatβl(βl)T is a square symmetric matrix. For any nonzero vector τ=(τ1,τ2,⋯,τL)T in RL, we get(68)τTβlβlTτ=βlTτTβlTτ=βlTτ2≥0.Consequently, the matrix D=βl(βl)T is positive semidefinite.Theorem 5.
LetHlE be a matrix. Then (HlE)THlE and HlE(HlE)T are all symmetric positive semidefinite matrices.Proof (apparent).
We can prove Theorem5 using the same idea as in Theorem 4.


## 3.1. Forced Positive-Definite Cholesky Factorization (FC)

The main difficulty of the MRELM algorithm is the calculation of the inverse matrix and the MP generalized inverse matrix involved in the training process, including the inverse calculation of symmetric positive semidefinite matrix. In this case, the improved MRELM mode based on the traditional Cholesky factorization could not be realized, because the Cholesky factorization of the symmetric positive semidefinite matrix might not exist. Even if such a factorization exists, the calculation process is generally numerically unstable for the elements of the matrix factorization factor may be unbounded. In order to overcome these difficulties, we put forward a modified approach based on the forced positive definite Cholesky factorization for the MRELM algorithm with multiple hidden layers, which is a numerical stability approach.When the forced positive definite strategy is adopted to improve the MRELM algorithm, the key problem is how to form the positive definite matrix from the modified Cholesky decomposition of the undetermined matrix. If the matrixG is not a positive definite matrix, the Cholesky factorization method, which forces the matrix to have positive definite property, is to find a unit lower triangular matrix L and a positive definite diagonal matrix D for the general symmetric matrix G, so that the matrix G-=LDLT is positively definite, and it is only one diagonal matrix E away from the matrix G.(17)G-=LDLT=G+EIn fact, the Cholesky factorization of symmetric positive definite matrix can be described as follows:(18a)djj=gjj-∑s=1j-1dssljs2(18b)lij=1djjgij-∑s=1j-1dssljslis,i≥j-1where gij represents the element of matrix G and djj denotes the main diagonal element of matrix D. Here, the Cholesky factorization factors L and D are required to satisfy two requirements: one is that all elements of D are strictly positive, and the other is that the elements of the factorization factor L are uniformly bounded. That is, for k=1,⋯,n and a positive number ρ, the formula (19) is required:(19)dkk>δ,rik≤ρ,i>kwhere the auxiliary quantity rik=likdkk, δ is a given small positive number. The matrix satisfying the above conditions is said to be sufficiently positive definite, where E is a zero matrix.Next, we describe thej-th step of this factorization. Suppose the j-1 column of the forced positive-definite Cholesky factorization has been calculated. For k=1,⋯,j-1, equation (19) holds. First calculate(20)γj=ξj-∑s=1j-1dssljs2,where ξj is taken as gjj, and the test value d- is defined as(21)d-=max⁡γj,δwhere δ is a small positive number. In order to determine whether d- can accept as the j-th element of D, we check whether rij=lijd- satisfies the formula (19). If so, let djj=d-, and get the j-th column of L from lij=rij/djj. Otherwise,(22)djj=ξj-∑s=1j-1dssljs2let ξj=gjj+ejj, select positive number ejj to make max⁡rij=ρ, and produce the j-th column of L.If the above process is completed, we obtain the Cholesky factorization formula (17) of the positive definite matrix G-, where E is a nonnegative diagonal matrix and the diagonal element is ejj. For the given matrix G, this nonnegative diagonal matrix E depends on ρ. If n>1, (23)Eρ∞≤ξρ+n-1ρ2+2γ+n-1ρ2δwhere ξ is the maximum norm of the nondiagonal elements of G and γ is the maximum norm of the diagonal elements of G. If ρ2=ξ/n2-1, the upper bound is minimized. So, let ρ satisfy formula (24):(24)ρ2=max⁡γ,ξn2-1,εMwhere εM represents the machine precision. We increase εM to prevent G from being small.Finally, we present the forced positive-definite Cholesky factorization algorithm, where the auxiliary quantitycis=lisdss, s=1,⋯,j, i=j,⋯,n. These values need not be stored separately; they can be stored in the matrix G.Algorithm 1.
Forced Positive-Definite Cholesky Factorization (FC)
S tep 1. Calculate the bounds of the elements of the factorization factor. Let ρ2=max⁡γ,ξ/v,εM, where v=max⁡1,n2-1 and γ and ξ are the maximum norm of the diagonal and nondiagonal elements of G, respectively.
S tep 2. Initialization. Let j=1, cii=gii, i=1,⋯,n.
S tep 3. Determine the minimum index q, so that cqq=maxj≤i≤n⁡cii, and exchange the information of q rows and j rows, q columns and j columns of G.
S tep4. Calculate the j-th row of L, and solve the maximum norm of lijdjj. Let ljs=cjs/dss, s=1,⋯,j-1, calculate cij=gij-∑s=1j-1ljscis, i=j+1,⋯,n, and let θj=maxj+1≤i≤n⁡cij. If j=n, let θj=0.
S tep 5. Calculate the j-th diagonal element djj=max⁡δ,cjj,θj2/ρ2 of D. The diagonal element of E is modified to Ej=djj-cjj. If j=n, stop.
S tep 6. Correct the diagonal elements and column index, and let cii=cii-cij2/djj, i=j+1,⋯,n and j=j+1; jump to Step 3.


## 3.2. Process of Matrices Decomposition for FC-IMRELM

According to the MRELM training process shown in equations (1) to (16), its essence is to solve the connection weight matrix βl between the hidden layer and the output layer. However, it can be seen from equations (5a), (5b), (10a), (10b), (15a), and (15b) that the solution method given in literature [18] involves matrix inversion, and the solution process of each hidden layer's parameters (Wl)HE (l is the number of hidden layers) involves the calculation of the MP generalized inverse of matrix. The problem of large amount of calculation reduces the modeling efficiency of the MRELM prediction model. In order to solve the above problems effectively, we propose a solution method of the weights matrix βl and hidden layer parameters (Wl)HE based on the forced positive definite Cholesky factorization.First, on the basis of equations (4a), (4b), and (4c), β1=(H1)TαT can be obtained from equation (4a), and ε=-αT/λ can be obtained from equation (4b), and then substituting β1 and ε into equation (4c), we can obtain(25)Aβ1k=bAk,k=1,2,⋯,m.where A=λ-1I+(H1)TH1, ((bA)(1),(bA)(2),⋯,(bA)(m))=(H1)TT, and β1=((β1)(1),(β1)(2),⋯,(β1)(m)), and(26)BUBk=Tk,k=1,2,⋯,mwhere B=λ-1I+H1(H1)T, T=(T(1),T(2),⋯,T(m)), and β1=(H1)T((UB)(1),(UB)(2),⋯,(UB)(m)).Therefore, the process of solvingβ1 according to equation (5a) can be transformed into solving m linear equations in the form of equation (25), and the process of solving β1 according to equation (5b) can be transformed into solving m linear equations in the form of equation (26); m is the dimension of observation vector.If the number of training samples is greater than the number of hidden nodes, that is,N>L, the solution process of β1 based on the Cholesky factorization is as follows. First, calculate the Cholesky factorization of matrix A:(27)A=SASATwhere SA is a lower triangular matrix with positive diagonal elements. The nonzero element (sA)ij in SA can be calculated by the element aij of A according to equation (28):(28)sAij=aii-∑n=1i-1sAin2i=j,aij-∑n=1j-1sAinsAjnsAjji>jwhere i=1,2,⋯,L, j=1,2,⋯,L. Substitute equation (27) into equation (25) and multiply both sides of the equation by (SA)-1: (29)SATβ1k=VAk,k=1,2,⋯,m,where (VA)(k)=(SA)-1(bA)(k). Solving β1=((β1)(1),(β1)(2),⋯,(β1)(m)) is equivalent to solving equation (29). Because (VA)(k)=(SA)-1(bA)(k) is equivalent to SA(VA)(k)=(bA)(k), the calculation formula for the elements (vA)i(k) of (VA)(k) can be obtained by comparing the elements on both sides of the equation (30)vAik=bAiksAiii=1,bAik-∑n=1i-1sAnivAnksAiii>1where i=1,2,⋯,L, (bA)i(k) is the element of (bA)(k). Finally, on the basis of obtaining SA and (VA)(k), the element (β1)i(k) of (β1)(k) can be calculated by using the elements of SA and (VA)(k): (31)β1ik=vAiksAiii=L,vAik-∑n=1L-isAi+n,iβ1i+nksAiii<Lwhere i=1,2,⋯,L, k=1,2,⋯,m, β1=((β1)(1),(β1)(2),⋯,(β1)(m)).If the number of hidden nodes is greater than the number of training samples, that is,N<L, the solution process of β1 based on the Cholesky factorization is as follows. First, calculate the Cholesky factorization of B:(32)B=SBSBTwhere SB is a lower triangular matrix with positive diagonal elements. The nonzero element (sB)ij in SB can be calculated by the element bij of B according to equation (33): (33)sBij=bii-∑n=1i-1sBin2i=j,bij-∑n=1j-1sBinsBjnsBjji>jwhere i=1,2,⋯,N, j=1,2,⋯,N. Substitute equation (32) into equation (26) and multiply both sides of the equation by (SB)-1: (34)SBTUBk=VBk,k=1,2,⋯,m,where (VB)(k)=(SB)-1T(k). Solving β1=(H1)T((UB)(1),(UB)(2),⋯,(UB)(m)) is equivalent to solving equation (34). Because (VB)(k)=(SB)-1T(k) is equivalent to SB(VB)(k)=T(k), the calculation formula for the element (vB)i(k) of (VB)(k) can be obtained by comparing the elements on both sides of the equation (35)vBik=TiksBiii=1,Tik-∑n=1i-1sBnivBnksBiii>1where i=1,2,⋯,N, Ti(k) is the element of T(k). Finally, on the basis of obtaining SB and (VB)(k), the element (uB)i(k) of (UB)(k) can be calculated by using the elements of SB and (VB)(k). (36)uBik=vBiksBiii=L,vBik-∑n=1L-isBi+n,iuBi+nksBiii<Lwhere i=1,2,⋯,N, k=1,2,⋯,m, β1=(H1)T((UB)(1),(UB)(2),⋯,(UB)(m)).At this point, we get the connection weight matrixβ1 between the first hidden layer and the output layer. The connection weight matrix βl between the other hidden layer and output layer can also be calculated by the above method. Compared with the solution method of the connecting weight matrix βl as shown in equations (5a), (5b), (10a), (10b), (15a), and (15b), the solution of βl based on Cholesky factorization does not involve the inverse operation of the matrix, and it can be achieved by using simple algebraic operations.The MRELM model contains multiple hidden layers, and the solving process of each hidden layer parameter(Wl)HE needs to calculate the MP generalized inverse of the corresponding matrices. However, it can be concluded from equations (7), (8), (11), and (13) that the solution method of (βl)+ and (HlE)+ using orthogonal projection method [5] involves the inverse calculation of symmetric semipositive definite matrixes (βl)Tβl, βl(βl)T, (HlE)THlE, and HlE(HlE)T, which has the problem of large computational cost and numerical instability, and if the condition number of the above matrices is too large, the calculation results of MP generalized inverse of matrices (βl)+ and (HlE)+ are usually unable to be obtained. This not only affects the modeling efficiency and prediction effect of MRELM model, but also may make the modeling process impossible to complete. However, the traditional Cholesky factorization method can only be used to solve the calculation of symmetric positive definite matrix. To effectively overcome the above difficulties, we use the forced positive definite Cholesky factorization to solve the MP generalized inverse of matrices (βl)+ and (HlE)+, and we then get the hidden layer parameters (Wl)HE.If(βl)Tβl is nonsingular, then take (βl)+=((βl)Tβl)-1(βl)T and substitute it into (7) and (11): (37)Cβl+k=βlTk,k=1,2,⋯,L,where C=(βl)Tβl, H(l+1)∗=T(((βl)+)(1),((βl)+)(2),⋯,((βl)+)(L)).Ifβl(βl)T is nonsingular, then take (βl)+=(βl)T(βl(βl)T)-1 and substitute it into (7): (38)DTβl+Tk=βlk,k=1,2,⋯,m,where D=βl(βl)T, H(l+1)∗=T((((βl)+)T)(1),(((βl)+)T)(2),⋯,(((βl)+)T)(m))T.Therefore, the process of solvingH(l+1)∗ can be transformed into solving L linear equations in the form of equation (37), or solving m linear equations in the form of equation (38), L is the number of hidden nodes, and m is the dimension of the observation vector.If(βl)Tβl is nonsingular, the solution process of H(l+1)∗ based on forced positive definite Cholesky factorization is as follows. First, calculate the modified Cholesky factorization result of matrix C:(39a)C=PCMCPCTwhere PC is the unit lower triangular matrix and MC is the positive definite diagonal matrix. (pC)ij is the nonzero element below the main diagonal in PC, and (mC)ii is the main diagonal element in MC; it can be calculated according to Algorithm 1 by using element cij of C, where i=1,2,⋯,m, j=1,2,⋯,m. If we set MC=(MC)1/2(MC)1/2, then (39a) can be written as(39b)C=QCQCTwhere QC=PC(MC)1/2. Substitute equation (39b) into equation (37) and multiply both sides of the equation by (QC)-1. (40)QCTβl+k=VCk,k=1,2,⋯,L,where (VC)(k)=(QC)-1((βl)T)(k). Solving H(l+1)∗=T(((βl)+)(1),((βl)+)(2),⋯,((βl)+)(L)) is equivalent to solving equation (40). Since (VC)(k)=(QC)-1((βl)T)(k) is equivalent to QC(VC)(k)=((βl)T)(k), by comparing the elements on both sides of the equation, the calculation formula for the element (vC)i(k) of (VC)(k) can be obtained: (41)vCik=βlTikqCiii=1,βlTik-∑n=1i-1qCnivCnkqCiii>1where ((βl)T)i(k) is the element of ((βl)T)(k), and (qC)ij is the element of QC, i=1,2,⋯,m, j=1,2,⋯,m. Finally, on the basis of obtaining QC and (VC)(k), the element ((βl)+)i(k) of ((βl)+)(k) can be calculated by using the elements of QC and (VC)(k). (42)βl+ik=vCikqCiii=L,vCik-∑n=1L-iqCi+n,iβl+i+nkqCiii<Lwhere i=1,2,⋯,m, k=1,2,⋯,L, (βl)+=(((βl)+)(1),((βl)+)(2),⋯,((βl)+)(L)); consequently, we can get H(l+1)∗=T(((βl)+)(1),((βl)+)(2),⋯,((βl)+)(L)).Ifβl(βl)T is nonsingular, the solution process of H(l+1)∗ based on forced positive definite Cholesky factorization is as follows. First, calculate the modified Cholesky factorization results of matrix D: (43a)D=PDMDPDTwhere PD is the unit lower triangular matrix and MD is the positive definite diagonal matrix. (pD)ij is the nonzero element below the main diagonal in PD, and (mD)ii is the main diagonal element in MD; it can be calculated according to Algorithm 1 by using the element dij of D, i=1,2,⋯,L, j=1,2,⋯,L. If we set MD=(MD)1/2(MD)1/2, then (43a) can be written as(43b)D=QDQDTwhere QD=PD(MD)1/2. Substitute equation (43b) into equation (38) and multiply both sides of the equation by (QD)-1. (44)QDTβl+Tk=VDk,k=1,2,⋯,m,where (VD)(k)=(QD)-1(βl)(k). Solving H(l+1)∗ = T((((βl)+)T)(1),(((βl)+)T)(2),⋯,(((βl)+)T)(m))T is equivalent to solving equation (44). Since (VD)(k)=(QD)-1(βl)(k) is equivalent to QD(VD)(k)=(βl)(k), by comparing the elements on both sides of the equation, the calculation formula for the element (vD)i(k) of (VD)(k) can be gotten as (45)
v
D
i
k
=
β
l
i
k
q
D
i
i
i
=
1
,
β
l
i
k
-
∑
n
=
1
i
-
1
q
D
n
i
v
D
n
k
q
D
i
i
i
>
1where (βl)i(k) is the element of (βl)(k) and (qD)ij is the element of QD, i=1,2,⋯,L, j=1,2,⋯,L. Finally, on the basis of obtaining QD and (VD)(k), the element (((βl)+)T)i(k) of (((βl)+)T)(k) can be obtained by using elements of QD and (VD)(k). (46)βl+Tik=vDikqDiii=L,vDik-∑n=1L-iqDi+n,iβl+Ti+nkqDiii<Lwhere i=1,2,⋯,L, k=1,2,⋯,m, (βl)+=((((βl)+)T)(1),(((βl)+)T)(2),⋯,(((βl)+)T)(m))T; therefore, we can get H(l+1)∗=T((((βl)+)T)(1),(((βl)+)T)(2),⋯,(((βl)+)T)(m))T.If(HlE)THlE is nonsingular, then let (HlE)+=((HlE)THlE)-1(HlE)T and substitute it into (8) and (13): (47)EUEk=g-1Hl∗Tk,k=1,2,⋯,L,where E=(HlE)THlE, (Wl)HE=((UE)(1),(UE)(2),⋯,(UE)(L))T(HlE)T.IfHlE(HlE)T is nonsingular, then let (HlE)+=(HlE)T(HlE(HlE)T)-1 and substitute it into (8) and (13): (48)FWlHETk=bFk,k=1,2,⋯,L,where (Wl)HE = ((((Wl)HE)T)(1),(((Wl)HE)T)(2), ⋯, (((Wl)HE)T)(L))T, ((bF)(1),(bF)(2),⋯(bF)(L))T = g-1(Hl∗)(HlE)T, F=HlE(HlE)T.Therefore, the process of solving(Wl)HE can be transformed into solving L linear equations in the form of equation (47), or solving L linear equations in the form of equation (48), where L is the number of hidden nodes.If(HlE)THlE is nonsingular, the solution process of (Wl)HE based on the forced positive definite Cholesky factorization is as follows. First, calculate modified Cholesky factorization results of matrix E: (49a)E=PEMEPETwhere PE is the unit lower triangular matrix and ME is the positive definite diagonal matrix. (pE)ij is the nonzero element below the main diagonal in PE, and (mE)ii is the main diagonal element in ME; it can be calculated according to Algorithm 1 by using element Eij of E, i=1,2,⋯,N, j=1,2,⋯,N. If we set ME=(ME)1/2(ME)1/2, then (49a) can be written as(49b)E=QEQETwhere QE=PE(ME)1/2. Substitute equation (49b) into equation (47) and multiply both sides of the equation by (QE)-1. (50)QETUEk=VEk,k=1,2,⋯,L,where (VE)(k)=(QE)-1((g-1(Hl∗))T)(k). Solving (Wl)HE=((UE)(1),(UE)(2),⋯,(UE)(L))T(HlE)T is equivalent to solving equation (46). For (VE)(k)=(QE)-1((g-1(Hl∗))T)(k) is equivalent to QE(VE)(k)=((g-1(Hl∗))T)(k), by comparing the elements on both sides of the equation, the calculation formula for the element (vE)i(k) of (VE)(k) can be calculated as (51)vEik=g-1Hl∗TikqEiii=1,g-1Hl∗Tik-∑n=1i-1qEnivEnkqEiii>1where ((g-1(Hl∗))T)i(k) is the element of ((g-1(Hl∗))T)(k) and (qE)ij is the element of QE, i=1,2,⋯,N, j=1,2,⋯,N. Finally, on the basis of obtaining QE and (VE)(k), the element (UE)i(k) of (UE)(k) can be calculated by using the elements of QE and (VE)(k). (52)UEik=vEikqEiii=L,vEik-∑n=1L-iqEi+n,iUEi+nkqEiii<Lwhere i=1,2,⋯,N, k=1,2,⋯,L, UE=((UE)(1),(UE)(2),⋯,(UE)(L)); accordingly, we can obtain the matrix (Wl)HE=((UE)(1),(UE)(2),⋯,(UE)(L))T(HlE)T.IfHlE(HlE)T is nonsingular, the solution process of (Wl)HE based on the forced positive definite Cholesky factorization is as follows. First, calculate modified Cholesky factorization results of matrix F:(53a)F=PFMFPFTwhere PF is the unit lower triangular matrix and MF is the positive definite diagonal matrix. (pF)ij is the nonzero element below the main diagonal in PF, and (mF)ii is the main diagonal element in MF; it can be calculated according to Algorithm 1 by using element Fij of F, i=1,2,⋯,L+1, j=1,2,⋯,L+1. If we set MF=(MF)1/2(MF)1/2, then (53a) can be written as(53b)F=QFQFTwhere QF=PF(MF)1/2. Substitute equation (53b) into equation (48) and multiply both sides of the equation by (QF)-1. (54)QFTWlHETk=VFk,k=1,2,⋯,L,where (VF)(k)=(QF)-1(bF)(k). Because solving (Wl)HE=((((Wl)HE)T)(1),(((Wl)HE)T)(2),⋯,(((Wl)HE)T)(L))T is equivalent to solving equation (54) and (VF)(k)=(QF)-1(bF)(k) is equivalent to QF(VF)(k)=(bF)(k), by comparing the elements on both sides of the equation, the calculation formula for the element (vF)i(k) of (VF)(k) can be calculated as (55)vFik=bFikqFiii=1,bFik-∑n=1i-1qFnivFnkqFiii>1where (bF)i(k) is the element of (bF)(k) and (qF)ij is the element of QF, i=1,2,⋯,L+1, j=1,2,⋯,L+1. Finally, on the basis of obtaining QF and (VF)(k), the element (((Wl)HE)T)i(k) of (((Wl)HE)T)(k) can be calculated by using the elements of QF and (VF)(k). (56)WlHETik=vFikqFiii=L,vFik-∑n=1L-iqFi+n,iWlHETi+nkqFiii<Lwhere i=1,2,⋯,L+1, k=1,2,⋯,L. So we can get (Wl)HE=((((Wl)HE)T)(1),(((Wl)HE)T)(2),⋯,(((Wl)HE)T)(L))T.So far we have obtained the parameters(Wl)HE of each hidden layer. Compared with the solving method of hidden layer parameters (Wl)HE shown in equations (7), (8), (11), and (13), the calculation approach of (Wl)HE based on the forced positive definite Cholesky factorization does not involve the operation of inverse matrix, which can be realized only by simple algebraic calculations. Meanwhile, the method of forcing the matrix to be positive definite matrix can guarantee the numerical stability of MRELM neural network training process.


## 3.3. FC-IMRELM Training Algorithm

Studies have shown that the number of hidden layers determines the learning accuracy and generalization ability of the MRELM model [25], and it is also a key factor that must be determined in advance when designing the MRELM network structure. Due to the complexity of various training samples applied in MRELM prediction model, it is difficult to accurately determine the optimal number of hidden layers by human experience, so that the MRELM prediction model has enough hidden layers to ensure its learning accuracy, while, at the same time, it has as few hidden layers as possible to maintain its contracted network structure. To avoid the disadvantages and difficulties of artificially selecting the number of hidden layers, we propose an incremental MRELM training algorithm based on forced positive Cholesky factorization, which can automatically determine the optimal number of hidden layers in MRELM, and the training process is as follows.Algorithm 2.
Incremental MRELM based on Forced Positive-Definite Cholesky Factorization (FC-IMRELM)
S tep 1. Suppose the training sample dataset is X,T=xi,ti(i=1,2,⋯,N), where X=[x1,x2,⋯,xN]T is the input samples, T=[t1,t2,⋯,tN]T is the labeled samples, each hidden layer contains the number of hidden nodes L, and the activation function is g(x).
S tep 2. Let the number of hidden layers in MRELM l=1, randomly initialize the input weights matrix W1 between the input layer and the first hidden layer and the bias vector B1 of the first hidden layer, and let (W1)HE=[B1W1], H1E=[1X]T.
S tep 3. Calculate the output matrix of the first hidden layer H1=g((W1)HEH1E), and calculate the matrices A, B, and (bA)(k)(k=1,2,⋯,m).
S tep 4. Calculate the connection weights matrix between the first hidden layer and the output layer β1=(λ-1I+(H1)TH1)-1(H1)TT or β1=(H1)T(λ-1I+H1(H1)T)-1T.
(1) IfN>L, the Cholesky factorization factor SA of matrix A is calculated according to the formula (28), and (VA)(k)(k=1,2,⋯,m) is calculated according to the formula (30) using matrices SA and (bA)(k)(k=1,2,⋯,m). Then matrices SA and (VA)(k)(k=1,2,⋯,m) are used to calculate (β1)(k) according to equation (31), and we can get β1=((β1)(1),(β1)(2),⋯,(β1)(m)).
(2) IfN<L, the Cholesky factorization factor SB of matrix B is calculated according to the formula (33), and (VB)(k)(k=1,2,⋯,m) is calculated according to the formula (35) using SB and T(k)(k=1,2,⋯,m). Then the matrices SB and (VB)(k)(k=1,2,⋯,m) are used to calculate (UB)(k)(k=1,2,⋯,m) according to the equation (36); subsequently, we can obtain β1=(H1)T((UB)(1),(UB)(2),⋯,(UB)(m)).
S tep 5. Calculate the expected output H(l+1)∗=T(βl)+ of the (l+1)-th hidden layer.
(1) If(βl)Tβl is nonsingular, the forced positive-definite Cholesky factorization factor QC of matrix C is calculated according to equations (39a) and (39b), and (VC)(k)(k=1,2,⋯,L) is calculated according to equation (41) using QC and ((βl)T)(k)(k=1,2,⋯,L). Then the matrices QC and (VC)(k)(k=1,2,⋯,L) are used to calculate ((βl)+)(k)(k=1,2,⋯,L) according to formula (42), so we can obtain H(l+1)∗=T(((βl)+)(1),((βl)+)(2),⋯,((βl)+)(L)).
(2) Ifβl(βl)T is nonsingular, the forced positive-definite Cholesky factorization factor QD of D is calculated according to equations (43a) and (43b), and (VD)(k)(k=1,2,⋯,m) is calculated according to equation (45) using QD and (βl)(k)(k=1,2,⋯,m). Then the matrices QD and (VD)(k)(k=1,2,⋯,m) are used to calculate (((βl)+)T)(k)(k=1,2,⋯,m) according to formula (46); thus we can get H(l+1)∗=T((((βl)+)T)(1),(((βl)+)T)(2),⋯,(((βl)+)T)(m))T.
S tep 6. Let l=l+1, calculate the parameter (Wl)HE=g-1(Hl∗)(HlE)+ of the l-th hidden layer, set (Wl)HE=[BlWl], Wl is the weights matrix between the (l-1)-th hidden layer and the l-th hidden layer, and Bl is the bias vector of the l-th hidden layer.
(1) If(HlE)THlE is nonsingular, the forced positive-definite Cholesky factorization factor QE of matrix E is calculated according to equations (49a) and (49b), and (VE)(k)(k=1,2,⋯,L) is calculated according to equation (51) using QE and ((g-1(Hl∗))T)(k)(k=1,2,⋯,L). Then the matrices QE and (VE)(k)(k=1,2,⋯,L)are used to calculate (UE)(k)(k=1,2,⋯,L) according to formula (52); therefore we can obtain (Wl)HE=((UE)(1),(UE)(2),⋯,(UE)(L))T(HlE)T.
(2) IfHlE(HlE)T is nonsingular, the forced positive-definite Cholesky factorization factor QF of matrix F is calculated according to equations (53a) and (53b), and (VF)(k)(k=1,2,⋯,L) is calculated according to equation (55) using QF and (bF)(k)(k=1,2,⋯,L). Then the matrices QF and (VF)(k)(k=1,2,⋯,L) are used to calculate (((Wl)HE)T)(k) according to formula (56); thereby we can get (Wl)HE=((((Wl)HE)T)(1),(((Wl)HE)T)(2),⋯,(((Wl)HE)T)(L))T.
S tep 7. Calculate the predicted output matrix of the l-th hidden layer according to the equation Hl=g(WlHl-1+Bl)=g((Wl)HEHlE).
S tep 8. Recalculate the connection weight matrix βl=(λ-1I+(Hl)THl)-1(Hl)TT or βl=(Hl)T(λ-1I+Hl(Hl)T)-1T between the l-th hidden layer and the output layer according to the principle shown in (28)-(36).
S tep 9. On the basis of βl, the MRELM prediction model with l hidden layers is established, and the final output result f(x)=Hlβl is calculated.
S tep 10. Calculate the sum of the empirical risk and the structural risk of the MRELM prediction model. (57)Rl=12βlTβl+λ2εTε.Then jump to Step 5. Judging the formula (58) from l=5,(58)Rl-i-Rl-i-1Rmax≤ξ,where ξ is the learning precision and Rmax is the maximum value of R1,⋯Rl, i=0,⋯,3. If the formula (58) is satisfied, the training process terminates, determine l as the optimal number of hidden layers, and establish the corresponding MRELM prediction model; otherwise continue to increase l until the condition is met.The number of hidden layers in MRELM increases successively from the initial value, and the expansion stops whenRl is no longer significantly reduced. At this time, even if the hidden layer is continued to be added, the Rl representing the learning accuracy and generalization ability of MRELM model will not be significantly improved but will lead to a large number of redundant hidden layers in MRELM. Therefore, the MRELM model at this time has the optimal number of hidden layers.


## 3.4. Proof of Positive Definiteness for FC-IMRELM

In the implementation of FC-IMRELM algorithm, the process of solving the connection weights matrixβl and the learning parameters (Wl)HE in the l-th hidden layer can be transformed into solving linear equations in the form of equations (25), (26), (37), (38), (47), and (48). The premise of applying standard Cholesky factorization to solving linear equations is that its coefficient matrix must be a symmetric positive definite matrix, and thus we need to prove that the matrices A=λ-1I+(H1)TH1 and B=λ-1I+H1(H1)T are symmetric positive definite matrices. The precondition of applying the forced positive definite Cholesky factorization to solving linear equations is that the coefficient matrix must be a symmetric matrix; therefore it is necessary to prove that the matrices C=(βl)Tβl, D=βl(βl)T, E=(HlE)THlE, and F=HlE(HlE)T are symmetric semipositive definite matrices. The following theorem shows that the matrices (Hl)THl and Hl(Hl)T are symmetric semipositive definite matrices, and λ-1I is symmetric positive definite matrix; hence the matrices λ-1I+(Hl)THl and λ-1I+Hl(Hl)T are symmetric positive definite matrices, too. Consequently, Eq. (25) and Eq. (26) can be solved by the standard Cholesky factorization, and the connection weights matrix βl in the l-th hidden layer can be gotten. In addition, the matrices (βl)Tβl, βl(βl)T, (HlE)THlE, and HlE(HlE)T are symmetric semipositive definite matrices, so Eq. (37), Eq. (38), Eq. (47), and Eq. (48) can be calculated by force positive-definite Cholesky factorization, and the learning parameters (Wl)HE in the l-th hidden layer can be obtained.Theorem 3.
LetHl be a matrix. Then λ-1I+(Hl)THl and λ-1I+Hl(Hl)T are all symmetric positive definite matrices.Proof.
Obviously, for the matrixA=λ-1I+(Hl)THl, then by properties of transpose we have (59)AT=λ-1I+HlTHlT=λ-1I+HlTHlTT=λ-1I+HlTHl=A.So, A is symmetric.
Similarly, for the matrixB=λ-1I+Hl(Hl)T, applying properties of transpose again, we get (60)BT=λ-1I+HlHlTT=λ-1I+HlTTHlT=λ-1I+HlHlT=B.Hence, B is symmetric.
On the other hand, letη=(η1,η2,⋯,ηL)T be a nonzero vector in RL. We can easily verify that (61)ηTλ-1Iη=λ-1∑i=1Lηi2>0.In addition,(62)ηTHlTHlη=HlηTHlη=Hlη2≥0.Thus, the matrix A=λ-1I+(Hl)THl is positive definite
Also, letμ=(μ1,μ2,⋯,μN)T be a nonzero vector in RN. Then we deduce that (63)μTλ-1Iμ=λ-1∑i=1Lμi2>0.Additionally,(64)μTHlHlTμ=HlTμTHlTμ=HlTη2≥0.Consequently, the matrix B=λ-1I+Hl(Hl)T is positive definite.Theorem 4.
Letβl be a matrix. Then (βl)Tβl and βl(βl)T are all symmetric positive semidefinite matrices.Proof.
First, for the matrixC=(βl)Tβl, it is not hard to see from properties of transpose that (65)CT=βlTβlT=βlTβlTT=βlTβl=C.Thus, C is symmetric.
Second, by the same method we have forD=βl(βl)T(66)DT=βlβlTT=βlTTβlT=βlβlT=D.Therefore, D is symmetric.
Finally, it is clear that(βl)Tβl is a square symmetric matrix. For any nonzero vector σ=(σ1,σ2,⋯σm)T in Rm, we have(67)σTβlTβlσ=βlσTβlσ=βlσ2≥0.Thus, the matrix C=(βl)Tβl is positive semidefinite.
It is also clear thatβl(βl)T is a square symmetric matrix. For any nonzero vector τ=(τ1,τ2,⋯,τL)T in RL, we get(68)τTβlβlTτ=βlTτTβlTτ=βlTτ2≥0.Consequently, the matrix D=βl(βl)T is positive semidefinite.Theorem 5.
LetHlE be a matrix. Then (HlE)THlE and HlE(HlE)T are all symmetric positive semidefinite matrices.Proof (apparent).
We can prove Theorem5 using the same idea as in Theorem 4.


## 4. Results and Discussion

In this section, experiments of our proposed FC-IMRELM are conducted on benchmark data sets for classification and regression problems. In order to investigate the improvement of learning accuracy of our methods, original ELM [5], TELM [17], and MRELM [18] are also evaluated. All the performance assessments are conducted in the MATLAB 2014b computational environment running on Windows 10 operating system with Intel Core™ i5-7200U CPU @2.7 GHz, 8 GB RAM, NVIDIA M150 Graphics card, and 2 GB GDDR5 video memory. Furthermore, to comprehensively compare resulting performances, the activation function of each algorithm used in the experiments is uniformly assigned as sigmoid function: g(x)=1/(1+exp⁡(-λx)), where λ is 1. The number of hidden neurons is set to 20. Moreover, each algorithm runs for 100 trials and the results obtained will be averaged as the final value.


### 4.1. Classification Problems

#### 4.1.1. Characteristics of Classification Datasets

The classification datasets are obtained from the UCI website [26] and the literature [27, 28]. To evaluate the robustness of our FC-IMRELM algorithm, we conducted the tests using simple benchmark datasets and real datasets collected from coal and iron ores industries. The characteristics of the datasets are shown in Table 1.Table 1
Properties of the classification datasets.
Data set Attributes Classes Training Testing Banknote 4 2 857 515 Blood 4 2 500 258 Diabetic 19 2 719 432 Image 19 7 1443 867 Wilt 5 2 4339 500 Coal spectral data 8 4 196 119 Iron spectral data 5 2 86 37


#### 4.1.2. Evaluation of Testing Accuracy on Classification Datasets

In order to make the performance evaluation more comprehensive, the real datasets that are related to complex industrial data were added to our performance evaluation. The original ELM, TELM, and MRELM algorithms are tested using the simple benchmark datasets and real datasets to validate the improvement of learning accuracy of our IMRELM algorithm. From Table2 and Figure 1, we can see that each algorithm has a good classification accuracy for Banknote dataset. For Blood, Diabetic, Wilt, Coal spectral, and Iron spectral datasets, the algorithms TELM, MRELM, and FC-IMRELM all outperform the ELM algorithm, and the FC-IMRELM algorithm has the highest classification accuracy. For Image dataset, the classification accuracy of the algorithms TELM, MRELM, and FC-IMRELM is much higher than that of the ELM algorithm, and the classification accuracy of the FC-IMRELM algorithm is still the highest, reaching 91.98%. The experimental results show that the average classification accuracy of our FC-IMRELM algorithm is significantly higher than that of the original ELM, TELM, and MRELM algorithms, and the computational experiments using Coal spectral and Iron spectral datasets also demonstrate that our FC-IMRELM algorithm can be easily extended to practical application.Table 2
Average testing classification correct percentage for the algorithms ELM, TELM, MRELM, and FC-IMRELM using classification datasets (%).
Data set ELM TELM MRELM FC-IMRELM Banknote 99.82 99.88 99.93 99.95 Blood 78.95 80.78 80.71 81.12 Diabetic 72.90 73.79 73.92 74.05 Image 86.61 91.09 91.29 91.98 Wilt 62.39 66.35 66.32 67.26 Coal spectral data (Coal) 89.31 91.22 91.02 94.12 Iron spectral data (Iron) 96.26 96.15 97.56 97.83Figure 1
The average testing classification accuracy for the algorithms ELM, TELM, MRELM, and FC-IMRELM using classification datasets.


### 4.2. Regression Problems

For the regression model, the root-mean-square error (RMSE) and the coefficient of determination (R2) [29] are used as the model performance evaluation indexes in this study to verify the effectiveness of the proposed FC-IMRELM algorithm. R2 and RMSE are expressed as follows:(69)RMSE=∑i=1Ntestyi-y^i2Ntest(70)R2=1-∑i=1Ntestyi-y^i2∑i=1Ntestyi-y-2where Ntest is the number of samples in the prediction set; yi is the actual value of the sample; y- is the average of the actual values; y^i is the predicted value calculated by the model. The value range of R2 is between (0, 1), the closer the value of R2 to 1 and the smaller the value of RMSE, the better the performance of the model.


#### 4.2.1. Characteristics of Regression Datasets

The regression datasets are obtained from the LIBSVM website [30]. The characteristics of the datasets are shown in Table 3.Table 3
Specifications for the regression datasets.
Data set Attributes Training Testing Bodyfat 14 212 40 Pyrim 27 64 10 Triazines 60 161 25


#### 4.2.2. Evaluation of Estimation Accuracy on Regression Datasets

As shown in Table4 and Figure 2, we can conclude that the algorithms TELM, MRELM, and FC-IMRELM all have good prediction results for the Bodyfat dataset, where RMSE of these algorithms are small, and R2 are above 0.98, while the prediction ability of the ELM algorithm is slightly worse, and R2 of ELM is equal to 0.8. For Pyrim and Triazines datasets, the prediction results utilizing the FC-IMRELM algorithm are better than those utilizing the algorithms ELM, TELM, and MRELM, among which RMSE of the FC-IMRELM algorithm is the smallest, and R2 is the highest. The above analysis of experimental results indicates that the advantage of the FC-IMRELM algorithm is that it can better extract data characteristics for the multiattribute data, so it has better predictive ability.Table 4
Performance comparison among the original ELM, TELM, MRELM, and FC-IMRELM on regression datasets.
Data set ELM TELM MRELM FC-IMRELM RMSE R2 RMSE R2 RMSE R2 RMSE R2 Bodyfat 0.0084 0.80 0.002 0.986 0.0021 0.984 0.002 0.99 Pyrim 0.12 0.34 0.08 0.47 0.09 0.61 0.06 0.67 Triazines 0.12 0.18 0.12 0.21 0.13 0.18 0.091 0.38Figure 2
The average RMSE (a) and R2 (b) for the algorithms ELM, TELM, MRELM, and FC-IMRELM using regression datasets.
(a) (b)


## 4.1. Classification Problems

### 4.1.1. Characteristics of Classification Datasets

The classification datasets are obtained from the UCI website [26] and the literature [27, 28]. To evaluate the robustness of our FC-IMRELM algorithm, we conducted the tests using simple benchmark datasets and real datasets collected from coal and iron ores industries. The characteristics of the datasets are shown in Table 1.Table 1
Properties of the classification datasets.
Data set Attributes Classes Training Testing Banknote 4 2 857 515 Blood 4 2 500 258 Diabetic 19 2 719 432 Image 19 7 1443 867 Wilt 5 2 4339 500 Coal spectral data 8 4 196 119 Iron spectral data 5 2 86 37


### 4.1.2. Evaluation of Testing Accuracy on Classification Datasets

In order to make the performance evaluation more comprehensive, the real datasets that are related to complex industrial data were added to our performance evaluation. The original ELM, TELM, and MRELM algorithms are tested using the simple benchmark datasets and real datasets to validate the improvement of learning accuracy of our IMRELM algorithm. From Table2 and Figure 1, we can see that each algorithm has a good classification accuracy for Banknote dataset. For Blood, Diabetic, Wilt, Coal spectral, and Iron spectral datasets, the algorithms TELM, MRELM, and FC-IMRELM all outperform the ELM algorithm, and the FC-IMRELM algorithm has the highest classification accuracy. For Image dataset, the classification accuracy of the algorithms TELM, MRELM, and FC-IMRELM is much higher than that of the ELM algorithm, and the classification accuracy of the FC-IMRELM algorithm is still the highest, reaching 91.98%. The experimental results show that the average classification accuracy of our FC-IMRELM algorithm is significantly higher than that of the original ELM, TELM, and MRELM algorithms, and the computational experiments using Coal spectral and Iron spectral datasets also demonstrate that our FC-IMRELM algorithm can be easily extended to practical application.Table 2
Average testing classification correct percentage for the algorithms ELM, TELM, MRELM, and FC-IMRELM using classification datasets (%).
Data set ELM TELM MRELM FC-IMRELM Banknote 99.82 99.88 99.93 99.95 Blood 78.95 80.78 80.71 81.12 Diabetic 72.90 73.79 73.92 74.05 Image 86.61 91.09 91.29 91.98 Wilt 62.39 66.35 66.32 67.26 Coal spectral data (Coal) 89.31 91.22 91.02 94.12 Iron spectral data (Iron) 96.26 96.15 97.56 97.83Figure 1
The average testing classification accuracy for the algorithms ELM, TELM, MRELM, and FC-IMRELM using classification datasets.


## 4.1.1. Characteristics of Classification Datasets

The classification datasets are obtained from the UCI website [26] and the literature [27, 28]. To evaluate the robustness of our FC-IMRELM algorithm, we conducted the tests using simple benchmark datasets and real datasets collected from coal and iron ores industries. The characteristics of the datasets are shown in Table 1.Table 1
Properties of the classification datasets.
Data set Attributes Classes Training Testing Banknote 4 2 857 515 Blood 4 2 500 258 Diabetic 19 2 719 432 Image 19 7 1443 867 Wilt 5 2 4339 500 Coal spectral data 8 4 196 119 Iron spectral data 5 2 86 37


## 4.1.2. Evaluation of Testing Accuracy on Classification Datasets

In order to make the performance evaluation more comprehensive, the real datasets that are related to complex industrial data were added to our performance evaluation. The original ELM, TELM, and MRELM algorithms are tested using the simple benchmark datasets and real datasets to validate the improvement of learning accuracy of our IMRELM algorithm. From Table2 and Figure 1, we can see that each algorithm has a good classification accuracy for Banknote dataset. For Blood, Diabetic, Wilt, Coal spectral, and Iron spectral datasets, the algorithms TELM, MRELM, and FC-IMRELM all outperform the ELM algorithm, and the FC-IMRELM algorithm has the highest classification accuracy. For Image dataset, the classification accuracy of the algorithms TELM, MRELM, and FC-IMRELM is much higher than that of the ELM algorithm, and the classification accuracy of the FC-IMRELM algorithm is still the highest, reaching 91.98%. The experimental results show that the average classification accuracy of our FC-IMRELM algorithm is significantly higher than that of the original ELM, TELM, and MRELM algorithms, and the computational experiments using Coal spectral and Iron spectral datasets also demonstrate that our FC-IMRELM algorithm can be easily extended to practical application.Table 2
Average testing classification correct percentage for the algorithms ELM, TELM, MRELM, and FC-IMRELM using classification datasets (%).
Data set ELM TELM MRELM FC-IMRELM Banknote 99.82 99.88 99.93 99.95 Blood 78.95 80.78 80.71 81.12 Diabetic 72.90 73.79 73.92 74.05 Image 86.61 91.09 91.29 91.98 Wilt 62.39 66.35 66.32 67.26 Coal spectral data (Coal) 89.31 91.22 91.02 94.12 Iron spectral data (Iron) 96.26 96.15 97.56 97.83Figure 1
The average testing classification accuracy for the algorithms ELM, TELM, MRELM, and FC-IMRELM using classification datasets.


## 4.2. Regression Problems

For the regression model, the root-mean-square error (RMSE) and the coefficient of determination (R2) [29] are used as the model performance evaluation indexes in this study to verify the effectiveness of the proposed FC-IMRELM algorithm. R2 and RMSE are expressed as follows:(69)RMSE=∑i=1Ntestyi-y^i2Ntest(70)R2=1-∑i=1Ntestyi-y^i2∑i=1Ntestyi-y-2where Ntest is the number of samples in the prediction set; yi is the actual value of the sample; y- is the average of the actual values; y^i is the predicted value calculated by the model. The value range of R2 is between (0, 1), the closer the value of R2 to 1 and the smaller the value of RMSE, the better the performance of the model.


### 4.2.1. Characteristics of Regression Datasets

The regression datasets are obtained from the LIBSVM website [30]. The characteristics of the datasets are shown in Table 3.Table 3
Specifications for the regression datasets.
Data set Attributes Training Testing Bodyfat 14 212 40 Pyrim 27 64 10 Triazines 60 161 25


### 4.2.2. Evaluation of Estimation Accuracy on Regression Datasets

As shown in Table4 and Figure 2, we can conclude that the algorithms TELM, MRELM, and FC-IMRELM all have good prediction results for the Bodyfat dataset, where RMSE of these algorithms are small, and R2 are above 0.98, while the prediction ability of the ELM algorithm is slightly worse, and R2 of ELM is equal to 0.8. For Pyrim and Triazines datasets, the prediction results utilizing the FC-IMRELM algorithm are better than those utilizing the algorithms ELM, TELM, and MRELM, among which RMSE of the FC-IMRELM algorithm is the smallest, and R2 is the highest. The above analysis of experimental results indicates that the advantage of the FC-IMRELM algorithm is that it can better extract data characteristics for the multiattribute data, so it has better predictive ability.Table 4
Performance comparison among the original ELM, TELM, MRELM, and FC-IMRELM on regression datasets.
Data set ELM TELM MRELM FC-IMRELM RMSE R2 RMSE R2 RMSE R2 RMSE R2 Bodyfat 0.0084 0.80 0.002 0.986 0.0021 0.984 0.002 0.99 Pyrim 0.12 0.34 0.08 0.47 0.09 0.61 0.06 0.67 Triazines 0.12 0.18 0.12 0.21 0.13 0.18 0.091 0.38Figure 2
The average RMSE (a) and R2 (b) for the algorithms ELM, TELM, MRELM, and FC-IMRELM using regression datasets.
(a) (b)


## 4.2.1. Characteristics of Regression Datasets

The regression datasets are obtained from the LIBSVM website [30]. The characteristics of the datasets are shown in Table 3.Table 3
Specifications for the regression datasets.
Data set Attributes Training Testing Bodyfat 14 212 40 Pyrim 27 64 10 Triazines 60 161 25


## 4.2.2. Evaluation of Estimation Accuracy on Regression Datasets

As shown in Table4 and Figure 2, we can conclude that the algorithms TELM, MRELM, and FC-IMRELM all have good prediction results for the Bodyfat dataset, where RMSE of these algorithms are small, and R2 are above 0.98, while the prediction ability of the ELM algorithm is slightly worse, and R2 of ELM is equal to 0.8. For Pyrim and Triazines datasets, the prediction results utilizing the FC-IMRELM algorithm are better than those utilizing the algorithms ELM, TELM, and MRELM, among which RMSE of the FC-IMRELM algorithm is the smallest, and R2 is the highest. The above analysis of experimental results indicates that the advantage of the FC-IMRELM algorithm is that it can better extract data characteristics for the multiattribute data, so it has better predictive ability.Table 4
Performance comparison among the original ELM, TELM, MRELM, and FC-IMRELM on regression datasets.
Data set ELM TELM MRELM FC-IMRELM RMSE R2 RMSE R2 RMSE R2 RMSE R2 Bodyfat 0.0084 0.80 0.002 0.986 0.0021 0.984 0.002 0.99 Pyrim 0.12 0.34 0.08 0.47 0.09 0.61 0.06 0.67 Triazines 0.12 0.18 0.12 0.21 0.13 0.18 0.091 0.38Figure 2
The average RMSE (a) and R2 (b) for the algorithms ELM, TELM, MRELM, and FC-IMRELM using regression datasets.
(a) (b)


## 5. Conclusions

(1) First of all, compared with MRELM, the FC-IMRELM algorithm proposed in this paper uses the idea of forced positive definite Cholesky factorization to determine the hidden layer parameters. The training process is more simplified, with low calculation amount and high numerical stability. In addition, the MRELM algorithm needs to set the network structure in advance, and the number of hidden layers remains unchanged in the training process, while the FC-IMRELM algorithm can automatically select the optimal number of hidden layers through the principle of structural risk and empirical risk minimization and adjust the network structure adaptively according to the training samples.(2) Secondly, compared with CF-FORELM, the FC-IMRELM algorithm proposed in this paper is designed for the semipositive definite matrices appearing in the parameter solving process of MRELM model. The condition number of the matrix is improved while forcing the matrix positive definite, thereby accelerating the convergence speed of the MRELM model and ensuring the numerical stability of the modeling process.(3) Finally, by introducing parameterλ to weigh the structural risk and empirical risk of ELM model, the FC-IMRELM algorithm has significantly improved its generalization ability compared with the traditional neural network. In addition, the forced positive definite Cholesky factorization is used to calculate its output weights, effectively reducing the computational cost brought by the increasing process of hidden layers. The prediction example shows that the FC-IMRELM algorithm can effectively avoid the numerical instability of MRELM model and has the advantages of high prediction accuracy and fast calculation speed, which can provide a novel and efficient solution to the prediction problem.



---

*Source: 6740523-2019-04-24.xml*