This work is motivated by a speciﬁc problem of a modern paper m anufacturing industry, in which boiler inlet water quality is a major conc ern. Boiler helps to produce power and steam which is used to cook wood chips (alon g with chemicals) to be supplied to a paper machine. If water treatment plant ca n’t produce water of desired quality as speciﬁed by the boiler, then it results in poor health of the boiler water tube and consequently aﬀects the quality of the paper. Variation in inlet water quality of the boiler is due to several crucial process parameters. We formulate this problem into a typical regression problem that involve s ﬁnding crucial process parameters from the set of available suspected causal varia bles. Then we build a novel nonparametric hybrid model to solve this regression p roblem and call it as optimal neural regression tree (NRT) model with the aim to so lve the issue of water quality. The major advantage of the proposed optimal NRT mod el is that it has very less tuning parameters and is easily interpretable as compa red to “black-box-like” complex machine learning models. Statistical properties i ncluding optimal values of the important model parameters of the proposed model are pro ved in this paper. We have also experimentally assessed the performance of the proposed model in comparison with the other state-of-the-art models. KEYWORDS Water quality; Demineralization process; Regression Tree ; Neural networks; Optimal neural regression tree 1. Introduction This work is motivated by a particular problem in a modern pap er industry which produces papers for multiple uses. Paper machines produce p apers by using pulp, ﬁbers, ﬁllers, chemical lubricants and a huge amount of wate r. Boiler produces steam for power generation purposes and also helps to make pulp for paper production. The steam produced in the boiler is used to cook wood chips (al ong with the cooking chemicals). Steam is also sent to dryer cans to remove water f rom the sheet produced which the drainage, vacuum, and mechanical pressing sectio ns of the paper machine can’t accomplish. The steam from the boilers is also used thr oughout the mill in heat exchanging, steam-traced piping, stock chests, etc. The bo iler stipulates the desired level of water quality to be received from the water treatmen t plant (WTP). In Corresponding author: Tanujit Chakraborty. Email: tanuji tr@isical.ac.in WTP, the process of demineralization (DM Process) is used fo r removal of dissolved solids by the ion exchange process (IEP). IEP involves two st ages of demineralization (Batchelder 1965). In the ﬁrst stage, it removes cations fro m water by the cation exchange process and then removes anions from water by the an ion exchange process in the second stage (Lhassani et al. 2001). DM process outlet pH happens to be the key performanceindicator (KPI) of WTP. It was foundthat WTP can’t producewater of desired quality speciﬁed by the boiler. By controlling th e water quality, we will be able to address the problem of variation in DM outlet water pH as well as improve boiler water tube health and ﬁnally make improvement in the p aper manufacturing process. An extensive preliminary analysis was conducted t o determine a set of possible causal variables that happen to be the key to water p H level variations. The aim of this paper is to determine the optimal levels of import ant process parameters to maintain water pH level within the desired standard as spe ciﬁed by the boiler. For this, we would develop a regression model which will also be useful for future prediction of water pH based on the important causal paramet ers. Controlling water pH will necessarily improve the water quality of the boiler w hich will be economical for the company as well. The formulation of the above-stated problem along with data collection plan is described in Section 2. Statistical and machine learning tools like factor analysi s, SPC techniques, decision trees, and artiﬁcial neural networks, etc have been applied to solve several problems in analyzing water quality of river (Singh et al. 2009; Ouyan g et al. 2006; Mahuli, Rhinehart, and Riggs 1993) and surface water planning (Gmar et al. 2017; Bhat- tacharya and Solomatine 2005). Multivariate statistical t ools and pattern recognition methods are found to be eﬀective to capture the temporal and sp atial variations in river water quality (Singh et al. 2004; Alberto et al. 2001 ) and manufacturing process eﬃciency (Chakraborty, Chakraborty, and Chattopa dhyay 2018). It was found that decision trees and neural networks (NN) have the abilit y to model arbitrary decision boundaries. Regression tree (RT) is found more rob ust when limited data are available, unlike NN. But decision trees are high variance e stimators and greedy in nature (Breiman 2017), whereas neural networks are popular for tackling prediction problems. Advanced neural networks have many free tuning pa rameters and hence there is a risk of over-ﬁtting in case of small and medium size d data sets. To utilize the positive aspects of these two powerful models, theoretical frameworks for combining both classiﬁers are often used jointly to make decisions. Ma pping tree-based models to neural networks allows exploiting the former to initiali ze the latter. The idea of mapping tree based models into NN are presented by several pa pers to solve many supervised learning problems (Sethi 1990; Sirat and Nadal 1 990). Similar idea in the interface area of decision trees and neural networks can be found in some recent literature as well (Chen et al. 2005; Balestriero 2017). But the major disadvantages of these algorithms are having many free tuning parameters and thus poor robustness. In spite of the use of neural tree models in practical problem s of classiﬁcation, regression, and forecasting, very little is known about the statistical properties of these models. In this paper, we obtained the upper-bound of n eural regression tree model parameters and designed an optimal NRT model to solve t he regression task of water quality problem in the paper manufacturing system. Motivated by the above discussion, we have proposed an optim al NRT model with optimal values of the model parameters. Harnessing the ense mble formulation, we try to exploit the strengths of RT and NN models and overcome thei r drawbacks. The 2 proposed model has the advantages of signiﬁcant accuracy, v ery less number of tuning parameters and easy interpretability as compared to more “b lack-box-like” advanced neural networks. Our model is very useful for small or medium sized complex data structures, like in the current water quality problem i n paper manufacturing industries. The theoretical bounds for model parameters ar e proposed. We practically recommend using the highest value of the model parameters as the optimal values when working with small and medium-sized data sets to achiev e ‘best’ performance of the model and this justiﬁes the name of the model as “optimal N RT model”. Finally, the robustness of the proposed algorithm is shown through ex perimental evaluation. The paper is structured as follows. Section 2 describes the m otivating problem and its formulation. In Section 3, we introduce optimal neural r egression tree model along with the upper bound of its parameters. Section 4 analyzes th e data by the proposed method and compares it with other state-of-the-art. Finall y, the concluding remarks are given in Section 5. 2. Motivating example In this section, we describe the motivating real-life probl em and the data collection plan for our study. The data on pH is taken on a regular basis fr om the Boiler lab. However, when the study was taken up, we considered only the m ost recent 12 months data. Measurement of water pH level at Boiler lab is done thro ugh a digital meter which is calibrated once a month using a standard solution. H ence, the normal gage R&R analysis was not necessary to be carried out to ensure cor rectness of the values observed using laboratory analysis. DM process outlet pH ha ppens to be the KPI of water treatment plant which gives water of desired quality s peciﬁed by the boiler. Water treatment plant comprises of two units: fresh water tr eatment and condensate water treatment. Freshwater treatment plant involves ﬁltr ation of Freshwater ob- tained from the tube well followed by treatment of water with chemicals, making the water suitable for Boiler and Turbine. Condensate treatmen t involves ﬁltration and treatment of condensate water obtained from Paper machines . Quality of condensate coming from Paper machines is the main constraint of the curr ent problem. Once the ﬁltration of water/condensate is done, it is sent to Demi neralization plant where water is made suitable to be used by the boiler. Details of the types of equipment used in DM plant processes are given (Vedelago and Millar 201 8) along with the process ﬂow diagram in Figure 1. •Strong acid cation (SAC) exchanger extracts the cations fro m the water by the cation resin and converts it into mineral acids. A vessel is p rovided therewith manually operated valves for controlling the ﬂow rates of th e inlet and outlet water. •Strong base anion (SBA) exchanger is an anion exchanger in wh ich mild-steel a vessel is internally lined with rubber to prevent corrosion . An external vessel is also provided with manually operated valves for controllin g the ﬂow rates of the inlet and outlet water. •The mixed bed (MB) unit comprises of a mild steel rubber-line d pressure vessel. Externally, the unit is provided with piping and valves to co ntrol the ﬂow of water during service and various regeneration stages. 3 Figure 1. Process Flow Diagram of DM Plant Process •Activated carbon ﬁlters (ACF) are used for the removal of chl orine and organic matter from water. Normally, the activated carbon ﬁlters ar e placed downstream of multi-grade ﬁlters and need to be regularly back-washed b y reversal of ﬂow so as to keep the surface of the carbon particles clean. Once the water is treated, it is stored in a tank called DM wate r storage tank (DMST). The water is pumped from DMST to Boiler via DM Transfe r Pump. Chem- icals like HCl and Caustic Soda are used for regeneration of t he resins in SAC, SBA and MB Exchanger. Morpholine is used for scaling the pH of DM W ater. Condensate water is also sent through Mixed Bed exchanger for treatment of water. Transferring DM water of pH in the range of 8.5-9.2 is essential to maintain drum water quality and the health of the water tubes running through Boiler. It was f ound that the variation in pH of DM water is large and maintaining proper pH level is es sential to maintain boiler water quality. The original problem can be considere d to be a two-stage prob- lem: 1) developing a prediction model for outlet water pH of D M process; 2) using the model and further actions to reduce variations in pH of DM out let water and maintain it close to the speciﬁed target (i.e., mid value of 8 .5−9.2). While developing a pre- diction model for outlet water pH of DM process, one needs to i dentify the important parameters of the process which together will be able to expl ain a signiﬁcant amount of variation present in the system. This work will also have a n impact on the amount of chemical usage, the health of the boiler tube and quality o f DM outlet water. 3. Model In this section, we will describe how a pre-trained RT can be r eformulated as a two-layered NN with similar types of predictive behavior. N eural networks and tree-based models are performing well in many real-life eng ineering problems. The idea of combining decision trees with neural networks was di scussed in previous literature but mainly for performing classiﬁcation tasks ( Sethi 1990; Sirat and Nadal 1990; Foresti and Micheloni 2002; Chen, Yang, and Abraham 20 07; Chakraborty, Chattopadhyay, and Chakraborty 2018). We extend the approa ch for regression problems and try to bridge the gap between theory and practic e by ﬁnding the upper bounds for the tuning parameters of the model. We present an e xample of a mapping from tree based model to a two-hidden layered neural network in Figure 2. 4 Figure 2. An example of neural regression tree model. A regression tree (left) and its corresponding NN structure (righ t). The circle nodes in the tree belong to split nodes and square nodes to leaf nodes. The p ath to the green shaded leaf (4) consists of all red nodes (0, 1, 3). Numbers in neurons cor respond to numbers in tree model nodes. The highlighted connections in the network are those relevant for the activity of the green neuron and its output value. We are given a training sample Dn={(X1,Y1),(X2,Y2),...,(Xn,Yn)}withn observations on pindependent variables. Neural regression tree (NRT) is a no npara- metric regression model in which pinput features X∈Cp= [0,1]pare observed and we try to predict a square integrable output function Y∈ /CA. A regression tree is a regression function estimate that uses a hierarchical axes -parallel split of the input space, where each tree node corresponds to one of the segment ation subsets in Cp. For simplicity and easy interpretability, let us consider o rdinary binary RT where a node has exactly two child nodes or zero child nodes (leaf no des). A regression tree consists of split nodes (for example, x(i)≥αfor some i∈ {1,2,...,p}and someα∈C) and leaf nodes. The feature space Cpis partitioned into axes-parallel hyper-rectangles. The standard splitting criteria MSE (me an squared error) has been used here to create RT. During prediction, the input is ﬁ rst passed into the tree root node. It is then iteratively transmitted to the lea f node that belongs to the subspace in which the input is located; this is repeated u ntil a leaf node is reached. If a leaf represents region R, then the natural regression function estimate takes the simple form tn(x) = (/summationtextn i=1Yi /BDXi∈R)/Nn(R),x∈R, whereNn(R) is the number of observations in cell R(where, by convention, 0 /0 = 0). In other words, the prediction for a query point xin leaf node Ris the average of the Yiof all training instances that fall into this region R. Let us assume for now that we have at hand a regression tree tn(whose construction eventually depends upon the data Dn), which takes constant values on each of k≥2 terminal nodes. It turns out that this estimate may be reinterpreted as two hidden layers neural networks, as summarized here. Let HL1 ={H1,...,H k−1}be the collection of all hyperplanes participating in the construction of tn. We note that each Hk′∈HL1 is of the form Hk′={x∈Cp:hk′(x) = 0}, where hk′(x) =x(ik′)−αik′for some (eventually data-dependent) ik′∈ {1,2,...,p}andαik′∈C. To reach the leaf of the query point x, we try to ﬁnd, for each hyperplane Hk′, the side on which xfalls (+1 for right and−1 for left). Using the above notations, the tree estimate tnis mapped to the neural network as discussed below. Also, an example is prese nted in Figure 2 with the description for a clear understanding of the neural regr ession tree model. 5 Designing the ﬁrst hidden layer ( HL1).The input layer (IL) supplies the features to the ﬁrst hidden layer of neurons which correspon ds tok−1 perceptrons, with the threshold activation function deﬁned as τ(hk′(x)) =τ(x(ik′)−αik′),where τ(u) = 2 /BDu≥0−1. So, for each split in the tree, there is a neuron in HL1 whose activity encodes the relative position of an input xwith respect to the concerned split. The output of the ﬁrst layer are ±1-vector ( τ(h1(x)),...,τ(hk−1(x))), which describes all decisions of the inner tree nodes (it also incl udes the nodes oﬀ the tree path ofx). Note that τ(hk′(x)) is +1 if xis on one side of the hyperplane Hk′,−1 ifxis on the other side of Hk′, and by convention +1 if x∈Hk′. It is important to remember that each neuron k′of this layer is connected to one and only one input x(ik′), and that the connection has weight 1 and oﬀset −αik′. Designing the second hidden layer ( HL2).HL1 outputs a ( k−1)-dimensional vector of ±1-bits that encodes the precise position of xin the leaves of the tree. The leaf node identity of xcan be extracted from the above-mentioned vector using a weighted combination of the bits along with an appropriate t hreshold function. The second hidden layer has kneurons, one for each leaf, and assigns a terminal cell to x. LetHL2 ={L1,...,L k}be the collection of all tree leaves, and let L(x) be the leaf containing x. We connect a unit k′fromHL1 to a unit k′′fromHL2 iﬀ the hyperplane Hk′is involved in the sequence of splits forming the path from th e root to the leaf Lk′′. The connection has weight +1 if the split by Hk′is from a node to a right child in that path, and −1 otherwise. Suppose we have ( u1(x),...,u k−1(x)) as the vector of±1-bits seen at the output of HL1, then the output vk′′(x)∈ {−1,1}of neuron k′′ isτ(/summationtext k′→k′′bk′,k′′uk′(x)+b0 k′′), where notation k′→k′′means that k′is connected to k′′andbk′,k′′=±1 is the corresponding weight. The oﬀset b0 k′′is set to b0 k′′=−l(k′′)+1/2, (1) wherel(k′′) is the length of the path from the root to Lk′′. To understand the intuition behindthe choice (1), we can see that there are exactly l(k′′) connections starting from the ﬁrst layer and pointing to k′′, and that /braceleftbigg/summationtext k′→k′′bk′,k′′uk′(x)−l(k′′)+1/2 = 1/2 ifx∈Lk′′/summationtext k′→k′′bk′,k′′uk′(x)−l(k′′)+1/2≤ −1/2 ifx/ne}ationslash∈Lk′′(2) Using (1), the argument of the threshold function is 1/2 if x∈Lk′′and is smaller than −1/2 otherwise. Hence vk′′(x) = 1 iﬀ the terminal cell of xisLk′′. To summarize, HL2 outputs a vector of ±1-bits (v1(x),...,vk(x)) whose components equal −1 except the one corresponding to the leaf L(x), which is +1. Output layer ( OL).Let(v1(x),...,vk(x))betheoutputofthe HL2.Ifvk′′(x) = 1, then the output layer computes the average ¯Yk′′of theYicorresponding to Xifalling inLk′′. This is equivalent to taking tn(x) =k/summationdisplay k′′=1wk′′vk′′(x)+bout, (3) wherewk′′=¯Yk′′/2 for all k′′∈ {1,...,k}, andbout=1 2/summationtextk k′′=1¯Yk′′. 6 3.1.Optimal Neural Regression Tree We have already described how a given RT can be mapped to a two h idden layered (2HL) NN model. Now we are going to critically summarize the f unctionality of the neural regression tree (NRT) model (see Figure 2) as follows : •An RT is built having ( kn−1) split nodes and knleaf nodes which is further mapped into a 2HL NN model having ( kn−1) andknhidden neurons in HL1 andHL2 respectively. Since the value of kdepends on n, so we are going to use kninstead of kin the rest of the paper. •InHL1, the neurons compute all the tree split decisions and indic ate the split directions for the inputs. Further, HL1 passes the information to HL2. •The probabilistic interpretation of the network output can be obtained by inter- preting the activation functions in the HL. The tree estimate tn, depending on Dn, can be seen as a neural network estimate. The architecture of this network is ﬁxed, and so are the weigh ts and oﬀsets of the three layers. A natural idea is then to keep the structure of t he network intact and let the parameters vary in a subsequent network training procedure with backpropagation training. In other words, once the connect ions between the neurons have been designed by the tree-to-network mapping, we can le arn network parameters in a better way by minimizing the empirical MSE for this netwo rk over the sample Dn. We propose an optimal NRT model in which we theoretically ﬁnd the upper bound of these free parameters which are nothing but the optimal va lues of the parameters for small and medium sample sized data sets. In order to incre ase the generalization capabilities of the NRT model, we replace the original relay -type activation function τ(u) with a hyperbolictangent activation function σ(u) := tanh( u) which has a chosen range from −1 to 1. More precisely, we use in optimal NRT model, σ1(u) =σ(β1u) at every neuron of the ﬁrst hidden layer and σ2(u) =σ(β2u) at every neuron of the second hidden layer. Here, β1andβ2are positive hyper-parameters that determine the contrast of the hyperbolic tangent activation: larger t he parameters β1andβ2, the sharper is the transition from −1 to 1. As β1andβ2approach to inﬁnity, the continu- ous functions σ1andσ2converge to the threshold function. We call our model optima l NRT in the sense that the optimal values (upper bounds for thi s case) of the tuning parameters kn,β1andβ2are derived (see Section 3.2). Besides eventually providin g better generalization, the hyperbolic tangent activation functions favor smoother de- cision boundaries and permit a relaxation of crisp tree node membership. Lastly, they allow operations with a smooth approximation of the discont inuous step activation function. This makes the loss function of the network diﬀeren tiable with respect to the parameters almost everywhere, and gradients can be back propagated to train the model. 3.2.Upper bound of model parameters Let us consider the tree in the ensemble and denote by G1≡G1(Dn), the bipartite graph which models the connections between the input vector sx= (x(1),...,x(p)) and the kn−1 hidden neurons of HL1. Similarly, let G2≡G2(Dn) be the bipartite graph represents the connections between the ﬁrst layer and theknhidden neurons of HL2. LetM(G1) be the set of p×(kn−1) matrices A= (aij) such that aij= 0 if (i,j)/∈G1. Also let M(G2) be the ( kn−1)×knmatrices B= (bij) such that bij= 0 7 if (i,j)/∈G2. The parameters that specify the ﬁrst hidden units are conta ined in a matrixAofM(G1) identiﬁed by the weights over the edges of G1and a column vector of biases b1and it is of size kn−1. Similarly, the parameters of the second hidden units are represented by a matrix BofM(G2) of weights over G2and by a column vectorb2of oﬀsets having size kn. Let us take the output weights and oﬀset to be Wout= (w1,...,w kn)⊤∈ /CAknandbout∈ /CA, respectively. And the parameters that specify the NRT model are represented by a “vector” as shown b elow: λ= (A,b1,B,b2,Wout,bout)∈M(G1)× /CAkn−1×M(G2)× /CAkn× /CAkn× /CA. We further assume that there exists a positive constant c1such that /bardblB/bardbl∞+/bardblb2/bardbl∞+/bardblWout/bardbl1+|bout| ≤c1kn, (4) where/bardbl·/bardbl∞is the supremum norm of a matrix and /bardbl·/bardbl1is theL1-norm of a vector. The rationale behind this assumption (4) is that the weights and oﬀsets are taken by the computation units of the second layer and the output laye r. We note that this condition is satisﬁed by the original random tree estimates as soon as Yis assumed to be bounded. Finally, we can assume that the absolute value of/bardblY/bardbl∞≤L <∞ almost surely, for some L. Therefore, letting Λ n=/braceleftbig λ= (A,b1,B,b2,Wout,bout) : (4) is satisﬁed/bracerightbig ,we see that the neural network implements functions of this particular form fλ(x) =W⊤ outσ2/parenleftBig B⊤σ1(A⊤x+b1)+b2/parenrightBig +bout, x∈Rp, whereλ∈Λn. Our aim is to tune the parameters λusing the data Dnsuch that the function realized by the obtained network becomes a ‘goo d’ estimate that can minimize the empirical error. Let Fn=/braceleftbig fλ:λ∈Λn/bracerightbig , where Fnbe the class of neural networks and mnbe the network that minimizes the empirical L2error which is deﬁned as Jn(f) =1 nn/summationdisplay i=1/vextendsingle/vextendsingleYi−f(Xi)/vextendsingle/vextendsingle2 over all functions f∈Fn, i.e.,Jn(mn)≤Jn(f), where Fnis a rich class of functions, including additive functions, polynomial functions havin g coeﬃcients of the same sign, products of continuous functions and etc. Thus, we will try to ﬁnd out an estimate mn:Cp→Rof the regression function m(X) =E[Y|X=x]. We say mnis consistent if E[mn(X)−m(X)]2tends to 0 as n→ ∞(the expectation is evaluated over Xand the training sample Dn). We can write using (Gy¨ orﬁ et al. 2006, Lemma 10.1)/BX/vextendsingle/vextendsinglemn(X)−m(X)/vextendsingle/vextendsingle2≤2 /BX/bracketleftbigg sup f∈Fn)/vextendsingle/vextendsingle/vextendsingle/vextendsingle1 nn/summationdisplay i=1/vextendsingle/vextendsingleYi−f(Xi)/vextendsingle/vextendsingle2− /BX/vextendsingle/vextendsingleY−f(X)/vextendsingle/vextendsingle2/vextendsingle/vextendsingle/vextendsingle/vextendsingle/bracketrightbigg + /BX/bracketleftbigg inf f∈Fn/integraldisplay Cd/vextendsingle/vextendsinglef(x)−m(x)/vextendsingle/vextendsingle2µ(dx)/bracketrightbigg ,(5) 8 whereµdenotes the distribution of X. To ﬁnd the upper-bound of the model parameters, we need to show the estimation error (ﬁrst term in the R.H.S. of Eqn. 5) and approximation error (second term in the R.H.S. of Eqn. 5) tend to 0. The former can be proved using non-asymptotic uniform devia tion inequalities and covering numbers corresponding to Fn, as shown in Theorem 3.1. Approximation error can be handled using a pseudo-estimate similar to RT ge nerated tnand application of Lipschitz property on the activation functi on of the NRT model, as shown in Theorem 3.2. Throughout the paper, we assumed Xis uniformly distributed inCpand/bardblY/bardbl∞≤L <∞almost surely, for some L. The next two theorems state that with certain restrictions i mposed on the number knof terminal nodes and with the parameters β1,β2being properly regulated as functions of n, the empirical L2risk-minimization provides upper bounds for the NRT model parameters. Larger the value of kn,β1andβ2, the better is the model. That is why, we call the upper bound of the model parameters as the opt imal values of NRT model. It is important to note that the depth of the tree and β1,β2are controlled for risk-minimization. Theorem 3.1. Ifkn,β2satisfykn,β2→ ∞such that kn=o(n1/6)andβ2= o(log(n)), then estimation error converges to 0asn→ ∞. Proof.For proof, refer to Appendix 1. Theorem 3.2. Ifkn,β1satisfykn,β1andβ2→ ∞such that kn=o(eβ2),kn=o(β1/4 1) andβ1=o(n2/3), then approximation error converges to 0asn→ ∞. Proof.For proof, refer to Appendix 2. Hence, the optimal choice of parameters obtained using empi rical risk minimization for optimal NRT model are: (a) kn=o(n1/6) (b)β1=o(n2/3) and (c) β2=o(log(n)) for practical use in regression problems having small and me dium-sized datasets. 4. Application 4.1.Data To identify the causal parameters causing the variations in the DM water outlet pH, at ﬁrst a brainstorming session was held with the process exp erts. Since not all the pa- rameters are controllable by the usersof the processes, a li st of controllable parameters is prepared and accordingly data are collected for the two pr ocesses (as shown in Fig- ure 1). The data set collected from the process over a year con sists of the observations from the following causal variables: Inlet Flow, Water Pres sure(water inlet pressureto the exchanger), Air Pressure, MB stroke and Amount of Morpho line/Chemical dosing (Liter per hr). The values of the response variable (DM water pH) varies between 7 to 10 (refer to Figure 3 for a graphical summary). Sample data sets for DMST-1 and DMST-2 are given in Tables 1 and 2. These data sets will be used for ﬁnding a pre- diction model that can help the company to forecast future wa ter pH level and take further steps for the reduction in water pH variation. 9 Table 1. Sample data set for DMST-1 Sl. No. Inlet Flow Water Pressure Air Pressure MB stroke Amount of DM water chemical outlet pH 1 1980 5.8 5.0 70 9.96 9.276 2 2150 7.0 7.0 60 8.69 9.094 3 1780 6.0 5.0 45 7.73 8.594 4 2808 5.2 6.4 50 6.54 8.738 5 1590 6.2 5.7 40 5.56 8.592 6 2995 6.0 6.0 50 7.23 9.099 . . . . . . . . . . . . . . . . . . . . . . . . . . . . Table 2. Sample data set for DMST-2 Sl. No. Inlet Flow Water Pressure Air Pressure MB stroke Amount of DM water chemical outlet pH 1 1489 5.4 6.0 80 10.04 9.246 2 1139 5.8 5.0 65 10.66 9.033 3 1448 5.4 5.0 45 7.69 8.483 4 1703 6.2 5.8 40 7.54 8.594 5 1258 6.2 5.7 35 5.99 7.589 6 1139 6.0 5.2 50 6.97 9.021 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2.Performance metrics The metrics used in this study to evaluate the performance of diﬀerent regression models (including the proposed model) are Root mean square e rror (RMSE), Mean absolute percentage error (MAPE), Co-eﬃcient of multiple d etermination ( R2) and Adjusted R2(AdjR2): RMSE=/radicalBig 1 n/summationtextn i=1(yi−/hatwideyi)2;MAPE=1 n/summationtextn i=1/vextendsingle/vextendsingle/vextendsingle/vextendsingleyi−/hatwideyi yi/vextendsingle/vextendsingle/vextendsingle/vextendsingle; R2= 1−/bracketleftbigg/summationtextn i=1(yi−/hatwideyi)2 /summationtextn i=1(yi−y)2/bracketrightbigg ;AdjR2= 1−/bracketleftbigg (1−R2)(n−1) n−k−1/bracketrightbigg ; where,yi,y,/hatwideyidenote the actual value, average value and predicted value o f the dependent variable, respectively for the ithinstant. Here nandkdenote the number ofdata pointsandindependentvariables usedfor performan ceevaluation, respectively. The lower the value of RMSE and MAPE and the higher the value of R2and Adj R2, the better is the model is. 4.3.Preliminary data analysis We started the analysis of the data with Anderson-Darling no rmality test on DM water outlet pH and it conﬁrms that the dependent variable do esn’t follow normal distribution (see Figure 3). An absence of normality in the D M water outlet pH data actually removed the possibility of application of convent ional parametric regression methods. This leads us to think about the nonparametric regr ession approaches. The stability of the process was checked using X−Rcontrol chart (see Figure 4). FromX−Rchart, it can be easily concluded that DM water outlet pH has h igh variation and it also indicates the presence of some assigna ble causes in the process. The dataset contains only numerical features and no missing entries, so no data cleaning task was performed.Weonly transformedthefeatur evalues into standardized form to achieve consistent results while using NN and NRT mod els. The standardiza- tion is done using the following equation: 10 Figure 3. Summary statistics for DM water outlet pH Xnew i=Xi−min(Xi) max(Xi)−min(Xi); whereXnew iis the standardized transformed value of the feature vector and lies within 0 and 1. Astheoutputofthenetworkswillalsoliebetween0and1, “lo gsig” isusedasoutput transfer function to bring it back in the original form at the end of the modeling by using, Ypred i=Zi{max(Xi)−min(Xi)}+min(Xi); hereZiis the standardized output within 0 and 1, and Ypred iis the predicted output of the model. 4.4.Results Sincethedatasetsdon’t satisfythebasicassumptionsof pa rametricregressionmodels, we decided to use nonparametric regression models. We furth er shuﬄed the observa- tions of the water pH dataset randomly and split it into train ing and testing data sets in a ratio of 70 : 30. Each experiment is repeated 5 times with d iﬀerent randomly as- signed training and test sets and we will ﬁnally report the av erages of the performance metrics observed over 5 times validations in Table 3. A few po pular nonparametric regression models such as k-Nearest neighbor (kNN), suppor tvector regression (SVR), Regression splines (B-splines), multiple adaptive regres sion spline (MARS), RT and NN with 2 hidden layers (NN with 2HL) were applied to the data s ets and the re- sults were recorded. Then we started experimenting with our proposed model. The training procedure for the optimal NRT model is as follows. R T is ﬁrst built using the scikit-learn implementation (Pedregosa et al. 2011) for tree designing. From the tree, we extracted the set of all split directions and split positi ons and used them to build neural network initialization parameters. The NRT models a re then trained using the TensorFlow framework (Abadi et al. 2016). The optimization with the net work model is done by minimizing the MSE on the training set. It is achiev ed by employing an iterative gradient-descent optimization technique. We ha ve used the default functions available in TensorFlow for this. Each neural network (including NRT) model was trained for 100 epochs. From the theoretical results on the u pper bounds of the model parameters, we can easily ﬁnd out the values of the initial co ntrast parameters of the activation functions of the ﬁrst and second HL of the model. T he value of β2will be 11 Figure 4. Control chart for DM water outlet pH lower than β1which is conﬁrmed from the theoretical bounds presented in S ection 3.2. This is practically very much signiﬁcant, since for a re latively small β2, the tran- sition in the activation function from −1 to +1 is smoother and a stronger gradient signal reaches the ﬁrst hidden layer in backpropagation tra ining. Similarly a converse explanation can be given for β1. In all experiments, we have used β1=o(n2/3) and β2=o(logn), where nis the number of training samples. Another property of the da ta set is that it is small sized and here tree based models perfor m better as compared to NN, since the latter typically need plentiful training da ta to perform well. For two data sets given in Section 4.1 standard NN (with 2HL) do not ac hieve comparable performance as compared to proposed optimal NRT model. In a s imilar way neural decision trees (NDT) were also implemented to asses its perf ormance on the water pH data sets. Since it has many free tuning parameters, the trai ning time and memory requirement are quite high as compared to the proposed model . Our proposed opti- mal NRT model is faster, especially when trained on a GPU. In T able 3, we present the results of diﬀerent regression model on the water pH data s ets and best results are displayed in bold fonts. From the experimental evaluati on of diﬀerent regression models, it can be concluded that our proposed optimal NRT mod el outperforms the state-of-the-arts. 4.5.Recommendation The proposed optimal NRT model will be useful for forecastin g the future values of water pH given the values of the process variables. This wi ll be helpful for the management and engineers to take preventive actions in the f uture. However, our objective of the work is not only to develop a competitive pre diction model for water pH level in DM process but also to ﬁnd out optimal level of proc ess parameters (in other words, causal variables) using statistical modeling . To keep DM water pH in the range of 8 .5−9.2 as speciﬁed by the boiler manual, the recommendations of ou r model is to maintain MB strokes, water pressure and chemical consumption within a speciﬁed range as shown in Table 4. Inlet ﬂow can’t be control led by the user of the process and Air pressure need not be controlled as far as the r ecommendation of the 12 Table 3. Quantitative measures of performance for diﬀerent regress ion models on test data set (average values of the metrics after 5-fold cross validations) Regression Models Data set RMSE MAPE R2AdjR2 kNNDMST 1 5.56 6.6760.65 56.21 DMST 2 6.44 7.9556.70 50.83 SVRDMST 1 4.56 5.5068.56 62.18 DMST 2 4.98 6.2063.70 57.83 B-splinesDMST 1 4.32 5.4069.85 63.30 DMST 2 6.94 7.2156.70 49.78 MARSDMST 1 4.29 5.2665.95 58.90 DMST 2 6.74 7.9357.53 47.05 RTDMST 1 3.44 4.1280.52 75.56 DMST 2 3.89 4.7876.56 71.23 NN with 2HLDMST 1 3.86 4.8076.95 70.03 DMST 2 4.12 5.9170.10 64.73 NDTDMST 1 3.18 3.5784.70 80.96 DMST 2 3.40 3.9282.18 78.14 Optimal NRTDMST 1 2.90 3.2586.80 83.55 DMST 2 3.15 3.6584.56 82.10 proposed model goes. Table 4 gives the optimum range of contr ollable parameters for DMST 1 and DMST 2 based on the recommendation of regression an alysis performed by optimal NRT model. Table 4. Optimal range of causal variables for achieving desired pH l evel Process Range for Range for Range for Expected range Water chemical for DM water Pressure MB stroke consumption outlet pH DMST 1 5.0-6.0 45-60 6.5-7.5 8.5-9.0 DMST 2 5.0-6.0 40-55 7.5-9.5 8.5-9.1 Table4depictsarangeforall theimportantprocessvariabl esforwhichtheexpected DM water outlet pHwill bewithintherequiredspeciﬁcation, i.e., 8.5−9.2. However, in ordertoﬁndouttheexact values oftheprocessparameters wi thintherangeprescribed in Table 4, a design of experiment would be necessary, which w as subsequently carried out to solve the problem. Though we have discussed only the pr oposed model and its accuracy level, our model also helped the manufacturing pro cess industry to improve its water quality level and to gain monetary beneﬁts due to a r eduction in the chemical consumption. 5. Conclusions This paper proposes a novel machine learning paradigm to sol ve a process eﬃciency problem in a paper manufacturing industry. The purpose of th is article is to develop a model for reducing the variation of water pH level in the DM p rocess in a manu- facturing industry. Our study proposed an optimal NRT model that maps RT into a 2-HL NN model with the theoretically obtained upper bound of the parameters of the model.Consequently, theproposedensemble modelsuccessf ully demonstratedthebest performanceand oﬀered a practical solution to the problem of ﬁnding optimal level for 13 the tuning parameters to improve water quality in the DM proc ess. The major advan- tage of the proposed optimal NRT model is that it has very less tuning parameters, easily interpretable as compared to “black-box-like” adva nced NN models and theo- retically sound as compared to many similar types of algorit hms. It was theoretically shown that there is some gain in considering 2HL in NN, but it i s not really necessary to go beyond 2HL in NN (Devroye, Gy¨ orﬁ, and Lugosi 2013). The proposed ensemble model will have an edge while working with limited data sets w here the beneﬁts of neural optimization can usually not be exploited, but with a speciﬁc inductive bias of the optimal NRT, it becomes possible. This also answers th e question one may ask about the need for a two-step pipeline (like optimal NRT mode l) over advanced NN models. Appendix 1. Proof of Theorem 3.1 The idea of the proof is based on the Theorem 16.1 of (Gy¨ orﬁ et al. 2006). The set Fn contains all neural networks constrained by Eqn. (4) of Sect ion 3.2 with inputs in Rp, two hidden layers of respective size kn−1 andkn, and one output unit. We note that Fnis deterministic, in the sense that it does not depend on Dn, but only on the size nof the original data set. We can derive a useful exponential i nequality to prove the Theorem 3.1. using uniformly bounded classes of functions. We have assumed that for each f∈Fn,fsatisﬁes/bardblf/bardbl∞≤c1knandYis also bounded ( /bardblY/bardbl∞≤L <∞). Letzn 1= (z1,...,zn) be a vector of nﬁxed points in /CApand let Hbe a set of functions from /CAp→ /CA. For every ε >0, we let N1(ε,H,zn 1) be the L1ε-covering numberof Hwith respect to z1,...,zn.N1(ε,H,zn 1) is deﬁnedas the smallest integer Nsuchthat thereexist functions h1,...,h N: /CAp→ /CAwiththepropertythat forevery h∈H, there is a j=j(h)∈ {1,...,N}such that 1 nn/summationdisplay i=1/vextendsingle/vextendsingleh(zi)−hj(zi)/vextendsingle/vextendsingle< ε. Note that if Zn 1= (Z1,...,Z n) is a sequence of i.i.d. random variables, then N1(ε,H,Zn 1) is a random variable as well. Now, let Z= (X,Y),Z1= (X1,Y1),...,Z n= (Xn,Yn), andCp= [0,1]p, we write Hn=/braceleftbigg h(x,y) :=|y−f(x)/vextendsingle/vextendsingle2: (x,y)∈Cp×[−L,L] andf∈Fn/bracerightbigg . Thefunctions in Hnwill satisfy the following: 0 ≤h(x,y)≤2c2 1k2 n+2L2≤4c2 1k2 n.This assumption holds for large nsuch that c1kn≥Lis satisﬁed. Using Pollard’s inequality 14 (Gy¨ orﬁ et al. 2006), we have, for arbitrary ε >0,/C8/braceleftbigg sup f∈Fn/vextendsingle/vextendsingle/vextendsingle1 nn/summationdisplay i=1/vextendsingle/vextendsingleYi−f(Xi)/vextendsingle/vextendsingle2− /BX/vextendsingle/vextendsingleY−f(X)/vextendsingle/vextendsingle2/vextendsingle/vextendsingle/vextendsingle> ε/bracerightbigg = /C8/braceleftbigg sup h∈Hn/vextendsingle/vextendsingle/vextendsingle1 nn/summationdisplay i=1h(Zi)− /BX(h(Z))/vextendsingle/vextendsingle/vextendsingle> ε/bracerightbigg ≤8E/bracketleftbigg N1/parenleftbigε 8,Hn,Zn 1/parenrightbig/bracketrightbigg exp/parenleftbigg −nε2 128(4c2 1k2n)2/parenrightbigg . (1) Next we try to bound the covering number/parenleftbig N1(ε 8,Hn,Zn 1/parenrightbig/parenrightbig . Let us consider two functions hi(x,y) =|y−fi(x)|2ofHnfor some fi∈Fnandi= 1,2. We get 1 nn/summationdisplay i=1/vextendsingle/vextendsingleh1(Zi)−h2(Zi)/vextendsingle/vextendsingle=1 nn/summationdisplay i=1/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingleYi−f1(Xi)/vextendsingle/vextendsingle2−/vextendsingle/vextendsingleYi−f2(Xi)/vextendsingle/vextendsingle2/vextendsingle/vextendsingle/vextendsingle =1 nn/summationdisplay i=1/vextendsingle/vextendsinglef1(Xi)−f2(Xi)/vextendsingle/vextendsingle×/vextendsingle/vextendsinglef1(Xi)−Yi+f2(Xi)−Yi/vextendsingle/vextendsingle ≤4c1kn nn/summationdisplay i=1/vextendsingle/vextendsinglef1(Xi)−f2(Xi)/vextendsingle/vextendsingle. Thus, N1/parenleftBigε 8,Hn,Zn 1/parenrightBig ≤N1/parenleftBigε 64c1kn,Fn,Xn 1/parenrightBig . (2) The covering number N1(ε 64c1kn,Fn,Xn 1) can be upper bounded independently of Xn 1 by extending the arguments of (Lugosi and Zeger 1995) from a n etwork with one hidden layer to a network with two hidden layers. We will now a pply Theorem 9.4, Lemma 16.4, and Lemma 16.5 repeatedly (Gy¨ orﬁ et al. 2006). L et the neurons of the ﬁrst hidden layer output belong to the class G1={σ1(a⊤x+a0) :a∈ /CAp,a0∈ /CA}, and for 0 < ε <1/4, N1(ε,G1,Xn 1)≤32/parenleftbigg4e εlog2×3e ε/parenrightbiggp+2 = 9/parenleftBig6e ε/parenrightBig4p+8 . Next, letting G2={bg:g∈G1,b∈[−c1kn,c1kn]}we get N1(ε,G2,Xn 1)≤4c1kn εN1/parenleftBigε 2c1kn,G1,Xn 1/parenrightBig ≤/parenleftBig36ec1kn ε/parenrightBig4p+9 . The second units compute the functions of the collection G3=/braceleftBig σ2/parenleftBigkn−1/summationdisplay i=1gi+b0/parenrightBig :gi∈G2,b0∈[−c1kn,c1kn]/bracerightBig . 15 Note that σ2satisﬁes the Lipschitz property |σ2(u)−σ2(v)| ≤β2|u−v|for all (u,v)∈/CA2. Thus, N1(ε,G3,Xn 1)≤2c1β2k2 n εN1/parenleftBigε 2β2kn,G2,Xn 1/parenrightBigkn−1 ≤/parenleftBig72ec1β2k2 n ε/parenrightBig(4p+9)kn+1 . Also, letting G4={wg:g∈G3,w∈[−c1kn,c1kn]}, By assuming without loss of generality c1,β2≥1, we get N1(ε,G4,Xn 1)≤4c1kn εN1/parenleftBigε 2c1kn,G3,Xn 1/parenrightBig ≤/parenleftBig144ec2 1β2k3 n ε/parenrightBig(4p+9)kn+2 . Finally, we can write Fn=/braceleftBigkn/summationdisplay i=1gi+bout:gi∈G4,bout∈[−c1kn,−c1kn]/bracerightBig , We conclude N1(ε,Fn,Xn 1)≤2c1kn(kn+1) ε/bracketleftbigg N1/parenleftBigε kn+1,G4,Xn 1/parenrightBig/bracketrightbiggkn ≤/parenleftBig144ec2 1β2(kn+1)4 ε/parenrightBig(4p+9)k2 n+2kn+1 . (3) Combining (1)-(3) together, we obtain/C8/braceleftbigg sup f∈Fn/vextendsingle/vextendsingle/vextendsingle1 nn/summationdisplay i=1/vextendsingle/vextendsingleYi−f(Xi)/vextendsingle/vextendsingle2− /BX/vextendsingle/vextendsingleY−f(X)/vextendsingle/vextendsingle2/vextendsingle/vextendsingle/vextendsingle> ε/bracerightbigg ≤8/parenleftBig9216ec3 1β2(kn+1)5 ε/parenrightBig(4p+9)k2 n+2kn+1 exp/parenleftbigg −nε2 2048c4 1k4n/parenrightbigg . Now let /tildewideZbe a non-negative random variable. Using this for any 0 < ε <1/4 and largen, we can write/BX/bracketleftbig/tildewideZ/bracketrightbig =/integraldisplay∞ 0 /C8[/tildewideZ > t]dt≤ε+/integraldisplay∞ ε /C8[/tildewideZ > t]dt. (4) 16 Using (4) for any 0 < ε <1/4 and large n, we can write/BX/bracketleftbigg sup f∈Fn/vextendsingle/vextendsingle/vextendsingle1 nn/summationdisplay i=1/vextendsingle/vextendsingleYi−f(Xi)/vextendsingle/vextendsingle2− /BX/vextendsingle/vextendsingleY−f(X)/vextendsingle/vextendsingle2/vextendsingle/vextendsingle/vextendsingle/bracketrightbigg ≤ε+8/integraldisplay∞ ε/parenleftBig9216ec3 1β2(kn+1)5 t/parenrightBig(4p+9)k2 n+2kn+1 exp/parenleftbigg −nt2 2048c4 1k4n/parenrightbigg dt ≤ε+8/parenleftBig2048c4 1k4 n nε/parenrightBig exp/bracketleftbigg/parenleftbig (4p+9)k2 n+2kn+1/parenrightbig log/parenleftBig9216ec3 1β2(kn+1)5 ε/parenrightBig −nε2 2048c4 1k4n/bracketrightbigg →2ε(n→ ∞),if the conditions of Theorem 3.1 holds. Asε→0, we get the estimation error converges to 0 and the proof is c omplete. Appendix 2. Proof of Theorem 3.2 To prove Theorem 3.2, let us consider a piecewise constant fu nction (pseudo-estimate) similar to the RT-tree tn, with only one diﬀerence: the function computes the true conditional expectation E[Y|X∈Lk′′] in each leaf Lk′′, but not the empiri- cal one, viz. ¯Yk′′. In another way, we can write ( W⋆ out)k′′=E[Y|X∈Lk′′]/2 and b⋆ out=/summationtextkn k′′=1E[Y|X∈Lk′′]/2 in Eqn. (3) of Section 3. This tree-type pseudo-estimate has the form tλ⋆(x) =W⋆⊤ outτ/parenleftBig B⋆⊤τ(A⋆⊤x+b⋆ 1)+b⋆ 2/parenrightBig +b⋆ out, x∈Rp, for some λ⋆= (A⋆,b⋆ 1,B⋆,b⋆ 2,W⋆ out,b⋆ out). We have/BX/bracketleftbigg inf f∈Fn/integraldisplay Cp/vextendsingle/vextendsinglef(x)−m(x)/vextendsingle/vextendsingle2µ(dx)/bracketrightbigg = /BX/bracketleftbigg inf f∈Fn/integraldisplay Cp/vextendsingle/vextendsinglef(x)−m(x)/vextendsingle/vextendsingle2µ(dx)/bracketrightbigg −2 /BX/bracketleftbigg/integraldisplay Cp/vextendsingle/vextendsingletλ⋆(x)−m(x)/vextendsingle/vextendsingle2µ(dx)/bracketrightbigg +2 /BX/bracketleftbigg/integraldisplay Cp/vextendsingle/vextendsingletλ⋆(x)−m(x)/vextendsingle/vextendsingle2µ(dx)/bracketrightbigg ≤2 /BX/bracketleftbigg/integraldisplay Cp/vextendsingle/vextendsinglefλ⋆(x)−tλ⋆(x)/vextendsingle/vextendsingle2µ(dx)/bracketrightbigg +2 /BX/bracketleftbigg/integraldisplay Cp/vextendsingle/vextendsingletλ⋆(x)−m(x)/vextendsingle/vextendsingle2µ(dx)/bracketrightbigg .(5) We are going to show that the two expectation terms in the R.H. S. of Eqn. (5) tends to zero under certain restrictions on kn,β1andβ2. The second term is easy to understand since m∈Fnand forn→ ∞it will converge to 0. But to deal with the ﬁrst term we take recourse to the following steps. By deﬁnition, tλ⋆(x) =W⋆⊤ out/bracketleftbigg τ/parenleftBig B⋆⊤τ(A⋆⊤x+b⋆ 1)+b⋆ 2/parenrightBig/bracketrightbigg +b⋆ out and fλ⋆(x) =W⋆⊤ out/bracketleftbigg σ2/parenleftBig B⋆⊤σ1(A⋆⊤x+b⋆ 1)+b⋆ 2/parenrightBig/bracketrightbigg +b⋆ out. Using the fact |a−b|2≤ |a|2− |b|2,/bardblY/bardbl∞≤L <∞and/bardblW⋆ out/bardbl ≤knL2 4; we can 17 write for all x∈Rp /vextendsingle/vextendsinglefλ⋆(x)−tλ⋆(x)/vextendsingle/vextendsingle2 ≤knL2 4×/vextenddouble/vextenddouble/vextenddoubleσ2/parenleftBig B⋆⊤σ1(A⋆⊤x+b⋆ 1)+b⋆ 2/parenrightBig −τ/parenleftBig B⋆⊤τ(A⋆⊤x+b⋆ 1)+b⋆ 2/parenrightBig/vextenddouble/vextenddouble/vextenddouble2 ≤2knL2 4/bracketleftbigg/vextenddouble/vextenddouble/vextenddoubleσ2/parenleftBig B⋆⊤τ(A⋆⊤x+b⋆ 1)+b⋆ 2/parenrightBig −τ/parenleftBig B⋆⊤τ(A⋆⊤x+b⋆ 1)+b⋆ 2/parenrightBig/vextenddouble/vextenddouble/vextenddouble2 +/vextenddouble/vextenddouble/vextenddoubleσ2/parenleftBig B⋆⊤σ1(A⋆⊤x+b⋆ 1)+b⋆ 2/parenrightBig −σ2/parenleftBig B⋆⊤τ(A⋆⊤x+b⋆ 1)+b⋆ 2/parenrightBig/vextenddouble/vextenddouble/vextenddouble2/bracketrightbigg =knL2 2/bracketleftbigg I1+I2/bracketrightbigg . (6) To ﬁndthe upperboundsof I1andI2, we will usethe following properties of functions: •Recall that σiistanhyperbolicactivation functionand τisathresholdactivation function, then wecan writefor all u∈ /CA,|σi(u)−τ(u)| ≤2e−2βi|u|for alli= 1,2. •Also,σisatisﬁes the Lipschitz property for continuous functions |σi(u)−σi(v)| ≤ βi|u−v|for all (u,v)∈ /CA2andi= 1,2. Using the above properties of functions, we can write I1≤4kn/summationdisplay j=1exp/bracketleftBig −4β2/vextendsingle/vextendsingle/parenleftbig B⋆⊤τ(A⋆⊤x+b⋆ 1)+b⋆ 2/parenrightbig j/vextendsingle/vextendsingle/bracketrightBig ≤4kne−2β2, since for every j,/vextendsingle/vextendsingle/parenleftbig B⋆⊤τ(A⋆⊤x+b⋆ 1) +b⋆ 2/parenrightbig j/vextendsingle/vextendsingle≥1/2 from the deﬁnition of (A⋆,b⋆ 1,B⋆,b⋆ 2) (see Section 3 of the paper). Using Lipschitz property, Cauchy-Schwarz inequality and t he deﬁnition of B⋆, we can ﬁnd the upper-bound for I2as follows: I2≤kn/summationdisplay j=1β2 2/vextendsingle/vextendsingle/vextendsingle/parenleftBig B⋆⊤/parenleftbig σ1(A⋆⊤x+b⋆ 1)−τ(A⋆⊤x+b⋆ 1)/parenrightBig j/vextendsingle/vextendsingle/vextendsingle2 ≤β2 2k2 n/vextenddouble/vextenddoubleσ1(A⋆⊤x+b⋆ 1)−τ(A⋆⊤x+b⋆ 1)/vextenddouble/vextenddouble2 ≤4β2 2k2 nkn−1/summationdisplay j=1exp/parenleftbigg −4β1|(A⋆⊤x+b⋆ 1)j|/parenrightbigg ≤4β2 2k3 nexp/parenleftbig −4β1ε/parenrightbig (for some ﬁxed jand arbitrary ε >0) For allnlarge enough, choosing ε=log(β1) 4β1, to getI2≤4β2 2k3 n β1. Putting these upper bounds of I1andI2together and using Eqn. (6) we get the following: E/integraldisplay Cp/vextendsingle/vextendsinglefλ⋆(x)−tλ⋆(x)/vextendsingle/vextendsingle2µ(dx)≤/bracketleftbigg 2L2k2 ne−2β2+2β2 2k4 nL2 β1/bracketrightbigg . (7) R.H.S. of (7) tends to 0 if the conditions of Theorem 3.2 holds and hence the proof. 18