on the Utilization of an ensemble of Meta-Heuristics for Simulating energy Consumption in Buildings

Predicting energy consumption has been a substantial topic because of its ability to lessen energy wastage and establish an acceptable overall operational efficiency. Thus, this research aims at creating a meta-heuristic-based method for autonomous simulation of heating and cooling loads of buildings. The developed method is envisioned on two tiers, whereas the first tier encompasses the use of a set of meta-heuristic algorithms to amplify the exploration and exploitation of Elman neural network through both parametric and structural learning. In this regard, 10 meta-heuristic were utilized, namely differential evolution, particle swarm optimization, invasive weed optimization, teaching-learning optimization, ant colony optimization, grey wolf optimization, grasshopper optimization, moth-flame optimization, antlion optimization, and arithmetic optimization. The second tier is designated for evaluating the meta-heuristic-based models through performance evaluation and statistical comparisons. An integrative ranking of the models is achieved using average ranking algorithm.

countries (Eom et al., 2020;Conevska et al., 2020). The amount of fuel consumed by this sector is equivalent to 2 billion tons of oil equivalent (TOE), representing 31% of fuel for worldwide energy use (International Energy Agency, 2015). More specifically, in terms of electricity and heating, about 0.84 billion TOE is consumed by the building sector. Moreover, the building sector consumes 12% of fresh water and contributes to 40% of global solid waste and 40% of CO2 emissions. Approximately 20 -25% and 30 -40% have been reported as consumed energy for developing and developed countries, respectively (Akande et al., 2015;Al-Sakkaf, et al., 2019). In the United States and European Union, up to 40% energy consumption is attributed to the building sector (Cao et al., 2016).
Based on the aforementioned points, improved sustainability in buildings is essential to reduce environmental impacts and improve the wellbeing of individuals. This is also dependent on the efficiency of several elements in buildings. For example, consumed energy in the heating and cooling processes could be managed by heating, ventilation, and air conditioning (HVAC) systems to ensure a comfortable indoor environment in buildings. Therefore, proper design of HVAC systems, based on climate and building attributes, is critical to the energy efficiency of buildings (McQuiston & Parker, 1982;Castaldo et al., 2018;Bui et al., 2019). In this regard, energy prediction models are required for facility managers to better understand sustainability in buildings (Ürge-Vorsatz et al., 2007;Egan et al., 2018;Fanti et al., 2018;Al-Sakkaf et al., 2019). Forward and inverse models are typically applied to evaluate the energy performance of buildings (Zhao et al., 2012). For forward modeling, the building attributes are utilized despite being not easily determined, decreasing the computation accuracy, increasing the computation time, and hindering the application for occupied buildings (Yezioro et al., 2008;Park et al., 2016). Inverse modeling, on the other hand, is a machine learning model that can be applied as an alternative to identify the correlation between energy consumption and building attributes or factors (Catalina et al., 2008;Yu et al., 2010;O'Neill & O'Neill, 2016). This modeling technique is easy and has a fast computation speed .
The application of machine learning has emerged as a powerful tool in building energy management . For example, Alvarez et al. (2018) predicted the energy performance of a house using artificial neural network (ANN) models. The models were developed using a dataset of different buildings located in Spain. The inputs to these models were the construction year and surface area while the model output comprised the U-opaque value. It was found that the constructed area factor was associated with the highest prediction power. The modeled U-Opaque values were compared against the real values for the same buildings, yielding a correlation coefficient of 0.967. The results confirmed the ability to provide an accurate estimate of the building's energy efficiency using some modeled building characteristics. Jihad & Tahiri (2018) applied different ANN models to predict the heating and cooling loads of residential buildings in Morocco. The forecasting model included six input variables; orientation, glazing rate, relative compactness, wall surface area, and building height and surface area. The training samples were generated using parametric analysis by changing the construction mode and building usage. The developed model yielded an accuracy of 98.7% and 97.6% for the training and testing data, respectively. This model could be useful to determine the energy consumption of a new building without the need for using simulation software nor calculating a thermodynamic balance. Moayedi et al. (2019) estimated the heating load in energy-efficient buildings using six machine learning techniques. These techniques were the multi-layer perceptron regressor, random forest, lazy locally weighted learning, alternating model tree, ElasticNet, and radial basis function regression. The calculated outcomes of the above-mentioned models were analyzed using five statistical indexes; root relative squared error (RRSE), root mean squared error (RMSE), mean absolute error (MAE), Pearson correlation coefficient (PCC), and relative absolute error (RAE). The random forest model yielded PCC, MAE, RMSE, RAE, and RRSE for the training dataset of 0.99, 0.19, 0.24, 2.08, and 2.38, respectively. Besides, it determined the PCC, MAE, RMSE, RAE, and RRSE for the testing dataset to be 0.99, 0.34, 0.46, 3.68, and 4.59, respectively. Therefore, it was nominated as the most powerful predictive network. Amber et al. (2018) compared the performances of five intelligent Yuan et al. (2018) adopted an ANN to simulate the seasonal hourly electricity consumption for three areas of the university compass in Japan. Six input parameters were used namely, day of the week, hour of the day, hourly dry-bulb temperature, hourly relative humidity, hourly global irradiance, and recorded hourly electricity consumption at the same hour and day of the week in the previous weeks. The Levenberg-Marquardt algorithm was used to train the developed feed-forward ANN model. The developed model accomplished a correlation coefficient that ranged from 96% to 99% at the training phase and from 95% to 99% at the testing phase. Zeng et al. (2020) presented a Gaussian regression process-based model for simulating electricity consumption of different types of buildings. It was derived that the developed model obtained higher prediction accuracies in office buildings than shopping malls and hotels attaining average deviations below 15%. It was also highlighted that the developed model could achieve a short testing time of 0.02 seconds per prediction.
In view of the above, it can be noticed that most of previous studies relied on classical machine learning models such as artificial neural network, multilayer perceptron and support vector machines to forecast building energy consumption. However, these conventional machine learning models are vulnerable to local minima entrapment and over fitting issues (Khan et al., 2020;Li et al., 2019). In addition, manual tuning-based prediction models suffer from inferior prediction capacity and long computational time originated from the need to calibrate the hyper parameters of the machine learning models (Jiang et al., 2020;Al-Allaf, 2011). Artificial neural network is the most widely utilized supervised learning algorithm that is usually trained using gradient descent back propagation algorithm. Its architecture is composed of input layer, output layer and hidden layers. The latter is composed of neurons which are responsible of analyzing the pattern of the fed dataset (Wang et al., 2021;Han et al., 2019). Multilayer perceptron is a feed forward artificial neural network that is composed of input layer, output layer and one or more hidden layers that are interposed between input and output layers, and it is normally trained using gradient descent algorithm (Hicham et al., 2017;Demirci et al., 2015). Support vector machines is one of the supervised learning algorithms that was proposed by Cortes and Vapnik (1995) and it can be used for either classification or regression problems. It is empowered based on mapping the original data into higher dimensional feature space through kernel functions, whereas the type of kernel function and its respective parameters are the most influential factors on the performance of support vector machines (Yetilmezsoy et al., 2021;Tezel & Buyukyildiz, 2016). In the light of foregoing, the main objective of the present study is to build a meta-heuristicbased Elman neural network model that an accurately simulate the heating and cooling loads of buildings. The developed model deploys utilizes a set of ten meta-heuristic algorithms for boosting the exploration-exploitation search capabilities of the ENN through optimizing both of its parameters and hyper parameters. The ten different types of meta-heuristics are: differential evolution (DE), particle swarm optimization (PSO), invasive weed optimization (IWO), teaching learning-based optimization (TLBO), ant colony optimization (ACO), grey wolf optimization (GWO), grasshopper optimization (GO), moth-flame optimization (MFO), antlion optimization (ALO), and arithmetic optimization (AO) algorithms. The performance of hybrid machine learning models is evaluated based on split validation using a set of performance metrics. The average ranking algorithm is adopted to create an integrative ranking of the hybrid models according to their performance with respect to several indicators. Besides, a correlation matrix is built in an attempt to examine the degree of correlation of the explanatory input variables with each other, and between the input and output variables.

PRoPoSed ModeL
The main objective of the present study is to build optimized Elman neural network model that can accurately estimate building energy consumption. As illustrated in Figure 1, the framework of the developed model is composed of two main components, namely hybrid machine learning model and model assessment. The first step of the first component is to construct the energy efficiency dataset. In this regard, the developed model is established based on a dataset published by the UCI machine learning repository (Asuncion & Newman, 2007) for a work carried out by Tsanas & Xifara (2012). The utilized dataset encompasses 768 data instances generated using the Ecotect energy analysis platform. The input explanatory features of the developed model incorporate a set of residential building characteristics including glazing area, glazing area distribution, orientation, overall height, roof area, wall area, surface area, and relative compactness. The output predicted variables of the developed model comprise the amounts of heating and cooling loads. The simulated heating and cooling loads are produced from different scenarios and combinations of the input variables. The numbers of possible values of glazing area, glazing area distribution, orientation, overall height, roof area, wall area, surface area, and relative compactness are 4, 6, 4, 2, 4, 7, 12, and 12, respectively.
The first component is designated for designing meta-heuristic-based Elman neural network model that can predict the amounts of heating and cooling loads. The developed model employs a set of meta-heuristic algorithms for amplifying the search capabilities of the ENN through both parametric and structural learning process. The learning process involves optimizing both the architecture and weighed connection of the ENN. This is expected to address the explorationexploitation trade-off dilemma which usually hinders the training process of artificial neural networks (Chong et al., 2021;Chiroma et al., 2016). The ENN is trained based on formulating a variable length single-objective optimization function designated for minimizing the mean absolute percentage error of heating and cooling loads. Elman neural network is one of the recurrent neural networks that was first proposed by the scientist Elman (1990). Its basic configuration comprises input layer, hidden layer, output layer and context layer. Elman neural network is augmented by additional context layers and feedback loop of unit delay element (Z -1 ) which aid in memorizing the computations of the previous time step. In this regard, preserving states between succeeding time steps creates a non-linear dynamic behavior of Elman neural network which improves its learning capability over feed forward neural network (Bianchi et al., 2017). The computations in the conventional feed forward neural networks depends solely on the current time step and don't account for preceding events which diminishes its learning abilities (kurach & Pawlowski, 2016;Ghasemi & Rasekhi, 2016). Meta-heuristics proved themselves as an efficient global search algorithms to improve the prediction performances of machine learning models in diverse civil engineering applications such as daily stream flow simulation (Malik et al., 2020), estimation of vertical settlement of raft-pile foundation , conceptual cost estimation of power plants (Hashemi et al., 2019) and prediction of structural failure of multistore buildings (Chatterjee et al., 2017). The developed model investigates the implementation of ten meta-heuristics such as DE, PSO, IWO, TLBO, ACO, GWO, GO, MFO, ALO and AO. In this regard, the applicability of ten different types of meta-heuristics are investigated as a result of the case dependency nature of meta-heuristics such that there is no absolute rule on how these algorithms can perform in a designated optimization problems. In this context, some meta-heuristics can provide better solutions than other implemented algorithms in a certain optimization problem. However, it may yield inferior solutions in other optimization problems (Teimouri et al., 2016;Martinez-Álvarez et al., 2014). There are two types of outputs from the first component of the model, whereas the first output includes the optimum transfer function, optimum numbers of hidden and context layers, optimum numbers of hidden and context neurons, and optimum weights. The second type involves the minimum mean absolute percentage error accomplished by the trained models.
Selection of the best performing heating and cooling load prediction model is a vital and complicated process (Shahsavar et al., 2020;Dutta et al., 2019). Hence, the second component of the developed model aims at evaluating the accuracy levels achieved by the hybrid meta-heuristicbased Elman neural network model as per two folds of assessment. In the first fold, the developed models are assessed based on split validation using a set of performance metrics. The utilized performance indicators are MAPE, RMSE, PCC, Nash-Sutcliffe coefficient of efficiency (NSCE), and Willmott's index of agreement (WI). The average ranking algorithm is adopted to create an integrative ranking of the hybrid machine learning models according to their performance with respect to several indicators. The second fold involves building a correlation matrix for the heating and cooling loads in an attempt to examine the degree of correlation and association of the input variables with each other, and between the input and output variables. All the algorithms and mathematical operations in this research paper are implemented using Matlab R2017a.

ModeL deVeLoPMeNT
This section describes the background of the utilized meta-heuristics and the designed automated training optimization function.

Differential Evolution
This algorithm starts by initiating different individuals (i.e., candidate solutions) in the population. These old solutions are updated to produce new solutions only one time in each iteration. Each solution is updated using three different random solutions and the mutation factor. One of the random solutions searches around the space, while the two other solutions create a distance from the first one. There must be a mixture between the previous and current solutions to generate new solutions in the next iteration. The selection operation phase is performed to maintain better solutions in the current iteration. Finally, the solution associated with the lowest fitness function is chosen as the best solution (Storn & Price, 1997).

Particle Swarm Optimization
This algorithm, which was introduced by Eberhart & Kennedy (1995), imitates the movement of a school of fish or flocks of birds searching for their food. This algorithm commences by forming the particles, whereas each particle in the swarm is characterized by its current position, best position, and velocity. The position of a particle is updated by considering randomized values in some directions. Meanwhile, the velocity of any particle relies on particle best position, global best position, and the random function (Shi & Eberhart, 1998). This process is repeated until reaching a stopping criterion (e.g., satisfactory solution, maximum number of iterations, and constant fitness for a certain number of iterations). Finally, the current position of the best particle in the last iteration is defined as the global best solution .

Invasive Weed Optimization
This algorithm, which was proposed by Mehrabian & Lucas (2006), simulates the invasive behavior of weeds in finding the most suitable place for growth and reproduction (Elshaboury et al., 2021). The randomness, resistance, and adaptability of weeds are imitated in this algorithm. By definition, weeds are powerful herbs that grow unintentionally and cause a serious threat to crops. The basic features of this algorithm are considered as follows: (1) population initialization where the seeds are distributed in the search space, (2) reproduction each flowering plant produces seeds based on its fitness value, (3) spectral spread that eliminates inappropriate plants, (4) competitive deprivation that ensures removing the grass with worst fitness from the colony to keep a constant number of herbs in the colony, (5) repeating the previous steps until reaching the maximum number of iterations, and (6) storing the minimum cost function of the grasses (Misaghi & Yaghoobi, 2019).

Teaching Learning-Based Optimization
It is a recent meta-heuristics algorithm which is inspired by the behaviors of teachers and learners in a classroom. The teacher is determined to be the learner with the best grade in the population. The teacher is responsible for training the learners and improving the mean grade of the class (known as the teacher phase). Meanwhile, each learner randomly interacts with another learner (known as the learner phase). In this algorithm, a group of students (i.e., learner) is considered as the population. Besides, the different subjects resemble the optimization design variables. The fitness function of the optimization problem is determined based on the results of the learner. This process is repeated iteratively until the termination condition is met and the best solution in the entire population is considered as the teacher (Rao et al., 2011).

Ant Colony Optimization
This algorithm, which was proposed by Dorigo et al. (1999), was inspired by the behavior of ant colonies searching for their food sources. The ants communicate with each other using the pheromone, such that the ant that finds a food source deposits pheromone along the path to be followed by other ants. In this way, ants could find the shortest distance between their nests and the food source (Tran et al., 2019). The procedures of this algorithm are described as follows: (1) deriving solutions to the problem by using greedy information and pheromone values and (2) updating the pheromone values using the best solutions. Therefore, this optimization technique relies on learning from positive examples (i.e., positive learning) (Nurcahyadi & Blum, 2021).

Grey Wolf Optimization
This algorithm, which was introduced by Mirjalili et al. (2014), is inspired by the leadership and hunting process of grey wolves (Panda & Das, 2019). The leadership hierarchy is simulated using four types of grey wolves; alpha, beta, delta, and omega. The alphas are responsible for deciding the hunting time and resting place for the whole group. The beta wolves maintain discipline in the group and advise the alphas in their decisions. The delta wolves are authorized to dominate omegas while following the orders of alphas and betas. The last group of wolves follows the orders of all other dominant wolves. The hunting mechanism of grey wolves comprises tracking and chasing the prey, harassing the target, and attacking the prey target (Jitkongchuen et al., 2016).

Grasshopper Optimization
This algorithm, which was developed by Saremi et al. (2017), mimics the social behavior of grasshoppers. Grasshoppers are considered harmful insects and cause damage to crops. The lifecycle of grasshopper consists of three stages; egg, nymph, and adult. The nymph grasshopper is born when the eggs hatch. These young insects start rolling and eat everything that comes in their way. Then, they become adults and start swarming in the air (Gad et al., 2020). Each insect in the swarm is characterized by two important motions; cooperation between grasshoppers in larvae and insect phases and foraging for the food movement. A swarm is formed of randomly generated grasshopper agents. The fitness value evaluation determines the best search agent, and the search agents start moving toward the best grasshopper search agent (Ullah et al., 2020).

Moth-Flame Optimization
This algorithm, which was introduced by Mirjalili (2015a), is inspired by the navigational mechanisms of moths during the night. The moth maintains a fixed angle of the moonlight to be able to travel long distances in a straight path. Since the moon is far away from the moth, it flies in a straight line by using the near-parallel light near the surface. Despite the effectiveness of lateral orientation, moths often fly spirally around lights. The same case applies when the moths preserve a stable angle with the light source in the case of human-made artificial light. This algorithm commences by generating random moths within the search space, then evaluating the position (i.e., fitness value) of each moth, and determining the best position by flame. The moths' positions are updated based on a spiral movement function, and the previous processes are updated until the termination criteria are met.

Antlion Optimization
This algorithm mimics the hunting mechanism of antlions (Mirjalili, 2015b). Antlions act as search agents that hunt preys by building traps, entrapping ants, catching preys, and rebuilding traps (Tian et al., 2018). This algorithm is characterized by a good exploration capability with the help of random walk and roulette wheel to build traps. Besides, the exploitation efficiency of this algorithm is highlighted by the time-varying boundary shrinking mechanism and elitism. The major advantages of this algorithm are as follows: ease of implementation, high precision, avoidance of local optima, and reduced need for parameter adjustment (Horng & Lee, 2021).

Arithmetic Optimization
This algorithm is motivated by the application of arithmetic operators (i.e., multiplication, division, subtraction, and addition) in solving arithmetic problems. This optimization algorithm starts by generating a set of candidate solutions. The arithmetic operators determine the feasible positions of the near-optimal solution. The position of each solution is updated based on the exploration and exploitation mechanisms. Eventually, the algorithm stops by satisfying the end criterion (Abualigah et al., 2021).

Automated Training of elman Neural Network
The developed meta-heuristic-based Elman neural network model is characterized by its both parametric and structural learning nature. It aims at the autonomous optimization of the numbers of hidden layers, hidden neurons, context neurons, context layers, type of transfer function, and values of weighted connections in the Elman neural network. The developed model explores the implementation of eight activation functions, namely radial basis function, normalized radial basis, triangular basis function, linear function, positive linear function, log-sigmoid function, hyperbolic tangent sigmoid function and Elliot symmetric sigmoid function. As a result of the structural learning features of the developed model, the number of weighted connections varies adaptively during each training epoch according to the numbers of the hidden layer, hidden neurons, and context layers. Hence, the number of weighted connections should be computed during each training epoch to design an efficient variable-length optimization model. The developed estimator is created for the sake of handling the variability in the length of the optimization model, and it can be mathematically expressed as follows.
Where; I . and N denote the numbers of input neurons and hidden neurons, respectively. C represents the number of context neurons. P denotes the number of hidden and context layers. O stands for the number of output neurons. It is with mentioning that architecture of Elman neural network is triggered by the type of simulation problem (heating or cooling), and the search exploration ability of the employed meta-heuristic optimizer.

ModeL ASSeSSMeNT
This section describes the scoring metrics used to evaluate the developed prediction models and the average ranking algorithm.

Performance Comparison
The present research study utilizes five performance indicators to evaluate the developed models. In this context, the evaluation metrics of MAPE, RMSE, NSCE, WI, and PCC can be computed using Equations 2-6, respectively (Altunkaynak, 2019;Guo et al., 2019;Natarajan & Nachimuthu, 2020).
where; P i and O i indicate the predicted and observed heating or cooling load, respectively. O and P are the average observed and predicted heating or cooling load, respectively. k stands for the size of data observations. A smaller value of MAPE and RMSE indicates a smaller deviation of the predicted values from the actual values. The values of NSCE and WI range from 0 to 1 while the values of PCC vary from -1 to 1. In this regard, the closer values of NSCE , WI, and PCC to 1 imply an increasing agreement between the predicted and actual heating or cooling loads.

Unified Ranking
Each one of the developed meta-heuristic-based Elman neural network models behaves differently according to the tackled type of the performance indicator. As such, the average ranking algorithm is utilized to generate a final sorting of the meta-heuristic-based Elman neural network models according to their scores with respect to the utilized performance evaluation metrics. The mean and standard deviation of the rankings can be obtained using Equations 7 and 8, respectively (Yu et al., 2018;Mohammed Abdelkader et al., 2021).
where; ¼ a and Ã a indicate the mean and standard deviation of the rankings of the meta-heuristicbased Elman neural netwoek models, respectively. B denotes the number of performance metrics. A lower value of ¼ a and Ã a imply a high-performing and more robust prediction model.

ReSULTS ANd dISCUSSIoN
The utilized dataset is composed of 768 observations, whereas 614 (80%) and 154 (20%) data instances are randomly selected and deployed for training and testing purposes. Figures 2, 3, 4 and 5 describe the relationship between each of the influential independent variables and the output dependent variables. Table 1 summarizes a set of statistical indicators for the input and output variables. The statistical analysis includes minimum, maximum, average, median, standard deviation, mode, skewness, sample variance, range, kurtosis, and coefficient of variation. For instance, the highest value of standard deviation is for glazing area distribution (88.03) followed by overall height (45.14) while the glazing area has the least standard deviation (0.11). With regards to the coefficient of variation, the surface area (56.8%) is associated with the highest coefficient of variation and then relative compactness (55.11%). The least coefficients of variation are for glazing area (13.83%), glazing area distribution (13.11%) and orientation (13.69%).

Figure 2. Spatial distributions of glazing area and glazing area distribution against heating and cooling loads
Different initializations of parameters' configuration were experimented in order to search for their optimum setting. In order to create a fair comparison between the ten meta-heuristics, the number of iterations and population size are assumed in all meta-heuristics to be 200 and 50, respectively. In the differential evolution algorithm, the mutation rate is assumed to be uniform distribution ranging between 0.2 and 0.8. In addition, the crossover probability rate is set as 0.2. In the particle swarm optimization algorithm, the social parameter and cognitive learning parameter are assumed as two, and inertia weight is assumed as 0.5. In the invasive weed optimization algorithm, the non-linear modulation index is assumed as two, and the initial and final standard deviations are 0.5 and 0.001, respectively. The maximum and minimum numbers of seeds are equal to five and zero, respectively. With regards to the ant colony optimization algorithm, the sample size, intensification factor and deviation-distance ratio are assumed as 40, 0.5 and 1, respectively. For the grey wolf optimization algorithm, the trade-off parameter which manages the balance between exploitation and exploration is assumed to be linearly varying from two to zero. In the grasshopper optimization, the maximum and minimum deceleration values of grasshoppers when approaching food source and eating it are assumed 1 and 0.00004, respectively. In the moth-flame optimization algorithm, the convergence constant is assumed to be decreasing from -1 and -2, and the logarithmic spiral motion constant is set as 1. The search control parameter and exploitation control parameter in the arithmetic optimization algorithm are assumed 0.5 and 5, respectively. Each meta-heuristic was then run for five times in an independent manner to circumvent unstable solutions and the best solution associated with the run of the least prediction error was appended to be used for the forthcoming performance comparison analysis.
As mentioned earlier, ten hybrid meta-heuristic-based Elma neural network models are developed for the sake of predicting the heating and cooling loads. The convergence of the models for predicting heating loads is demonstrated in Figure 6. ENN-PSO accomplished the least MAPE and ENN-GWO The optimum structures of the developed meta-heuristic-based Elman neural network models for the prediction of heating and cooling loads are reported in Tables 2 and 3, respectively. It can be noticed that different architectures are obtained for either the prediction of the amount of heating or cooling loads, which imply different optimal solutions are generated by the ten meta-heuristic Elman neural

Figure 6. Convergence of the meta-heuristic-based Elman neural network models for forecasting heating loads
network models. This exemplifies that the complexity of optimization problem and the explorationexploitation search capabilities play a principal role in the quality of generated optimal solutions, which evinces the need for studying the need for studying the implementation of five different types of meta-heuristics. In the prediction of heating loads, the optimum structure of ENN-PSO consists of 3 hidden and context layers, 4 hidden and context neurons, and hyperbolic tangent sigmoid as the optimum transfer function. The optimum structure of ENN-GO is composed of 4 hidden and context layers, 2 hidden and context neurons, and linear function as the optimum transfer function. At the level of predicting cooling loads, the optimum structure of ENN-GWO comprises 3 hidden and context layers, 6 hidden and context neurons, and hyperbolic tangent sigmoid as the optimum transfer function. The optimum topology of ENN-GO encompasses 2 hidden and context layers, 8 hidden and context neurons, and log-sigmoid as the optimum transfer function. In addition to that, it is found that ENN-DE obtained varying topologies in the interpretation of heating and loading patterns. For instance, ENN-DE obtained one hidden layer, one context layer, two hidden and context neurons and hyperbolic tangent sigmoid function in heating loads assessment. In cooling loads evaluation, ENN-DE achieved an optimum topology of two hidden layers, two context layers, seven hidden and context neurons, and Elliot symmetric sigmoid activation function.
Illustrations of the performances of ENN-PSO, ENN-GWO, ENN-ACO, and ENN-ALO in the prediction of heating loads are depicted in Figures 8 and 9. X-axis represents a certain data instance

Figure 7. Convergence of the meta-heuristic-based Elman neural network models for forecasting cooling loads
which corresponds to a specific scenario for a combination of the input variables. Y-axis denotes the corresponding heating load, whereas five hundred data instances that were randomly sampled from training and testing partitions are plotted in these figures. The black and blue curves denote the actual and predicted heating values that are obtained by the meta-heuristic-based Elman neural network models It can be inferred that ENN-PSO, ENN-GWO have successfully predicted the heating loads pattern, whereas they produce very close predicted values to the actual heating loads. However, ENN-ACO and ENN-ALO failed to predict the actual heating loads efficiently providing considerable   Figures 10 and 11 show the performances of ENN-PSO, ENN-GWO, ENN-ACO, and ENN-ALO in the prediction of cooling loads. Y-axis represents the cooling load in these figures, ao five hundred data instances were randomly picked from the training and testing partitions, and visualized. The black and blue curves depict the actual and predicted amounts of cooling loads generated from the meta-heuristic-based Elman neural network models. Graphical comparisons reveal that the cooling pattern is better recognized by ENN-PSO, ENN-GWO than by ENN-ACO and ENN-ALO. In this regard, the prediction trend of the cooling loads provided by the models of ENN-PSO, ENN-GWO models resembled the measured actual values. However, explicit variations are encountered between the actual and predicted values by the models of ENN-ACO and ENN-ALO. Figure 12 describes the correlation between the expected and measured heating loads for all the training and testing instances based on ENN-PSO. The regression chart shows that the determination coefficient ( R 2 ) is 92.6%, which demonstrates a very good agreement between the predicted and measured heating loads. A graphical comparison between the predicted and actual cooling loads for all data instances based on ENN-GWO is shown in Figure 13. It can be found that an acceptable consistency ( R % 2 88 3 = . ) is obtained between the predicted and actual cooling loads. A performance comparative analysis between the meta-heuristic-based Elman neural network models for predicting heating loads is reported in Table 4. It can be concluded that the ENN-PSO model  Table 5. It can be noticed that ENN-GWO achieved the highest prediction performance, whereas it provided MAPE, RMSE, NSCE, WI, and PCC of 8.154%, 3.343, 0.876, 0.956, and 0.884, respectively. The developed models of ENN-PSO and ENN-IWO yielded the second highest prediction performance. In addition, it is observed that ENN-AO attained the lowest prediction performance, such that MAPE, RMSE, NSCE, WI, and PCC are equal to 17.342%, 6.156, 0.581, 0.844, and 0.781, respectively. As such, it can be interpreted that the developed ENN-PSO and ENN-GWO provided lower error rate than ENN-AO by 37.46% and 35.18% in the prediction of amounts of heating and cooling loads, respectively. At the level of both heating and cooling loads, it can be argued that the developed modes of ENN-PSO, ENN-GWO, ENN-IWO and ENN-TLBO demonstrated satisfactory good prediction performance. This can be explained by the high exploration and local search behaviors of these models which allowed them to efficiently investigate the full range of entire search space meanwhile avoiding local minima stagnation. On the The average ranking algorithm is utilized to blend the performances of the ten meta-heuristicbased Elman neural network models with respect to the five evaluation metrics. The results generated from the average ranking algorithm for the prediction of heating and cooling loads are presented in Tables 6 and 7, respectively. In the prediction of heating loads, it is found that ENN-PSO accomplished the first ranking and the most stable performance followed by ENN-TLBO while ENN-AO provided the tenth ranking. In this regard, µ a and Ã a of ENN-PSO are equal to 1 and 0, respectively. For ENN-TLBO, µ a and Ã a are equal to 2.4 and 0.49, respectively. For cooling loads, ENN-GWO accomplished the first ranking followed by ENN-PSO and then ENN-IWO. ENN-GWO achieved µ a and Ã a of 1.6 and 1.2, respectively. It is also shown that µ a and Ã a of ENN-IWO are equal to 2.4 and 0.49, respectively. Furthermore, it is worth noting that the performances of meta-heuristics vary from one optimization problem due to their case dependent nature, whereas some meta-heuristics can outperform their counterparts in some sort of problem. However, they may underperform in other optimization problems. For instance, ENN-ACO accomplished higher rank than ENN-ALO in heating loads prediction. However, it is outranked by it in cooling loads prediction. Moreover, ENN-TLBO Further comparative performance evaluations of meta-heuristic-based Elman neural network models against classical machine learning models are reported in Table 8. In the back propagation artificial neural network (BPANN) model, the numbers of hidden layers and hidden neurons are assumed 5 and 2, respectively. The momentum coefficient and learning rate are assumed 0.8 and 0.01, respectively. The activation function is log sigmoid and gradient descent algorithm is used for  training the network. In the support vector machines (SVM) model, Gaussian radial basis function is utilized with sigma and scaling factor equal to 1. In both heating and cooling loads, the developed ENN-PSO and ENN-GWO performed better than BPANN and SVM. In this context, ENN-PSO and ENN-GWO provided higher prediction accuracies than BPANN and SVM in heating and cooling,  Computational time is an important aspect when comparing the prediction models. Table 9 reports the average training time and average testing time of the meta-heuristic-based Elman neural network models. The used laptop is of specifications: Intel Core i7 CPU, 2.21 GHz and 16 GB of memory. It can be observed that ENN-GWO and ENN-MFO are the most time efficient prediction models, whereas they consumed lesser training times than the remainder of other prediction models. In this regard, the training times of ENN-GWO and ENN-MFO are 10.88 and 14.26 minutes, respectively. It can be also interpreted that ENN-IWO and ENN-AO have relatively short training times of 16.34 and 19.42 minutes, respectively. On the other hand, ENN-GO has the longest training time with 94.49 minutes. In addition, it can be noticed that most of the developed meta-heuristic-based Elman neural network models nearly consume the same testing times, whereas all of their testing times vary from 1.25 to 2.45 seconds. Table 10 describes the correlation between the input and output variables. As can be seen, the glazing area has a negative correlation with orientation and overall height. The roof area has a positive strong correlation with heating load and cooling load. The overall height has a negative strong correlation with heating load and cooling load. In this regard, the correlation coefficients of the roof area with heating load and cooling load are 0.889 and 0.896, respectively. Additionally, the correlation coefficients of overall height with heating and cooling loads are -0.862 and -0.863, respectively. It is also revealed that wall area and relative compactness manifested low levels of dependencies with heating and cooling loads. It is also evinced that the amount of heating load is highly correlated with cooling load (correlation coefficient is 0.976).

CoNCLUSIoN
Analyzing energy performance is crucial to maintain sustainable energy consumption in residential buildings. As such, this research paper proposed a novel hybrid meta-heuristic-based Elman neural network model for autonomous forecasting of amounts of heating and cooling loads in residential buildings. Performance comparison demonstrated that the complexity of the designated optimization problem and search capability of meta-heuristic significantly implicate the performance of meta- heuristics. In this context, different meta-heuristics obtained different optimum architectures in the same optimization problem. In addition, same meta-heuristic yielded different optimum topologies in different optimization problems (heating and cooling assessment). Comparative analysis also illustrated that the developed ENN-PSO and ENN-GWO accomplished the highest performance accuracies among the meta-heuristic-based Elman neural network models in heating and cooling loads prediction, respectively. In this regard, they managed to reduce the prediction error of ENN-AO by 37.46% and 35.18% in heating and cooling loads, respectively. Against the ENN-DE model, the developed ENN-PSO and ENN-GWO managed to improve the prediction accuracies by 13.44% and 19.55% in heating and cooling loads, respectively. The outcome of the average ranking algorithm manifested that ENN-PSO, ENN-AO, and ENN-ALO models exhibited the most stable ranking in the heating loads prediction ( a Ã = 0) meanwhile the ENN-TLBO, ENN-MFO, ENN-AO, and ENN-ALO models accomplished the lowest standard deviation of rankings ( 0.4 a Ã = ) in cooling loads. In addition to that, the average ranking algorithm exemplified that the relative importance (i.e., ranking) varied from heating to cooling loads prediction. For example, EN-PSO obtained a higher ranking (first) than ENN-GWO (third) in heating loads simulation. Nevertheless, ENN-GWO (first) outranked ENN-PSO (second) in cooling loads simulation.
When compared against the conventional machine learning models of BPANN and SVM, the developed ENN-PSO succeeded in diminishing the forecasting error of heating loads by 29.75% and 27.68%, respectively. Additionally, the developed ENN-GWO was able to reduce the prediction error of cooling loads estimation by 40.84% and 14.98%, respectively. Analysis of computation times of the meta-heuristic-based Elman neural network models revealed that ENN-GWO and ENN-MFO had the shortest training times. Besides, ENN-IWO and ENN-AO consumed comparatively short training times. It was also observed that ENN-GO requires substantially long training time although it produces high prediction error (eighth ranking in both heating and cooling). The constructed correlation matrices illustrated that the highest level of correlation were exhibited between the input factors of roof area and overall height, and the output variables of heating and cooling loads. It was also noticed that wall area and relative compactness are lowly correlated with heating and cooling loads. It is expected that the developed model can aid architectural designers through providing them with a platform to study the implications of several design parameters on buildings' energy efficiency.


Stationary Behavior of Constant Stepsize SGD Type Algorithms: An Asymptotic Characterization

Stochastic approximation (SA) and stochastic gradient descent (SGD) algorithms are work-horses for modern machine learning algorithms. Their constant stepsize variants are preferred in practice due to fast convergence behavior. However, constant step stochastic iterative algorithms do not converge asymptotically to the optimal solution, but instead have a stationary distribution, which in general cannot be analytically characterized. In this work, we study the asymptotic behavior of the appropriately scaled stationary distribution, in the limit when the constant stepsize goes to zero. Specifically, we consider the following three settings: (1) SGD algorithms with smooth and strongly convex objective, (2) linear SA algorithms involving a Hurwitz matrix, and (3) nonlinear SA algorithms involving a contractive operator. When the iterate is scaled by $1/\sqrt{\alpha}$, where $\alpha$ is the constant stepsize, we show that the limiting scaled stationary distribution is a solution of an integral equation. Under a uniqueness assumption (which can be removed in certain settings) on this equation, we further characterize the limiting distribution as a Gaussian distribution whose covariance matrix is the unique solution of a suitable Lyapunov equation. For SA algorithms beyond these cases, our numerical experiments suggest that unlike central limit theorem type results: (1) the scaling factor need not be $1/\sqrt{\alpha}$, and (2) the limiting distribution need not be Gaussian. Based on the numerical study, we come up with a formula to determine the right scaling factor, and make insightful connection to the Euler-Maruyama discretization scheme for approximating stochastic differential equations.

1. Introduction. Stochastic approximation (SA) algorithms are the major work-horses for solving large-scale optimization and machine learning problems. Theoretically, to achieve asymptotic convergence, we should use diminishing stepsizes (learning rates) with proper decay rate (Nemirovski et al., 2009). However, constant stepsize SA algorithms are preferred in practice due to their faster convergence (Goodfellow et al., 2016). In that case, instead of converging asymptotically to the desired solution, the iterates of constant stepsize SA have a stationary distribution. Although in many cases such weak convergence to a stationary distribution was established in the literature, it is not possible to fully characterize the limiting distribution. The reason is that, when constant stepsize is used, the distribution of the noise sequence within the SA algorithm plays an important role in the stationary distribution of the iterates. Since the distribution of the noise is in general unknown, the stationary distribution cannot be analytically characterized. In this work, building upon the works on stationary distribution of constant stepsize SA algorithms, we aim at understanding the limiting behavior of the properly scaled stationary distribution as the constant stepsize goes to zero.
More formally, with initialization X (α) 0 ∈ R d , consider the SA algorithm where F : R d → R d is a general nonlinear operator, α is the constant stepsize, and {w k } is the noise sequence. Observe that SA algorithm (1) can be viewed as an iterative algorithm for solving the equation F (x) = 0 in the presence of noise (Robbins and Monro, 1951). A typical example is when F (x) = −c∇f (x) (where c > 0 is a constant) for some objective function f (·), in this case Algorithm (1) becomes the popular stochastic gradient descent (SGD) algorithm for minimizing f (·) (Lan, 2020; Bottou, Curtis and A version of this work was submitted to a conference on 28th May, 2021 Equal contribution between Zaiwei Chen and Shancong Mou Nocedal, 2018). Another example lies in the context of reinforcement learning, where F (x) = T (x) − x, and T (·) is the Bellman operator (Sutton and Barto, 2018). In this case, Algorithm (1) captures popular reinforcement learning algorithms such as TD-learning (Sutton, 1988) and Q-learning (Watkins and Dayan, 1992). Under some mild conditions on the operator F (·), it was shown in the literature that the sequence {X (α) k } converges weakly to some random variable X (α) (Dieuleveut, Durmus and Bach, 2017;Bianchi, Hachem and Schechtman, 2020;Yu et al., 2020;Durmus et al., 2021). However, for a fixed α, it is not possible to fully characterize the distribution of X (α) because it depends on the distribution of the noise sequence {w k }, which is usually unknown. In this work, we further consider letting α go to zero, and study the distribution of a properly centered and scaled iterate. Specifically, let Y (α) k := (X (α) k − x * )/g (α), where x * is the solution of F (x) = 0, and g : R → R is a properly chosen scaling function 1 . When k goes to infinity, we expect that Y (α) k converges weakly to some random variable Y (α) . Then we let α go to zero, and our goal is to further characterize the weak limit of Y (α) . Notice that proper scaling of the iterates is essential for raveling its fine grade behavior because otherwise the limiting distribution of the un-scaled iterates will converge to a singleton as the stepsize α goes to zero, which is analogous to the almost sure convergence result for using diminishing stepsizes in SA literature.
To summarize, we want to find a suitable scaling function g(·) and characterize the following two-step weak convergence of the centered scaled iterate Y where we use the notation ⇒ for weak convergence (or convergence in distribution).

Main Contributions. Our main contributions are twofold.
Characterizing the Distribution of Y . We propose a general framework for characterizing the distribution of Y in the following 3 cases: (1) SGD with a smooth and strongly convex objective, (2) linear SA with a Hurwitz matrix, and (3) SA involving a contractive operator. In particular, we show that in all three cases above the correct scaling function is g(α) = √ α, and the distribution of Y is Gaussian with mean zero and covariance matrix being the unique solution of an appropriate Lyapunov equation. Our proof is to use the characteristic function as a test function to obtain an integral equation of the distribution of Y , and then show that the desired Gaussian distribution solves the integral equation.
Determining the Suitable Scaling Function. For more general SA algorithms, we show empirically that the scaling function need not be g(α) = √ α and the distribution of Y need not be Gaussian. Inspired by this observation, we propose a method to find the the correct scaling function for general SA algorithms. In particular, our results indicate that the scaling function g(α) should be chosen such that (1) lim α→0 1.2. An Illustrative Example. In this section, we provide an example to further illustrate the problem we are going to study. Consider SA algorithm (1). Suppose that F (x) = −x is a scalar valued function, and {w k } is a sequence of i.i.d. standard normal random variables. We make such noise assumption here only for ease of exposition, and it will be relaxed in later sections. In this case, the SA algorithm (1) becomes This algorithm has the following two simple interpretations: (1) it can be viewed as the SGD algorithm for minimizing a quadratic objective function f (x) = x 2 /2, which has a unique minimizer x * = 0, (2) it 1 The scaling function is unique up to a numerical factor can also be viewed as an SA algorithm for solving the fixed-point equation T (x) = x with T (x) being identically equal to zero, therefore x * = 0 is the unique fixed-point.
and V[ · ] denote the mean and variance, respectively, we have converges weakly to a random variable Y (α) , whose distribution is N (0, 1 2−α ). In this case, we are able to analytically characterize the distribution of Y (α) for a fixed α because of the simplicity of the algorithm (3) and the noise {w k } being i.i.d. Gaussian. For the general SA algorithm (1) with limited information on the noise sequence {w k }, it is in general not possible to fully characterize the distribution of Y (α) . Now consider the second weak convergence in Eq.
(2). Since we have already shown that Y (α) ∼ N (0, 1 2−α ). As α goes to zero, we see that Y (α) converges weakly to a random variable Y , whose distribution is N (0, 1 2 ). As opposed to the first weak convergence in Eq.
(2), where the distribution of Y (α) in general cannot be fully characterized, we are able to characterize (in later sections) the distribution of Y for the more general algorithm (1), and for more general noise assumptions. Intuitively, the reason is that as the constant stepsize decreases, the effect of the entire distribution of the noise {w k } on the distribution of Y α is weakened. In the limit, only the mean and the variance of w k play roles in determining the distribution of Y . This is analogous to central limit theorem type of results.
To summarize, we have shown in the special case of Algorithm (3) that the correct scaling function is g(α) = √ α, and the distribution of the limiting random variable Y is a Gaussian distribution with mean zero and variance 1/ √ 2. In Section 2, we extend this result to the more general algorithm (1) with weaker noise assumptions.
1.3. Related Literature. Since proposed in (Robbins and Monro, 1951), SA has been popular for solving large scale optimization in modern machine learning applications. Although require using diminishing stepsizes to achieve convergence (Hu et al., 2017;Xie, Wu and Ward, 2020;Mertikopoulos et al., 2020;Shamir and Zhang, 2013;Li and Orabona, 2019;Fehrman, Gess and Jentzen, 2020;Gower, Sebbouh and Loizou, 2021), constant stepsize is used in practice (Goodfellow et al., 2016). In contrast to the success in machine learning practice, there is little discussion about the stationary distribution of constant step size SGD (Dieuleveut, Durmus and Bach, 2017;Bianchi, Hachem and Schechtman, 2020;Yu et al., 2020;Durmus et al., 2021). Among them, Dieuleveut, Durmus and Bach (2017) bridges Markov chain theory and the constant step size SGD algorithm. They provided an explicit asymptotic expansion of the moments of the averaged SGD iterates. Bianchi, Hachem and Schechtman (2020) studied the asymptotic behavior of constant step size SGD for a nonconvex nonsmooth, locally Lipchitz objective function. It was shown that in a small step size regime, the interpolated trajectory of the algorithm converges in probability towards the solutions of the differential inclusionẋ = ∂F (x) and the invariant distribution of the corresponding Markov chain converges weakly to the set of invariant distribution of the differential inclusion. Yu et al. (2020) established an asymptotic normality result for the constant step size SGD algorithm for a non-convex and non-smooth objective function satisfying a dissipativity property. Durmus et al. (2021) first established non-asymptotic performance bounds under Lyapunov conditions and then proved that for any step size, the corresponding Markov chain admits a unique stationary distribution.
The work mentioned before established the stationary distribution for almost strongly convex and smooth functions for a fixed constant stepsize. Since the SGD iterates will converge to a singleton as the constant step size goes to zero, none of the previously mentioned literates can be applied to study the limiting behavior of SGD in this regime. To understand such behavior, we propose to study the properly centered and scaled iterates. Although not directly related, it shares a similar flavor when studying the limiting behavior of the stationary distribution of the stochastic gradient Langevin dynamics (SGLD) iterates as step size goes to zero.
Another set of related literature is on the diffusion approximation of SGD (Li, Tai and Weinan, 2017;Feng, Li and Liu, 2017;Yang, Hu and Li, 2021;Sirignano and Spiliopoulos, 2020;Latz, 2021). Authors aim to approximate the trajectory of SGD by a diffusion process which solves an SDE. Notice that they also study the scaled version of the diffusion limit of SGD. However, different from our approach, their scale is in temporal domain and cannot be applied in our research.
The Markov chain perspective of studying SGD iterates when step size goes to zero (Dieuleveut, Durmus and Bach, 2017) is related to the heavy traffic analysis in queuing theory (Eryilmaz and Srikant, 2012). It has been studied in the literature using fluid and diffusion limits (Gamarnik and Zeevi, 2006;Harrison, 1988Harrison, , 1998Harrison and López, 1999;Stolyar et al., 2004;Williams, 1998) where the interchange of limit is usually problematic (Eryilmaz and Srikant, 2012). An alternative approach in stochastic networks is based on drift arguments introduced by (Eryilmaz and Srikant, 2012) and further generalized by (Maguluri and Srikant, 2016;Maguluri, Burle and Srikant, 2018;Hurtado-Lange, Varma and Maguluri, 2020;Mou and Maguluri, 2020). We adopt similar techniques in quantifying the limiting distribution of the scaled SGD iterates. Notice that in stochastic networks, people mainly focus on finite state space Markov chains. However, when it comes to SGD iterates, the state space is continuous and thus more challenging.
The rest of this paper is organized as follows. In Section 2, we characterize the distribution of Y (cf. Eq (2)) in the following cases: (1) F (·) is the negative gradient for some smooth and strongly convex function is a contraction operator. In all three cases above, the scaling function is chosen to be g(α) = √ α. Then in Section 3, we first empirically show that for more general SA algorithms beyond these cases, the scaling function need not be g(α) = √ α and the distribution of Y need not be Gaussian. Then we present a method to determine the scaling function for more general SA algorithms and make connection to the Euler-Maruyama discretization scheme for approximating SDE. Finally, we conclude this paper in Section 4.
2. Characterizing the Asymptotic Stationary Distribution. Through out this section, we make the following assumption regarding the noise sequence {w k }.
ASSUMPTION 2.1. The sequence {w k } is independent and identically distributed with mean zero and a positive definite covariance matrix Σ ∈ R d×d .
Note that Assumption 2.1 is much weaker than the assumption used in Section 1.2, where the noise is assumed to obey the Gaussian distribution. Nevertheless, extending our results to the more general noise setting (e.g. martingale difference noise, Markovian noise, etc) is one of our future directions.

SGD for Minimizing a Smooth and Strongly
is an objective function. Then the SA algorithm becomes which is the well-known SGD algorithm for minimizing f (·). To proceed, we require the following definition.
DEFINITION 2.1. A convex differentiable function h : R d → R is said to be L-smooth and σ-strongly convex with respect to the Euclidean norm · 2 if and only if To characterize the asymptotic behavior of Algorithm (4), we make the following assumption.
ASSUMPTION 2.2. The function f : R d → R is twice differentiable, and is L-smooth and σ-strongly convex.
Under Assumption 2.2, the function f (x) has a unique minimizer (or F (x) = 0 has a unique solution), which we have denoted by x * . To proceed, let Y which is obtained by first subtracting both sides of Eq. (4) by x * , and then dividing by √ α. We next characterize the two-step weak convergence (cf. Eq. (2)) in the following theorem. Let H f ∈ R d×d be the Hessian matrix of the objective function f (·) evaluated at x * , which is well-defined because f (·) is twice differentiable (cf. Assumption 2.2). THEOREM 2.1. Consider the iterates {Y (α) k } generated by Algorithm (5). Suppose that Assumptions 2.1 and 2.2 are satisfied, then the following statements hold.
(1) There exists a thresholdᾱ > 0 such that for all α ∈ (0,ᾱ), the sequence of random variables {Y converges weakly to a random variable Y , which satisfies the following integral equation for all t ∈ R d . In addition, suppose that Eq. (6) has a unique solution (in terms of the distribution of Y ), then the distribution of Y is a Multivariate normal distribution with mean zero and covariance matrix Σ Y being the unique solution of the Lyapunov equation REMARK. To establish Theorem 2.1 (2), we require Eq. (6) to have a unique solution in terms of the distribution of Y . Such assumption will be relaxed to some extent in Section 2.4.
Since Σ is positive definite (cf. Assumption 2.1), and H f is also positive definite under Assumption 2.2, it is well established in the literature that the Lyapunov equation H f Σ Y + Σ Y H f = Σ has a unique solution, which is explicitly given by Consider the special case where f (x) = x 2 /2. In this case we have H f = 1, and hence Σ Y = 1 2 Σ by the Lyapunov equation. As a result, the distribution of the limiting random variable Y is Gaussian with mean zero, and covariance matrix being 1 2 Σ. This agrees with the illustrative example (which is for scalar valued iterates) presented in Section 1.2.
From Theorem 2.1, we see that the distribution of Y only depends on the Hessian of f (·) at x * . This makes intuitive sense because we are studying the asymptotic behavior of SA algorithm (4), and only the properties of f (·) around x * should play a role in determining the stationary distribution.

Stochastic Approximation for Solving Linear Systems of Equations
which aims at iteratively solving the linear equation Ax + b = 0. Note that since the matrix A is not necessarily symmetric, F (x) need not be a gradient of any objective function. Such linear SA algorithm arises in many realistic applications. One typical example is the TD-learning algorithm for solving the policy evaluation problem in reinforcement learning, where the goal is to solve a linear Bellman equation. See Bertsekas and Tsitsiklis (1996); Srikant and Ying (2019) for more details.
To study the asymptotic behavior of Algorithm (7), we make the following assumption regarding the matrix A.
ASSUMPTION 2.3. The matrix A is Hurwitz, i.e., all the eigenvalues of the matrix A have strictly negative real part.
REMARK. Since A being Hurwitz implies A being non-singular, Assumption 2.3 implies that the target equation Ax + b = 0 has a unique solution, which we denote by x * . Assumption 2.3 is standard in studying linear SA algorithms. In particular, it was shown in the literature that under Assumption 2.3 and some mild conditions on the noise {w k }, Algorithm (7) converges in the mean square sense to a neighborhood around x * (Bertsekas and Tsitsiklis, 1996).
To study the asymptotic distribution, for a fixed stepsize α, we define the centered scaled iterate Y The full characterization of the two-step weak convergence (cf. Eq. (2) (1) There exists a thresholdᾱ > 0 such that for all α ∈ (0,ᾱ ), the sequence of random variables {Y (2) For any positive sequence {α i } satisfying α i ∈ (0,ᾱ ) for all i and lim i→∞ α i = 0, the sequence of random variables {Y (αi) } converges weakly to a random variable Y , which satisfies the following equation In addition, suppose that Eq. (9) has a unique solution, then the distribution of Y is a Multivariate normal distribution with mean zero and covariance matrix being the unique solution of the Lyapunov equation AΣ Y + Σ Y A + Σ = 0.
Since the matrix A is Hurwitz, and Σ is positive definite, the existence and uniqueness of a positive definition solution to the Lyapunov equation AΣ Y + Σ Y A + Σ = 0 are guaranteed (Haddad and Chellaboina, 2011). Lyapunov equations were used extensively in studying the stability of linear ordinary differential equations (ODE). Interestingly, it also shows up in determining the limit distribution of centered scaled iterates of discrete linear SA algorithms.

Stochastic Approximation under Contraction
which can be interpreted as an SA algorithm for finding the fixed-point of the operator T (·). These type of algorithms arise in the context of reinforcement learning, where T (·) is the Bellman operator. To proceed, we need the following definition.
Next, we state our assumption regarding the operator T (·).
Assumption 2.4 essentially states that the operator T (·) is a contraction mapping with respect to the weighted 2 -norm · µ . By Banach fixed-point theorem (Banach, 1922), the operator T (·) has a unique fixed-point, which we denote by x * .
We next write down the update equation of the centered scaled iterate Y In the next theorem, we characterize the distribution of the limiting random vector Y (cf. Eq. (2)). Let J ∈ R d×d be the Jacobian matrix of T (·) evaluated at x * . THEOREM 2.3. Consider the iterates {Y (α) k } generated by Algorithm (11). Suppose that Assumptions 2.1 and 2.4 are satisfied, then the following statements hold.
(2) For any positive sequence {α i } satisfying α i ∈ (0,ᾱ ) for all i and lim i→∞ α i = 0, the sequence of random variables {Y (αi) } converges weakly to a random variable Y , which satisfies the following equation In addition, suppose that Eq. (12) (12)), to conclude that the distribution of Y is Gaussian, we need to assume that the equation has a unique solution. In this section, we show that such uniqueness assumption can be relaxed to some extend.
where the interchange of integral and differentiation is justified (Flanders, 1973). Now Eq. (13) is an ODE, which has solutions of the form where C is a constant. Since φ Y (t) as a characteristic function satisfies φ Y (0) = 1, we have C = 1 and hence φ Y (t) = exp(− a 4b t 2 ), which is characteristic function for a Gaussian random variable with mean zero and covariance a/(2b). Therefore, the uniqueness assumption can be removed in the uni-dimensional setting.

Multi-Dimensional Setting.
Moving to the multi-dimensional setting, consider Eq. (6) of Theorem 2.1 as a representative example. To show the same result of Theorem 2.1 (2) without imposing the uniqueness assumption, we consider the case where (1) the Hessian matrix H f of the objective function f (·) evaluated at x * is the identity matrix, and (2) the covariance matrix of the noise w k is also an identity matrix. Extending the result to the more general setting where H f and Σ can be any positive definite matrices is a future research direction.
Similarly let φ Y (t) = E[e it Y ] be the characteristic function of the random vector Y . Then in this case Eq. (6) becomes t tφ Y (t) + 2t ∇φ Y (t) = 0, which is equivalent to where ψ Y (t) := log(φ Y (t)). To solve the partial differential equation (PDE) (14), we will first convert the PDE from Cartesian coordinates to spherical coordinates, which then is directly solvable.
2.5. Proof of Theorem 2.1. In this section, we present our proof for Theorem 2.1. The proofs for Theorems 2.2 and 2.3 are similar and hence are omitted.
Before going into details, we first highlight the main ideas for the proof. For Theorem 2.1 (1), to show that {Y (αi) } converges weakly to a some random variable Y , we show that any sub-sequence of {Y (αi) } further contains a weakly convergent sub-sequence, with a common weak limit. We do it by first showing that the sequence of random variables {Y (αi) } is tight. In particular, we show that there exists a constant B independent of α such that E[ Y (α) 2 2 ] ≤ B for all small enough α. This result in conjunction with the Markov inequality implies the tightness of {Y (αi) }. As a result of tightness, {Y (αi) } contains a weakly convergent sub-sequence. Now consider Theorem 2.1 (2). For any positive sequence {α i } such that lim i→∞ α i = 0, since the family of random variables {Y (αi) } is tight, there is a weakly convergent subsequence {Y αi k }. We further show that the weak limit Y of the subsequence {Y (αi k ) } solves Eq. (6). In this case, under the assumption that Eq. (6) has a unique solution, the random variable Y is a Gaussian random variable with mean zero, and covariance matrix Σ Y being the unique solution of the Lyapunov equation H f Σ Y + Σ Y H f = Σ. Since for every sequence {Y (αi) }, there is a weakly convergent subsequence {Y (αi k ) } with a common weak limit, the sequence of random variables {Y (αi) } also converge weakly to the same random variable Y .
2.5.1. Proof of Theorem 2.1 (1). To prove the result, we will apply Proposition 2.1 in Yu et al. (2020). For completeness, we first state this proposition using our notation in the following.
(a) Since the objective function f (·) is assumed to be L-smooth, we have for any (b) Since the objective function is assumed to be σ-strongly convex, we have for any x ∈ R d : (c) This is immediately implied by Assumption 2.1, with L = Trace(Σ) 1/2 . Now apply Proposition 2.1, when the stepsize α satisfies α < σ 2(3L 2 +Trace(Σ)) , the SGD iterates {X (α) k } converge weakly to some random variable X (α) , which is distributed according to the unique stationary distribution π α . In addition, we have k } converges weakly to some random variable Y (α) and E[ Y (α) 2 2 ] < ∞.
2.5.2. Proof of Theorem 2.1 (2). Following the road map described in the beginning of this section, we present and prove a sequence of lemmas in the following. Together they imply the desired result.
LEMMA 2.1. Let α 0 = σ/L 2 . the family of random variables {Y (α) } 0<α≤α0 is tight. PROOF OF LEMMA 2.1. We first show that that there exists an absolute constant C > 0 such that E[ Y (α) 2 ] ≤ C for any α ∈ (0, α 0 ]. Using the update equation (1) The existence and uniqueness of a stationary distribution Y (α) is proved in Part (1) of this theorem. We next show that the family of random variables {Y (α) } 0≤α≤α0 is tight. Using the equation and we have By smoothness, we have By strong convexity, we have . When α ∈ (0, α 0 ], we have from the previous inequality that Hence, for any α > 0, let M = Trace(Σ)/σα, then we have for any α ∈ (0, α 0 ]. It follows that the family of random variables {Y (α) } 0<α≤α0 is tight.
LEMMA 2.2. Let {α i } be a positive sequence of real numbers such that lim i→∞ α i = 0. Suppose that {Y αi } converges weakly to some random variable Y . Then Y satisfies the following equation PROOF OF LEMMA 2.2. For any i ≥ 0, we have which implies for any t ∈ R d : Using the Taylor's theorem and we have Using Theorem 3.3.20 from Durrett (2019) and we have Using the previous two inequalities in Eq. (18) and we have Simplify the above equality and we obtain We next let i go to infinity on both sides of the previous inequality and evaluate the limit of the terms Since {Y (αi) } converges weakly to some random variable Y , we have by continuity theorem (Theorem 3.3.17 in Durrett (2019)) that For the term T 6 , we have by bounded convergence theorem that lim αi→0 T 6 = 0. To evaluate the terms T 2 to T 5 , the following definition and result from Van der Vaart (2000) (2000)). Let f : R d → R be measurable and continuous at every point in a set C. Let X n ⇒ X, where X takes its values in C.
if and only if the sequence of random variables f (X n ) is asymptotically uniformly integrable.

Now consider the term T
which goes to zero as M → ∞, we have by Theorem 2.4 that Using the same line of analysis, we have lim i→∞ T 3 = lim i→∞ T 4 = lim i→∞ T 5 = 0. It follows that Rearranging terms and we obtain the resulting equation.
LEMMA 2.3. Suppose Eq. (17) admits a unique solution. Then the random variable Y given in Lemma 2.2 follows a Gaussian distribution with mean zero, and covariance matrix Σ Y being the unique solution of the Lyapunov equation PROOF OF LEMMA 2.3. Suppose that Eq. (6) has a unique solution, we only need to verify that the multinormal distribution with mean zero and covariance matrix being the solution to the Lyapunov equation (6).
3. Identifying the Suitable Scaling Function for More General SA Algorithms. In the previous section, we have shown that for several particular SA algorithms (e.g. SGD, linear SA, and contractive SA), the scaling function is g(α) = √ α and distribution of the limiting random variable Y is Gaussian. In this section, we consider more general SA algorithms. We first show impirically in the following section that in general the scaling function need not be g(α) = √ α, and the distribution of Y need not be Gaussian. (1) is the SGD algorithm for minimizing the scalar objective f (x) = x 4 /4. That is:

Numerical Experiments. Suppose that Algorithm
Note that f (·) in this case is neither smooth nor strongly convex. It is clear that the unique minimizer of f (·) is zero. Let the centered scaled iterate Y k /g(α). We next use numerical simulation to show that the correct scaling function in this case is g(α) = α 1/4 instead of g(α) = √ α. In Figures 1 and 2, we plot the empirical density function of Y (α) for different α. For the right scaling function, we expect the density function to converge as α decreases, while for the wrong scaling function, we expect the density function to change drastically for order-wise different α. As we see, it is clear that g(α) = √ α is not suitable in this case, and g(α) = α 1/4 seems to be the right scaling. To further verify this result, we plot the logarithmic empirical density function as a function of y 4 in Figure 3. We observe linear growth in Figure 3. This indicates that the density function p Y (y) is proportional to e βy 4 , where β is some numerical constant. Therefore, numerical experiments suggest that the distribution of Y is not Gaussian but super Gaussian in this problem.

3.2.
A Method to Determine the Suitable Scaling Function. Inspired by the numerical simulations provided in the previous section, we here provide a method to determine the correct scaling function for general SA algorithms.
To gain intuition, we consider the centered scaled iterates Y Notably, the factor in terms of the stepsize α in front of the term (Y (α) k ) 3 is α 3/2 , which is equal to the square of the factor α 3/4 in front of the noise term w k . Now for general SA algorithm (1), by rewriting the update equation (1) in terms of the centered scaled iterate Y In view of the previous equation and the empirical observations in the previous section, we see that we need to choose a scaling function g(α) such that the following condition is satisfied.
We next verify the choice of scaling functions in Section 2 using our proposed Condition 3.1. For SGD with a smooth and strong convex objective, since In view of the previous inequality and Condition 3.1, it is clear that the only possible choice of g(α) is g(α) = √ α. For linear SA algorithm studied in Section 2.2, one can also easily show using Condition 3.1 that g(α) = √ α. As for contractive SA studied in Section 2.3, using the contraction property and we have Since all norms are "equivalent" in finite dimensional spaces, the previous inequality implies that we must choose g(α) = √ α. Now to further verify the correctness of the scaling function suggested by Condition 3.1, consider the SGD algorithm with the following two choices of objective functions: (1) f (x) = e x 2 , and (2) f (x) = x 4 4 + sin 2 (x) 2 . Note that in these two cases the function f (·) is not smooth and strongly convex.
Case 1. In the first case where f (x) = e x 2 , since when choosing g(α) = √ α, we haveF (y) = lim α→0 g(α) 2 α 2ye (yg(α)) 2 = 2y. One interesting implication of this example is the following. In the three SA algorithms considered in Section 2, it seems that it is the function F (·) that appears in the algorithm determines the distribution of Y . However, the above example suggests that it is the functionF (·) of Condition 3.1 instead of F (·) that directly impacts the distribution Y . In SGD with a smooth and strongly convex objective, linear SA, and contractive SA, F (·) andF (·) happen to coincides, but this is in general not the case. In fact, since we have x 2i−1 by Taylor series,F (·) in this example is exactly the dominant terms appears in the Taylor series. In addition, this suggests that the distribution of the limiting random variable Y has a density function proportional to e β x 2 , where β is a numerical constant.
We next verify this choice of g(α) and the distribution of Y (α) for small enough α using numerical simulation in the following. We see from Figure 4 that with the scaling function g(α) = √ α, the empirical density function of the random variable Y (α) seems to converge. and Figure 5 further justifies this result by showing that the density function p Y (y) of the distribution of Y in this case is proportional to e β x 2 , where β is a numerical constant.
Since lim x→0 sin(x) x = 1, the only possible choice of the scaling function g(α) to satisfy Condition 3.1 (2) is g(α) = √ α. In this case, we haveF (y) = lim α→0 1 √ α y 3 α 3/2 + sin(y √ α)cos(y √ α) = y by L'Hôpital's rule. This is another example where F (·) =F (·). In fact, since x 4 is dominated by sin 2 (x) as x approaches x * (which is 0), the scaling function and the functionF (·) are determined only by the dominant term. . Similarly, we verify this choice of scaling function via numerical experiments. In Figures 6 and 7, we plot the empirical density function of the random variable Y (α) for different stepsize α, and see if the density function converges as α goes to zero. The results suggest that g(α) = α 1/2 seems to be the correct scaling. To further verify this result, we plot the logarithmic function of the empirical density of Y (α) as a function of y 2 and observe straight lines. Therefore, the distribution of Y (α) is proportional to e β x 2 , where β is a numerical constants.
3.3. Connection to Euler-Maruyama Discretization Scheme for Approximating SDE. The choice of the scaling function suggested by Condition 3.1 has an insightful connection to the Euler-Maruyama discretization scheme for approximate the solution of an SDE, as elaborated below. Let (B t ) t≥0 be a Brownian motion. Consider the following SDE: with initial condition X 0 . The Euler-Maruyama discretization {X k } to the solution (X t ) of SDE (21) is defined as follows. Let ∆t be the discretization accuracy. SetX 0 = X 0 , and recursively defineX k according toX k+1 =X k + ∆tF (X k ) + (B (k+1)∆t − B k∆t ).
Since (B t ) t≥0 is a Brownian motion, we have (B (k+1)∆t − B k∆t ) ∼ N (0, ∆t). Therefore, by letting {Z k } be an i.i.d. sequence of standard normal random variables, we can rewrite the previous equation aŝ The approximation property of the Euler-Maruyama discretization (22) to its corresponding SDE (21) has been studied in the literature, see Wenlong et al. (2019). Specifically, it was shown that under some mild conditions on F (·), the Euler-Maruyama scheme is known to have the first-order accuracy of the SDE (21). As a consequence, intuitively, when (X t ) t≥0 has a stationary distribution µ, the limiting distribution µ ∆t of {X k } as a function of the discretization accuracy ∆t should converge weakly to µ as ∆t tends to zero. If we view the discretization accuracy ∆t as the stepsize in Eq. (22). In order for µ ∆t to converge to some nontrivial distribution µ as ∆t tends to zero, it is important to notice that the scaling factor of the noise Z k in terms of ∆t must be order-wise equal to the square root of the scaling factor of F (X k ). This observation coincides with Eq. (20) in the previous section, which eventually leads to our Condition 3.1.

Conclusion.
In this paper, we characterize the asymptotic stationary distribution of properly centered scaled iterate of SA algorithms. In particular, we show that for (1) SGD with smooth and strongly convex objective, (2) linear SA, and (3) contractive SA, the scaling function is g(α) = √ α and the corresponding stationary distribution is Gaussian. For SA beyond these cases, we empirical show that the stationary distribution need not be Gaussian, and provide a method for determine the suitable scaling function. Since our paper is the first study for this problem, it might open a door for research in this direction.