Inconsistent Cluster Analysis With Disease Feature Enhancement ( ICADFE ) For American Cotton Leaf Disease Recognition

1497 Published By: Blue Eyes Intelligence Engineering & Sciences Publication Retrieval Number:E12200585C19/2019©BEIESP DOI: 10.35940/ijeat.E1220.0585C19 ABSTRACT--The broad leaves of cotton plant carry various visible disease symptoms. The ability of visual analysis by experts motivated the development of the plant disease recognition model. There are several visual feature descriptors, which can be primarily distinguished on the basis of pattern, texture or color. This system has been developed for the convenience of the farmers, who can avail the benefit by submitting the pictures of infected cotton leaves on the interface and the plant disease recognition system will return type of disease. In this paper, a dynamic feature descriptor is designed with inconsistent cluster analysis (ICA) and disease feature enhancement (DFE), which are combined as hybrid descriptor known as ICADFE for the recognition of the cotton plant disease. The ICADFE is found to improve the detection accuracy (approx. 80%), precision (approx 95%) and f1-measure (approx. 88%) on average in comparison with traditional shape and texture based feature descriptors such as scale invariant feature transform (SIFT), speeded up robust features (SURF) and fast retina keypoints (FREAK) with multicategory SVM (mSVM) for disease recognition.

I. INTRODUCTION
According to the statistics of the Ministry of Public Affairs, nearly 50-60% people are completely or partially involved in the agricultural practices in India, whereas nearly 60-70% of total population is entirely dependable over agriculture for their daily needs [2]. The cotton crop falls in the broad leaf type crops, where the disease symptoms are clearly visible and can be spotted over the leaf with the image processing models, which inspired research in this field. Methods have been developed to determine the type and stage of the leaf-borne cotton crop disease after evaluating the shape and texture of the given disease sample in order to recognize the crop disease, which can be proved to be efficient and helpful for the novice or un-skilled farmer and laborers to understand the crop disease in order to help them apply the control measures in the earlier stages. The proposed model is based upon the difference of Gaussian (DoG) & difference of hessians (DoH) based shape analysis from the given cotton crop samples, whereas the texture features are extracted in the form of histogram of oriented gradients (HoG) and Statistical Analysis of Structural Information (SASI) features for the texture feature description. In the proposed model, the overall working methodology and results are based upon the accuracy of statistical measures for the assessment of the overall performance of this model in comparison with the existing models.

II. RELATED WORK
Ole Mathis Opstad Kruse et al. [1] have proposed the visual feature extraction model based upon the color-based pixel and pixel group analysis. S. Arivazhagan et al. [2] have worked upon the robust machine learning and flexible classification method for the recognition of the leaf-borne disease. A.H.Kulkarni et al. [3] have worked towards the extraction of the wide variety of the features from the given samples, which includes the shape, vein structure and color based features from the given objects in the images. Qinghai He et al. [4] team has worked on the utilization of the flexible and high contrast enhancement of the features as well as the de-noising for the filtering of the Gaussian, salt & peers and spike noise from the obtained samples in order to improve the overall accuracy for the cotton leaf images. Prof. Sanjay B. Dhaygude et al. [7] has designed and implemented the texture based feature description for the detection and classification of the leaf-borne diseases with the probabilistic classification based upon the support vector machine (SVM). Haiguang Wang et al. [8] studied the various models for the robust and flexible classification methods for the leaf-borne diseases, where the powerful feature extraction based upon the color, pattern and texture have been deeply studied to understand their impact in detail.

III. DESIGN & IMPLEMENTATION
The proposed model has been designed by combining the shape and texture based feature descriptors along with the probabilistic classification based upon the SVM classification algorithm. The proposed model utilizes the shape and texture based features for the localization and extraction of the disease regions, which is followed by the disease classification in the given samples. The shape and texture oriented feature are designed upon the basis of different of Gaussian (DoG), difference of hessians (DoH), inconsistent cluster analysis, computes the combination of directional derivatives on X and Y axes using Gaussian Kernel, non-maxima suppression to eliminate the outliers and dependence of state (known as hysteresis) along with the multi-category enabled support vector machine (SVM) for the disease feature analysis of the affected cotton leaves. In this paper, the low level feature descriptors are employed for the automatic recognition of the disease in the cotton leaves to evaluate their performance. These texture-based feature descriptors are used to extract the low level features from the target images. The low-level features are known to extract the strongest features in the given image, hence, their performance is being evaluated in comparison with the proposed feature descriptors. The SIFT feature descriptor utilizes the different of Gaussians (DoG) based technique to describe the minima matrix of the input image. The Gaussian kernel, known as scale-space kernel, is used for the SIFT feature descriptor in this scenario. This Gaussian function is defined with L(x,y,σ), whereas the variable-scale Gaussian is computed with equation G(x,y,σ) over the image matrix I(x,y). Hence, the final equation can be described as following:

Shape based Features:
In this paper, the low level feature descriptors are employed for the automatic recognition of the disease in the cotton leaves to evaluate their performance. These texturebased feature descriptors are used to extract the low level features from the target images. The low-level features are known to extract the strongest features in the given image, hence, their performance is being evaluated in comparison with the proposed feature descriptors. The SIFT feature descriptor utilizes the different of Gaussians (DoG) based technique to describe the minima matrix of the input image. The Gaussian kernel, known as scale-space kernel, is used for the SIFT feature descriptor in this scenario. This Gaussian function is defined with L(x,y,σ), whereas the variable-scale Gaussian is computed with equation G(x,y,σ) over the image matrix I(x,y). Hence, the final equation can be described as following: L(x,y,σ) = G(x,y,σ) * I(x,y) ....Eq. 1 Where x and y give the location of pixels in the image matrix, denotes the intensity of the Gaussian filter and * gives the convolution procedure in the pixel (x, y), and G(x,y,σ)= e -(x*x + y*y) /2σ* σ ÷ (2πσ*σ) ....Eq. 2 Where σ denotes the intensity of the Gaussian filter, x and y denotes the given pixel position and value in the image matrix. The scale space extrema is computed with symbol e, in the above equation.
The SURF descriptor is the feature descriptor with increased execution speech, and uses the difference of Hessian matrix instead of Gaussian. The fast computation under SURF is achieved with box-type convolution filter.
The sum of the image I on the specific location (x,y) T represents the sum of all intensities to compute the integral image I s (x). The equation to compute the integral image can be given as: The integral image computation utilizes the sequence of the pixels, their intensities and summation procedure. The x and y give the number of rows and columns respectively, whereas i and j denote the current pixel position. I(i,j) denotes the value of current pixel, the sum of which is computed along with its neighbors. After the computation of the integral image, the hessian matrix is computed over the integral image matrix. The determinant of Hessian is observed by using the scale selection method. Each of the given point (x,y) in the integral image I is used to compute the Hessian matrix H(x,σ), where denotes the scale of the Hessian. The equation can be given as following: .....Eq. 4 Where L xx (x,σ), L xy (x,σ), L yx (x,σ) and L yy (x,σ) denote the different convolutions of second order derivative of Hessian, which are computed with (σ 2 /σx 2 )*g(σ). The FREAK feature descriptor is computed in the form of binary mask over the Hessians matrix computed under the SURF algorithm. The binary mask is prepared after computing the threshold value to describe the pixel with 0 or 1. The binary mask is further used as the feature descriptor rather than the Hessian matrix prepared by the SURF algorithm.

K-means analysis:
Image clustering is applied to the image data under the proposed disease recognition model. The k-means based image clustering is learnt to segment the images multiple clusters by fragmentation of pixel frequencies in the target image. The k-means algorithm is applied over the set of observations (x 1 ,x 2 ,x 3 …..x n ), which are pixels (Gray or RGB) in the target image, to partition the n no. of observations to k clusters, where (k! >n) and (k <=n) are satisfied.
The figure 1 shows the results of different executions of k-means algorithm on the same image. The cluster sequence is noticeable in all of the rounds. The results of first and second rounds are similar, whereas second and third clusters shuffled in 3 rd run. The 4 th and 5 th run also produced the results similar to 3 rd run, which again got shuffled in 6 th run. This problem is resolved by using the inconsistent cluster analysis method explained in the following section.

Figure 1: Proving the inconsistency of the K-means algorithm
The k-means clustering method aims to select the sets or clusters with minimum variance or difference in withincluster sum of squares (WCSS). Technically, the objective of k-means algorithm is to discover the following (Eq. 5): Where denotes the mean value of data points (or pixels), denoted , in image matrix. k gives the number of clusters, n represents the number of entities (or pixels) in cluster k. But there is problem of inconsistency associated with k-means clusters, which prevents the crop disease recognition model from detecting the correct cluster containing the disease symptoms. Hence, layer analysis model is used to select the required cluster.

Inconsistent Cluster Analysis:
The inconsistent cluster analysis is the method to determine the correct cluster containing the disease region in the target image. The adaptive inconsistent cluster analysis is applied to determine the target cluster containing the disease symptoms. The figure 2 explains the procedure used to determine the target cluster from the clusters produced by k-means algorithm. The inconsistent cluster analysis is performed with the supervision of training data prepared for selective clusters, which is used for classifier learning.  in figure 3. The inconsistency of the clusters can be noticeably verified according to the results shown by the inconsistent cluster selection approach in the proposed model. The result has been demonstrated on 6 step executions, where the k-means clustering approach is showing the inconsistent behavior, which is naturalized by inconsistent cluster selection as shown in figure 3.

Disease Feature Enhancement:
The image features learnt using inconsistent feature analysis method, are sometimes insignificant due to the use of unsupervised clustering, which must be transformed to produce consistent features. A reliable edge detection model has created to determine the regions containing the features in target cluster, which computes the combination of directional derivatives on X and Y axes using Gaussian Kernel, non-maxima suppression to eliminate the outliers and dependence of state (known as hysteresis). The X and Y derivatives of the target image cluster are computed by using the following kernel function (Eq. 1) and its derivative (Eq. 6).
K=(6 * σ) + 1 ....Eq. 6 Where K is the kernel function and denotes sigma (or frequency) to transform the input value. The following equations compute the first bi-directional derivative on X (Eq. 8) and Y (Eq. 7) axes respectively. Where R and C represent the rows and columns respectively, K denotes the kernel, and i and j gives the current row and column respectively. X2(i,j)=X1(i,j)+I(R+K1-1,C+K2-1)*X1(K,1) ....Eq. 11 Y2(i,j)=Y1(i,j)+I(R+K1-1,C+K2-1)*Y1(K,1) .....Eq. 12 Where X2 and Y2 denote the second derivatives of the image pixels, K1 and K2 denote the steps in range of kernel from 0 to K (Kernel value), X1 & Y1 represent the first derivatives and the i & j denote the current location of pixels in the terms of rows and columns respectively. Now, compute the magnitude of the second derivatives using following equations. The magnitude is computed between the Kernel ranges of to vertically and to horizontally, where H and W are number of rows and columns in the image matrix.
....Eq. 13 Where G denotes magnitude, X2 and Y2 are second derivatives computed using Eq. 11 and 12 respectively. To compute the non-maxima suppression, the kernel ranges similar to magnitude applied over following method: Where T is term frequency, and computed over X2 and Y2 derivatives of target image. The non-maxima suppression is created over the magnitude matrix, where the values falling in certain range are transformed to 0 on the basis of decision taken term frequency as shown in following equation: NS=G ….Eq. 15 .....Eq. 16 Where NS denotes non-maxima suppression matrix, and certain set of ranges (T) are applied to determine the outliers, which are converted to absolute zero in order to compute the non-maxima suppression matrix. Further, the hysteresis method is applied to determine the significance of each pixel in the feature matrix selected through the inconsistent cluster analysis method. The following method is used to determine the feature significance:  In figure 4, the disease feature enhancement is explained step by step. The visualization begins with original feature image before enhancement in the first placeholder. The multi layered derivatives are described in the first row (columns 2-4) and 2 nd row (column 1). In 2 nd row (column 2) the non maxima derivative is extracted on the basis of X and Y derivatives. Afterwards, the hysteresis size is corrected from the hysteresis matrix, which is computed for multiple Kernel values.
Then, the edge localization is performed on the hysteresis image after size correction. After edge localization, the neighboring pixel selection is performed, which is followed by gap filling method using morphological operations. Finally, the enhanced disease feature is extracted from the target image w.r.t. the binary mask.In the proposed model, we have combined the SIFT, SURF and FREAK to eliminate the false recognition cases. These feature descriptors work to determine the shape and texture of the maxima points in the given image matrix. The shape and texture features describe the strong visual pattern in the given image in the form of shape and texture in order to recognize the cotton disease. The proposed model uses the hybrid feature descriptor for combining the shape and texture based features has been developed as the improved classification model for plant disease detection from the infected leaves. The shape and texture features are considered to be the best features in the case of estimation and localization of the disease infected area in the given set of images. The SIFT feature descriptor model utilizes the difference of Gaussians for the creation of minima matrix, whereas the SURF descriptor utilizes the difference of Hessians for the similar purpose. The support vector machine (SVM) model has been incorporated for the purpose of classification of the plant disease in the given sample over the extracted features with SIFT, SURF and FREAK extractors. The following algorithm describes the workflow of the proposed model in detail: Proposed Algorithm 1.
Load the testing image matrix into runtime memory 2.
Compute the clusters over the testing image using K-means algorithm 3.
Return the clusters in similar sized image matrices 4.
Load the training data for inconsistent cluster analysis (ICA) 5.
Configure the mSVM (multi support vector machine) classification

6.
Train mSVM classification model with ICA training data matrix and label vector 7.
Classify the clusters using the mSVM trained with ICA data 8. Select cluster according to the classification decision, Set sigma, maximum & minimum threshold values 9.
Compute the kernel function (Eq. 6), Compute the first derivatives horizontally and vertically (Eq. 7 and 8) 10. Compute the second derivatives horizontally and vertically (Eq. 11 and 12) 11. Compute the magnitude matrix over second derivative (Eq. 13) 12. Store the magnitude value matrix as term frequency matrix (Eq. 14) 13. Normalize the term frequency for each magnitude value to determine the outliers to create the non-maxima suppression by filtering out the outliers (Eq. 15 and 16) 14. Compute the hysteresis matrix over the nonmaxima suppression matrix (Eq. 18) 15. Load training data for disease classification with disease labels, Train mSVM classifier with training data 16. Evaluate the testing feature computed on step 16 to determine the plant disease 17. Return the detected type of disease The multilayered decision logic has been incorporated in this algorithm, where all the features including Hysteresis, inconsistent cluster analysis, etc are analyzed under the layered design in order to finalize the observation. The voting based classification has been used over the observations provided by SVM with all visual features collectively.

IV. EXPERIMENTAL CLASSIFICATION RESULTS AND ANALYSIS
For the purpose of result analysis, a number of samples are tested under this research, which involves the different grades of disease with variable disease covered regions. The image data for crop diseases also includes the following image ( Figure 5), which carries the minimum disease covered regions, which can be analyzed with difficulty. This classification model is the proposed case is capable to detect the infected leaves even with the minimum disease covered regions.

Figure 5: Dataset image after clustering (a) Infected leaves (clusters) (b) Healthy leaves (clusters)
The results have been evaluated in the form of the statistical parameters and the performance measures in order to understand the effectiveness of the proposed model. The proposed model has been tested with 10% ratio for cross validation and remaining 90% data is used for the classifier learning. The results of the Inconsistent Cluster Analysis based Disease Feature Enhancement (ICADFE) model are collected and evaluated against the existing models for disease recognitions, as per shown in following table. The statistical parameter of accuracy, precision, recall and f1measure are computed for each of the feature descriptor, as shown in the following table: The table 1 contains the performance measures of the cross validation testing with selected samples from the given dataset. The proposed ICADFE descriptor has been recorded with 79.55% of accuracy, which outperformed SIFT (65%), SURF (45%) and FREAK (40%). Similarly, the ICADFE descriptor has outperformed using precision (96.67%) in comparison with SIFT (62.50%), SURF (45%) and FREAK (41.67%).

V. CONCLUSION
The proposed model is made capable of working with the variations in the variety of the diseases defined among the given image dataset of disease affected leaves of cotton, which are determined from the local or global image features such as texture, pattern or color. Also the visual orientations along the size of disease covered regions, occlusion and several color illumination variations are used for the purpose of image classification. The proposed model has been designed to produce the rigid and improved results than the existing models, which has been analyzed from the obtained results. The proposed model has been found improved in the terms of the all of the performance parameters, which includes the precision, recall, f1-measure and accuracy. The proposed model has been found highly accurate at approx 80%, against the 67.50%, 56.50 and 50% percent readings obtained from the other disease feature descriptors.
In the future, the proposed model can be improved by using the Self-learning AI techniques deep learning with the help of the neural network along with the shape and curve based features, which can further improve the performance of the script recognition models.

VI. ACKNOWLEDGEMENT
We are highly thankful and grateful to the Department of Research, Innovation & Consultancy of IK Gujral Punjab Technical University, Kapurthala Road, Jalandhar for giving us an opportunity and providing us all the required assistance to write this Research Paper. The constant guidance and encouragement received from Dr. A.P.Singh, Dean RIC, I.K Gujral Punjab Technical University has been of great help in carrying out the research work and is acknowledged with reverential thanks. The authors would like to express a deep sense of gratitude and thanks profusely to Dr.Pankaj Rathore, Director, Regional Research Station, Punjab Agricultural University, Faridkot, Punjab and Dr.Paramjit Singh, Director, Regional Research Station, Punjab Agricultural University, Bathinda, Punjab. Without his wise counsel and able guidance, it would not have been possible to complete the research work in this manner.


Traffic flow models with looking ahead-behind dynamics

Motivated by the traffic flow model with Arrhenius look-ahead relaxation dynamics introduced in [A. Sopasakis and M.A. Katsoulakis, SIAM J. Appl. Math., 66 (2006), p. 921--944], this paper proposes a traffic flow model with look ahead relaxation-behind intensification by inserting look behind intensification dynamics to the flux. Finite time shock formation conditions in the proposed model with various types of interaction potentials are identified. Several numerical experiments are performed in order to demonstrate the performance of the modified model. It is observed that, compare to other well-known traffic flow models, the model equipped with look ahead relaxation-behind intensification has both enhanced dispersive and smoothing effects

Introduction
In this paper, we are concerned with the shock formation phenomena -bounded solutions with unbounded derivatives -for a class of nonlocal conservation laws, where u is the unknown, F is a given smooth function, andū is given by (1.2)ū(t, x) = (K * u)(t, x) = R K(x − y)u(t, y) dy, where K to be chosen later. The nonlinear advection couples both local and nonlocal mechanism. This class of conservation laws, identified in [10], appears in several applications including traffic flows [20,8], the collective motion of biological cells [2,4], dispersive water waves [22,6,14], high-frequency waves in relaxing medium [7,16] and the kinematic sedimentation model [9,23,1].
Along with the above nonlocal models, the model that motivates the present work is the traffic flow model with looking ahead relaxation, introduced by Sopasakis and Katsoulakis: • A traffic flow model with Arrhenius look-ahead dynamics [20] (1. Here, u(t, x) represents a vehicle density normalized in the interval [0, 1], γ a is a positive constant proportional to the look-ahead distance and K 0 is a positive interaction strength. This model takes into account interactions of every vehicle with other vehicles ahead within the look ahead distance γ a . Some careful numerical study of the above traffic flow model is carried out in [8]. In addition to this, an improved interaction potential . This linear interaction potential is intended to take into account the fact that a vehicle's speed is affected more by nearby vehicles than distant ones. In the case of a good visibility(large γ a ), the numerical examples in [8] suggest that (1.5) with the linear potential yields solutions that seem to better correspond to reality. Several finite time shock formation scenarios of solution to (1.5) with (1.6) were presented in [11]. The authors in [10] identified threshold conditions for the finite time shock formation of (1.5) subject to two different potentials above. The sub-thresholds for finite time shock formation conditions in [10] are consistent with the numerical result in [8].
We set K 0 = 1 through out this paper, since this parameter is not essential in our blow-up analysis. Then we can rewrite nonlocal termū = K * u associated with (1.6) and (1.7) as follows, respectively: and (1.9)ū(t, x) = 2 γ a x+γa x 1 + x − y γ a u(t, y) dy.
In this paper, we extend (1.5) by considering a look behind intensification. That is, we consider the traffic flow model with look ahead relaxation and look behind intensification: Here,ū is given in (1.8) and (1.9). For theũ = K b * u, we shall consider constant and linear interaction potentials. i.e., Here, γ b is a nonnegative constant proportional to the the look-behind distance. By setting K 0 = 1 again, we can rewrite the nonlocal termũ associated with the above kernels as follows, respectively: This look behind intensification model is intended to take into account the driving behavior of some drivers who actively react to vehicle distributions of one's ahead and behind. More precisely, considering the flux in (1.10), a vehicle's velocity is determined by (1 − u)e −ū+ũ . Here, the local traffic density u at one's location plays an major roles. In addition to this, averaged ahead and behind traffic densities of each driver adjust the velocity via relaxation and intensification effects. The drivers equipped with this strategy monitor/compare the densities ahead and behind within assigned distances, and prefer to choose accelerate(decelerate) when one has a relatively high density behind(ahead). It is natural to assume that γ a ≥ γ b > 0.
The motivation of the behind intensification may be explained as follows. On one lane highway, if there is a relatively high density behind a driver, the driver may feel that he/she is holding up the traffic and want to i) escape from pressure by driving faster if conditions allow, ii) accelerate to improve the flow of entire convoy of cars.
The objective of this article is twofold: i) We identify threshold conditions for the finite time shock formation of the traffic flow model (1.10) subject to two different potentials (1.8)-(1.13) and (1.9)-(1.14). As is known that the typical well-posedness result(e.g., [10]) asserts that either a solution of a hyperbolic partial differential equation exists for all time or else there is a finite time such that some norm of the solution becomes unbounded as the life span is approached. The natural question is whether there is a critical threshold for the initial data such that the persistence of the solution regularity depends only on crossing such a critical threshold. This concept of critical threshold and associated methodology is originated and developed in a series of papers by Engelberg, Liu and Tadmor [5,15,21] for a class of Euler-Poisson equations.
ii) We investigate performance of the proposed model (1.10) via numerical examples in comparison with (1.5) and the Lighthill -Whitham -Richards(LWR) [13,17] model We are interested in both dispersive and smoothing effects of the proposed global flux. It is well known that the u x of the LWR model blows-up if there is a point such that u x (t, x) t=0 > 0. Indeed, the derivative of d = u x along the characteristic satisfies the Riccatti equationḋ = d 2 that leads to the blow-up of d unless d t=0 ≤ 0. For some circumstances, this result is unrealistic because no shock formation is observed in the free flow(i.e., sparse traffic).
On the other hand, in [8], it is observed the global flux in look-ahead model (1.5) has some smoothing effect that seems to be able to prevent the shock formation. But simultaneously, the relaxation makes waves less dispersive. In our numerical examples in Section 3, it is observed that the proposed model (1.10) greatly improves waves' dispersion and smoothing phenomena. Now, the finite time shock formation results are collectively stated as follows. In the theorems, we assume that γ a ≥ γ b > 0.
then u x must blow-up at some finite time. .
then u x must blow-up at some finite time.
Regarding these results some remarks are in order.
i) The condition (1.16) reflects some balance between sup x∈R [u 0 (x)] and inf x∈R [u 0 (x)]. It seems the nonpositive term inf x∈R [u 0 (x)] is more negative, then sup x∈R [u 0 (x)] needs to be large for the finite time shock formation. It can be interpreted that not only the car density behind the traffic jam but also the car density ahead of the traffic jam contribute to the formation of shock.
ii) There is no direct comparison between (1.16) and look-ahead only model's((1.5)-(1.8)) blow-up condition which is obtained in [10]. Roughly speaking, if inf x∈R [u 0 (x)] is not too negative, the threshold in (1.16) is higher than the one in (1.18). Also, for example, if γ a = 1 and γ b = 0.5, then the right hand side of (1.16) is bigger than that of (1.18) for any initial data u 0 (x). The same thing holds for majorities of γ a and γ b values. This is interesting because, even though we insert behind intensificationũ which in turn increases the waves speed, the blow-up thresholds actually increased. These blow-up results are consistent with numerical experiments in Section 3. iii) In contrast to (1.16), the blow-up condition in (1.17) depends only on γ a , γ b and the size of the initial slope. It is interesting to observe that a little difference in the interaction kernels leads to different types of blow-up conditions; one involves two global terms sup x∈R [u 0 (x)] and inf x∈R [u 0 (x)], the other involves local u (x) term only. Furthermore, the blow-up threshold in (1.17) is higher than the one in (1.16) for not too negative inf x∈R [u 0 (x)]. Indeed, Example 3.4 in Section 3 shows that the model with linear interaction potentials has less steep wave(behind the traffic jam) than the model with constant potentials for all time.
We now summarize the main arguments in our proofs. We want identify some threshold condition for the shock formation of solutions to (1.10). The local existence and blow-up alternative results for (1.1) or (1.10) can be found in [10]. From Corollary 1 in [10], it suffices to track the dynamics of u x . The idea is based on tracing M (t) : . The existence and differentiability(in almost everywhere sense) of M (t) and N (t) are proved by Constantin and Escher [3], which we summarize in the following. (t) = u tx (t, η(t)) a.e. on (0, T ).
We also state a useful result, which is obtained in [12].
) Consider the following differential inequality, We remark that the above lemma remains valid even if (1.19) holds almost everywhere. We now conclude this section by outlining the rest of the paper. In Section 2, we prove Theorems 1.1 and 1.2. In Section 3, we demonstrate the performance of the proposed model via several numerical experiments.

Proof of theorems
Proof of Theorem 1.1: Let d := u x and apply ∂ x to the first equation of (1.10), Then, along (t, ξ(t)), using (1.8) and (1.13), we have and (2.1) can be written as, a.e.

Substituting this lower bound into (2.4), we obtaiṅ
Rewriting of this inequality gives where M 2 (≥ M 1 ) is given by We claim that M 2 has an uniform upper bound,

By Lemma 1.4 (i), if
then M (t) will blow up a finite time. This is exactly the threshold condition as stated in Theorem 1.1.
To complete our proof we still need to verify both claims (2.9) and (2.7).
To verify (2.9), we set From 0 ≤ u(t) ≤ 1 it follows that −1 ≤ v ≤ 1. If suffices to find upper bound for M 2 over the set Here, we use max (u,v)∈Ω {−2(1 − 2u)v + u(1 − u)} = 2 which can be verified easily since the underlying function is linear in v and quadratic in u. For the next one, By rearranging, (2.10) It follows that where the maximum value is achieved at ∂Ω. This gives Since is arbitrary, we choose = 1 2 to get Q 2 ≤ 1, hence Q ≥ −1, which gives (2.7). This completes the proof of Theorem 1.1.
Proof of Theorem 1.2: Let d := u x and apply ∂ x to the first equation of (1.10), (2.12) Using (1.9) and (1.14), we display the derivatives of the global terms(we omit the t−dependence). Also, We reorder d(x) terms in (−ū xx +ũ xx ) and obtaiṅ (2.15) We should point out that in contrast to the model with constant potentials(see (2.3)), the above equation does not have u 3) ). Later this difference leads to the structurally different blow up conditions. Rewriting of (2.15) gives where D 2 (≥ D 1 ) is given by To find an upper bound of D 2 , we set Then, one can see that and 0 ≤ u(x) ≤ 1 we obtain the desired bound. Now consider where v and C are defined above.
As in the proof of Theorem 1.1 we let Then it holds Here, we use max (u,v)∈Ω {(1 − 2u)v − u(1 − u)/2} = 1 which can be verified easily since the underlying function is linear in v and quadratic in u.
To find an upper bound of γaγ b 2(γa+γ b ) 2 C, we expand and obtain The last inequality holds because of (2.18) Therefore, (2.20) leads to Applying part (i) of Lemma 1.4, we complete the proof of Theorem 1.2. .

Numerical examples.
In this section, we demonstrate the performance of the proposed look ahead relaxationbehind intensification model (1.10), in comparison with look ahead relaxation only model (1.5) and LWR(non-global flux) model (1.15).
In all numerical examples, Lax-Friedrichs numerical scheme is applied, and we take large enough computational domain so that no waves touch its boundary within the final computational time. The Courant-Friedrichs-Lewy(CFL) numbers are chosen as 0.5 or 0.25, and all solutions are computed on uniform grids with ∆x = 1/100 or ∆x = 1/50. In order to compute the nonlocal termū(t, x) = K * u, we used the trapezoidal method.
We re-display the models of considerations here for readers' convenience: • (LWR) Lighthill -Whitham -Richards model: • (Look-A) Look ahead relaxation model: • (Look-AB) Look-ahead relaxation-behind intensification model: Here, the global terms associated with constant interaction potentials((1.6) and (1.11)) are given by Also, the linear interaction potentials((1.7) and (1.12)) lead to (3.5) In all examples, we should point out that our look ahead-behind model is better in the sense of dispersive and smoothing effects. with γ a = 1 and γ b = 0.5 subject to the following initial data: The situation corresponds to highly congested traffic(−4 < x < 2) is preceded by free traffic and followed by less congested traffic. See Figure 1.
In this example, on can clearly see the effects of the global fluxes. Look-A model's waves(Red) lag behind that of the LWR waves(Blue). The waves in Look-AB model have faster pace than the ones in Look-A. This can be explained by the awareness of traffic behind. While the waves with global fluxes develop no shock, the LWR wave(Blue) develop a shock discontinuity at t = 2. It is interesting to observe that the valley between two initial plateaus is filled quickly in Look-AB model, while Look-A model has persistent valley until t = 2.5. Overall, we should point out that Look-AB model has both enhanced dispersive and smoothing effects.  This initial condition corresponds to a situation in which the red traffic light is located at x = −2 and it is turned to green at the initial time instance. See Figure 2.
In this example, it is observed that waves in the Look-AB model have the fastest pace and the most prominent dispersive effect. Furthermore, if we zoom in near the steepest gradients(lower right), one can clearly see that the Look-AB wave has less steep slope. It looks like the behind intensification dynamics in the Look-AB model not only increase the flux, but also contribute to the dispersive effect through global interactions.  In this example, we observe that the LWR model's solution(blue) develops the shock. In contrast to this, the Look-AB solutions with linear potentials(magenta) and constant potentials(green) have no shock. Furthermore, the valleys between three plateaus are filled relatively quickly in Look-AB models.  With this very steep initial wave, one can again observe that the waves in Look-AB model have enhance dispersive effects compare to the LWR wave(blue). One can clearly see the shock formation in the LWR wave. Furthermore, it is observed that the Look-AB model with linear interaction potentials(magenta) has less steep wave than the model with constant potentials(green); see the lower right figure in Figure 4. This observation is consistent with remark iii) in the introduction.


s2fieldsofstudy:
Computer Science
Mathematics

text:
Multigrid method for pricing European options under the CGMY process

: We propose a fast multigrid method for solving the discrete partial integro-di ﬀ erential equations (PIDEs) arising from pricing European options when the underlying asset is driven by an inﬁnite activity L´evy process. We consider the CGMY model whose kernel singularity gets worse when the parameter Y approaches two. Due to the integral term, the discretization matrix is dense. In order to obtain an e ﬃ cient multigrid method, we apply a ﬁxed point iteration as a smoother for multigrid. In each smoothing step, we only need to solve a sparse matrix corresponding to the di ﬀ erential operator and compute a matrix-vector product involving the integral operator by a fast Fourier transform (FFT). We prove that the ﬁxed point iteration smoother is e ﬀ ective reducing the high frequency components. Moreover, we also prove a two-grid convergence of the multigrid method by a local mode analysis. We demonstrate the e ﬀ ectiveness of the multigrid method by solving the option pricing equation under the CGMY model with ﬁnite and inﬁnite variation processes.

Introduction
The standard Black-Scholes model for pricing options assumes the underlying asset price satisfies the geometric Brownian motion model. However, it has been observed in practice that market prices often have large jumps that do not follow the geometric Brownian motion. Jump diffusion models have been proposed to more accurately capture the actual market behavior; see e.g., [1]. Recently, models based on the Lévy process have become popular in the financial literature [1][2][3][4][5][6][7]. Option pricing under the exponential Lévy process with finite activity [5,[8][9][10][11][12] and infinite activity [13][14][15][16][17] has been extensively studied. Under this model, pricing of the European option requires solving a partial integrodifferential equation (PIDE), which consists of differential and integral terms in the equation.
Numerical solution of PIDEs is challenging since the discretization of the integral term will give rise to a dense matrix. In addition, the singularity of the Lévy measure inside the jump integral may lead to a less accurate numerical approximate solution. (More discussion about the singularity of the jump integral will be given in the next section.) In the literature, the jump process can be classified into finite activity, infinite activity with finite variation, and infinite activity with infinite variation, where the singularity increases from the first to the later cases. In [14], the jump integral was split into local and nonlocal parts with the local term computed using implicit time stepping and the nonlocal term computed using explicit time stepping. In [13,18], an integration by parts technique was used to transform the integral part into a weakly singular Volterra equation. A collocation method was then applied to the CGMY (Carr-Geman-Madan-Yor) process [2] in the case of infinite activity and finite variation. In [16], a wavelet Galerkin finite element technique was applied to an infinite activity process.
In this paper, we are concerned with the numerical solution method proposed by Wang et al. [19]. This approach is able to achieve second order accuracy under an infinite activity and finite variation process, and a better than first order accurate for an infinite variation process. Moreover, it is applicable to different types of options such as European and American, non-constant coefficient PIDEs, nonconstant volatility, and various types of initial and boundary conditions. In each time step, a fixed point iteration is performed to compute the option value. In this approach, an FFT is used to efficiently compute the jump integral and more importantly, it avoids inverting a full dense matrix arising from the discretization of the jump integral. However, the convergence of the fixed point iteration can be very slow. In fact, it can be proved that the convergence rate deteriorates as the number of grid points increases (Y larger than 1) or the singularity of the Lévy measure increases. It will be made more precise in the next section.
Recently, it has been applied to solving American option equations [36]. Multigrid methods have been developed for option pricing in the literature, for instance in [37][38][39]. However these methods do not consider the CGMY or any jump models.
The success of a multigrid technique depends on smoothing and its interaction with the coarse grid correction. Our idea is to apply the fixed point iteration as a smoother in the multigrid framework, rather than to use it as a solver, for solving the PIDE arising from option pricing. We will perform a smoothing analysis and show that the fixed point iteration can effectively reduce high frequency errors. Moreover, using a local mode analysis, we provide a two-grid convergence analysis of the method. Direct discretization is used for the coarse grid matrices. Linear interpolation and full weighting are used for inter-grid transfer.
We remark that the proposed multigrid method can be used with the penalty method for pricing American options. The discretized linear system is almost the same as that of the European option with some of the diagonal entries modified as a result of the early exercise constraint. Due to the similarity, we will focus only on European option in this paper.
In Section 2, we describe the PIDE for option pricing under the Lévy process and derive the linear system resulting from a finite difference method. In Section 3, we present a multigrid method for solving the linear system. Moreover, we provide a smoothing analysis for the proposed smoother and a convergence analysis for the two-grid method. Finally, numerical results are given in Section 4 to demonstrate the fast convergence of the proposed multigrid method by solving PIDEs under a Lévy process with finite activity, infinite activity with finite variation, and infinite activity with infinite variation. Compared with the fixed point iteration which may take tens to hundreds or even thousands of iterations to converge, our proposed multigrid method converges within 10 iterations as the Y parameter of the CGMY model increases from 0 (mildly singular) to close to 2 (highly singular).

Model equation
In Black-Scholes modelling with jumps, the value of the European option price, V, satisfies the following partial integro-differential equation (PIDE) [1]: where r is the risk-free interest rate, q the continuous dividend yield, and σ the volatility. The parameter τ is the time from expiry; i.e., τ = 0 denotes the expiry time T and τ = T denotes the start of the option. The evolution of the underlying risky asset price S is driven by a Lévy process whose Lévy measure ν satisfies |y|<1 y 2 ν(dy) < ∞, |y|≥1 ν(dy) < ∞.
Under the CGMY model [2], ν is defined as follows are the indicator variables. The parameter C > 0 is the measure of the overall level of activity, G ≥ 0 and M ≥ 0 control the rate of exponential decay on the left and right of the Lévy density, and Y < 2 describes the behavior of the Lévy density in the neighborhood of zero where the density tends to infinity. The parameter Y is of interest mathematically and numerically since its value determines the regularity of ν. If Y < 0, the measure ν integrates to a finite value yielding a process of finite activity. If Y ∈ [0, 1], the process displays infinite activity but finite variation since |y|<1 y 2 ν(dy) < ∞. If Y ∈ (1, 2), the process is said to have infinite activity and infinite variation.

Discretization
There have been a number of discretization methods for solving PIDEs [13,14,16,18]. Here, we will consider the method proposed in [19], which is easy to implement and applicable to non-constant coefficients. To be self contained, we present the major steps here and refer the readers to [19] for details. Meanwhile, notations are introduced which will be used in the smoothing analysis in Section 3.1.
The jump integral in (2.1) defined on [−∞, ∞] is first approximated by a finite domain [y min , y max ] with appropriately chosen y min and y max so that the error of approximation is negligibly small. The interval is divided into N subintervals [y j − ∆y/2, y j + ∆y/2], j = 0, 1, . . . , N − 1, where ∆y = (y max − y min )/N. Note that the Lévy measure, ν(y) in (2.3) goes to infinity as y goes to 0, and it is only second moment integrable in general; see (2.2). The computation of the approximate integral is splitted into three parts: I Ω 0 , I Ω 1 and I Ω 2 , corresponding to the three regions , , Notice that ν(y) is smooth on Ω 2 but exhibits singular behaviours on Ω 0 and Ω 1 . Using the discretization techniques introduced in [19] to handle the singularity, the integral can be approximated bȳ andγ(y j ) is the value of the following integral γ(y j ) computed by a quadrature rule: j y j +∆y/2 y j −∆y/2 y 2 ν(y)dy; if y j ∈ Ω 1 y j +∆y/2 y j −∆y/2 ν(y)dy; if y j ∈ Ω 2 0; if y j ∈ Ω 0 .
As mentioned above, ν(y) is smooth on Ω 2 . The quadrature rule can be just any standard numerical integration method such as the composite trapezoidal rule. On Ω 1 , a special quadrature rule is needed to capture the singular behaviour of ν(y). The idea is to choose quadrature weights such that the integral is exact for polynomials of up to degree 3. Thus the accuracy of the quadrature rule is fourth order which is needed to compensate for the 1/y 2 j factor when y j is near 0. Combining (2.1) and (2.4), the equation to be solved can be written as We see that the complex integral is simplified to a discrete sum. A semi-Lagrangian method is used to discretize the first order term. Let [S 1 , . . ., S m ] be a set of grid points for S and τ n be the n-th time step. Also let V n ≈ V(S , τ n ) be the approximate value at (S , τ n ). The semi-Lagrangian method traces node S at τ = τ n+1 along the characteristic trajectory back to τ = τ n at a point S * , which in general is not a grid node S j . Let V n * be the option value at S * , which is approximated by interpolating the nearby values of V n . In order to achieve second order accuracy, a quadratic interpolation Φ n+1 is used and the approximate value is denoted by Φ n+1 V n . Note that the superscript n + 1 for Φ n+1 emphasizes that the interpolation operator depends on time τ n+1 since the location to be interpolated, S * , is derived from S at τ n+1 .
Applying a fully-implicit timestepping with step size ∆τ along the characteristic direction gives: We note that the Crank-Nicolson scheme can also be used to obtain second order accuracy in time.
However, since there is no difference in applying the multigrid method for solving the discrete linear system arising from the two timestepping schemes, we will focus on the fully-implicit discretization. The computation of (QV n+1 ) can be expensive in general. However, it can be accurately approximated by an FFT method [12], and the approximate value is denoted by (BV n+1 ) . Applying central finite differencing to discretize (LV n+1 ) , the fully discrete system can be written as: where α and β are given by , for = 2, . . . , m − 1. When = 1, we set α = β = 0 at S 1 = 0; and when = m, we set V n+1 m equal to the relevant Dirichlet boundary condition.
In matrix form, the discrete system can be rewritten as: where I is the identity matrix, L is defined such that and B is a dense matrix whose entries satisfy the following properties [12] The matrix-vector product BV n+1 can be computed efficiently by using an FFT.

Fixed point iteration
Each time step requires solving the linear system (2.7), where the coefficient matrix, denoted by A, is dense since B is dense. In [19], the authors proposed a fixed point iteration which splits A into The solution is computed by solving iteratively where V n+1,k is an approximate solution to V n+1 after k iterations. This method only involves inverting a sparse matrix A 1 and the matrix-vector product with the dense matrix B, which can be computed efficiently.
It can be proved that the fixed point iteration converges and the convergence rate, ρ, is bounded by [19]: Although the fixed point iteration is convergent, the convergence rate can be very slow. Note that λ is a function of the density function ν(y). It can be estimated by for Y > 0. As a result, as ∆y approaches 0.

Multigrid method
As shown in (2.11), the convergence of the fixed point iteration can be very slow when Y > 0 and ∆y is small. This arises in the infinite activity case for the Lévy process. The situation gets worse when Y approaches 2 in the infinite variation case.
We propose an efficient multigrid method for solving the linear system (2.7). The idea is to apply the fixed point iteration as a smoother in the multigrid context. While the convergence rate of the fixed point iteration can be slow, especially when Y gets close to 2, it turns out that the fixed point iteration is an effective smoother which damps the high frequency.
A V-cycle multigrid method is used to solve the PIDE (2.6). We remark that other cycles such as W-cycle, F-cycle, or full multigrid can also be used. But they are generally more expensive. We find that V-cycle is the more effective option in this case.
In a two-grid setting, the algorithm consists of three major steps. We start with an initial guess for the solution,V k h ≈ V n+1 , on the fine grid. The first step is pre-smoothing. A small number of the fixed point iterations (typically ν 1 = 1 or 2) is applied to compute an approximate solution of (2.7), where S denotes the fixed point iteration given by (2.9), and A h and f h are the fine grid matrix and right-hand side from (2.7), respectively. The role of the smoother is to damp away the high frequency components of the error, thus resulting in a smooth error. The smoothing property is analyzed in Section 3.1. The second step is coarse grid correction. Note that the error, e h , on the fine grid satisfies the equation, where r h = f h − A hṼh is the residual on the fine grid. Since the error is smooth after pre-smoothing, we can solve it accurately and efficiently on the coarse grid with grid size H, where A H is the discretization matrix on the coarse grid and r H = Rr h is the restriction of r h given by a restriction operator R. Here we use full weighting for R; i.e., Note that Dirichlet conditions are used for the boundary points (see [19]). Thus the residual values are zero on the boundary points which will be used in the restriction near the boundary.
After e H is obtained, we update the fine grid solution by where P is a prolongation operator which interpolates the coarse grid error e H onto the fine grid. Linear interpolation is used for P; i.e., The last step is to apply smoothing again, where ν 2 is usually 1 or 2.
In a multigrid setting, the coarse grid equation (3.1) is solved by applying the two-grid algorithm recursively.
The success of the multigrid method hinges on the effectiveness of the smoother. In the next section, we perform a smoothing analysis for the fixed point iteration.

Smoothing analysis
To simplify the analysis, we assume that an evenly spaced grid along the log S coordinate is used together with periodic boundary conditions. By a change of variable x = log S , we define V(x, τ) = V(e x , τ). Then (2.6) can be written as the following PIDE with constant coefficients, where V ⊗γ denotes a correlation product; i.e. (V ⊗γ) i = N−1 j=0V (x i , y j )γ(y j ). The PIDE is discretized similarly as described in Section 2.1 where the unknowns are now V n+1 . More precisely, we have where the index j * is defined such that when applying the semi-Lagrangian discretization, the grid point x j at time τ n+1 is traced back to x j * at time τ n . The value of V n j * ≡ V(x j * , τ n ) is computed by an upwind quadratic interpolation. After rearranging terms, (3.3) becomes where α = (σ 2 +σ)/(2(∆x) 2 ). Then, the corresponding fixed point iteration (2.9) is given in component form by is an approximate solution to V n+1 after k fixed point iterations. Note that the choice of the grid points for the upwind quadratic interpolation depends on the sign of the coefficient for V x in (3.2). Without loss of generality, assuming the coefficient is negative, we have where ψ's are the Lagrangian basis functions and x j * lies between x i−p−1 and x i−p for some integer p. Substituting (3.6) into (3.5), we obtain V n+1,k+1 j The fixed point iteration is used as smoother in our multigrid method. We will estimate the smoothing factor of the fixed point iteration by applying the local Fourier analysis [20] to (3.7).
be the error of V n+1,k . Let C k µ and G µ be the symbols of E k andγ, respectively, where µ is the Fourier mode; i.e. C k µ and G µ are the Fourier transform of E k andγ, respectively [20]. Also, let h = 1/N be the mesh size and i = √ −1. Then (3.7) becomes C k+1 Note the negative subscript of G −µ is a result of the correlation operator. The left-hand side of (3.8) can be simplified as Similarly, the right-hand side of (3.8) can be written as where η j = j−p m= j−p−2 ψ j,m e −( j−m)µπhi . LetŜ µ be the symbol of the smoother. Then one can easily show thatŜ µ = C k+1 µ /C k µ ; i.e., It has been proved that |η j | ≤ 1 [40, Section 5.2]. Furthermore, since G µ is the symbol ofγ, by (2.5), we have |G −µ | ≤ λ. Hence, we can boundŜ µ as |Ŝ µ | = |∆τG −µ ||η j | |1 + (r + λ + 4α sin 2 (µπh/2))∆τ| ≤ λ∆τ 1 + (r + λ + 4α sin 2 (µπh/2))∆τ . (3.11) Let ζ be the smoothing factor of the fixed point iteration, which is defined as Considering the high frequency Fourier modes; i.e., N/2 ≤ µ ≤ N, it follows from (3.11) that To bring more insight into the formula, we estimate the upper bound for ζ by noting that Considering the limit where ∆x ≈ 0 and keeping ∆τ fixed, we have In the Lévy process under the CGMY model, we always have Y < 2. Hence for small ∆x, the smoothing factor for the fixed point iteration is small, which is fundamental to achieve an effective multigrid method.

Two-grid convergence analysis
In the previous section, we have shown by a local Fourier analysis that the fixed point iteration smoother is effective damping high frequency errors. In this section, we will perform a two-grid convergence analysis based on the local Fourier analysis. For easy exposition, we assume there is no post-smoothing.
The iteration matrix of the two-grid iteration can be written as where L h and L H are the fine and coarse grid operator, respectively, P is the prolongation, R is the restriction, and S is the smoother. For two-grid analysis [21], it is customary to pair up the low and high frequency modes such that the Fourier transformed matrix,M, is a block diagonal matrix with each block a 2 × 2 matrix,M µ , which is given bŷ Here c µ = cos(µπh)/2, s µ = sin(µπh)/2, and h = ∆x. By (3.10), we haveŜ (3.14) whereσ = (σ 2 +σ)/2. Also, we have assumed the interest rate r to be zero, just to simplify the algebraic expression.
A similar calculation yieldŝ where η j is defined similarly as η j on the coarse grid. Details can be found in Appendix A.

Numerical results
In this section, we present results of numerical experiments to demonstrate the effectiveness of the proposed multigrid method for pricing European options under the CGMY model. The PIDE is discretized by the semi-Lagrangian method in the spatial domain and implicit Euler's method for time. The discrete linear system of equations is solved by the proposed multigrid method. We also include numerical results using a fixed point iteration for comparison.
In the multigrid procedure, a V-cycle is used with one pre-and one post-smoothing where the smoother is a fixed point iteration. The coarse grid operators are obtained from direct discretization with linear interpolation and full-weighting restriction for intergrid transfer. The stopping criterion is when the relative residual l 2 -norm was less than 10 −8 .
As explained in Section 2, the fixed point iteration (cf. Section 2.1) used in the literature converges slowly when the grid size is small and becomes worse for Y ≈ 2. We will show the multigrid convergence results where Y changes from 0 to 1.98.
Example 1: In this example, Y = 0. In this case, the CGMY process is also known as variance gamma. The input parameters for the PIDE and CGMY model are given in Table 1.  Table 2 presents the convergence of the proposed multigrid method (MG) for different grid sizes. In this example, the total number of time steps is 500, and so ∆τ = 10 −3 , which is typical for pricing option values with reasonable accuracy. Multigrid converges very fast in only 2 iterations. The number of fixed point iterations (FP) is relatively constant which shows that the convergence is independent of the grid size under a Variance Gamma process, as explained in (2.11). Example 2: In this example, Y = 0.6442. When 0 < Y < 1, the Lévy process has infinite activity but finite variation. The input parameters are given in Table 3. The convergence results are shown in Table 4. In this case, the number of fixed point iteration increases for smaller grid size, supported by the formula in (2.11) with Y > 0. The multigrid convergence remains constant, independent of the grid size. In this example, we also use a larger ∆τ to test its effect on the convergence. In practice, a large time step size may be necessary when pricing options with long expiry dates; e.g., T = 10. However, a larger time step size will generally result in a more ill-condition discrete linear system. Table 4 shows that the fixed point iteration numbers increase dramatically, by a factor between 5 to 7. The number of multigrid iteration only increases mildly.
Example 3: When 1 < Y < 2, this corresponds to a Lévy process with infinite activity and infinite variation. Table 6 shows the convergence results for pricing a call option for Y = 1.4 and Y = 1.98, with input model parameters given in Table 5. When Y = 1.4, the discrete problem becomes much more difficult to solve, as can be seen from the number of fixed point iteration which increases from hundreds to thousands as the grid size decreases. The multigrid iteration remains very efficient with mesh independent convergence rate. The case where Y = 1.98 is numerically challenging since the kernel function of the Lévy process becomes very singular. The fixed point iteration converges very slowly. The multigrid iteration deteriorates slightly from 4 to 9, but still converges very quickly.
Example 4: In this example, we illustrate the convergence behaviour of multigrid, fixed point iteration and BiCGstab with the fixed point iteration as preconditioner [19]. The results are shown in Figure 1.
Here we consider the case where Y = 1.4. The convergence behaviour is similar for other values of Y and they are not shown here. We can see that multigrid converges quickly with a constant rate and the convergence is essentially unchanged for different grid sizes. Fixed point iteration converges at a much slower rate and the rate of convergence deteriorates as the grid size h decreases. The convergence of preconditioned BiCGstab is somewhere between the two. It is much better than fixed point iteration but not as good as multigrid. Also, the convergence is not as smooth and uniform as the other two, and it is mesh dependent. Table 7 shows the rates of convergence for multigrid and fixed point iteration, which are computed by taking the average rates of the last 3 iterations. The convergence of BiCGstab is not very uniform and we do not include its rate of convergence here.
In Sections 3.1 and 3.2, we presented a local Fourier analysis to analyze the smoother and the two-grid convergence. It would be interesting to show how the theoretical analysis results compare with the actual convergence. However, both the smoothing factor and two-grid convergence factor include a quantity η j (3.9) which comes from the interpolation operator during the semi-Lagrangian discretization. This quantity is difficult to compute and it varies at different grid points, grid sizes and time steps. As such, we were unable to compute the smoothing factor and the spectral radius ofM. Instead, we compute an upper bound for the smoothing factor and the infinity norm ofM as described in (3.11), (3.19) and (3.20). They are shown in columns 3 and 4 in Table 7. Since they are upper bounds, they do not match very well with the actual convergence of multigrid. Nevertheless, we can see that the smoothing factor decreases with h as discussed in Section 3.1 and the two-grid factor is bounded by 0.75 as shown in Theorem 3.2. We also remark that the discrepancy may be due to the fact that the underlying PIDE is a nonlocal operator where the local mode analysis may not be the right tool for analysis.

Conclusion
In this paper, we have proposed a fast multigrid method for solving the discrete equations from a finite difference discretization of the partial integro-differential equations arising from pricing European options under the CGMY process. The fixed point iteration originally proposed in [19] converges well when the density kernel is not singular (Y ≈ 0) but very slowly when the kernel becomes more singular (Y ≈ 2). While the fixed point iteration is not an efficient solver by itself, we have proved that it is effective removing high frequency error components. By using the fixed point iteration as a smoother for multigrid, we have obtained a fast multigrid solver which converges within 10 iterations for different values of Y ranging from 0 to 1.98. In addition, we have provided a two-grid convergence analysis for the multigrid method.
The current paper focuses primarily on European options. For American options under the CGMY process, the corresponding PIDE is nonlinear with a linear complimentarity condition. The proposed multigrid method would not be directly applicable since it is designed for solving linear systems. An interesting idea is to combine the techniques in this paper with a nonlinear multigrid method such as Full Approximation Scheme (FAS), which will be a future work. Another future work is to perform better local mode analysis of the method, such as three-grid analysis, in order to obtain improved quantitative results.


extfieldsofstudy:
Computer Science
Medicine

text:
Stacking Ensemble Technique for Classifying Breast Cancer

Objectives Breast cancer is the second most common cancer among Korean women. Because breast cancer is strongly associated with negative emotional and physical changes, early detection and treatment of breast cancer are very important. As a supporting tool for classifying breast cancer, we tried to identify the best meta-learner model in a stacking ensemble when the same machine learning models for the base learner and meta-learner are used. Methods We used machine learning models, such as the gradient boosted model, distributed random forest, generalized linear model, and deep neural network in a stacking ensemble. These models were used to construct a base learner, and each of them was used as a meta-learner again. Then, we compared the performance of machine learning models in the meta-learner to determine the best meta-learner model in the stacking ensemble. Results Experimental results showed that using the GBM as a meta-learner led to higher accuracy than that achieved with any other model for breast cancer data and using the GLM as a meta learner led to low root-mean-squared error for both sets of breast cancer data. Conclusions We compared the performance of every meta-learner model in a stacking ensemble as a supporting tool for classifying breast cancer. The study showed that using specific models as a metalearner resulted in better performance than single classifiers, and using GBM and GLM as a meta-learner is appropriate as a supporting tool for classifying breast cancer data.

there is a difficulty in accurate diagnosis since fine-cell aspiration is affected by various external variables, such as the small size of cancer, noninvasive breast cancer, and low stage [7]. In addition, diagnosis can be inaccurate due to various factors, such as the skills of the examiner and the similarity between organisms [8,9]. To overcome these limitations, it is necessary to provide an auxiliary tool to diagnose breast cancer accurately.
Machine learning, which is a field of artificial intelligence technology that makes machine judgments through data learning, has been widely used in industries such as medicine, manufacturing, and finance [10]. Machine learning is largely divided into supervised learning and unsupervised learning. Supervised learning is a method that enables a model to deduce correct labels for new data by learning data with labels. On the other hand, unsupervised learning refers to a learning technique in which data without labels is learned, and data are clustered according to the similarity of their characteristics. The common characteristics that supervised and unsupervised learning have are the importance of the data to use and the importance of choosing a prediction model for learning. There is a problem that if the label of the data is biased, or the selected model is over-fitted with the corresponding data, the model cannot show the correct performance. The machine running technology that emerged to solve these problems is ensemble techniques. Ensemble techniques are largely divided into three categories, namely, bagging, boosting, and stacking. Bagging a method of allowing redundancy when randomly classifying training data to a prediction model in the process of learning. Boosting is similar to bagging, but it reuses the data that the prior model cannot classify well in the learning process by providing weight . Stacking is an ensemble technique that is also known as super learning. It uses various models, such as random forest, k-nearest neighbors, and support vector machine as base learners to generate new data based on predicted results. Then this new data is used for another predictive model, which is called meta-learner, to finally derive the predicted value [11]. Studies on ensemble techniques have been steadily progressing, and many of them have been limited to the techniques of bagging and boosting ensembles [12][13][14]. If bagging and boosting ensemble techniques are used to generate strong predictive models using the same incomplete and weak prediction models, then staging ensemble techniques are available to combine various predictive models [15]. Many previous studies have introduced new stacking ensemble techniques for breast cancer prediction and to increase the performance. As prior research on stack-ing ensembles, research to determine the best meta-learner model is needed as when the same machine learning models are used in the base learner and meta-learner. Therefore, we compared the performance of meta-learners and identified the best meta-learner model in a stacking ensemble as a predictive aid for breast cancer patients.

Data Source
The Wisconsin Breast Cancer Diagnostic (WBCD) data provided by the University of Wisconsin consisted of data from 569 consecutive fine-needle aspiration tests, with 10 variables representing the characteristics of the nucleus, each with mean, standard deviation, and ideal values [16]. For example, each variable radius contains three variables, namely, the mean radius (mean radius of cell), radius SD (standard deviation of cell), and worst radius. Therefore, the WBCD consists of 30 variables, except for diagnostic results that indicate benign and malignant.
The Wisconsin Breast Cancer -Original (WBCO) data contained nine features and one diagnostic value for 699 breast biopsies, like the WBCD data.
For the WBCD data, each parameter value has different numerical measurement values. For this study, other numerical measurement values except for diagnosis results and ID were normalized by using the following equation: Unlike WBCD data, WBCO data has the same range of numerical measurement values and is not normalized. Only 10 rows with missing values are removed. In this study, the WBCO data and the WBCD data were randomly extracted at a ratio of 8:2 to train data and test data, so that the values were not duplicated. Training data was used to train the model, and the performance was evaluated using the test data.

Model Development
The algorithm used for the stacking ensemble in this study consists of four models. www.e-hir.org

Stacking Ensemble for Breast Cancer
that compensates for the error of the previous tree and is expressed as follows: Here, d n (x) is generated by minimizing the error of d n-1 (x) with each sequentially generated tree.
2) Distributed Random Forest Distributed Random Forest (DRF) is an ensemble model used to reduce over-fitting while maintaining predictive performance by combining decision trees with slightly different variable values. DRF extracts features randomly with the bootstrap technique and generates several decision trees with corresponding values. The average of the predicted results in each decision tree is used as the final prediction result, and DRF has excellent performance even when the parameter is used as the default value.

3) Generalized Linear Model
Generalized Linear Model (GLM) is mainly used when the dependent variable is a categorical variable. The dependent variables used in this study are categorical variables (benign and malignant) and binomial variables at the same time, so it is appropriate to use the model that is expressed as follows: where E is the expected value, and β 0 ,β 1 ··· β n are the regression coefficients. Here, y is the probability that an observation belongs to a group, x is the probability of not belonging to the group, and the higher the E value, the higher the probability of belonging to the group.

4) Deep Neural Network
Deep Neural Network (DNN) is a neural network with multiple hidden layers between the input layer and the output layer. DNN is a neural network structure that can produce different results depending on the number of layers and the number of nodes in each layer.

Experiment Environment and Model Architecture
This study was implemented using the H2O package of the R program (ver.3.4.3). The H2O package is an open-source library from H2O.ai, which provides a variety of supervised and unsupervised learning algorithms [17]. It also supports various languages, such as R, Python, and Java. Most of the prediction models in this study were implemented using the basic parameter values of the H2O package, and 5-fold cross-validation was performed to ensure the reliability of the data. As shown in Figure 1, base learner models were implemented using h2o.gbm, h2o.randomForest, h2o.glm, and h2o.deeplearning functions, and a true value is given to determine whether to keep cross validation as the parameter value of the corresponding functions.
First, h2o.gbm has a default number of 50 trees with a maximum depth of 5. h2o.randomForest has a default number of 50 trees with a default depth of 20. h2o.glm has a family parameter Gaussian as a default value, but in this study, a binomial was used because the dependent variable is binomial distribution. By default, h2o.deeplearning has two hidden layers and 200 nodes, the epoch is 10, and the activation function is the ReLU function. The meta-learner was implemented using the h2o.stackedEnsemble function, and the functions used in the base-learner model were used once again.

III. Results
In this study, we compared the accuracy and root-meansquare error (RMSE) of each prediction model for performance evaluation. In performance evaluation, accuracy indicates the accuracy of the prediction model, and the RMSE is a value calculated by comparing the predicted value with the observed value to determine how much error there is. Therefore, the lower the RMSE is, the better the performance is. To guarantee the reliability of the results, we randomly extracted the data and distributed it to the training data and the test data by randomly repeating this process five times and compared the average values of accuracy and RMSE. Figures 2 and 3 compare the accuracy and RMSE values of each model with the accuracy and RMSE values obtained when a stacking ensemble was used. In Figures 2, 3 and Tables 1, 2, 'single classifier' means a prediction model using single GBM, DRF, GLM and DNN. 'Ensemble classifier' means a prediction model using the GBM, DRF, GLM, and DNN models as the base learner, and each model is used as a meta-learner again. Figure 2 shows the results of each prediction model for WBCD data. Figure 2 shows that the ensemble model for WBCD data achieved better prediction accuracy than the single models used. Also, the RMSE value of the stacking ensemble using DRF, GLM, and DNN was lower than that of the conventional single model. Table 1 shows the average prediction accuracy and average RMSE of each model used to draw the graph in Figure  2. Table 1 shows that the ensemble model achieved higher accuracy than the single model. In particular, the ensemble model using GLM as a meta-learner showed the highest average prediction accuracy of 97.37%. Also, the ensemble model using GLM as the meta-learner showed the lowest RMSE value, 0.1873. In contrast, the single GLM model

Stacking Ensemble for Breast Cancer
showed the highest average RMSE value of 0.2206. According to Figure 3, the WBCO data also shows that each model achieved high performance overall, but other ensemble models, except GBM and DNN, showed lower average predictive accuracy rates than those using a single model. The RMSE for the ensemble model using GBM and GLM as the meta-learner was lower, while the RMSE value for the ensemble model using DRF and DNN as the metalearner was higher than that for the single model. Table 2 shows the average prediction accuracy and average RMSE of each model used to draw the graph, as shown in Table 1. As seen in Table 2, only the ensemble model using GBM and DNN as the meta-learner showed an increase in average accuracy. On the contrary, the model with single GLM had the highest average prediction accuracy of 98.10%. In the case of average RMSE, the RMSE value was lower only when GBM and GLM were used as a meta-learner, and the model using GLM as a meta-learner showed the lowest RMSE value.

IV. Discussion
Unlike in the West, breast cancer in Korea has a very high incidence in people under 50 years old, rather than after 50 [18]. This means that the age which breast cancer can occur is lower than in the West. Therefore, early diagnosis and accurate prediction of breast cancer is needed through convergence with various fields. Recently, various studies have been conducted to improve the accuracy of cancer prediction using machine learning techniques [19]. In particular, research has focused on ensemble techniques that improve performance by combining various prediction models rather than one of the machine learning techniques [20]. In this study, we used a stacking ensemble technique that can combine various prediction models among ensemble techniques to generate prediction models and verified their performance. In this study, a stacking ensemble was constructed using GBM, DRF, GLM, and DNN, which are powerful prediction models, and each prediction model was used again as a meta-learner. Since GBM, DRF, GLM, and DNN were predictive models that achieved high performance even if they were used alone, the same performance improvement could not be obtained in all models when the stacking technique was used. However, in this study, the stacking ensemble using GBM as a meta-learner showed better prediction accuracy than the existing single model for both sets of breast cancer data. In addition, the RMSE value was confirmed to be the metalearner, and the stacking ensemble using the GLM decreased in different data. The results of this study demonstrated that a stacking ensemble using GBM and GLM as a meta-learner is suitable as an auxiliary tool for predicting breast cancer.
This research identified the best performance meta-learner model when the same models were used in the base learner and meta-learner. The result of this research could be a valuable reference for choosing machine learning models for further stacking ensemble research. However, this study had a limitation that it was based on the default parameter values provided by H2O; we expect that it would be possible to develop a stacking ensemble with higher performance by


s2fieldsofstudy:
Computer Science
Physics

text:
Computing 3 point correlation function randoms counts without the randoms catalogue

As we move towards future galaxy surveys, the three point statistics will be increasingly leveraged to enhance the constraining power of the data on cosmological parameters. An essential part of the three point function estimation is performing triplet counts of synthetic data points in random catalogues. Since triplet counting algorithms scale at best as $\mathcal{O}(N^2\log N)$ with the number of particles and the random catalogues are typically at least 50 times denser than the data, this tends to be by far the most time consuming part of the measurements. Here we present a simple method of computing the necessary triplet counts involving uniform random distributions by means of simple one dimensional integrals. This speeds up the computation of the three point function by orders of magnitude, eliminating the need for random catalogues, with the triplet -- and pair -- counting of the data points alone being sufficient.

INTRODUCTION
Current and past galaxy redshifts surveys have heavily relied upon the analysis of two point statistics in order to constrain cosmological parameters down to the percent level (see e.g. Anderson et al. 2012;Chuang et al. 2016;Cuesta et al. 2016;Gil-Marín et al. 2016a,b;Alam et al. 2017). As upcoming galaxy redshift surveys seek to push constraints on cosmological parameters to the sub-percent level, three points statistics -the three point correlation function (3PCF) and the bispectrum -will begin to play a bigger role in analyses. Recent works have shown that the baryon acoustic oscillation features can be detected in both the 3PCF (Slepian & Eisenstein 2015;Slepian et al. 2017a,b) and the bispectrum (Pearson & Samushia 2018, hinting at the possibility of increased constraining power on cosmological parameters via the inclusion of three point statistics. The large scale bispectrum has also been used to supplement the redshift-space distortion measurements from the power spectrum (Gil-Marín et al. 2015a,b, 2017. Turning to small scales, the Halo Occupation Distribution (HOD) model (Scherrer & Bertschinger 1991;Jing et al. 1998;Peacock & Smith 2000;Seljak 2000;Scoccimarro et al. 2001;Berlind & Weinberg 2002;Zheng et al. 2009, and references therein) is a popular ⋆ E-mail: dpearson@phys.ksu.edu method of linking the galaxy and dark matter distributions (see e.g. Tinker et al. 2005;Guo et al. 2014Guo et al. , 2016Rodríguez-Torres et al. 2016). Yuan, Eisenstein & Garrison (2018) have shown that the squeezed 3PCF can help tightly constrain the parameters of the HOD, making it very likely that the 3 point statistics will become important to the development of mock catalogs.
To measure the 3PCF from data, one must count the triangles of specific shapes and sizes from the data -e.g. the galaxy catalogue -as well as a set of unclustered random points -the random catalogue -a process which is naively O(N 3 ) in time complexity. The counts from the data are then compared to the expected mean numbers from the unclustered random points to get an estimate of the 3PCF (Peebles & Groth 1975;Peebles 1980), as the unclustered mean triplet counts are very sensitive to the geometry of the survey and its number density variations. The most popular estimator used for the 3PCF is that of Szapudi & Szalay (1998) due to its superior edge effect corrections (Kayo et al. 2004), [DDD(r 1 , r 2 , r 3 ) − 3DDR(r 1 , r 2 , r 3 ) where the combinations of Ds and Rs tell you how many vertices of the triangle come from either the data or the random points, respectively. The random catalogue usually contains 50 or more times as many objects as the data to make the shot-noise in this Monte-Carlo estimation subdominant to the variance in the data. Combine this with the O(N 3 ) complexity and the counts for DDR, DRR, and RRR can consume a large amount of time.
While the computational complexity of the three point statistics can be somewhat mitigated (Baldauf et al. 2015;Slepian & Eisenstein 2015;Pearson & Samushia 2018), their calculation still tends to take a significant amount of CPU time. Additionally, the studies which stand to benefit from the inclusion of the three point statistics may need to measure them tens of thousands of times, making any potential reduction in the computational complexity a welcome improvement.
In this letter we present a simple method of obtaining the counts involving the random points that does not require actually counting triangles from random catalogues. Since triplet-counting from random catalogues takes most of the computational time, this results in a significant reduction of required CPU hours. The method is restricted to measurements of 3PCF from uniform periodic cubes, which means that, unfortunately, it cannot be applied to the measurements from real survey data. There are, however, many stages in the cosmological analysis of the 3PCF that can be performed on uniform periodic cubes and do not need to account for survey geometry. These include the validation of theoretical templates -comparing model predictions to measurements from simulations and quantifying biasesand the HOD parameter fitting. For these applications our method allows one to set the number density of randoms arbitrarily high with no performance degradation.

Pair counts
We start by reminding the reader why random pair counting is not really necessary when computing the two-point statistics from a preiodic cube. The number of unclustered pairs separated by a distance r ±∆r/2 can be easily computed analytically. There are on average n R V box particles in a uniform periodic cube with volume V box and a number density of points n R . The volume between two concentric spheres of radii r − ∆r/2 and r + ∆r/2 -i.e. a spherical shell -is and the total number of pairs is since The cube is periodic there are no edge effects. We will now generalize this idea for triplet counts.

RRR counts
We need to estimate number of triplets separated by (r 1 ± ∆r/2, r 2 ± ∆r/2, r 3 ± ∆r/2) in a uniform periodic cube with a number density of n R and volume V box . We will start by computing this number for (r 1 , r 2 ± ∆r/2, r 3 ± ∆r/2). This can be achieved by finding the volume of the two overlapping In essence, we want to predict how many points, p 3 , there will be, which is simply the volume of the overlap region of the two spherical shells with central radii r 2 and r 3 , multiplied by the number density.
spherical shells whose cross-sections are show in Figure 1.
Multiplying that by the number density will then tell us the average number of p 3 points given p 1 and p 2 . We will make use of the equation for the volume of overlap of two spheres with radii r and R separated by a distance d (Weisstein 1999), This equation is valid only from the point where the two spheres just touch, up to the point where the smaller sphere completely falls within, but still touches the larger sphere, i.e. R − r d R + r. Because of this, we need to explicitly define the volume outside of those bounds, From careful study of Figure 1, we can find that the volume of interest is given by first finding the overlap volume of the two outer spherical surfaces, then subtracting the overlap volume of the outer spherical surface of one shell with the inner spherical surface of the other and vice versa. However, this ends up removing the overlap volume of the two inner spherical surfaces twice, so we have to add one back. Mathematically, this can be expressed as V cs (r 1 , r 2 , r 3 , ∆r) = V(r 1 , r 2 + ∆r/2, r 3 + ∆r/2) − V(r 1 , r 2 + ∆r/2, r 3 − ∆r/2) − V(r 1 , r 2 − ∆r/2, r 3 + ∆r/2) a.) b.) Figure 2. Examples of the r 2 -r 3 spherical shell overlap region in two special cases which are automatically accounted for by equation (5) but would not work with (4). Panel a shows a case where r 1 is less than ∆r. Panel b shows the cross-section for nearly co-linear triangles.
Here, care must be taken that the r of equation (5) is actually the smaller of the two spherical shell radii. For example, when considering isosceles or equilateral triangles, it is possible that r 2 +∆r/2 is larger than r 3 −∆r/2, since r 2 may equal r 3 . Exercising this caution, along with the special considerations of equation (5) will yield the correct volume even in special cases where the overlap volume looks quite different than in Figure 1 -see e.g. Figure 2. This volume times n R is the number of p 3 points falling into the overlap region. For a finite bin width, p 2 can be anywhere inside a spherical shell r 1 ± ∆r/2 around p 1 so we will have to integrate equation (6) with respect to 4πn R r 2 1 dr 1 . Finally, we have to account for the fact that there are on average n R V box points p 1 . Combining all of this results in an exact expression for the expected number of RRR counts without any shot-noise, expressed as a simple one dimensional integral, Here s is the number of unique permutations of the side lengths, with s = 1, 3, and 6 for equilateral, isosceles, and general triangles, respectively. We note that for most triangles this reduces to a surprisingly simple expression, However, this expression will break down in special cases such as those shown in Figure 2. For this reason, we recommend simply evaluating the integral in equation (7) numerically. 1

DRR counts
For DRR counts we have to replace the first vertex by a data point. Since data and random points are uncorrelated 2 between each other the formula for the DRR counts coincides 1 The terms in equation (6) will at most be polynomials of degree 5 in r ′ 1 , meaning a simple 3 point Gaussian quadrature rule is all that is needed. 2 This would not be the case for nontrivial survey geometries where the exact placement of data points with respect to the boundary makes a difference.
with the expression for RRR in equation (7) with one minor change due to the potentially differing number densities of data and random points where n D is the average number density of data points -i.e. n D = N D /V box .

DDR counts
The predictions of DDR(r 1 , r 2 , r 3 ) are almost identical to the procedure outlined in section 2.2, except for the last step where we integrate over 4πn R r 2 1 dr 1 . Since the data points are clustered, the distribution of p 2 around p 1 is not uniform and the number density depends on the distance r 1 resulting in Here n D (r) is the nonuniform data number density, which can easily be found from the data pair counts by first computing the number density in each spherical shell and then linearly interpolating 3 to the value at the particular r needed for the numerical integration. Note that we can no longer use a simple integer multiple to account for the unique permutations of the side lengths as the value of n D (r) will change as we permute. Given this, there will be 1, 3 or 6 terms between the braces for equilateral, isosceles or general triangles, respectively. We note that when summing these terms it is necessary to store the number of triangles as a floating point number to ensure accuracy. After the summing of all the necessary terms in the braces, the number can be rounded to an integer if desired.
Equations (7), (9), and (10) express all terms in the estimator of equation (1) that involve random points without any shot-noise and in terms of simple one dimensional integrals. DDR counts do require the pair counts of data points but these can be precomputed along with the DDD counts with no additional computational time. We have publicly released our code for doing the predictive calculations for use by the community. 4 . Differences in the average actual and predicted counts divided by their standard deviations added in quadrature for the cases of 1× and 10× randoms. Note the difference in vertical scale between the left and right hand panels. The triangle index range was truncated from 3128 to 1452, though we note that the remainder of the range looks very similar to what is seen after a triangle index of ∼600. It is clear that the RRR (top) and DRR (middle) predictions agree quiet well with the predictions for both the 1× and 10× cases. The DDR (bottom) counts agree quiet well for the vast majority of triangle indices, though there are notable deviations at lower triangle indices. We note that these effects are not systematic and decrease with increased number density of randoms, as well as with increased grid resolution in the generation of the lognormal mocks.

RESULTS
To test equations (7), (9) and (10) we compared their predictions with actual triplet counts. We generated 100 lognormal mock catalogues (Coles & Jones 1991) -N D ∼ 500000, V box = (1024 Mpc) 3 -to serve as a proxy for data points 5 and 100 random distributions with the same number density, and 100 with 10× the number density of data. We then performed the direct counts in bins of width ∆r = 1 Mpc, from r min = 0 Mpc to r max = 32 Mpc, and calculated our predictions. We then averaged together all 100 counts, and separately all of the predictions in order to make our comparisons.
To perform the actual counts we used a relatively straightforward, GPU accelerated, O(N 3 ) counting algorithm with periodic boundary conditions to remove any edge effects along with a few simple optimizations. 6 The first optimization employed is to bin the particles to a grid with spacing equal to r max , the maximum separation considered 5 While lognormal mocks do not adequately reproduce the threepoint clustering statistics observed in simulations or the real Universe, for the purposes of these tests we simply needed data with some clustering and are not concerned with the recovered clustering statistics. As we already had many lognormal mocks at our disposal, we utilized some of them out of convenience.

6
While we do not release our exact code used here publicly, we utilized the same algorithm in a library for use in a different project. For those who may be interested, you can view the GPU implementation at https://github.com/dpearson1983/ganpcf/blob/master/source/ganpcf.cu.
for any of the sides of the triangle. This allows us to check for triangles with points that are within r max of our first vertex by only checking points in the grid cell containing the vertex, and those immediately adjacent to that cell. The other, more minor optimization is to simply check that the distance between the first and second vertex is less than r max before checking for the third vertex, and if not, proceed to the next point for consideration as the second vertex. While it would usually be faster to not repeat count the same triangle by simply relabeling which is the first, second and third vertex, due to the peculiarities of GPU memory accesses, we found that it was actually more efficient to do the repeat counting. Figure 3 shows the difference of the average RRR (top), DRR (middle) and DDR (bottom) counts with the predictions from equations (7), (9) and (10) divided by their standard deviations added in quadrature, with the 1× and 10× randoms cases on the left and right, respectively. For clarity, the triangle index range was truncated. We note that in all cases the remainder of the range looks nearly identical to what is seen after a triangle index of ∼600. The predictions match the actual counts remarkably well, especially as the number of random points increases.
The DDR counts tend to show larger deviations between the direct counts and our predictions at small triangle indices -i.e. small scales. We verified that these deviations are not systematic -i.e. they will fluctuate up and down with equal probability for independent realisations averaging to zero -and that they tend to zero in the limit of high number density. They look systematic because at small triangular indexes the measurements share the same small number  (top) show that the grid spacing is partly to blame for the deviations. The Gaussian haloes (bottom) do not use a grid for their creation, and have an exactly known number density profile, instead of one estimated from the data-data pair counts. We see that there are no seemingly systematic deviations at small triangle indices. We note that there are some slight deviations at high indices, particularly if a finer binning is used, suggesting that the deviations are have more to do with the discrete nature random catalogues and grid effects.
of points and are strongly correlated -or anti-correlated depending on the shapes. Additionally, the deviations all occur when one side of the triangle is smaller than 5 Mpc, or 3 Mpc for the 1× and 10× randoms cases, respectively. The lognormal mocks were created using a grid spacing of 2 Mpc, where a number of galaxies was determined for each grid cell, then placed inside that cell in a uniform random fashion. We ran tests on mocks with a grid spacing of 0.5 Mpc but the same number density, and note that the deviations were reduced -see Figure 4. To further test that these small scale deviations were not due to a failing of the algorithm, but rather some combination of low number densities and grid effects in the mocks, we generated haloes with an exactly known, Gaussian number density profile, eliminating the need to use the data-data pair counts, making the predictions from equation (10) far more accurate. The random catalogues where set to the same number density -e.g. N D ∼ N R . We then only counted triangles which used the central particle as one of the vertices. The number density profile was such that the maximum particle radius was ∼25 Mpc with very few particles in the outer parts of the haloes. This led to seeing a similar -though significantly smaller -seemingly systematic effect for high triangle indices.
The results of the using a higher resolution grid for creating the lognormal mocks and the Gaussian haloes are shown in Figure 4. These results suggest that the seemingly systematic deviations in the DDR counts versus predictions have more to do with the discrete nature of the random catalogues, and grid based effects than a failing of the predictions.

CONCLUSIONS
We present a method for predicting the counts involving a random catalogue for the 3PCF analysis of simulated or mock data that is free of shot-noise and does not actually require random catalogues. We have shown that the predictions from equations (7), (9), and (10) agree remarkably well with the actual counts, while keeping the same computational complexity for arbitrarily high number densities of randoms.
The method only works for uniform periodic cubes which may lead the reader to believe that the method is not useful for the analysis of real survey data. This, however, is not true as there are many stages in the survey data analysis that require computing the 3PCF from a large number of periodic cubes.
The most obvious example is the validation of theoretical templates on N-body simulations. To make sure that theoretical predictions of the 3PCF are sufficiently accurate -and to calibrate any systematic effects in case they are not -they will have to be compared to the measurements from N-body simulations with a known cosmology, or a hidden cosmology in the case of data challenges. For these validation tests to be meaningful they need to be performed on periodic cubes to cleanly separate theoretical systematics from other observational effects. Since the cumulative volume of N-body simulations used for this purpose needs to be much larger than the size of the data, the method presented in this paper can save a significant amount of computational time.