A Method for Fault Detection in Wireless Sensor Network Based on Pearson’s Correlation Coefficient and Support Vector Machine Classification

Sensor nodes are tiny low-cost devices prone to various faults. So, it is imperative to detect those faults. This paper presents a sensor measurement fault detection algorithm based on Pearson’s correlation coefficient and the Support Vector Machine(SVM) algorithm. As environmental phenomena are spatially and temporally correlated but faults are somewhat uncorrelated, Pearson’s correlation coefficient is used to measure correlation. Then SVM was used to classify faulty readings from normal readings. After classification, faulty readings are discarded. Here each sensor nodes periodically collects environmental features and sends them to their associated cluster heads. Each cluster head analyze collected data using the classification algorithm to detect whether any fault is present or not. Network simulator NS-2.35 and Matlab are used for evaluation of our proposed method. The fault detection algorithm was evaluated using performance metrics, namely, Accuracy, Precision, Sensitivity, Specificity, Recall, F1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F_1$$\end{document} Score, Geometric Mean(G_mean), Receiver Operating Characteristics (ROC), and Area Under Curve(AUC). Performance evaluation shows, the proposed method performs well for high fault percentages.

Introduction
Wireless Sensor Network(WSN) consists of a large number of sensor nodes scattered over a large spacial region with one or more base stations(BS). The basic functionality of the WSN node is to collect information and send it to BS for analysis [10,19]. These sensor nodes are tiny, low-cost devices equipped with one or more sensors. Sensor nodes are prone to various faults [8]. There are many factors involved that can cause a fault, such as 1 3 hardware failure, software fault, communication error, etc. [21] These faults may persist for a long time or may be instantaneous. Since faults are unavoidable, discovering faulty and fault-free nodes is crucial in the field of WSN [4]. The fault detection algorithm needs to be designed in such a way that it can detect such inconsistent behavior that can disrupt the normal functionality of the WSN. Sensors are deployed to monitor the region, and if any suspicious event occurs, report the same to the base station. Moreover, these cheap nodes are prone to various faults, including measurement faults. This implies that sensor measurements are not reliable. So the event can't be concluded based on a single measurement instance. Measurements need to be correlated spatially and temporally to decide whether a measurement is faulty, normal, or event. It is shown that environmental phenomena are spatially and temporally correlated [1]. In this work, this property was exploited to detect faulty reading. The Person's correlation coefficient [20] between two-time slots was utilized to detect the presence of a fault. Then SVM [2] was used to classify between faulty reading and normal reading.
The contribution in this work is as follows, • Pearson's correlation coefficient ( ) was used to measure the correlation between twotime slots. is invariant of scaling and self normalizing, which makes it suitable for spatially correlated environmental features. • SVM was used to classify faulty and normal reading. SVM performs well in higher dimensional feature space. Moreover, using the kernel trick of SVM, nonlinear data can also be solved. As the used dataset is nonlinear and has two features, namely humidity, and temperature, SVM is suitable for the proposed method. • NS2 and Matlab are used to simulate and analyze our proposed method, respectively.
The proposed method was also tested over a real-world dataset.
The rest of the paper is organized as follows: related works are discussed in Sect. 2. Problem definition is given in Sect. 3. Preliminaries are described in Sect. 4. Our proposed method is presented in Sect. 5. Simulation results and performance evaluation are shown in Sect. 6. The conclusion is drawn in Sect. 7.

Related Works
In [1], authors proposed General anomaly detection(GAD), a fully distributed scheme for practical large scale networked industrial sensing systems(NISSs). Real-time detection, Distributed solution, and General solution are three properties of GAD. Their method, distributed matching-based grouping algorithm (DMGA), divides all sensing components into small, strongly correlated groups in a fully distributed way. Spatial correlation is a natural property in various physical phenomena. "Since physical phenomena are continuous, these spatial correlations should be temporally correlated to previous mappings." They assumed measurement errors follow Gaussian distributions. They evaluated their method using a successful detection rate (SDR) and false-positive detection rate (FPDR).
In [21], authors used SVM with Gaussian kernel to detect faults in WSN. They classified Data faults as Offset fault, Gain fault, Stuck-at fault, Out of bounds. In their proposed method, received sensor data is classified using SVM with Gaussian kernel function. Their learning phase is performed in BS. Then, the decision function was transmitted to each cluster head to classify new measured data. In their proposed method, every time a new data is measured, an observation vector composed of last three data measurement of two sensors ( V t , V t−1 , V t−2 ) is constructed by the data preparation block although SVM is capable of multidimensional classification. After that, this observation vector is classified by SVM using the decision function. If SVM output is positive new data is normal else data is faulty. Their labeled dataset is based on an existing dataset published by the researchers of the University of North Carolina at Greensboro [12].
A fault diagnosis protocol based on gradient descent and evolutionary approach was proposed by Swain et al. in [13]. Their method detects the faulty nodes and isolates them from the network. Their proposed protocol comprises four phases, namely clustering phase, communication phase, fault detection and classification phase, and isolation phase. They used a genetic algorithm and neural network for fault detection. Their method classifies faults into four types according to the rate of faults. Then faults are isolated. Their comparison showed that evolutionary approach outperforms gradient descent for their sensor data.
In [14], the authors proposed a fault diagnosis protocol that detects heterogeneous faults in WSN. Their protocol is capable of detecting soft permanent, hard permanent, transient, and intermittent faults, as claimed by authors. In their method, hard permanent faults are identified by a time out status register mechanism. Soft permanent, intermittent, and transient faults are detected by a statistical test, namely analysis of variance (ANOVA) test. In this phase, m number of measurements of n nodes in a cluster was used for the ANOVA test. This test was repeated for r times to detect faults. They also used the feed-forward probabilistic neural network (PNN) to classify these heterogeneous faults. However, as the protocol relies on the coordinator node, any inconsistent behavior such as node failure or erroneous results of the coordinate node will lead to the degradation of the performance of the protocol.
In [11], Saeedi et al. proposed a density-based spatial clustering of applications with noise (DBSCAN) algorithm for detecting anomalies. Their proposed method used three features from eight features of the Intel Berkeley Research lab(IRLB) dataset. A significant portion of the data is normal, and they are stored in CH with high Density. Authors utilized this idea in their work. DBSCAN algorithm identifies low-density regions and labels them as an anomaly. Then they trained the SVM classifier using these labeled data.
In [5], authors objective was to design a fault detection and diagnosis system using existing ML techniques. They used three machine-learning algorithms, namely, Fuzzy Deep Neural Network (FDNN), Support vector machine (SVM), and Neural Network(NN), for fault detection in their work.
In [15], authors proposed a complete fault diagnosis methodology to detect faulty sensors in WSN. Their proposed method contains four phases, namely, initialization phase, fault detection phase, fault classification phase, and fault tolerance phase. They used checksum and Fletcher's checksum method to detect hard fault and link fault. They also used Mann-Whitney U statistical test for soft fault detection. Then they utilized Gaussian transformation function to classify the soft faults. They also used a stepwise regressional method to tolerate fault in their method. The Detection phase is performed on each cluster head with its own and member sensor measurements. In Mann-Whitney U statistical test, P-value decides the status of the sensor measurement value. They used a threshold value based on application and situation of the sensor network. If P-value is less than the threshold, their method declares it as soft fault.
In [17], the authors proposed DODS (Distributed Outlier Detection Scheme), where outliers are detected locally by each node. They used four data types, i.e., temperature, voltage, humidity, and light. The main idea is to clean sensed data (measurements) from outlier (incorrect data). The scheme operates in nodes that made the sensing operation and does not require any neighbor's communication. The solution exploits the temporal correlations existing in the sensed data (current and history sensed data) of the same node and its remaining energy level. Outlier detection is performed using Bayes' classifier for each type of data. Only nodes belong to an interesting region (IR) participate in the outlier detection process. To learn the prior probability and to compute all conditional probability, they used a supervised off-line method. They considered a different set of classes (small, medium, large) for different data types and used the maximum a posteriori (MAP) concept in order to determine optimal class. But their method does not differentiate between faults and events. No spatial correlation is considered in their method.

Problem Description
Let, N number of nodes are deployed in a region. Each node is equipped with k sensors to sense k(k > 1 ) number of environmental features. These nodes periodically sense environmental features and send them to their associated Cluster Head (CH). Upon receiving sensed data, CH starts analyzing them to detect if there is any fault.
Suppose, sensor nodes periodically sense environmental data with time interval T. Then total time is divide into time slots t (t = 0, 1, 2 … ) depending on T as illustrated in Fig. 1.
where m is member nodes of the cluster head. As sensor nodes are faulty, this Z t is composed of both faulty and normal reading. CH detects these faults and eliminates them. In this work, CH computes Pearson's correlation coefficient( ) between Z t and Z t−1 i.e., (Z t−1 , Z t ) . Then, CH classifies (Z t−1 , Z t ) as faulty or normal using Support Vector Machine Classifier.

Preliminaries
Spatiotemporal correlation is the nature of various physical phenomena such as temperature, humidity, illumination, and many more. Spatial correlation implies a correlation mapping of measures between two neighboring nodes at time t. As physical phenomena are continuous, there should also be a mapping between current measure (t) and previous measure ( t − 1 ). This is called temporal correlation. Several possible correlation measures are there, among which Pearson's correlation coefficient is most popular.

Definition 1
The (Pearson) correlation coefficient between two random variable X and Y is defined as where cov[X, Y] is the covariance of X and Y. Whereas var[X] is the variance of random variable X. The correlation coefficient can be viewed as a degree of linearity between X and Y. [20]

Support Vector Machine
Let a training sample set of length of length k is given with two separable classes P and N: where y k ∈ {1, −1} labels X k to belong to either of the two classes. It is required to find a hyper-plane in terms of weight vector(w) and bias term(b), that separates the two classes.
In case of linear classification, SVM classifier computes a decision hyper-plane i.e.
This implies for both i ∈ P and i ∈ N has to satisfy (3).
From all possible hyperplanes, the objective is to find out the optimal hyperplane that satisfies the above condition. It is needed need to place the optimal hyperplane in a position such that distance from the hyperplane to the closest point of either side is the same. Now the problem of finding the optimal decision plane in terms of and b can be formulated as: Its solution gives us the optimal margin classifier. The solution can be obtained using the Lagrange multipliers method, as shown in (4). where i is Lagrange multiplier and i is support vector.

Definition 2
The vectors i residing on either side of separation hyperplane for which holds are called support vector(sv). SVM only depends on the support vectors. Other sample vectors are not important. [2] Substituting w from (3) we get, For the optimal weight vector and optimal bias b, we have:

SVM for nonlinear classification/Kernel Mapping
The method mentioned above is called linear SVM, which converges in case of linearly separable data. However, using kernel trick, SVM also works with nonlinear data set. With the help of Kernel trick, a sample x is mapped into a higher dimensional feature space where sample x is linearly separable as: Moreover, the decision function can be rewritten for the new space as: where and b are the parameters of the decision plane in the new space.
Furthermore, the classification function in new space becomes: From (5) and (6), we can see that vector x j appears only in inner products of both decision function and learning law. This implies that it is not needed to explicitly specify mapping function (X) . Just need the inner product of the vectors in new space .

Definition 3
In Machine Learning, a kernel refers to kernel trick, which is used to solve a non-linear problem using a linear classifier. A kernel function takes i and j vectors as arguments and returns the inner product of their images ( i ) and ( j ) : [9] A kernel function only returns the inner product of two vectors. So the dimension of kernel space is not so important. Kernel K( 1 , 2 ) needs to be positive semidefinite to fulfill the criteria of Reproducing Kernel Hilbert Space(RKHS) where optimization problem has a finite-dimensional solution that converges to an optimal one. Now by replacing T with ( ) T of (3), the new separation hyperplane in kernel space is obtained as shown in (7). We used this equation for classification of non linear problem.
The bias term (b) can be computed from any of the support vectors x i as shown in (8).

Proposed Scheme
In this work, Pearson's correlation coefficient( ) and Support Vector Machine(SVM) was applied to determine faults. Pearson's is suitable for fault detection because of some of its intrinsic properties such as is invariant to scaling i.e. (x, y) = (x, ax + b) , where a and b are constants. This implies normal reading, as well as event reading, will show a high correlation, whereas the presence of faulty reading will show a lesser correlation [7]. SVM is capable of classifying higher dimensional data, and using kernel trick SVM also can classify nonlinear data. As our data set in nonlinear and has two features, SVM classifier was used in this work. Figure 2 shows the learning process(training and Classification).

Training Phase
Support Vector Machine(SVM) was used to classify between faulty reading and normal reading. SVM is trained at Base Station with labeled data ( Z i,t , y ), where y ∈ {−1, 1} . Here label y = 1 represents not faulty and y = −1 represents faulty class.
Our labeled data Z i,t is composed of humidity and temperature. At first, Z i,t is separated in two vectors H Z i,t and T Z i,t for humidity and temperature respectively. Then two new vector is constructed as H Z t = ( H Z 1,t , H Z 2,t , … , H Z m,t ) and T Z t = ( T Z 1,t , T Z 2,t , … , T Z m,t ) . These are done in data preparation phase.
Then Pearson's correlation coefficient( ) is computed between two time slot t − 1 and t for all H Z t and T Z t as H and T respectively. After that vector is set as = [ H , T ] . Then SVM is trained with ( , y) using procedure given in Section (4.1). Then w and b of classification hyperplane is send to all cluster head nodes.

Fig. 2 Learning Process
In this method, Gaussian kernel function(K) was used as (10).
where x and x j are two vectors and is a free parameter.

Classification Phase
At each time slot(t), every sensor node collects environmental phenomena(e.g., humidity and temperature) and sends them to its associated cluster-head. Cluster head(CH) collects the reading of the neighboring node. CH computes the correlation between time-slot t and t − 1 of each feature of all the neighboring nodes.
At each time slot, Cluster Head (CH) stores measurements of its member nodes in a vector Z t . Then CH computes correlation coefficient between Z t and Z t−1 using (1). After that, decision function (5) is applied to this . If belongs to positive class then Z t is not faulty otherwise Z t is faulty.
The fault detection algorithm was presented in Algorithm 1. As described in the algorithm, like any other supervised learning, this method is also divided into two phases. First, the training phase, which is performed at BS using training data. The training phase produces a decision function in terms of w and b. Second, Classification phase, which is performed at each CH using the decision function.

Simulation Scenario and Performance Evaluation
To evaluate the proposed method, NS2.35 [3] and Matlab is used. NS2 is used to simulate the network scenario and generate measurement data. Then, this generated data is analyzed using Matlab. Performance evaluation is also done using Matlab.
The simulation data is generated using NS2.35. Here Gaussian distribution(N) with mean( ) and variance( ) is used of normal and event reading of the dataset [12]. This dataset contains only measurements of 4 sensor nodes. However, actual WSN consists of hundreds of sensor nodes. In the simulation for simplicity, only one cluster head with 100 member nodes is simulated. Results for actual dataset by computing temporal correlation also shown in Section 6.4.
In the network scenario, 100 sensors are randomly deployed in a 300m × 300m with CH at the center. The transmission range of each node is 60m.
Normal sensor readings for temperature are drawn from N( 1t , 2 1t ) and event readings from N( 2t , 2 2t ) . Where 1t = 28.1273, 2t = 29.3112 and 1t = 1.0952, 2t = 4.5588. Normal sensor readings for humidity are drawn from N( 1h , 2 1h ) and event readings from N( 2h , 2 2h ) . Where 1h = 59.6504, 2h = 78.4943 and 1h = 9.7391, 2h = 11.3831. All the faulty readings for temperature is also drawn from N( 2t , 2 2t ) . Where 2t = 29.3112 and 2t = 4.5588 . Moreover all the faulty readings for humidity is also drawn from N( 2h , 2 2h ) . Where 2h = 78.4943 and 2h = 11.3831. In the simulation, for analyzing the proposed method, the first 20% data are drawn from normal reading, and the last 20% data are drawn from event reading. Rest is a mix of normal and fault reading. Fault readings are mixed according to fault percentage. For example, in the case of 10% fault, 10% data is faulty, and the rest 90% data is normal. We run our simulation for 4690 time slots for 100 nodes for fault percentages 50%, 40%, 30%, 20%, 10%, and 5%. So for each fault percentage length of total data is 469000, 60% of which is used for training purpose and the rest 40% is used for testing. Simulation parameters are given in Table 1.

Computational Complexity Analysis
In this work, Pearson's correlation coefficient and SVM are used. The computational complexity of SVM training and testing is O(n 3 ) and O(n) [18]. The complexity of Pearson's correlation coefficient is O(n) . So in this work, computational complexity for training phase is O(n 3 ) , and for testing phase, it is O(n).

Confusion Matrix
In the case of binary classification, only four possible outcomes may occur. They are, Positive samples tested as positive (TP), Negative sample tested as negative (TN), Positive sample tested as negative (FN) and, Negative sample tested as positive (FP). In binary classification, data samples labeled as positive or negative can be predicted as positive or negative or vice versa. So, only four possible Test outcomes can occur, and they can be represented in the form of a confusion matrix [16], as shown in Table 2.
The proposed method is evaluated using performance matrices namely accuracy, precision, recall, F 1 Score, Gmean, sensitivity, and specificity [6]. These matrices are defined as follows.
Accuracy Measure: "Accuracy is ratio of observations that are correctly predicted to total observations."   Gmean: "Geometric mean (Gmean) is the square root of true positive and true negative." Sensitivity: "Sensitivity is the ratio between all the positive samples to true positive(TP) samples." Specificity: "Specificity is the ratio between actual negative samples to true negative(TN) samples." Accuracy, Precision, Sensitivity, Specificity, Recall, F 1 , Score, G_mean for faults 50%, 40%, 30%, 20%, 10%, and 5% is shown in Table 3. From the results presented in Table 3, it can be seen clearly that with increasing fault percentages, performances increases. This is because, with increasing fault percentages, the correlation between faulty and normal data decreases. However, the normal reading correlation remains high, making classification hyperplane more and more precise.
ROC Analysis: Receiver Operating Characteristics (ROC) analysis studies the sensitivity and the Specificity of the classifier. A ROC curve is a plot in which x-axis is the Specificity, and y-axis is the sensitivity of the classifier.
Area Under curve: The total area under the ROC curve is abbreviated as AUC. ROC curves are generally used to evaluate the performance of various machine learning algorithms as it gives a comprehensive and visual method of summering the accuracy of an algorithm [6]. Therefore, in this paper, ROC and AUC analysis is used as one of the performance metrics to evaluate our proposed method. From the ROC curve, AUC is calculated by calculating the size of the area under curve. The higher the area, the Better the performance of the method. Figure 3 shows ROC for various fault percentages. AUC is given in Table 3.

Comparison with existing work
Simulation results is also compared with the existing works of Zidi et al. [21] and Jan et al. [5]. Figure 4a shows comparative results of Detection Accuracy for our method and   Table 4). Figure 4a shows that our method works better with higher fault percentages. Our approach is based on the correlation coefficient. As we know, if the fault percentages increase, the correlation coefficient among data measurements decreases. For this reason, our method works better in high fault percentages. Figure 4b shows comparative results of Average False Positive Rate (FPR) for our method and the method of Zidi et al. Figure 4b shows that our method shows 14% improvement in terms of Average False Positive Rate compared to Zidi et al. Correlation among measurements of non-faulty node is high. So, a non-faulty node measurement detected as faulty is less likely to happen. That is why our method has better FPR.

Performance Evaluation For dataset
The proposed method also evaluated using a real-world environmental dataset of [12]. From this dataset, data of temperature and humidity for the indoor sensor node of the multi-hop scenario with anomalies is taken. In this dataset, the anomaly is achieved on a sensor node by a hot water kettle, which increases the temperature and the humidity simultaneously.
AUC and ROC are given in Fig. 5.

Conclusion
In this work, a fault detection algorithm based on Pearson's correlation coefficient and Support Vector Machine classification is presented. Simulation results show that the method works better in high fault percentages. This is because, with increasing fault percentage, correlation decreases. The proposed method is also evaluated using a dataset [12]. Simulation results, as well as results from real-world data, show that the proposed method successfully classifies fault and normal sensor reading. A possible future directive could be working with more than one machine learning algorithm to detect faults. Applicability of the machine learning methods over some complex networks and data with more event parameters may also be studied in the future.


extfieldsofstudy:
Computer Science
Mathematics

text:
Frozen Gaussian Sampling: A Mesh-free Monte Carlo Method For Approximating Semiclassical Schr\"odinger Equations

In this paper, we develop a Monte Carlo algorithm named the Frozen Gaussian Sampling (FGS) to solve the semiclassical Schr\"odinger equation based on the frozen Gaussian approximation. Due to the highly oscillatory structure of the wave function, traditional mesh-based algorithms suffer from"the curse of dimensionality", which gives rise to more severe computational burden when the semiclassical parameter \(\ep\) is small. The Frozen Gaussian sampling outperforms the existing algorithms in that it is mesh-free in computing the physical observables and is suitable for high dimensional problems. In this work, we provide detailed procedures to implement the FGS for both Gaussian and WKB initial data cases, where the sampling strategies on the phase space balance the need of variance reduction and sampling convenience. Moreover, we rigorously prove that, to reach a certain accuracy, the number of samples needed for the FGS is independent of the scaling parameter \(\ep\). Furthermore, the complexity of the FGS algorithm is of a sublinear scaling with respect to the microscopic degrees of freedom and, in particular, is insensitive to the dimension number. The performance of the FGS is validated through several typical numerical experiments, including simulating scattering by the barrier potential, formation of the caustics and computing the high-dimensional physical observables without mesh.

Introduction
The Schrödinger equation, which depicts the evolution of the wave function of a quantum system, is given as follows: where m o denotes the mass of the particles and denotes the reduced Planck constant. By adopting the nondimensionalisation techniques (see Appendix A), one obtains the following semiclassical Schrödinger equation in the dimensionless form: where the parameter ε, named the semiclassical parameter, describes the microscopic macroscopic scale ratio and is small within the semiclassical regime (i.e. 0 < ε 1). The quantum mechanical wave function u(t, x), valued in C, often serves as an auxiliary quantity to compute the physical quantities of interest. The physical quantities, including position, momentum and energy, etc., are referred to as "observables", which are of the form of quantum averages with respect to the wave functions of certain selfadjoint operators: Here,Ô denotes the operator related to a physical observable. For examples,Ô q = x and O p = −iε∇ are the operators corresponding to the position and momentum observables, respectively. In the semiclassical regime, namely ε 1, the solution to equation (1.2) is highly oscillatory both in space and time on the scale O(ε) and doesn't converge in strong sense as ε → 0, which gives rise to significant computational burdens [13,17]. Nevertheless, various effective numerical algorithms have been proposed so far, succeeding in capturing the highly oscillatory wave functions accurately with affordable computational complexity. The existing algorithms can be divided into two classes: direct discretization methods and asymptotic methods. Among the direct discretization methods, one of the most prevailing choices is the time splitting spectral method proposed in [1], where the meshing strategy ∆t = O(ε) and ∆x = O(ε) is practically sufficient for accurate approximation of the wave function. The time step size of the time splitting spectral method can be even relaxed to O(1), in other words, independently of ε, if one is only interested in capturing the correct physical observables, whereas one still needs to resolve the spatial oscillations. The asymptotic methods, including the WKB methods [5,6] and various wave packets based methods [29,14,18,19,20,27,12,24], provide accurate numerical approximation to the wave function with a relatively less computational cost with respect to the parameter ε compared to the direct discretization methods. The wave packets based methods have been extended to more complicated quantum dynamics, such as the vector potential cases and systems with random inputs [30,27,12,24], providing broad application prospects in many applications.
Although the existing algorithms have made great progress in efficiently simulating the wave functions and the physical observables in the semiclassical regime, they still suffer from the "the curse of dimensionality" when considering high dimensional cases, since those methods necessarily reply on an Eulerian mesh either in the initial decomposition or in the time evolution. When the dimension number m 1, the computation and memory amounts of the both classes of methods increase exponentially as m increases, which often leads to more infeasible computational cost when ε 1. For example, for the time splitting spectral method, the degree of freedom to approximate the wave function on an m−dimensional mesh is O(ε −m ), not to mention the increased computational complexity to evolve the wave function on this high dimensional mesh. Therefore, a novel numerical method for the Schrödinger equation in semiclassical regime, whose cost is insensitive to ε and the dimension number m, is highly desirable.
To the best of our knowledge, there has been very few existing methods, which are feasible or have the potential to be suitable for high dimensional simulations. From this perspective, the path integral sampling method proposed recently in [22,21] can be viewed as a pioneering work, which is constructed on the so-called Frozen Gaussian approximation (FGA). The FGA represents the wave function with a superposition of Gaussian functions in the phase space with fixed widths and is often wrote in the form of an integral on the phase space: A(t, q, p)e iΘ(t,x,q,p) ε dqdp, (1.4) where the variables A and Θ are can be approximately obtained through evolving a system of Lagrangian ODEs, thus the full quantum dynamics is reduced to an ensemble of wave packets dynamics. There have been quite a few works related to the FGA, including the theoretical studies in the chemistry community [8,9,15,16], the rigorous convergence analysis [28] and the FGA based algorithms for solving hyperbolic systems which reply on the direct phase space discretization [18,19,20]. However, it is not until [22,21] that people find another distinct prospective, giving the FGA ansatz (1.4) a probabilistic interpretation, and show that the Monte-Carlo sampling method based on the FGA gives a feasible algorithm which is mesh-free in computing the physical observables. The main spirit of [22,21] can be summarized as follows: firstly, as an integral on the phase space, the FGA ansatz can be viewed as an expectation of random variables on the phase space, which can be calculated by solving the parameter ODEs; in addition, when the initial data is Gaussian, the module function |A(0, q, p)| is also Gaussian and naturally serves as an ideal choice of a sampling density function on the phase space. However, what remains to be done is that there have been no rigorous numerical analysis on the sampling error nor discussion on how to extend the sampling method to general cases when the initial data is not Gaussian. Our current work aims to fill the vacancy on these issues.
In this work, we develop a Monte Carlo method which relies on the probabilistic interpretation on the FGA ansatz, and we shall name it the Frozen Gaussian Sampling (FGS). The essence of the FGS is to properly choose a reference density function on the phase space, and reformulate the FGA representation (1.4) as an expectation of a complex-valued functional with respect to the Frozen Gaussian trajectories. In this view, the FGS is a path-integral approximation of the wave function, and in the dynamical picture, the evolution of the quantum state is replaced by propagating an ensemble of Frozen Gaussian wave packets, where the initial phase space coordinates are sampled from the reference density function. The design of an efficient FGS algorithm boils down to proposing a suitable reference sampling measure, which is closely related the form of the initial condition of the wave function, and the quality of an FGS algorithm is measured by the variance of the functional to be sampled. Hence, the FGS provides a subtle shift in the computation paradigm and the analysis framework, which, as we shall elaborate, turns to show great promises as a simulation tool for high dimensional problems.
We focus on two of the most important types of initial quantum states: the Gaussian wave packet and the WKB data, and study the respective strategies in designing the density functions that balance the need of variance reduction and sampling convenience. For the Gaussian initial condition scenario, we rigorously prove that, with the feasible choice of the sampling density, the variance of the functional to be sampled is independent of the scaling parameter ε . Furthermore, the number of samples to reach a certain sampling error δ is in the order O c m δ −2 , where c is a O(1) constant and m is the dimension number. In particular, if one aims to halve the sampling error, only four times of the samples are needed, which is independent of the dimension m or the semiclassical parameter ε. Different from all the existing methods, this result shows that the number of wave packets does not to "resolve" any spatial oscillations, and is of a sublinear scaling of the microscopic degrees of freedom. Besides, the computation cost grows rather gently as m increases, and the whole FGS algorithm is parallelizable, which makes high dimensional simulations rather affordable. For WKB initial data cases, we propose a sequential sampling strategy leveraging the mono-kinetic structure of the highly oscillatory wave function. We numerically verify that the sample size to reach a certain error is also insensitive to ε, which still outperforms the existing mesh-based algorithms whose computational cost grows rapidly as ε → 0.
It is worthy emphasizing that, we only require that the initial condition to be in cer-tain forms to carry out efficient initial sampling, but we do not require that the solutions remain in such ansatz forms. This is a natural consequence of the construction of the FGS algorithm, while most asymptotic methods break down when the solution ansatz is no longer valid. And we have carried a set of 2D numerical tests to illustrate such a property, where the initial Gaussian wave packet is scattered into various profiles by different potential barriers, and the FGS method is able to capture the correct evolution behavior with a reasonably small sample size.
Another main highlight of the FGS lies in that it is mesh-free in computing the physical observables. By taking advantage of path integral representation, the FGS succeeds in computing a wide class of physical observables without reconstructing the wave functions. Therefore, if one is only interested in physical observables rather than the detained behavior of the wave function, the FGS serves as an efficient algorithm to simulate high dimensional Schrödinger equation, since it is a mesh-free algorithm of a sublinear scaling complexity and is naturally compatible with parallel computing. To demonstrate the performance of the FGS in computing observables, we implement the FGS to compute the physical observables in high dimensional quantum systems with a harmonic potential and it is shown numerically that the computational cost to reach a certain error is rather affordable in various dimensions. With the sample size quadrupling when the dimension number increases by one, we observe that the errors in the observables is uniformly bounded in all dimensions tested.
The rest of this paper is organized as follows: in Section 2, we review the basic setup the FGA and give the detailed procedures to compute the wave functions and physical observables with the FGS. We also present the error analysis framework of the FGS and elaborate the difficulties to be overcome when choosing the density function for the FGS. In Section 3, we focus on Equation (1.2) with Gaussian initial data and prove that the quantitative estimates of the FGS sampling error . In Section 4, we propose a feasible sampling strategy to apply the FGS on Equation (1.2) with WKB initial data. In Section 5, the numerical performance of the FGS is validated through five carefully designed numerical experiments, including simulating the Gaussian barrier potential and the caustics and computing the high-dimensional position and momentum observables without mesh.
we present the frozen Gaussian sampling (FGS) algorithm as a general computation methodology for semiclassical quantum dynamics and discuss the practical implementation of computing physical observables based on the FGS algorithm.

The frozen Gaussian representation for the semiclassical Schrödinger equation
The frozen Gaussian approximation is an asymptotic approximation to the solution of Equation (1.2) with O(ε) error [28]. It is based on an integral representation on the phase space: Here the phase function Θ is given by For a given initial condition q and p, the evolution of Q and P is governed by the classical Hamitonian H(q, p) = E(q) + 1 2 |p| 2 and defines a trajectory on phase space R 2m called an FGA tracjectory, which also determines the value of the amplitude A and action S. We name the variables Q, P, S, A as FGA variables. The evolution of the FGA variables is determined by the ODE system as follows: where we use the short hand notations Moreover, the initial condition is given by In addition, we remark that the phase space integral (2.1) can also be viewed as an integral of the FGA trajectory. Denote the FGA tracjectory by z t = (Q(t), P (t)), where we hide the dependence of the FGA variables on (q, p) for simplicity (we use the similar short hand notation in the rest of the paper). Then z t satisfies the Hamiltonian flow and defines a one-to-one map z 0 → z t on the phase space. In other words, the FGA trajectory z t is specified once the initial condition (q, p) is fixed and the FGA variables S and A are also specified with a given trajectory z t . The idea of viewing the FGA ansatz (2.1) through FGA trajectory is important in the introduction of the frozen Gaussian sampling.
To conclude this subsection, we directly quote two theorems from [18,28]. The first theorem shows the FGA ansatz integration (2.1) reproduces the initial condition: where the initial amplitude A(0, q, p) given in Equation (2.10) is determined by u in .
The second theorem shows that the L 2 error of u FGA with respect to the exact solution is O(ε): Theorem 2.2. Consider the semiclassical Schrödinger equation (1.2) with the exact solution u(t, x). Assume that the potential function E(x) ∈ C ∞ (R m ). For any t > 0, we have the following L 2 error estimate result for the FGA ansatz (2.1): where C A (t, m) is a parameter independent of the scaling parameter ε.

The frozen Gaussian sampling algorithm
We have seen in Section 2.1 that the FGA ansatz (2.1) is an integral on the phase space. With the FGA, the task of solving a PDE is converted into solving an ensemble of ODEs (2.3)-(2.6). However, difficulties still remain in computing the integral (2.1), which can be of high dimensions. Direct mesh method of numerical integration may suffer from "the curse of dimensionality" and does not lead to a practical algorithm. To deal with the high-dimensional integral (2.1), we seek a stochastic method for calculation. To this end, we write the FGA ansatz in the form of an expectation and approximate the integral (2.1) numerically with a Monte Carlo method. The resulting stochastic algorithm is called the frozen Gaussian sampling (FGS), which also gives a meshfree implement to compute physical observables. In this subsection, we elaborate the FGS algorithm and discuss its advantages and difficulties to be further investigated.
To write the FGA ansatz (2.1) into an expectation, we introduce a real-valued probability density function π(·) on the phase space R 2m . The density π then defines a probability measure P 0 on R 2m : where P 0 (R 2m ) = R 2m π(z 0 )dz 0 = 1. With the probability measure, we can rewrite the FGA ansatz (2.1) as (2.14) Thus, we may use a Monte Carlo sampling for u FGA (t, x): are M independent identically distributed samples with the probability density function π. We refer to u FGS as the frozen Gaussian sampling (FGS) wave function. We also note that Λ (t, x, z 0 ) can be treated as a random variable related to z 0 with the expectation Algorithmically, once we have sampled z 0 , the FGA trajectory and FGA variables are available through evolving the ODE system Equation (2.3)-(2.6) up to time t. With enough samples of the FGA trajectories, we approximates the FGS wave function u FGS numerically by Equation (2.15). From the discussion above, we give the pseudo-code of the frozen Gaussian sampling (FGS) in Algorithm 1. For simplicity, in the rest of this paper, we use the short hand notations to the FGA trajectories and FGA variables: We note that analogous short hand notation may be also applied to other variables without confusion. However, there are certain difficulties when actually executing the FGS algorithm, especially in high dimensional cases. First, to compute the initial amplitude, one needs to calculate the integral (2.10). For a Gaussian wave packet initial condition u in , the initial amplitude A(0, q, p) can be calculated explicitly in Equation (3.2), which has been pointed out in [22]. However, for a more complex initial condition u in such as a WKB initial condition, the explicit expression of A(0, q, p) is unavailable. To compute the initial amplitude (2.10), the computational cost of numerical integration is huge due to the high dimension and high frequency oscillation structure of the integrand.
Next, let us discuss the error analysis framework of the FGS wave function u FGS . Obviously, the L 2 error function E 0 = u FGS − u L 2 of the FGS wave function with respect to the exact solution can be divided it into two parts: Here E FGA = u FGA − u L 2 denotes the asymptotic error of the FGA ansatz and E S = u FGS − u FGA L 2 denotes the Monte Carlo sampling error. Thanks to Throrem 2.2, we arrive at where C A (t, m) is a constant depending on time t and dimension m and independent of the rescaled parameter ε. However, with unreasonable choice of density function π, the Monte Carlo sampling error E S may increase dramatically as the dimension number m increases, thus formidably many samples are needed. Therefore, to apply the FGS in high dimensional cases, variance reduction techniques such as importance sampling are necessary to reduce the number of samples. Ideally speaking, we aim to choose a reference density π which not only leads to a minimized sampling variance, but also easy to sample numerically. If we solely consider the variance reduction prospective, an obvious choice of sampling density can be given as follows: For a Gaussian wave packet initial condition, the density function given by Equation (2.20) actually defines a multivariate normal distribution and the resulting FGS algorithm is proven numerically practical and effective [22,21]. For a more complex WKB initial condition when the explicit expression of |A(0, z 0 )| in Equation (2.20) is unavailable, how to choose π(z 0 ) is still unclear. As far as we know, no rigorous numerical analysis on sampling error E S has been carried out prior to this work. The aim of the rest of this paper is to investigate and overcome the above difficulties of the FGS. In Section 3, we prove rigorously that the sampling error E S of FGS with sampling density (2.20) is independent of ε with a Gaussian initial condition. In Section 4, we introduce an asymptotic expression approximating the initial amplitude (2.10) when the initial data is of a WKB form and present a new initial sampling method of the FGS algorithm based on this approximation. In Section 5, we show various numerical examples to test the practical performance of the FGS with these new techniques, especially in high dimensional cases.

Computing observables without mesh
For many applications, the goal is not to approximate the wave function in semiclassical Schrödinger equation (1.2) itself, but rather computing certain observables. To calculate an observable, reconstructing the wave function on a mesh is often neither feasible nor necessary, especially for high dimensional cases. In this subsection, we introduce a mesh-free method to compute a wide class of observables, which only use the FGA trajectories sampled by the FGS and doesn't require the reconstruction of the FGS wave function. The method exhibits great advantages on computing high dimensional problems.
Let O = u|Ô|u whereÔ denotes a self-adjoint operator asscoiated with a physical observable. Then we compute O based on the FGS wave function u FGS as follows: We observe that the only dependence on x in Equation (2.21) is in the term g j,k , which is an integration over the Gaussian wave packets. For a wide class of observablesÔ, the term g j,k (t, x) can be carried out analytically through Gaussian integration formula, hence the observables O = u|Ô|u can be computed without mesh through computing the double summation in Equation (2.21). This class of observables include functions of the position variable and the momentum variable, which is in a polynomial form or can be approximated by a polynomial. Now we take the position observable and the momentum observable as examples to show the analytical expression of g j,k . For the position observableÔ = x, we obtain For the momentum observableÔ = −iε∇ x , we obtain To sum up, given the FGA trajectories and FGA variables, we can reconstruct the observables in a mesh-free fashion without computing the FGS wave function. Clearly, the complexity for the physical observable evaluation (2.21) is O(M 2 ). In the following sections, we shall further explore the dependence of the ensemble size M on the scaling parameter ε and the dimension number m and show great potential of the FGS in computing physical observables in high-dimensional cases.

Frozen Gaussian sampling with Gaussian initial conditions
In Section 3, we consider a m−dimension normalized semiclassical Gaussian wave packet centering atq with momentump: The semiclassical Gaussian wave packet is one of the simplest yet most important wave function in quantum dynamics, whose semiclassical limit arrives at a delta function on the phase space. With the Gaussian integration formula, we are able to calculate the initial amplitude A(0, q, p) (see Equation (2.10)) explicitly and the resulting density function π as shown in Equation (2.20) is Gaussian. Finally, we rigorously prove the sampling error E S of the FGS with density function is independent of ε and discuss its dependence on the dimension number m.

Initial sampling and initial amplitude for Gaussian initial condition
Plug the Gaussian initial condition (3.1) into Equation (2.10), we calculate A(0, q, p) analytically with the Gaussian integration formula: and As we have stated, the probability density function π(z 0 ) given by Equation (2.20) is a feasible choice of initial sampling. For Gaussian initial condition (3.1), the resulting probability distribution for the initial position and momentum z 0 = (q, p) is a multivariate normal distribution centered at q =q and p =p with probability density function: Algorithmically, we can directly sample z 0 = (q, p) from Gaussian density function in Equation (3.4), then the initial amplitude A(0, q, p) can be computed analytically.

Main result of sampling error
In this subsection, we aim to establish a rigorous error estimate for the FGS algorithm. We justify that the Monte Carlo sampling error E S = u FGS − u FGA L 2 is independent of the scaling parameter ε and study the relationship between E S and the dimension number m. With the error estimate, we are able to discuss the proper sample size for the FGS for given ε and m and illustrates the advantages of the FGS algorithm, especially in high-dimensional cases.
Recall that u FGA denotes the FGA ansatz defined in Equation (2.1) and u FGS denotes the FGS wave function given in Equation (2.15). Now that u FGS is determined by a series of random variables z (j) 0 on the phase space, we define the quadratic expectation of sampling error E S and total error E 0 as follows: Then plug Equation (2.15) into the expression of E S : where is regarded as the sampling variance of the FGS wave function. The derivation of Equation (3.7) uses the fact that z (j) 0 are independent identically distributed with density π. We aim to estimate an upper bound for the integral I FGA .
Before we present the main result of error estimate, we introduce the definition of the subquadratic condition. We say the potential function E(x) satisfies the subquadratic condition if there exists a constant C E such that The subquadratic condition is commonly used in analysis of hyperbolic systems [28,19,20,21] and we also assume the subquadratic condition in our theorem. With all of the preparation, let us show the main result of this section: 2) up to time t with the normalized Gaussian initial condition (3.1). Suppose that the potential function E(x) ∈ C ∞ (R m ) and satisfies the subquadratic condition (3.9). When applying the frozen Gaussian sampling (FGS) as in Equation (2.15) with the initial sampling density π as in Equation (2.20), there exists a constant C 0 (m, t) independent of the scaling parameter ε and the initial data u in such that the sampling variance given in Equation Moreover, if the assumptions in Lemma 3.6 are valid, the constant C 0 (t, m) in Equation where γ 0 (t) > 1 is a constant independent of ε and m.
We leave the proof of Theorem 3.1 to Section 3.3 and discuss on the theorem first. As a corollary, the following estimate on the quadratic expectation E S and E 0 defined in Equation (3.12) and (3.6) is obvious: Under the same assumptions in Theorem 3.1, we have the following estimates to E S and E 0 : and To conclude this subsection, let us demonstrate the advantages of the FGS over the time-splitting spectral methods (see [1]) through the error estimate in Corollary 3.2. For simplicity, we assume a j = 1 for j = 1, 2, · · · , m and ignore the asymptotic error (3.14) where is another constant independent of ε and m. In order to guarantee that E 0 ≤ δ, the sample size of the FGS should be more than In contrast, if we use a mesh-based algorithm such as time-splitting spectral methods in [1] to numerically approximate the semiclassical Schrödinger equation (1.2), the spatial step length of the mesh should be at least of order O(ε) to capture the highly oscillatory wave function. For an m−dimensional semiclassical Schrödinger equation, the spatial degree of freedom of the mesh-based method should be at least Recall that the sample size N FGA is the m−th power of an O(1) parameter. Therefore, when scaling parameter ε 1 or the dimension number m is large, the computational cost to solve a given semiclassical Schrödinger equation with the FGS is much smaller than those with the grid based methods, i.e. N FGA N TSP .
Remark 3.3. We remark that though the sampling size N FGS increases exponentially as m increases, we examine from numerical examples that the increasing rate of N FGS is much smaller compared to ε −1 . For more information of this issue, see Example 5 in Section 5.3.

Proof of Theorem 3.1
Before we prove Theorem 3.1, let us first collect some properties of the Hamiltonian flow on the phase space. In this subsection, we assume the potential function E(x) ∈ C ∞ (R m ) and satisfies the subquadratic condition (3.9) without special instructions.
For simplicity, we denote the map on the phase space from initial time 0 to time t as: t, q, p), P (t, q, p)) =: (Q κt (q, p), P κt (q, p)) . For the transformation κ on the phase space, we denote its Jacobian matrix as Note that the transformation κ is canonical, i.e. its Jacobian matrix J κt satisfies for any (q, p) ∈ R 2m . Here I m denotes m × m identity matrix.
The following propositions give uniform bounds of the Jacobian matrix J κt (q, p) and its derivatives. where | · | denotes the matrix norm.
Note that Proposition 3.4 was proven in Proposition 1 of [28]. For the canonical transformation κ, we define where ∂ z = ∂ q − i∂ p . The invertibility and boundness estimate of Z κt can be proven by mimicking the proof of Lemma 5.1 in [20]. We summarize them into the following proposition.
Proposition 3.5. For given time t > 0 and k ∈ N, there exists a constant K 0 (t, k) such that sup
With the properties of the matrix J κt and (Z κt ) −1 , we finally establish the estimate of the FGA variable A(t, q, p): Lemma 3.6. Consider the FGA variable A(t, q, p) that solves Equation (2.6) on the phase space. For given time t > 0, there exists a constant K 1 (t, m) such that Moreover, if we assume the constant K 0 (t, 1) in Proposition 3.4 and 3.5 is independent of the dimension number m, there exists a parameter γ * (t) > 1 such that The proof is straightforward through Gronwall's inequality.
Proof. Based on Equation (2.6), we arrive at Note that the constants in Proposition 3.4 and 3.5 are independent of the dimension number m, we use the constant K * (t) = K 0 (t, 1) to estimate Equation (3.25): Then apply Gronwall's inequality to Equation (3.25): where γ * (t) = e 2(K * (t)) 2 is independent of ε and m.
Remark 3.7. We remark that the importance of the three propositions and lemmas lie in the fact that the estimate to the FGA trajectories and FGA variables is uniform with respect to the initial phase z 0 = (q, p). Hence we can make global estimate on the integration I FGA . Validity of the propositions and lemmas rely on the subquadratic condition of the potential function E(x).
After the preparation above, we finally move to the proof of Theorem 3.1: Proof. According to Equation (3.8), we have (3.28) and it suffices that we focus on the following secondary moment: Now recall that the sampling density π is given by Equation (2.20), then we can simplify Equation (3.30): Next we integrate M (2) (x) in Equation (3.31) with respect to the position variable x: Now we apply Lemma 3.6 in Equation (3.32): (3.33) Here K 1 (t, m) is a constant independent to ε. According to Equation (3.3), we can calculate the integration of |A(0, z 0 )| on the phase space: (3.34) Finally we use Equation (3.34) in Equation (3.33): Hence the right hand side of Equation (3.35) is a constant independent to ε. Let C 0 (t, m) = 2 3m 2 (K 1 (t, m)) 2 , we finally arrive at Equation (3.10). Through the expression of C 0 (t, m), it is obvious that Equation (3.11) is valid provided that Equation (3.24) in Lemma 3.6 is satisfied.

Frozen Gaussian sampling with WKB initial conditions
In this section, we study the semiclassical Schrödinger equation (1.2) with the following normalized WKB initial condition: where a in and S in are two real-valued functions on R m , which are both independent of ε. The initial amplitude a in follows With the celebrated WKB approximation (see Appendix B for a brief introduction), the exact solution to Equation (1.2) has an asymptotic description of WKB type before caustics onset, provided that a in is smooth and decays rapidly [13]. Practically, we mainly consider the cases with the following Gaussian initial amplitude: where a j > 0 are sufficiently big andx j ∈ R (j = 1, 2, ·, m) are given parameters. We denote the mean of a in asx = (x 1 , x 2 , · · · , x n ) T .

(4.4)
The cases with general initial amplitude a in may be treated in a similar fashion and we refer to Remark 4.4 for discussion on this issue. Compared to Gaussian wave packet initial data cases where the explicit expression of A(0, q, p) is available, the difficulty of applying the FGS in WKB initial data cases lies in that there is no natural choice for the reference density function π(z 0 ) which achieves the balance between variance reduction and convenience of sampling. Moreover, the WKB approximation to Equation (1.2) may break down after caustics onset, which brings additional difficulties to the numerical approximation of Equation (1.2). In this section, we put forward a new initial sampling algorithm based on the stationary phase approximation to the initial amplitude (2.10). Then we validate the numerical performance of the FGS associating with the new initial sampling strategy through various typical examples in Section 5.

Description of the algorithm
As we have stated, one of the difficulties of the FGS algorithm with a WKB initial condition is that A(0, q, p) (given in Equation (2.10)) can not be calculated explicitly, hence there is no obvious way to do initial sampling. In this section, we apply the stationary phase method to approximate the initial amplitude A(0, q, p) in Equation (2.10) and give a practical sampling method based on the analytical approximation to A(0, q, p). Before we start, let us summarize the stationary phase method into the following lemma: Lemma 4.1. Let α(x) and Ψ(x) be two smooth functions on R m . Suppose that α(x) ∈ C ∞ 0 (R m ) and Ψ(x) ∈ C ∞ (R m ). Assume that ∇Ψ(x) has a finite number of zero points {y k } N k=1 on the support of α and det ∇ 2 Ψ(y k ) = 0 for each zero point y k . We consider the highly-oscillatory integration As ε → 0, I ε has the following asymptotic behavior: where the term sgn ∇ 2 Ψ (y k ) denotes the number of positive eigenvalues of the matrix ∇ 2 Ψ (y k ) minus the number of negative eigenvalues.
A proof to Lemma 4.1 can be found in [3].
According to Equation (4.6), the initial amplitude A(0, q, p) in Equation (2.10) can be calculated asymptotically as follows: In order to illustrate the spirit of sampling strategy, we make some reasonable assumptions which simplify the asymptotic description of A(0, q, p) in Equation (4.7). The cases when the assumptions break down will be discussed in Remark 4.3. We will show the Assumption 1. For any p ∈ R m , there exists at most one y ∈ R m such that ∇S in (y) = p, i.e. ∇S in is an injective mapping. The inverse mapping of ∇S in is given by T : p → y. The definition domain of T may not be R n and we denote it by D(T ), i.e.
T : D(T ) → R m , p → y. (4.8) Recall that we assume the amplitude of the WKB initial data a in (x) decays exponentially as x → ∞, the support of a in can be viewed as a small interval. Hence it is reasonable to assume that for any p ∈ R m , there exists at most one y ∈ R m such that ∇S in (y) = p and a in (y) = 0. With this assumption, the summation in Equation (4.7) can be reduced to only one term. Besides, we will show later in Remark 4.3 that the difficulties caused by multivalueness of the mapping T are not essential.
For simplicity, we use the variable y to denote the unique image of the mapping T defined in Equation (4.8) at p ∈ D(T ) in the rest of this section without special instructions.
The next assumption is a technique assumption without lose of generality: Assumption 2. The inverse mapping T defined in Equation (4.8) is smooth and satisfies det ∇ 2 S in (y) = 0, ∀p ∈ D(T ). (4.9) Based on the two assumptions, the asymptotic approximation (4.7) of the initial amplitude A(0, q, p) can be simplified as follows: where χ D(T ) (p) stands for the indicative function on D(T ). Then the modulus of Equation (4.10) is given by For Gaussian initial condition cases, we directly use the real-valued Gaussian function |A(0, q, p)| on the phase space as the sampling density π in Equation (2.14), which results in the minimum sampling variance. However, for the WKB initial condition cases, it is impossible to sample towards |A sp (q, p)| given by Equation (4.11) directly since the non-Gaussian function |A sp (q, p)| is hard to sample. To balance between variance reduction and convenience of sampling, we choose the following sampling density function: where Z is a normalization parameter to ensure π a probability density, i.e. R 2m π(q, p)dqdp = 1.
We note that the only difference between |A sp (q, p)| in Equation (4.11) and the probability density π in Equation (4.12) lies in the exponent of the term det ∇ 2 S in (y) in the denominator. Hence the choice of π in Equation (4.12) is rather acceptable from the prospective of variance reduction. The following proposition shows that the probability density π in Equation (4.12) is also easy to sample, provided that the the amplitude a in is given by Equation (4.3): Proposition 4.2. Consider the WKB initial data (4.1) where the amplitude function a in of is given by Equation (4.3). Assume that Assumption 1 and Assumption 2 are valid. Suppose that the random variables Q and P are sampled with respect to the probability density π(q, p) defined in Equation (4.12), then Q and P follows 1. The marginal distribution of the random variable Q is a normal distribution with expectation E π q =x, wherex as in Equation (4.4) denotes the mean of a in .
2. Denote Y = T (P), then the conditional distribution of the random variable Y with respect to Q is a normal distribution.
3. The normalization parameter Z in Equation (4.12) is given as follows: (4.13) Proof. Recall that the mapping T as in Equation (4.8) is the inverse mapping of the gradient function ∇S in . We refer to J T (p) as the Jacobian matrix of the mapping T at p. Based on Assumption 2, the determinant of J T (p) follows: , ∀p ∈ D(T ). (4.14) Then we apply Equation (4.14) and the integral substitution rule to compute the marginal density π Q (q) of the random variable Q: Then the marginal distribution of the random variable Q is a normal distribution N (x, Σ 1 ), where the covariance matrix Σ is a diagonal matrix as follows: , · · · , ε + 1 a n . (4.16) To calculate the expression of the normalization parameter Z given in Equation (4.13), we just need to integrate π Q (q) by q to ensure that π Q (q)dq = 1, i.e.
For fixed Q = q, we denote the conditional density of the random variables P and Y with respect to Q as π P|Q (p) and π Y|Q (y). According to Equation (4.15), we compute π P|Q (p) as follows: Then we can compute π Y|Q (y): Note that the exponent term of the conditional density function π Y|Q (y) is quadratic of y, therefore the conditional distribution of Y with fixed Q is a normal distribution. For fixed Q = q, the conditional distribution of Y follows N (µ 2 , Σ 2 ), where the mean vector and covariance matrix is given by and , · · · , ε εa n + 1 . Motivated by Proposition 4.2, one can sample the random variables (q, p) obeying the probability density π defined in Equation (4.12) on the phase space by sampling two Algorithm 2 Initial sampling of (q, p) on the phase space for WKB initial data (4.1) with amplitude (4.3) 1: Initial position sampling: Sample q ∼ N (x, Σ 1 ), wherex and Σ 1 is given by Equation (4.4) and (4.16).
normally distributed random variable q and y, where p = ∇S in (y). The initial sampling method bases on Proposition 4.2 can be summarized into the following algorithm: To conclude this subsection, we make two remarks on the initial sampling algorithm 2. The first remark elaborate what happens when Assumption 1 is not satisfied. The second remark discuss the cases of WKB initial data (4.1) with general amplitude a in and generalize Proposition 4.2 to these cases. Remark 4.3. Now We discuss the possible treatments when Assumption 1 is violated. If there exists more than one y ∈ R n such that ∇S in (y) = p for some p ∈ R n , the asymptotic approximation A sp (q, p) of the initial amplitude A(0, q, p)(see Equation (4.7)) can't be reduced to only one term as in Equation (4.11), i.e.
The sampling density associating with A sp (q, p) in Equation (4.22) can be given through mimicking the density function (4.12): Note that the density function π sp determines a mixture distribution and each component of the summation in the density function π sp can be sampled through algorithm 2. Then it is natural to adopt the method of sampling the mixture model (see Section 4.2 in [4]) to sample π sp . Therefore, the difficulties resulted from the invalidity of Assumption 1 is not essential.
Remark 4.4. We discuss the possible treatments when a in is not Gaussian. Since a in is of a general profile, Proposition 4.2 no longer holds and there is no direct way to sample the probability density π given in Equation (4.12). To overcome these problems, we may approximate a in with a Gaussian amplitude a S in and sample the following probability with Algorithm 2: where we have replaced a in in Equation (4.12) with the approximated Gaussian amplitude a S in . Ideally, the density π S balances the need of variance reduction and convenience of sampling such that it serves as the sampling density of the FGS when a in is not Gaussian.

Numerical examples
In the previous sections, we have elaborated the FGS algorithm with Gaussian initial data and WKB initial data and studied the sampling error with respect to the sampling size M , the scaled Planck parameter ε and the dimension number m. In this section, we present five typical examples using the FGS to numerically approximate the semiclassical Schrödinger equation (1.2) to illustrate the theoretical results on the sampling error and manifest the practical potential of the FGS in solving Equation (1.2), especially in high-dimensional cases with small ε. The numerical examples include exploring the relationship between the sampling error and parameter ε, applying the FGS in challenging cases such as caustics and computing the physical observable with the FGS in high-dimensional cases. The rest of this section is organized as follows. In Section 5.1 and 5.2, we investigate the cases with Gaussian initial data and WKB initial data respectively. In section 5.3, we compute the physical observables for a high-dimensional problem with the FGS and compare the numerical results with the exact value obtained through the Hagedorn wave packet method.

Example 1: 1D Gaussian initial data examples
In this example, we aim to justify the theoretical results on the sampling error E S in Section 3 through numerical examples. We implement the FGS to approximate Equation (1.2) at time t = 0.5 and study the relationship among the sampling error E S , the scaled Planck parameter ε and ensemble size M .
In the Gaussian initial condition (3.1), let m = 1, a 1 = 2,p = − 1 2 andq = 1 2 . Then we obtain the following Gaussian wave packet initial data: Here we consider two different potential functions respectively: the harmonic potential E 1 (x) = x 2 2 and the torsion potential E 2 (x) = 1 − cos x. Recall that the L 2 sampling error E S is defined as follows: where u FGA denotes the FGA ansatz (2.1), which deviates from the exact solution to Equation (1.2) by the asymptotic error of the FGA ansatz. Note that u FGA is numerically unavailable, we choose the FGS wave function with a sufficiently large ensemble size M 0 as an accurate approximation of the FGA ansatz where M 0 = 1.5 × 10 5 is chosen large enough.
To compute the FGS wave function u FGS and u FGA , we use the fourth order Runge-Kutta method to solve the ODE system (2.3)-(2.6) in the time interval [0, 0.5] with time step length ∆t = 0.01 and reconstruct the wave function u FGS at t = 0.5. We emphasize that ∆t is chosen sufficiently small, hence the error due to the approximation of the ODE system is negligible. The wave function reconstruction is carried on the spatial interval [−π, π]. Also note that the FGS wave function u FGS with a fixed ensemble size M may vary with different initial samples z (j) 0 M j=1 , the sampling error E S is actually a random variable. Therefore, when approximating the FGS wave function u FGS , we average over numerical results obtained through 30 sets of initial samples for each simulation of E S (M, ε).
In Table 1 and Figure 1, we present the sampling error E S (M, ε) of the FGS with for different scaling Planck parameter ε = 0.0256, 0.0064, 0.0016, 0.0004 and two potential functions E 1 and E 2 respectively. One can easily notice that E S is independent of ε for both potential functions and the convergence rate of the sampling error E S with respect to the ensemble size M is of order O(M − 1 2 ), which is in accordance with the theoretical results in Theorem 3.1 and Corollary 3.2. We also underline that it only takes 400 samples for the sampling error E S to reach O(10 −2 ). As a reference, we have stated in Section 3 that the computational cost of the mesh-based algorithms to reach certain error should be at least O(ε −1 ), which indicates that the FGS is a practical and accurate algorithm to numerically approximate Equation (1.2) with 1D Gaussian initial data, especially for the cases when ε is small. Example 2: 2D Gaussian initial data example As we have stated, the FGS wave functions is constructed by an ensemble of Gaussian wave packets Λ z (j) 0 (see Equation (2.16) for definition) related to an ensemble of FGA trajectories with a certain probability distribution. In Example 1, the potential functions are so simple that the exact solutions to Equation (1.2) are still Gaussian wave packets at time t = 0.5, which lowers the difficulty of constructing the wave function from the ensemble of Gaussian wave packets Λ z (j) 0 . In this example, we implement the FGS algorithm on 2D semiclassical Schrödinger equation with more complicated potential functions, where the exact solutions to Equation (1.2) are not Gaussian wave packets as time evolves. Though it is much more difficult to approximate a wave function that is not Gaussian from an ensembles of Gaussian wave packets, the FGS method is able to capture the correct evolution behavior with a reasonably small simple size.
In the Gaussian initial condition (3.1), let m = 2, a 1 = a 2 = 2,p = (−1, −1) and q = (1, 1). The scaling parameter is fixed at ε = 1 96 . Then we obtain the following 2D Gaussian wave packet: 2 + 10e −5|x| 2 . We illustrate of the potential E 3 and E 4 in Figure 2. Intuitively, due to the effect of the exponential term, the shapes of the potential functions E 3 and E 4 show peaks at the center. The main difference between the two potentials lies in that the peak of E 4 is much higher than E 3 , which leads to different behaviors of the exact solutions to Equation (1.2).
In this example, we apply the Strang Splitting Spectral Method (SP2) [1] with sufficiently small step lengths to compute the exact solution to Equation (1.2). We reconstruct the FGS wave function on spatial region [−π, π]×[−π, π] through Algorithm 1 and compare the resulting FGS density function with the exact density function ρ(t, x) approximated by the SP2. Similar to the previous example, we use the fourth order Runge-Kutta method to solve the ODE system (2.3)-(2.6) in the time interval [0, 0.5] with time step length ∆t = 0.01 and reconstruct the FGS wave function at different times. We choose the ensemble size M = 2000. Note that M = 2000 is of lower order than O(ε −2 ), which is the degree of freedom of the mesh-based method SP2.
In Figure 3 and 4, we show the evolution of the FGS density functions with two different potential functions E 3 and E 4 and compare them with the density functions computed via the SP2 algorithm. The density functions under both potential functions are no longer Gaussian after certain time. For E 3 in Figure 3, the wave packets pass through the potential peak at the center and disperse into different directions, which makes the wave packet crescent. In contrast, the potential peak of E 4 is so high and steep that the wave packets in Figure 4 rebound against the potential peak and disperse into the opposite directions against the wave functions in Figure 3. For both cases, the profile of the FGS density functions conform with the exact solutions when the wave packets are no longer Gaussian, which ensure the accuracy of the FGS method in computing non-Gaussian wave packets. In future work, we are ready to apply the FGS to approximating more complicated wave packets arising in scientific problems.

Examples with WKB initial condition
In this subsection, we apply the FGS with the initial sampling algorithm put forward in Section 4 to two examples of Equation (1.2) with WKB initial data. Caustics forms in the former example and does not form in the latter, which makes the former example more challenging for the FGS. To avoid confusion, we would like to point out that caustics become relevant to Equation (1.2) only in the semiclassical limit, ε → 0, but it also causes numerical challenges for moderately small ε.
Example 3: WKB initial data example with caustics.
In the WKB initial condition We choose the constant potential E 5 (x) = 10. Due to the compressive initial velocity d dx S in , caustics will form. This example was used in [23,1,14] as a reference to test whether certain algorithms can capture caustics.
In this example, we plot the FGS wave function and the resulting FGS current density function and compared them with the reference solutions computed via the SP2 method with sufficiently small step lengths. The computation is carried out on spatial interval [−π, π]. Similarly, we use the fourth order Runge-Kutta method to solve the ODE system (2.3)-(2.6) in the time interval [0, 0.54] with time step length ∆t = 0.01 and construct the wave function and current density at t = 0.54. We fix the scaling Planck parameter as ε = 0.0016 and plot the FGS wave functions and current density functions with ensemble sizes M = 400, 800, 1600 in Figure 5. The wave functions and current density functions computed via the FGS can reproduce the same profile as the reference solutions, which ensures the sampling strategy proposed in Section 4 is reliable and of low computational cost to approximate Equation (1.2) with the WKB initial conditions even in the presence of caustics.

Example 4: WKB initial data example without caustics
Now that we have shown the reliable results of the FGS in computing complicated WKB initial data problems with caustics, we further quantitatively study the sampling we obtain the following WKB initial data: We consider the null potential function E 6 (x) = 0 and the quadratic potential function These examples were used in [14] and caustics does not form. Similar to Example 1 in Section 5.1, we aim to compute the sampling error E S defined in Equation (5.2). The FGA ansatz u FGA in Equation (5.2) is approximated by a FGS wave function with sufficiently large ensemble size M 0 = 1.5 × 10 5 as shown in Equation (5.3). We use the fourth order Runge-Kutta method to solve the ODE system (2.3)-(2.6) in the time interval [0, 0.5] and reconstruct the wave function u FGS at t = 0.5. The computation is carried out on spatial interval [−π, π].
In Table 2 and Figure 6, we present the sampling errors E S (M, ε) of the FGS with for different scaling Planck parameter ε = 0.0256, 0.0064, 0.0016, 0.0004 and two potential functions E 6 and E 7 respectively. Similar to Example 1, the convergence of the sampling error E S is of order O(M − 1 2 ) as M → +∞. Moreover, we observe that the FGS sample errors of the WKB initial data cases mildly increase as the scaling parameter ε → 0 for both potential functions. Nevertheless, we underline that the increase of E S is insensitive to ε, which is manifested by the results in Table 2: when ε becomes 1 64 times smaller, the errors grow at most a couple times larger. Compared to the mesh-based algorithm whose computational cost is of order O(ε −1 ), the increase of the sampling error of the FGS is rather acceptable.

Computing physical observables in high dimensions
Example 5: High dimensional example In this subsection, we aim to compute the position observable and momentum observable as follows: where u denotes the exact solution to semiclassical Schrödinger equation (1.2). According to Section 2.3, the FGS provides a mesh-free algorithm to compute the position observable and momentum observable without reconstructing the FGS wave function u FGS . The physical observables computed via the FGS are denoted as follows: where we refer to Equation (2.21) for the computing expressions of O q FGA and O p FGA . We compute the mean-square error of the quantities of the physical observables: where we view the observables O q FGS and O p FGS computed via the FGS as a random variables, whose probabilistic distribution is determined by the FGS random samples z (j) 0 m j=1 on the phase space. We aim to discuss the mean-square errors E q and E p with respect to the semiclassical parameter ε and the ensemble size M .
To approximate the mean-square errors, we compute FGS are sampled by independent phase space samples. We choose the ensemble size M = 60 in our tests. To compute the exact quantities O q and O p of the position and momentum observables which are analytically unavailable, we consider a special case when the quantities of O q and O p can be given explicitly by a semiclassical Gaussian wave packet named the Hagedorn wave packet [10,11,30,17], which can be approximated numerically through evolving ODEs. For this special case, the position and momentum observables involved to the wave packet can be derived explicitly through the wave packet parameters. We directly quote Proposition 3.18 in [17] as the theoretical foundation of the Hagedorn wave packet method: Theorem 5.1. Consider the semiclassical Schrödinger equation (1.2) with a quadratic potential E(x) = |x| 2 2 . Let vectors q H (t), p H (t) ∈ R m and matrices Q H (t), P H (t) ∈ C m×m be the solution to the following ODE system: ) ds be the corresponding action function. Then the following complex Gaussian wave packet Moreover, the initial amplitude (a 1 , a 2 , · · · , a m ) ∈ R m is given as follows a j = 1 + {0.2j}, j = 1, 2, · · · , m, (5.27) where {x} denotes the fractional part of a rel number x.
Remark 5.2. As shown in Equation (3.12), the sampling error of the FGS is relevant to the term m j=1 1+a j √ a j 1 2 , which is related to the initial amplitude. To avoid the influence resulting from the initial conditions, we choose a j to be restricted in the interval [1,2).
We aim to compute the mean-square error E q and E p given by Equation ( The computation is carried out on spatial interval [−π, π] and the scaling parameter is fixed to ε = 0.0256. In figure 7, we plot the position and momentum observable errors E q and E p with different ensemble sizes M and dimension number m = 3, 4, · · · , 7. For fixed m, the convergence rate of E q and E p is roughly O(M −2 ), whereas the convergence slows down when M is large for m = 3. This is due to the fact that, when M is sufficiently large, the Monte Carlo sampling errors becomes much smaller than the asymptotic errors of the FGA ansatz that dominate the observable errors E q and E p . Moreover, we observe that, with fixed ensemble size M , the errors of the physical observables mildly increase as the dimension number m increases. As is manifested in , in order to reached a uniformly bounded mean-square error, the sample size should quadruple as the dimension number increases by one. Therefore, it is recommended that the sample size to compute physical observable in a m-dimensional problem should be of the order O(4 m ).
Furthermore, we record in Table 4 the system time to implement the FGS to compute an m−dimensional physical observables O q FGA or O p FGA . Since the complexity to compute the position and momentum observables are almost the same, we only present the system time to compute O q FGA in Table 4. To reached a uniformly bounded mean-square error, the system time increases at five times as m increases by one.

Conclusion and discussion
In this work, we have developed a novel stochastic simulation methodology for the semiclassical Schrödinger equation that shows great promises for high dimensional simulations. Its numerical cost is insensitive to the semiclassical parameter, grows gently with respect to the dimension number and it provides an efficient way to compute the physical observables without reconstructing the wave function. The core component of the algorithm is to properly choose a density function, that achieves a balance between variance reduction and sampling convenience, and should be specified for different types of initial wave functions. Although we have already obtained satisfactory results for the Gaussian initial conditions, and the WKB initial conditions, more advanced sampling techniques, in theory, can be utilized in the framework of the FGS methodology, and hence further improve its performance. Besides, the idea of turning computing highly oscillatory waves into sampling and propagating wave packets may also be applied to other wave equations in the high frequency regime. Thus, this work initiates many directions that we shall explore in the future.
With the characteristic flow X t : y → x W (r, y), the ODE system (6.7) yields a phase function as follows: Note that the characteristic flow X t : y → x W (r, y) may not be one-by-one, hence the phase function S W given in Equation (6.8) may be smooth only in a small time interval when characteristic flow X t : y → x W (r, y) doesn't cross. The cases when X t ceases to be a diffeomorphism are called caustics phenomena. See Figure 2 in [7] for an illustration of the caustics characteristic. Before caustics onset, the WKB ansatz in Equation (6.4) is an asymptotic approximation to the exact solution to Equation (1.2) of order O(ε), provided that the initial considition is exact and the amplitude a in decays rapidly (see Theorem 2.3 in [13]). After caustics forms, the exact solution may no longer be approximated by the WKB ansatz (6.4), hence the WKB approximation breaks down. The readers may refer to [2] for detailed discussion on the WKB approximation and analytic techniques after caustics onset.