Yuhao Huang

Shitao Chen

Yu Chen

Zhiqiang Jian

Nanning Zheng
email: nnzheng@mail.xjtu.edu.cn

Spatial-Temproal Based Lane Detection Using Deep Learning

Keywords: Lane detection, Inverse perspective transform, Convolutional neural network, Spatial-temporal correlation

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

In recent years, autonomous driving has received widespread attention. Identifying traffic scenes around autonomous vehicles is a key step in autonomous driving. Besides, the identification of traffic markings in traffic scenes, especially lane boundaries, provides the necessary information for decision making in autonomous driving modules like path planning. We all know that human drivers rely mainly on vision when driving cars and the vision sensors are low-cost in the market. Therefore, the vision-based lane detection algorithm is the mainstream direction of lane detection. While there already exist plenty of vision-based lane detection algorithms, several tough challenges remain to be handled in practice. Firstly, poor lighting conditions makes it difficult to recognize lane markings at night or in other bad weather conditions like rain, snow, fog, etc. Secondly, complex traffic environments, including wear and tear of lane boundaries, congestion of urban roads and lane boundaries being partially obscured, have increased the difficulty of lane detection. The last is the boundary type recognition. Compared to Advanced Driver Assistance Systems (ADAS), autonomous driving requires not only accurate lane boundary location but also precise boundary type classification.

Lane boundary detection is generally divided into three parts [START_REF] Bar | Recent progress in road and lane detection: a survey[END_REF], [START_REF] Sandipann P Narote | A review of recent advances in lane detection and departure warning system[END_REF]: preprocessement of the image and segmentation of the region of interest (ROI), extraction of lane boundary features, and finally lane fitting. Traditionally, lane detection has been performed by using hand-crafted features like gradient information and intensity change information [START_REF] Labayrade | A reliable and robust lane detection system based on the parallel use of three algorithms for driving safety assistance[END_REF], [START_REF] Nieto | Robust multiple lane road modeling based on perspective analysis[END_REF], [START_REF] Samadzadegan | Automatic lane detection in image sequences for vision-based navigation purposes[END_REF]. These algorithms achieve high accuracy on highways with well-lit scenes while having a comparatively bad performance on finding lane markings when lane boundaries are obscured or the experiments are run in bad weather conditions . This is because that lane boundaries have different features in different traffic scenarios, while the hand-crafted features cannot meet the requirements of multiple scenes at the same time.

With the development of artificial intelligence, deep learning is widely used in various fields. AlexNet [START_REF] Krizhevsky | Imagenet classification with deep convolutional neural networks[END_REF], VggNet [START_REF] Simonyan | Very deep convolutional networks for large-scale image recognition[END_REF], GoogleNet [START_REF] Szegedy | Going deeper with convolutions[END_REF] and ResNet [START_REF] He | Deep residual learning for image recognition[END_REF] achieved good results in the ImageNet Challenge from 2012 to 2015 respectively. CNNbased object detection methods like Faster R-CNN [START_REF] Shaoqing Ren | Faster r-cnn: Towards real-time object detection with region proposal networks[END_REF], YOLO [START_REF] Redmon | You only look once: Unified, real-time object detection[END_REF] and SSD [START_REF] Liu | Ssd: Single shot multibox detector[END_REF] have gained high accuracy in object detection task. Some algorithms utilized CNN as an feature extractor to identify lane boundaries and also achieved high accuracy in data sets of different conditions [START_REF] Lee | Vpgnet: Vanishing point guided network for lane and road marking detection and recognition[END_REF], [START_REF] Pan | Spatial as deep: Spatial cnn for traffic scene understanding[END_REF], [START_REF] He | Accurate and robust lane detection based on dual-view convolutional neutral network[END_REF]. However, they are computationally complex and are difficult to be applied to driverless systems with limited computing resources.

To deal with the above challenges, we propose a real-time, robust and accurate lane detection approach. We use CNN to accurately extract the features of lane boundaries and obtain the lane boundary position in different traffic scenes, instead of the traditional lane detection algorithm which uses image processing. In order to reduce the computational complexity, we combine the constraints of spatial and temporal information of lanes to reduce the scale of CNN network structure and to ensure that our algorithm can run accurately and robustly on a self-driving platform in real time.

Related Work

Lane detection algorithms can be divided into two types, one includes traditional methods which use hand-crafted features like the gradient information and the color space for lane detection. Another one contains CNN-based methods that are more frequently used by researchers because of the rise of deep learning.

Since lane boundaries have color and shape information that are different from that of the road, there are many gradient-based approaches developed to extract lane boundary features by calculating the gradient information in images and also to determine lane boundary locations. Aly et al. [START_REF] Aly | Real time detection of lane markers in urban streets[END_REF] used Hough transform to detect lane boundaries. In [START_REF] Kumar Satzoda | Hierarchical additive hough transform for lane detection[END_REF], hierarchical Hough transform was used by authors to detect lane boundaries with the corresponding intensity variations. Because of CNN's powerful feature extraction capabilities, many researchers use CNN for lane boundary detection. He et al. [START_REF] He | Accurate and robust lane detection based on dual-view convolutional neutral network[END_REF] used lane boundary detection based on the front view and the top view. It first extracted lane boundary candidates using weighted hat-like filtering in the top view and then classified lane boundary candidates according to Dual-View Convolutional Neutral Network (DVCNN). Lee et al. [START_REF] Lee | Vpgnet: Vanishing point guided network for lane and road marking detection and recognition[END_REF] proposed VPGNet based on multi-task learning to accomplish lane detection and ground traffic identification. In the process of lane detection, VPGNet first divided the image into 8 × 8 grids and determined if each grid belongs to a lane boundary class, and then they used sub-sampling and clustering algorithm to extract lane boundaries from the VPGNet's output. Pan et al. [START_REF] Pan | Spatial as deep: Spatial cnn for traffic scene understanding[END_REF] proposed a novel network structure called SCNN, which introduced spatial information into the network structure and accurately identified lane boundaries, even in the case of lane boundary obstructions. However, the CNN network structures in these three methods are complicated and require a large amount of computing resources, making it difficult to use them in actual autonomous vehicle applications. Chen et al. [START_REF] Chen | Brain inspired cognitive model with attention for self-driving cars[END_REF] also used regression to do lane boundary detection. This article proposed an end-to-end method which used front view images as input and got the lane boundaries position directly from the CNN. However, it was greatly influenced by the environment.

Spatial and Temporal Based Lane Boundary Detection

The workflow of our algorithm is shown in Fig. 1: first, we use inverse perspective transform to get top view images and use the temporal and spatial relevance of a lane to get its estimated location. Then we cut the image into 10-20 sub-images centering on the estimated position of the lane boundary, each sub-image con-tains local lane boundary information. We then use CNN to classify and regress sub-images to get the exact location and category of local lane boundaries. Finally, we combine the local lane boundary location with the linear information and through the spline fitting, we can get the lane boundary information for the entire image.

Spatial-temporal Based Lane Boundary Position Estimation in Top View

Inverse Perspective Transform: It is common sense that images taken by a camera have perspective effects which will destroy the original geometry of lane boundaries. Therefore, the first step in our algorithm is to convert the front view to the top view by inverse perspective mapping. Inverse perspective transform can eliminate the distortion caused by perspective effects in the picture. Through this transform step, we can get the top-down projection of the road. We select the top view as the basis for lane detection. Then we utilize the inverse perspective transform method proposed in [START_REF] Hoiem | Putting objects in perspective[END_REF]: Coordinate Estimate: Both Lee et al. [START_REF] Lee | Vpgnet: Vanishing point guided network for lane and road marking detection and recognition[END_REF] and Pan et al. [START_REF] Pan | Spatial as deep: Spatial cnn for traffic scene understanding[END_REF] used the single-frame lane detection methods. The benefit is that these methods don't require any prior information in the process of lane detection task, which can be completed only by inputting the image of a single frame. However, single-frame algorithms have high time and space complexity and redundancy. In application, lane boundaries have strong spatial and temporal constraints, especially in the top view. A lane boundary is a continuous ground road marking for dividing different lanes. In practice, the camera typically samples the pavement information at 30 frame per second (FPS). By combining these two points, it is not difficult to conclude that the lane boundary position between upper and lower frame images does not change much, which means that the lane boundary in the next frame image is near the lane boundary of the previous frame. Generally speaking, the lane width is the same for multi-lane roads, and the space constraints are obvious in the top view.

  u v 1   = 1 z   f 0 u c 0 f v c 0 0 1     1 0 0 0 0 cosθ x -sinθ x y c 0 sinθ x cosθ x 0       x y z 1     (1)
Based on the above considerations, we propose to use the temporal and spatial constraints of lane boundaries to reduce the range of searching area (shown in Fig. 2(c)). The main difference between our detection algorithm and the algorithm in [START_REF] Pan | Spatial as deep: Spatial cnn for traffic scene understanding[END_REF] is that ours is based on multi-frame image lane detection. We use the historical lane boundary information of the previous frame image as a priori knowledge to estimate the lane boundary position of the current frame. We then use CNN to accurately detect the position of the boundary based on the estimated lane boundary position with the time constraints of the lane boundaries discussed above. The mean value of the former 10 frames' lane boundary positions is utilized as the estimated position of the current frame.

Sub-image Extraction:

In our algorithm, we use control points to represent lane boundaries and a lane boundary is usually represented by five control points. Typically, we select 10-20 control points according to the boundary type classification result. These control points represent the lane information of the current lane and the lane beside the current lane. After we obtain the estimated position of the lane boundary through the above steps, both the lane boundary regression and the boundary type classification are performed in the neighborhood of the estimated position. So at this step we will divide the image into 10-20 sub-images centered on the estimated location and pass them to CNN (as shown in Fig. 2(d)). We use the estimated position of the lane as a priori information and cut the sub-image centered on the estimated position. After obtaining the exact location of the lane boundary through CNN, we use the exact position we get as the control point for the current frame and also as a reference for the next frame. The sub-image size needs to be determined by the camera's internal and external parameters. We consider the following rules when determining the sub-image size: first, each sub-image can only contain one lane which needs to be as large as possible. Second, sub-images should be long enough in the vertical direction to avoid lane misdetection of the white dotted boundary. However, an excess of vertical length will result in a decline in the accuracy of the curve section. In experiments, the size of the sub-image in our algorithm is 220 × 300.

CNN for boundary type classification and lane boundary regression

Network Structure: Inspired by the network structures in [START_REF] Bojarski | End to end learning for self-driving cars[END_REF] and [START_REF] Shaoqing Ren | Faster r-cnn: Towards real-time object detection with region proposal networks[END_REF], we proposed a novel convolutional neural network based on multitasking, which can accurately track the lane boundary position while classifying lane boundary types. For each task of classification and regression, the network consists of 9 layers, including 5 convolutional layers and 4 fully connected layers. Before entering CNN, we normalized the input image. The feature of our network structure is that the sub-tasks share the 5 convolutional layers and a fully connected layer, thus compressing the parameter space. Experiments show that our CNN based on multitasking can improve the performance of line boundary classification and lane regression tasks during training and eventually converges well. Our CNN is structured as follows: The first 3 convolutional layers we use the 5 × 5 kernel with 2 × 2 stride, and it followed with 2 convolutional layers which have 3 × 3 kernel and non-stride. The last convolutional layer is connected with a fully connected (FC) layer which has 4096 neurons. Then the network is divided into two parts (for the two sub-tasks of linetype classification and lane boundary regression), for each contains the same scale of 3 FC layers. The 3 FC layers have 1000, 500, 100 neurons respectively. We added a softmax layer for the sub networks which we used as the linetype classifier.

In our algorithm, we use a multitasking learning mechanism to learn and predict simultaneously the lane boundary type and the lane boundary position. Our loss function is defined as follows:

loss = { L cls , boundarytype = N oneline αL cls + (1 -α)L loc , boundarytype ̸ = N oneline (2)
where α is a factor to balance the classification loss and regression loss. L cls andL loc are respectively the classification loss function and regression loss function. In this paper, we use the cross-entropy loss as the classification loss function and the mean square error loss as the regression loss function.

Linetype Classification and Lane Boundary Regression:

In autonomous vehicles, linetype classification is necessary since self-driving cars need to change lanes and sometimes even at a high frequency. Traditional methods are not good at boundary classification while CNN can realize it easily with its powerful classification capabilities. Based on lane categories, we classify lane boundaries into double yellow lines, white solid lines, white dashed lines, single yellow lines, Noneline, etc. The Noneline category is designed to improve the robustness of the system, as lane boundary information for lane boundaries, especially for adjacent lanes, tends to be obstructed by moving vehicles in a realistic traffic environment. In this case, the control point of the corresponding position will be set to Noneline type. We will not refer to the CNN coordinate regression value of this control point. Generally, the lane boundaries are occluded 1-2 control points, so we can still get the lane information accurately by using lane fitting.

Different from the methods which transform the regression problems into classification task, we predict the exact value of the lane boundary location directly from the CNN. There are two main benefits of using lane boundary regression: a high regression accuracy which can reach pixel level and a direct regression which requires low space and time costs with a small amount of computation. The CNN output is shown in Fig. 3(a) Lane Fitting: The algorithm of the previous step gives information about the control points of the current lane and its adjacent lane, including the boundary type and the exact location of the lane boundary. In this step, we fit the lane boundary through control points with the Catmull-Rom (CR) spline [START_REF] Twigg | Catmull-rom splines[END_REF]. The CR spline is also known as Overhauster spline, which is one kind of cubic interpolation spline. It is calculated as follows:

Lane Fitting and Optimization

P (s) = [ 1 u u 2 u 3 ]     0 1 0 0 -τ 0 τ 0 2τ τ -3 3 -2τ -τ -τ 2 -τ τ -2 τ         P i-2 P i-1 P i P i+1     (3)
where P i , P i-1 , P i+1 are the current point, the previous point and next point on the spline respectively. τ is the tension parameter which affects the spline's sharpness. The value of u is set in the interval (0, 1), which represents the points between P i and P i-1 .

Optimization: Searching for lane boundary information throughout the whole image is unnecessary. In our algorithm, we will verify the linetype of the adjacent right and left lanes. If the left lane or the right lane cannot be changed, we will cancel the lane boundary position estimation of the adjacent lane, and in the next frame If the current lane can be changed lanes (e.g. the left lane is white dashed boundary), we will initialize the corresponding adjacent lanes based on the lane constraints. We noticed that multi-lane roads have the same lane width in general. Therefore, by calculating the deviation between the current lane's boundaries on each side, we can estimate the lane's position relative to the adjacent lane. Then in the next frame's prediction we can detect lane information for the corresponding adjacent lane. In this way the amount of computation will also be minimized.

Experiments

In this section, we have performed a series of contrast experiments for verifying the effectiveness of our proposed algorithm. These experiments are mainly executed on three different open-source datasets containing various traffic scenarios as well as different data:

1) Road Vehicle Dataset (RVD). This dataset contains images taken in different traffic scenarios and weather conditions, which were collected by our autonomous vehicle research group about one year ago.

2) Caltech lane dataset. This dataset contains 1224 tabbed images taken in four different places with the weather condition, however, limited to sunny days.

3) Tusimple benchmark dataset. The Tusimple benchmark dataset contains 6408 labeled images, but again lacks lane boundary data for different traffic or weather conditions.

In the three above-mentioned datasets, RVD contains mainly data of highway or urban roads in China, while most of the other two datasets' data was collected on Americans roads. Except for evaluations of algorithm quantitative tests on these open-source datasets, we have also compared our experiment results with that of other state-of-the-art methods. Besides, we also implemented our proposed algorithm on embedded platforms (i.e., Nvidia TX2) to verify the capability of its real-time performance in real autonomous vehicles, along with algorithm migration and application on a real self-driving car named Pioneer in our research group. Our experiment platform is shown in Fig. 3(b).

Model Training

The input of our proposed model contains images with a resolution of 1280 x 1024 pixels. In order to speed up the operation and reduce the computational complexity, a set of pre-processing and transforms of original image input are needed before sending these images to the core CNN network. The final resolution of CNN input image patches is 110 x 150 pixels. Due to the special structure of our CNN network, there is no pre-training model available, which means that we need to train CNN from scratch. Considering the lane boundary classification task and the lane boundary position regression task share the common convolution layers, and the CNN classification task is easy to converge, we first train the lane boundary classifier. Next, we change the CNN learning rate, batch size and dropout rate for training the lane boundary regression as well as classification tasks in the same time. After that use the Adaptive Moment Estimation optimization (Adam) and adjust the proportion of regression loss to classification loss to 4 : 1. We then use the Rectifier neural networks (ReLU) activation function and add the softmax layer to the output of the classification network. Using the GTX 1080ti as a training platform, we trained 10 epochs for each of the classification and regression tasks, which took about 12 hours.

Evaluation

To estimate our model's effectiveness, we evaluated our algorithm by using precision and recall. Since we did the boundary type classification and lane boundary regression simultaneously, only samples with the correct classification and regression results will be set as True Positive (TP). As for the measurement of the regression result, we set the threshold to a value which is supposed to be accurate when it is larger than the deviation of the predicted lane boundary coordinates and the ground-true. In the experiment we set the threshold to 5 pixels. Table II shows the accuracy of our algorithm in various traffic scenarios.

Road Vehicle Dataset (RVD): Since RVD contains images collected from different weather conditions and traffic scenarios, we can evaluate our methods' robustness on it. There are some visualized results which are shown in Fig. 4. And also we measured precision and recall which are shown in Table . 1. The performance on RVD shows that our method can accurately recognize the lane boundaries in several environments especially in bad weather conditions.

Caltech lane dataset: We compared our methods with Lee's [START_REF] Lee | Vpgnet: Vanishing point guided network for lane and road marking detection and recognition[END_REF] on Caltech lane dataset. Table . 2 shows that our methods outperformed Lee's method in such scenario. Since Caltech lane dataset only cantains 1224 labeled images, this result represent that our method can extract the lane boundaries' features from a relative small scale of dataset efficiently.

Tusimple benchmark dataset: Tusimple benchmark dataset evaluates the methods performance on highway. We also tested our methods on it. Our methods only achieved about 90 % accuracy. Although the images in Tusimple dataset are collected from undulated roads, which seriously impact the quality of the top view image, our methods show the robustness in this situation. Some of the visualized results on Caltech lane dataset and Tusimple benchmark dataset shown in Fig. 4.

Implementation on Embedded System

Since self-driving platforms have limited computing resources, it is meaningful to implement the lane boundary detection algorithm on the embedded platform. We choose the GPU-based embedded platform Nvidia TX2, which has 256 CUDA cores and low power consumption as the experiment's platform. After implemented and optimized proposed algorithm on Nvidia TX2, the frame rate has achieved 30 FPS on embedded system.

Conclusion and Future Work

In this paper, we proposed a lane detection algorithm based on deep learning and spatio-temporal information. We used the spatial and temporal constraints of a lane to estimate the position of the lane and also applied CNN to predict the lane boundary type and the exact location of a lane. Based on the recognization of boundary types, we propose a lane detection framework that corresponds to human driving habits, which reasonably reduces the computational complexity of the algorithm and improves the robustness of it. Experiments showed that our algorithm can accurately and quickly detect lane boundary information. In the future work, we hope to combine laser sensor information with our algorithm to further improve the accuracy and robustness of lane boundary inspection.

Fig. 1 .

 1 Fig. 1. Our method's workflow consists of three parts. First we perform the inverse perspective transform and estimate the rough position of lane boundaries. Next we use CNN to classify the boundary type and regress the lane boundaries position. Finally, we optimize the output of CNN and perform lane fitting.

  (x, y, z) is the world coordinates and (u, v) is the pixel coordinates of the image. The camera tilt is θ x and the focal length is f . (u c , v c ) is the camera optical center and y c is camera height. Fig. 2(a) is an image of the real physical scale obtained by the transformation. For convenience, we have resized the image to 1280 × 1024. The image we used is shown in Fig. 2(b).

Fig. 2 .

 2 Fig. 2. (a) is the top view image of real physical scale. (b) is the resized image. (c) shows the searching area, where the rectangles are the searching area and the estimated lane boundary positions are in the center of the rectangles. (d) illustrates some sub-images extracted from the whole image.

Fig. 3 .

 3 Fig. 3. (a) The regression and classification outputs from CNN. The red points are the regression output and the rectangles refer to the boundary type classification output, where blue, yellow and green ones represent the white solid line, yellow line and white dashed line respectively. (b) Experiment platform.

Fig. 4 .Fig. 5 .

 45 Fig. 4. Visualization result on RVD

Table 1 .

 1 Evaluation on RVD Dataset.

	Highway	5066	97.067	95.389 96.221
	Rainy and Snowy	3802	93.127	94.864 93.987
	Night	3566	92.376	91.272 91.821
	Urban	3467	98.261	96.246 97.192
	Heavy Urban	3257	93.491	92.341 92.912
	Nonline	1237	95.912	94.366 95.133
	Curve	1520	96.491	92.231 94.312
	Complex Illumination	1382	94.821	93.412 94.112

Scenario Total Image Precision[%] Recall[%] F 1 [%]

Table 2 .

 2 Comparison with Lee's[START_REF] Lee | Vpgnet: Vanishing point guided network for lane and road marking detection and recognition[END_REF] work on Caltech dataset.

	Traffic	Methods Precision[%] Recall[%] F 1 [%]
	Cordova1	Ours Lee's	97.326 -	96.651 96.987 -88.4
	Washington1	Ours Lee's	98.067 -	97.129 97.596 -86.9

Acknowledge This research was partially supported by the National Natural Science Foundation of China (No. 61773312, 61790563), the Programme of Introducing Talents of Discipline to University (No. B13043). We are grateful to the reviewers for taking the time to read this article.


Adrian Cătălin Florea
email: acflorea@unitbv.ro

Răzvan Andonie
email: andonie@cwu.edu

A Dynamic Early Stopping Criterion for Random Search in SVM Hyperparameter Optimization

We introduce a dynamic early stopping condition for Random Search optimization algorithms. We test our algorithm for SVM hyperparameter optimization for classification tasks, on six commonly used datasets. According to the experimental results, we reduce significantly the number of trials used. Since each trial requires a re-training of the SVM model, our method accelerates the RS optimization. The code runs on a multi-core system and we analyze the achieved scalability for an increasing number of cores.

1

Introduction

Most Machine Learning (ML) models are described by two sets of parameters. The first set consists in regular parameters that are learned through training. The other set, called hyperparameters or meta-parameters, consists of parameters which are set before the learning starts. It is essential to identify the combination of hyperparameter values which produce the best (or closed to the best) generalization performance. This is done by re-training multiple models with different combinations of hyperparameter values and evaluating their performance. We call this re-training + evaluation for one set of hyperparameter values a trial. Since training a model can be very resource intensive, it is important to reduce the number of trials.

In the specific case of SVM classifiers, the algorithm performance depends on several parameters and it is quite sensitive to changes in any of those parameters [START_REF] Cortes | Support-vector networks[END_REF]. The choice of the kernel, for example, can have a dramatic influence on the classification performance [START_REF] Chapelle | Choosing multiple parameters for support vector machines[END_REF]. The cost parameter , controlling the trade-off between margin maximization and error minimization is also highly important as, for the nonseparable case, the algorithm must allow training errors. For a polynomial kernel, a wrong choice of the degree can easily lead to over-fitting [START_REF] Bishop | Pattern Recognition and Machine Learning (Information Science and Statistics)[END_REF].

The most commonly used hyperparameter optimization strategy is a combination of Grid Search (GS) and manual tuning1 [START_REF] Lecun | Efficient Backprop, ser[END_REF][START_REF] Hinton | A practical guide to training Restricted Boltzmann Machines[END_REF][START_REF] Smusz | Exploiting uncertainty measures in compounds activity prediction using support vector machines[END_REF]. More elaborate techniques are: Nelder-Mead [START_REF] Nelder | A Simplex Method for Function Minimization[END_REF], simulated annealing [START_REF] Kirkpatrick | Optimization by simulated annealing:Quantitative studies[END_REF], evolutionary algorithms [START_REF] Hansen | Reducing the time complexity of the derandomized evolution strategy with covariance matrix adaptation (cma-es)[END_REF], and Bayesian methods [START_REF] Thornton | Auto-WEKA: Combined selection and hyperparameter optimization of classification algorithms[END_REF].

Random Search (RS) is another standard technique for hyperparameter optimization. A nice feature of RS is the possibility of adaptive early stopping. The key is to define a good stopping criterion, representing a trade-off between accuracy and computation time. The rise of the randomized methods begun with the work of Bergstra and Bengio [START_REF] Bergstra | Algorithms for hyper-parameter optimization[END_REF][START_REF] Bergstra | Random Search for Hyper-Parameter Optimization[END_REF]. Using the same number of trials, RS generally yields better results than GS or more complicated hyperparameter optimization methods. Especially in higher dimensional spaces, the computation resources required by RS methods are significantly lower than for GS [START_REF] Lemley | Big holes in big data: A monte carlo algorithm for detecting large hyper-rectangles in high dimensional data[END_REF]. Also, RS methods are relatively simple and easy to implement on parallel computer architectures.

Several software libraries dedicated to hyperparameter optimization exist, some of them being autonomous, while others being built on top of existing ML software. LIBSVM [START_REF] Chang | LIBSVM: A library for support vector machines[END_REF] and scikit-learn [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF] come with their own implementation of GS, with scikit-learn also offering support for RS. Spearmint [START_REF] Snoek | Practical bayesian optimization of machine learning algorithms[END_REF] and Bayesopt [START_REF] Martinez-Cantin | Bayesopt: A bayesian optimization library for nonlinear optimization, experimental design and bandits[END_REF] are software packages dedicated to Bayesian optimization. Auto-WEKA [START_REF] Kotthoff | Auto-WEKA 2.0: Automatic model selection and hyperparameter optimization in WEKA[END_REF] is also able to perform Bayesian optimization but, unlike the previous two which are standalone libraries, it is built on top of Weka [START_REF] Hall | The WEKA data mining software: An update[END_REF]. Hyperopt [START_REF] Bergstra | Hyperopt: a Python library for model selection and hyperparameter optimization[END_REF] and Optunity [START_REF] Claesen | Easy Hyperparameter Search Using Optunity[END_REF] are currently two of the most advanced libraries for hyperparameter optimization.

Our contribution is an improved RS optimization technique, which reduces the number of trials, without a significant impact on the prediction performance. The key is a new dynamically calculated early stopping condition for RS. The method is implemented in parallel and achieves a good scalability. Our experiments are on the SVM classification problem applied to six commonly used datasets and five hyperparameters. According to them, our method accelerates the RS optimization.

The paper proceeds as follows. Section 2 describes our algorithm and the dynamic stopping condition, with an emphasis on the algorithm's parallel nature. Section 3 presents the experimental results and the paper is concluded with Section 4.

Proposed Algorithm and Probabilistic Properties

A highly simplified version of a hyperparameter optimization algorithm is characterized by an objective fitness function f and a generator of samples g . The fitness function returns a classification accuracy measure of the target model, computed either through cross-validation or on a separate validation set. The generator g is in charge of providing the next set of values that will be used to compute the model's fitness. A hasNext method implemented by the generator offers the possibility to terminate the algorithm before the maximum number of N evaluations is reached, if some convergence criteria is satisfied.

In the particular case of RS, the generator g simply draws samples from the specific distribution of each of the hyperparameters to be optimized. Our goal is to reduce the computational complexity of the RS method in terms of number of trials. In other words, we aim to compute less than N trials, without a significant impact on the value of the fitness function.

For this, we introduce a dynamic stopping criterion, included in a randomized Dynamic Early Stopping Criterion for Random Search in SVM Hyperparameter Optimization 3 optimization algorithm (Algorithm 1). The algorithm is a two step optimizer. First, it iterates for a predefined number of steps n , N n << , and finds the optimal combina- tion of hyperparameter values, temp_opt. Then, it searches for the first result better than temp_opt. The optimal result, opt, is either the first result better than temp_opt or

temp_opt if N is reached.
= f(rndPoint) if (f_rnd > tmp_opt) { index = i; p = rndPoint; tmp_opt = f_rnd } } for i := n+1; i < N; i++ { rndPoint := g.Next(); f_rnd, _ := f(rndPoint) if (f_rnd > tmp_opt) { index = i; p = rndPoint; opt = f_rnd break } } return index, p, opt }
The following problems arise: i) Can we determine a value for n that maxim- izes the probability of obtaining the best results?; and ii) Can the algorithm be parallelized without impacting the probability of obtaining an optimal value?

Sequential execution

Algorithm 1 finds the optimum under the assumption that opt is in any position i , n i > , and no result better than temp_opt is in the

range 1] 1, [   i n .
We denote by i E1 the event that opt is reached on the i -th trial, and by i E 2 the event that no value better than temp_opt is obtained between the n -th and the i - th trial. The probability of i

E1 is N i E P 1/ = ) 1 ( (1)
The probability that all values in the range 1

] 1, [   i n
are worse than temp_opt is the same as the probability that the best result among the first 1

 i at- tempts lays in the range ] [1, n : 1) /( = ) 2 (  i n i E P (2)
Since the two events are independent, the probability that we hit opt after n i > attempts is:

1

)) ( /( = ) 2 ( ) 1 ( =   i N n i E P i E P i P (3)
The event E of finding opt after at most

m trials, N m n < < (where /2 = N m is a reasonable target), has probability i m n i N n i m n i N n i P m n i E P 1 1 = = 1) ( 1 1 = = 1 = = ) (        (4)
Since i 1/ is monotonically decreasing, the right term of eq. ( 4) has a lower bound:

i m n i N n dx m n x N n 1 1 = 1     (5)
We differentiate the left term of eq. ( 5):

  1), ln ln ( 1 = ln ln          n m N n m N n dn d ( 6
)
equate to zero and solve for n obtaining: e m n / = (7)

Fig. 1. Lower bound heatmap of the probability to obtain best result from a target space of maximum 300 attempts while terminating faster, depending on the values of m (x axis) and n (y axis). Darker shades correspond to greater probability.

Choosing for

n a value larger than the optimal one increases the probability of finding the combination of values that yields the optimal result but with an increased risk of a greater number of trials. The result from eq. ( 7) can be used to implement an improved version of the Algorithm 1 that can automatically set the value of n to e N/ . For example, in order to maximize the chances to obtain the best value, after a target maximum of 150 attempts, we must set n to e 150/ ( 55  ). For a target maximum of 100 attempts, n should be 37, and so on. Fig. 1 shows the lower bound heatmap of the probability to obtain the best results while stopping earlier with respect to the values of m and n .

Parallel execution

We generate a parallel implementation of our method as follows:

-Split the work between W workers (can be anything from lightweight threads of execution, OS threads, CPU cores or even different servers). We Dynamic Early Stopping Criterion for Random Search in SVM Hyperparameter Optimization 5 decided to use the GOLANG [START_REF] Google | The Go programming language[END_REF] support for goroutines 2 , which are basically lightweight threads managed by the GO run-time.

-Each worker w executes W N/ trials using the same early stopping criteri- on. In this case, W n w n / = , signifying that on average, with /2 = N m , W workers will terminate after ) /(2W N trials, with W N/ being the worst case.

-The manager gathers the results from all workers and selects the best candidate.

Algorithm 2 implements the above steps. The random values are generated and distributed by the manager, or each worker generates its own random sequence. Any of the following parallel pseudo-random number generation strategies can be selected [START_REF] Quinn | Parallel Programming in C with MPI and OpenMP[END_REF]: Manager-Worker (MW), Sequence Splitting (SS), Leapfrog (LF), and Parametrization (P).

The inverse problem

Given a restricted computational budget, expressed by a target number of trials m , we obtained the optimal value for n . We are now interested in solving the reverse problem: Given an acceptable probability P to achieve the best result among the N trials, which is the optimal value for n ? For the RS algorithm without the dynamic stopping criterion, if all trials are independent, the required number of trials needed to identify the optimum with a probability is given by The problem becomes interesting in the context of our stopping criterion when we are willing to compromise, by accepting a lower probability , for a further reduction of the number of trials. In case of Algorithm 2, according to eq. ( 5), probability P has a lower bound:

n em N n dx x m n N n P ln = 1 1            (8)
This, together with eq. ( 7) gives: 9) represents the probability to identify the optimum regardless of the activation of the stopping criterion -might also be among the first trials in which case the algorithm will test all the possible combinations. The probability to find the optimum after a number of trials strictly lower than N has a lower bound given by relation [START_REF] Hinton | A practical guide to training Restricted Boltzmann Machines[END_REF], which translates to:

(10) The value of n in eq. ( 8) can be adjusted in the interval ] , / [ m e m to increase the probability of identifying the optimal value, but at the same time increase the computational cost (the number of trials).

Experiments

We use our method to optimize the following five hyperparameters of a SVM [START_REF] Cortes | Support-vector networks[END_REF] classifier: kernel type (RBF, Polynomial or Linear chosen with equal probability), 

(drawn from an exponential distribution with 10 =  ); cost( C , drawn from an expo- nential distribution with 10 =  ); degree (chosen with equal probability from the set {2,3,4,5}) and 0 coef on [0,1] ). We run our experiments on six of the most popular datasets from UCI Machine Learning Repository 3 : Adult (a1a), Adult (a6a), Breast Cancer, Diabetes, Iris and Wine. Adult (a1a) and Adult (a6a) are variations of the same dataset but with different number of samples; the second one is around six times larger. Details of the datasets are presented in Table 1. We apply ten fold cross-validation to evaluate the classification accuracy [START_REF] Sokolova | A systematic analysis of performance measures for classification tasks[END_REF] and compare the obtained results, both in terms of classification performance and number of trials. We use the following optimizers (all implemented in the Optunity library): GS, RS, Particle Swarm, and Nelder-Mead. We also use the Weka SVM, with its implicit hyperparameters.

We run the Algorithm 2 with 8 = W and 250 = N

, which leads to 92 = n . We also run the four optimizers in Optunity, for a maximum number of 250 trials.

Accuracy estimation

Table 2 presents the results of applying Algorithm 2 for four parallelization strategies, compared with the results obtained with Optunity (RS, GS, Particle Swarm and Nelder-Mead) as well as with the results obtained using Weka with each of the three kernels (RBF, Polynomial and Linear) and the implicit values for the other parameters ( 1

.0 = C , features of number _ _ 1/ =  , 3 = degree , 0.0 = 0 coef
). The best results are marked in bold. Since we compare multiple classifiers on multiple datasets, we have to use additional statistical tests for further investigation, as suggested in [START_REF] Demsar | Statistical comparisons of classifiers over multiple data sets[END_REF].

We calculate the Friedman [START_REF] Friedman | A comparison of alternative tests of significance for the problem of m rankings[END_REF] and the Iman-Davenport [START_REF] Iman | Approximations of the critical region of the Friedman statistic[END_REF] statistics using eq. ( 11), respectively eq. ( 12), with N being the number of datasets, k the number of algorithms and j R the average rank of algorithm j from Table 3, and obtain 6.04 = 32.826, = 2

F F F  .             4 2 1) ( 2 1 = 1) ( 12 = 2 k k j R k j k k N F  (11) 2 1) ( 2 1) ( = F k N F N F F      (12)
With 11 algorithms and six data sets,

F F is distributed according to the F distribution with 10 = 1 11  and 50 = 1) (6 × 1) (11  
degrees of freedom. The critical value of (10,50) F for 0.05 =  is 2.03, so we reject the null-hypothesis, which means the algorithms are not equivalent in terms of prediction performance.

Adrian Cătălin Florea and Răzvan Andonie

The critical difference [START_REF] Nemenyi | Distribution-free Multiple Comparisons[END_REF][START_REF] Demsar | Statistical comparisons of classifiers over multiple data sets[END_REF] is given by:

N k k q CD 6 1) ( =    (13)
where critical values  q are based on the Studentized range statistic divided by 2 . At significance level of 0.05 =  , the critical difference is . ). At significance level 0.1 =  , the critical difference is and we observe that Optunity -NM is significantly worse than GO -LF and GO -P and also that GO -P is significantly better than Optunity -NM and Weka -RBF.

Efficiency of the stopping condition

Based on the above results, we exclude Nelder-Mead (due to its significantly worse classification performance) and the Weka SVM (since is it does not perform a real hyperparameter optimization) from the analysis. Table 4 depicts the rank across all datasets in terms of number of trials.

We perform another Friedman test and, using formulas [START_REF] Bergstra | Algorithms for hyper-parameter optimization[END_REF] and [START_REF] Bergstra | Random Search for Hyper-Parameter Optimization[END_REF], and obtain:

28.554 = ) is 2.420. This means that we can rule out the null- hypothesis and state that the algorithms are not equivalent with respect to the number of trials. We compute the critical difference according to formula (13) and obtain 3.678 = 0.05

CD

. Table 5 shows the difference in the average rank values for each pair of algorithms. The values greater than 0.05

CD

are marked in bold font. We can identify two groups of algorithms, the first group (GO -SS and GO -P) performs significantly better than the second group (Optunity -GS, Optunity -RS and Optunity -PS). It is not clear to which of the two groups GO -MW and GO -LF belong to. One possible explanation for the better results obtained by GO -SS and GO -P may be related to the superior parallel implementation of the random generators. However, since the number of random values generated in our tests is relatively small, this difference in Dynamic Early Stopping Criterion for Random Search in SVM Hyperparameter Optimization 9 performance is most probably coincidental. The Holm test rejects all four hypotheses, since the corresponding p values are smaller than the adjusted  's, leading to the conclusion that all four versions of our algorithm are significantly more efficient in terms of number of trials than the standard RS implementation.

Scalability

Besides the accuracy and the number of runs we also measure the algorithm's speedup (the ratio of the sequential execution time to the parallel execution time) as a Adrian Cătălin Florea and Răzvan Andonie measure of its scalability. The values are depicted in Table 7.

Conclusions

We introduced a new dynamic stopping condition for RS based hyperparameter optimization, together with its parallel implementation. In the context of SVM classification, on six of the most commonly used datasets, we obtained on par accuracy values with the existing mainstream hyperparameter optimization techniques. With all four of the parallel random generators used, the algorithm terminates after a significantly reduced number of trials compared to the standard implementation of RS, which leads to an important decrease in the computational budget required for the optimization.

The present work opens further research directions in terms of optimizing the hyperparameters for other ML algorithms where the search space has a larger number of dimensions and the required computational budget is currently a major issue. The algorithm implementation is flexible enough to allow a gradient-free optimization of any function.

Algorithm 1 .

 1 Parametric stop optimizer func Maximize(f, g, n, N){ index = -1; tmp_opt = -math.MaxFloat64; opt = tmp_opt for i := 0; i <= n; i++ { rndPoint := g.Next(); f_rnd, _ :

Algorithm 2

 2 Parallel stop optimizer func PMaximize(f, g, N, W){ // channels used by workers to communicate their results resultsChans := make(chan fn.Sample, W) for w := 0; w < W; w++ { go func(w int) { // index, point, function value, global optim, k i, p, v, gv, k := Maximize(f, g, N/W, w ) resultsChans <-fn.Sample{i, p, v, gv, k} }(w) } // Collect results results := make([]fn.Sample, W) for i := 0; i < W; i++ { results[i] = <-resultsChans } }

Table 1 .

 1 Details on used datasets

	Dataset	Instances	Features	Classes #
	Adult (a1a)	1,605	123	2
	Adult (a6a)	11,220	123	2
	Breast Cancer	683	10	2
	Diabetes	768	8	2
	Iris	150	4	3
	Wine	178	13	3

Table 2 .

 2 Accuracy and number of trials for Algorithm 2 using different parallelization strategies (MW, SS, LF, P), compared with Optunity (RS, GS, Particle Swarm and Nelder-Mead) and Weka's SVM.

		GO		GO		GO		GO		Optunity Optunity Optunity	Optunity	Weka	Weka	Weka
		MW	SS		LF		P		RS	GS	PS	NM	RBF	Poly	Linear
	Dataset	Acc	Runs	Acc	Runs	Acc	Runs	Acc	Runs	Acc	Acc	Acc	Acc	Runs	Acc	Acc	Acc
	Adult (a1a)	0.837	164	0.836 194 0.837 194 0.837	148	0.837	0.835	0.838 0.833 108	0.828 0.839 0.754
	Adult (a6a)	0.844	169	0.844 133 0.844 203 0.845	138	0.844	0.844	0.845 0.843 142	0.838 0.760 0.760
	Cancer 0.975	172	0.975 167 0.975 214 0.975	187	0.975	0.975	0.975 0.968 17	0.969 0.974 0.969
	Diabetes 0.776	223	0.776 203 0.780 220 0.779	130	0.777	0.777	0.776 0.651	6	0.775 0.686 0.777
	Iris	0.973	224	0.980 156 0.980 189 0.973	158	0.980	0.987	0.967 0.940	6	0.953 0.727 0.960
	Wine	0.989	194	0.989 215 0.989 162 0.989	177	0.989	0.984	0.989 0.956 62	0.983 0.404 0.972
	Average 0.899	191	0.900 178 0.901 197 0.900 156.334 0.900	0.900	0.898 0.865 56.834 0.891 0.732 0.865
	STDEV 0.091 27.188 0.092 31.241 0.091 20.746 0.090 22.223 0.092	0.092	0.090 0.120 57.711 0.088 0.190 0.112

Table 3 .

 3 Algorithms' accuracy ranking on the used datasets.

		GO	GO	GO	GO Optunity Optunity Optunity Optunity	Weka Weka	Weka
	Dataset	MW	SS	LF	P	RS	GS	PS	NM	RBF	Poly Linear
	Adult (a1a)	6	7	3.5	3.5	5	8	2	9	10	1		11
	Adult (a6a)	3	4	7	1	5	6	2	8	9 10.5 10.5
	Cancer	5.5	5.5	5.5	5.5	2.5	2.5	1	11	9.5	8		9.5
	Diabetes	8	7	1	2	5	4	6	11	9	10		3
	Iris	5.5	3	3	5.5	3	1	7	10	9	11		8
	Wine	2.5	2.5	2.5	2.5	5.5	7	5.5	10	8	11		9
	Average	5.083 4.833 3.750 3.333 4.333 4.750 3.917 9.833 9.083 8.583 8.500
	This clearly rules out the Nelder-Mead algorithm, which is significantly worse
	than the Parametrization based implementation of our RS (	9.833 	3.333	>	6.163

Table 4 .

 4 Algorithms ranking in terms of number of runs.

		GO	GO	GO	GO Optunity Optunity Optunity
	Dataset	MW	SS	LF	P	RS	GS	PS
	Adult (a1a)	2	3.5	3.5	1	6	6	6
	Adult (a6a)	3	1	4	2	6	6	6
	Cancer	2	1	4	3	6	6	6
	Diabetes	4	2	3	1	6	6	6
	Iris	4	1	3	2	6	6	6
	Wine	3	4	1	2	6	6	6
	Average	3.000	2.083	3.083	1.833	6.000	6.000	6.000

Table 5 .

 5 Algorithms ranking difference in terms of number of runs.

		GO	GO		GO			GO	Optunity	Optunity	Optunity
		MW	SS		LF					P			RS	GS	PS
	GO MW	-	-0.917	0.084		-1.167				3	3	3
	GO SS		-		1		-0.25	3.917	3.917	3.917
	GO LF				-		-1.25	2.917	2.917	2.917
	GO P									-	4.167	4.167	4.167
	Optunity RS													-	0	0
	Optunity GS													-	0
	Optunity PS													-
	Finally, since the main goal of our work is to obtain an improved version of
	RS, we compare our method directly with Optunity RS, using the Holm [29] test. The
	standard error for our experiment is	SE	=	k	(	k		1)/(6	N	)	=	1.247	. Table 6 shows the
	results of the Holm rejection test.

Table 6 .

 6 Performance in terms of number of trials required for GO -MW, GO -P, GO -SS and GO -LF against Optunity -RS in terms of the Holm test.

		i	z	=	(	R	0		SE )/ Ri	p		/	i
	GO P	1							3.341	0.00084	0.0084
	GO SS	2							3.140	0.00043	0.01
	GO MW	3							2.405	0.00808	0.0125
	GO LF	4							2.339 0.009668	0.0167

Table 7 .

 7 Algorithm speedup with increasing number of cores.

	Dataset/Cores	2	3	4	6	8
	Adult(a1a)	1.37	2.96	3.48	4.13	4.52
	Adult (a6a)	1.97	2.71	3.02	3.35	3.70
	Cancer	1.98	2.91	3.53	3.91	4.09
	Diabetes	1.86	2.72	3.34	3.70	3.88
	Iris	1.94	2.75	3.10	3.43	3.54
	Wine	1.99	2.81	3.28	3.96	4.11
	Average	1.85	2.81	3.29	3.75	3.97

https://github.com/jaak-s/nips2014-survey -82 out of 86 optimization related papers presented at the NIPS

conference used GS

https://tour.golang.org/concurrency/1


Daniel Stamate

Wajdi Alghamdi

Daniel Stahl

Alexander Zamyatin

Robin Murray

Marta Di Forti

Can Artificial Neural Networks Predict Psychiatric Conditions Associated with Cannabis Use?

Keywords: Machine Learning, Neural Networks, Prediction modelling, ROC optimisation, Monte Carlo, Computational Psychiatry, Cannabis, Psychosis

This data-driven computational psychiatry research proposes a novel machine learning approach to developing predictive models for the onset of firstepisode psychosis, based on artificial neural networks. The performance capabilities of the predictive models are enhanced and evaluated by a methodology consisting of novel model optimisation and testing, which integrates a phase of model tuning, a phase of model post-processing with ROC optimisation based on maximum accuracy, Youden and top-left methods, and a model evaluation with the k-fold cross-testing methodology. We further extended our framework by investigating the cannabis use attributes' predictive power, and demonstrating statistically that their presence in the dataset enhances the prediction performance of the neural network models. Finally, the model stability is explored via simulations with 1000 repetitions of the model building and evaluation experiments. The results show that our best Neural Network model's average accuracy of predicting first-episode psychosis, which is evaluated with Monte Carlo, is above 80%.

Introduction

Policy changes and the legalisation of cannabis across countries, especially the Netherlands, Uruguay and some states in the USA, indicate that cannabis is gaining greater global acceptance. It has been reported that cannabis has been the most popular illicit drug in the world in the last decade, with an estimated 183 million annual users [1]. However, studies show that countries with higher cannabis consumption also have higher proportions of people seeking treatment for psychotic disorders, including firstepisode psychosis, schizophrenia, anxiety and substance use disorders [START_REF] Radhakrishnan | Gone to Pot: A Review of the Association between Cannabis and Psychosis[END_REF]. This makes it imperative to understand the likely consequences of cannabis use, even if these consequences affect only a minority of users. Therefore, the link between cannabis consumption and the risk of experiencing psychotic disorders must be thoroughly scrutinised.

Recent researches have attempted to understand whether specific patterns of cannabis use, such as potency or age, are associated with a higher risk of developing psychotic disorders. One study concluded that nearly a quarter of all new psychosis patients in South London (UK) could be associated with the use of high-potency, skunk-like cannabis [START_REF] Di Forti | Proportion of Patients in South London with First-Episode Psychosis Attributable to Use of High Potency Cannabis: a Case-Control Study[END_REF]. Another study [START_REF] Alghamdi | A Prediction Modelling and Pattern Detection Approach for the First-Episode Psychosis Associated to Cannabis Use[END_REF] estimated that a person who uses cannabis daily for more than six months has a 70% likelihood of suffering from psychotic disorders.

Few studies have used risk prediction modelling or advanced machine learning algorithms to establish a link between cannabis use and first-episode psychosis. In fact, apart from our recent work [START_REF] Alghamdi | A Prediction Modelling and Pattern Detection Approach for the First-Episode Psychosis Associated to Cannabis Use[END_REF], we are not aware of the existence of any such studies. Most prior studies have relied solely on explanatory research strategies and been based on various conventional statistical techniques, such as hypothesis formulation and verification via statistical tests, logistic regression modelling, etc. These methods are wellrecognised and widely used in medical research, but they rarely match the high potential of machine learning methods. The domain of machine learning has continued to develop for many years, and advanced predictive techniques have been expanded and improved constantly. These advanced predictive techniques, in turn, are prompting extensive usage of computers and artificial intelligence techniques in many domains, such as medicine [START_REF] Zhou | Machine Learning for Medical Applications[END_REF] [START_REF] Iniesta | Machine learning, Statistical Learning and the Future of Biological Research in Psychiatry, Psychological Medicine[END_REF]. Artificial neural networks have proven their remarkable ability to detect predictive patterns in different types of datasets of various complexities, and often showed their superiority compared with other machine learning techniques.

Artificial neural networks are computational models that simulate the way biological neural networks process information in the human brain [START_REF] Kuhn | Applied Predictive Modelling[END_REF]. They are usually composed of several highly interconnected groups of artificial neurons that work together to solve specific problems. These computational models are typically used to model complex relationships between inputs and outputs, such as those in prediction models. Just like learning in biological neural networks, artificial neural networks process information by 'tuning' existing connections among neurons in order to process information.

While techniques involving support vector machines are still popular within the machine learning community [START_REF] Alghamdi | A Prediction Modelling and Pattern Detection Approach for the First-Episode Psychosis Associated to Cannabis Use[END_REF] [START_REF] Katrinecz | Predicting Psychosis Using the Experience Sampling Method with Mobile Apps[END_REF], artificial neural networks are gaining considerable attention again. Recently, artificial neural networks have been successfully used in understanding the heterogeneous manifestations of asthma [START_REF] Belgrave | Predictive Modelling Strategies to Understand Heterogeneous Manifestations of Asthma in Early Life[END_REF], diagnosing tuberculosis [START_REF] Elveren | Tuberculosis Disease Diagnosis Using Artificial Neural Network Trained with Genetic Algorithm[END_REF], classifying leukaemia [START_REF] Adjouadi | Classification of Leukaemia Blood Samples Using Neural Networks[END_REF], detecting heart conditions in electrocardiogram (ECG) data [START_REF] Yan | A Restricted Boltzmann Machine Based Two-Lead Electrocardiography Classification[END_REF], etc. These studies show that neural networks are capable of handling complex medical data, such as ambiguous ECG signal data, and achieving outstanding results not yet produced by other methods.

In this study, we offer a novel machine learning approach that uses neural networks to develop predictive models for the onset of first-episode psychosis. The dataset on which we based our study was collected by psychiatry practitioners and has been used in previously conducted studies, such as [3][4]. It comprises an extensive set of variables, including demographic, drug-related and other variables, with specific information on participants' histories of cannabis use, as seen in Table 1. Our paper proposes a novel data-driven computational psychiatry and machine learning approach to developing predictive models for the onset of first-episode psychosis, based on feed-forward artificial neural networks. The performance capabilities of the predictive models are enhanced and evaluated by a methodology consisting of model optimisation and testing, which integrates a phase of model tuning, a phase of model post-processing with ROC optimisation based on maximum accuracy, Youden and top-left methods, and a model evaluation with the k-fold cross-testing methodology. We further extend our framework by investigating the cannabis use attributes' predictive power, and demonstrating statistically that their presence in the dataset enhances the prediction performance of the neural network models. Finally, the model stability is explored via simulations with 1000 repetitions of the model building and evaluation experiments. The results show that our best models' accuracies in predicting first-episode psychosis in intensive Monte Carlo simulation fall between 75.03% and 85.13%, with an average of about 81%.

The rest of the paper is organised as follows. Section 2 presents our methodology for predicting the first-episode psychosis, based on experimenting with artificial neural networks, and our novel methodology for model optimisation and post-processing, and evaluation with optimized cut-off point selection on the ROC curve. The section also investigates the outcomes of the extensive Monte Carlo simulations in order to study the variation of the models' performance. In Section 3, we build optimised prediction models without the cannabis attributes to study if there is a statistically significant difference with respect to the performances of the models using the cannabis attributes. Finally, the conclusion and the directions for future work are presented in Section 4.

2

Building Prediction models

Data preparation

The data we used to build our predictive models were a part of a case-control study [START_REF] Di Forti | Proportion of Patients in South London with First-Episode Psychosis Attributable to Use of High Potency Cannabis: a Case-Control Study[END_REF]. The clinical data comprise 1106 records divided into 489 patients, 370 controls and 247 unlabelled records. The patients were individuals who presented with first-episode psychosis to the inpatient units of the South London & Maudsley Mental Health National Health Service (NHS) Foundation Trust. The controls were healthy people recruited from the same area served by the Trust. The control samples were similar to the patient samples in age, gender, ethnicity, educational qualifications and employment status.

Each record in the data refers to a participant in the study and has 255 possible attributes divided into four groups. The first group consists of demographic attributes, which represent general features like gender, race and level of education. The second group of drug-related attributes contains information on the use of non-cannabis drugs, such as tobacco, stimulants and alcohol. The third group contains genetic attributes. These were removed from the analysis for the purposes of this study. The final group contains cannabis-related attributes, such as the duration of use, initial date of use, frequency, cannabis type, etc. (see Table 1).

The goal of this stage is to perform a high-level simplification of the dataset and to prepare the dataset for use in our novel approach to predict first-episode psychosis. This stage involved several steps. First, records that were missing critical data were removed from the dataset. This included both records with missing labels and records with missing values on all cannabis-related variables. Second, certain variables were removed from the dataset. This primarily included variables that were deemed to be irrelevant to the study (e.g. those related to the individual IDs of the study participants) and variables that fell outside the scope of the study (e.g. certain gene-related variables). In addition, any numeric predictors with zero or near-zero variance were dropped. Third, we sought to standardise the encoding of missing values across the dataset. Prior to this step, values including 66, 99 and -99 all represented cases with missing values; thus, we replaced all such indicators with a consistent missing value indicator: NA. Fourth, some variables were re-labelled to provide more intuitive descriptions of the data they contained. In multiple situations, some variables had similar meanings but also records with missing values. Therefore, we conducted an imputation process to effectively combine the information from all the related variables into one. For example, two variables described alcohol use, but were inconsistently present and presented missing values across the records. These were combined in a way that created one single variable with consistent values that were as complete as possible. This process was used to generate value-reacher and value-consistent variables related to alcohol use, tobacco use, employment history and age.

Finally, any attribute that had more than 50% missing values was removed from the study. We then dropped any record for which more than 70% of the remaining attributes contained missing values. The resulting dataset, after the transformations above, contained 783 records and 78 attributes. The records are divided into 451 patients and 332 controls. A summary of some of these fields-specifically, those that relate to cannabis use, such as type, the age of first use and duration-can be seen in Table 1.

Missing values treatment

Although the data set was pre-processed and attributes with more than 50% missing values were removed, the final dataset still contained several missing values. Of the 783 records, only 22.8% were complete cases. This volume of missing information makes modelling more challenging but is often the reality of medical and social research.

The predictive power of the data may depend significantly on the way that missing values are treated. Some machine learning algorithms, such as decision trees [START_REF] Kuhn | Applied Predictive Modelling[END_REF], have the capability to handle missing data outright. However, most of machine learning algorithms do not have the capability to handle missing data. In many situations, missing values are imputed using a supervised learning technique, such as k-Nearest Neighbour (KNN. These imputation techniques do not have theoretical formulations but are often applied in practice [START_REF] Alghamdi | A Prediction Modelling and Pattern Detection Approach for the First-Episode Psychosis Associated to Cannabis Use[END_REF] [START_REF] Iniesta | Machine learning, Statistical Learning and the Future of Biological Research in Psychiatry, Psychological Medicine[END_REF]. Several imputations techniques, such as the KNN imputation, the tree bagging imputation from the caret package [START_REF] Kuhn | Applied Predictive Modelling[END_REF] and the random forest imputation from the randomForest package [START_REF] Liaw | Classification and Regression by randomForest[END_REF] were considered in this work. The last method, although it was the most computationally expensive, produced the best results regarding the performance of the final predictive models.

Training and tuning feed-forward artificial neural networks

To develop optimised predictive models for first-episode psychosis, we controlled the values of the parameters for each of the considered algorithms using chosen grids. Predictive models have been fitted in a five-fold cross-validation procedure, on each training set after pre-processing techniques were applied on the same training set, and have been tested on each test set. Models based on single-layer and multi-hidden-layer neural networks were optimised (tuned) to maximise AUC, the area under the ROC curve.

To avoid overfitting, the single-layer neural networks were tuned over 10 values for the size (i.e. the number of hidden units) and 10 values for the decay (i.e. the weight decay), which is the parameter in the penalisation method for model regularisation. This approach is like the penalisation method in ridge regression and is based on the L2 norm [START_REF] Kuhn | Applied Predictive Modelling[END_REF]. The optimal values were 17 and 0.01, respectively. Multi-layer neural networks were tuned over 10 values for each of the three hidden layers (i.e. 10 values for the number of hidden units in each layer) and over 10 values for the decay. The optimal values were 10, 10 and 10 for the three layers and 0.001 for the decay.

Treating unbalanced classes

When there is a priori knowledge of a class imbalance, one direct method to reduce the imbalance's influence on model training is to select training set samples with roughly equal event rates [START_REF] Kuhn | Applied Predictive Modelling[END_REF]. Treating data imbalances usually leads to better prediction models and a better trade-off between sensitivity and specificity.

In this study, we considered three sampling approaches to sub-sample the training data in a manner that mitigated the imbalance problem. The first approach was downsampling, in which we sampled (without replacement) the majority class to be the same size as the minority class. The second method was up-sampling, in which we sampled (with replacement) the minority class to be the same size as the majority class. The last approach was the synthetic minority over-sampling technique (SMOTE) [START_REF] Qazi | Effect of Feature Selection, SMOTE and Under Sampling on Class Imbalance Classification[END_REF]. SMOTE selects a data point randomly from the minority class, determines the K nearest neighbours to that point and then uses these neighbours to generate new synthetic data points using slight alterations. Our analysis used five neighbours. The results show that the up-sampling procedure yielded no real improvement in the AUC or the accuracy performances. A simple down-sampling of the data also had no positive effect on the model performances. However, SMOTE with neural networks models led to an increase in both the AUC and the accuracy.

As mentioned, data balancing supports a good trade-off between sensitivity and specificity. Another method that helps to balance sensitivity and specificity, or a good trade-off between the two performances, is model post-processing through the determination of new cut-off points on the ROC curves [START_REF] Kuhn | Applied Predictive Modelling[END_REF]. Our framework used three such methods, which can be seen as post-processing optimisations of the models. The first method found the point on the ROC curve closest to the top-left corner of the ROC plot, which represents the perfect model (100% sensitivity and 100% specificity). The second method is Youden's J index [START_REF] Bohning | Revisiting Youden's index as a useful measure of the misclassification error in a meta-analysis of diagnostic studies[END_REF], which corresponds to the point on the ROC curve farthest from the main diagonal of the ROC plot. The third method, which is "maximum accuracy" found the cut-off, which is the point with the highest accuracy.

In order to further improve the model performance, a specially designed post-processing procedure and model evaluation were adapted in our modelling procedure. First, the dataset is stratified splatted randomly into in 60% training data and 40% evaluation data. Then, the training data is used for training and for optimising the model, as explained in Subsection 2.3, in a cross-validation fashion, with AUC as optimisation criterion, with and without class balancing. Different pre-processing methods such as missing values imputation and sampling methods that we have explained above were appropriately integrated into the cross-validation. The optimal model obtained on the training data was then applied to the evaluation dataset in a specially designed postprocessing procedure, which we call the k-fold cross-testing method.

In the k-fold cross-testing method, we produce k post-processed model variants of the original optimised model. First, we create k stratified folds of the evaluation dataset. Then, k-1 folds are used to find an alternative probability cut-off on the ROC curve with one of the three specific methods presented above (top-left, Youden, and largest accuracy), obtaining a post-processed model variant. The remaining one-fold is scored with the post-processed model variant based on the newly found cut-off point. Finally, the whole procedure is repeated until all folds are used for scoring at their turn, then the predictions are integrated, and the model performance is measured on the whole evaluation dataset. We note here as an essential remark that in each such iteration of the procedure, the ROC optimisation data (the k-1 folds) and the scored data (the remaining fold) are always distinct, so the data for model post-processing and the data for scoring are always distinct.

Increasing model performance via optimized cut-off point selection on the ROC curve

The ROC curve is a graphical technique for evaluating the ability of a prediction model to discriminate between two classes, such as patients and controls. ROC curves allow visual analyses of the trade-offs between a predictive model's sensitivity and specificity regarding various probability cut-offs. The curve is obtained by measuring the sensitivity and specificity of the predictive model at every cutting point and plotting the sensitivity against 1-specificity. The left image in Fig. 1 shows the ROC curves obtained for both the single-layer neural networks and the multi-layer neural networks.

The curve shows that multi-layer neural network performs better regarding the evaluation dataset. to the ROC curve as shown in the right image in Fig. 1. This is the point with the largest value for the Youden index which is defined as (sensitivity + specificity -1) [START_REF] Bohning | Revisiting Youden's index as a useful measure of the misclassification error in a meta-analysis of diagnostic studies[END_REF]. These are the two most popular methods for establishing the optimal cut-off [START_REF] Perkins | The inconsistency of "optimal" cutpoints obtained using two criteria based on the receiver operating characteristic curve[END_REF] [START_REF] Kuhn | Applied Predictive Modelling[END_REF]. We used both of these methods, as well as the maximum accuracy approach, which determines the point on the ROC curve corresponding to the greatest accuracy (the red point in Fig. 1, right). In our analysis, the optimal cutting point was derived from independent sets, rather than from the training set or the evaluation sets, as shown previously. This is particularly important, especially, for small datasets.

Monte Carlo and models' stability

Due to the uncertainties introduced by the missing values in the data and due to expected variations of the predictive models' performance, depending on the datasets that were chosen for training and testing, we perform extensive Monte Carlo simulations to study the performances' variations and the models' stability. The simulations for each NN consisted of 1,000 iterations of the proposed procedure. The models' performances concerning accuracy, sensitivity, specificity and kappa were evaluated for each iteration on separate a testing dataset. The aggregation of all iterations yielded various distributions of the concerned performance measures. These distributions were then visualised using box plots in Fig. 2 to capture the models' performance capability and stability. The subfigures in Fig. 2 were grouped by their performances' measures into four subfigures. Each subfigure contained sex box plots for single-layer and multilayer neural networks with several optimized cut-off points on the ROC curve such as top-left, Youden index and the maximum accuracy. Also, estimations of the predictive neural networks' performances regarding means and standard deviation (SD) are shown in Table 2. The results as shown in Table 2 are regarding the models' performances when applied with the ROC optimisation techniques. On the one hand, single-layer neural networks with ROC optimisations based on Youden or top-left scored a mean accuracy of 0.76. The performance of the single-layer neural networks slightly improved when a ROC optimisation based on maximum accuracy was applied, resulting in a mean accuracy of 0.79 (95% CI [0.75, 0.83]) and a mean sensitivity of 0.83 (95% CI [0.76, 0.9]). On the other hand, the multi-layer neural networks with an ROC optimisation based on top-left cutting point achieved the best results with a mean accuracy of 0.81 (95% CI [0.77, 0.86]) and a mean sensitivity of 0.84 (95% CI [0.77, 0.88]), comparable to the results achieved by multi-layer neural networks with Youden and maximum accuracy.

This procedure is very computationally costly; therefore, a robust framework was essential. Parallel processing was performed on a data analytics cluster of 11 servers with Xeon processors and 832GB fast RAM. We used the R software with some packages, including caret, pROC, e1071, randomForest, ggplot2, plyr, DMwR, AppliedPre-dictiveModeling and doParallel [START_REF] Cran | The Comprehensive R Archive Network[END_REF]. In general, we detect that the proposed models have good predictive power and stability, based on an acceptable level of variation in their performance measures evaluated across extensive Monte Carlo experiments. However, the results indicate that the performance differences between the two types of neural networks with different methods for selecting the ROC cutting points are not significant regarding the 4 performances.

3

Cannabis use attributes' predictive power

After performing the Monte Carlo simulations, we further investigated the predictive models in order to better comprehend the predictive power of the cannabis-related attributes over first-episode psychosis. Moreover, we investigated the association between cannabis-related attributes and first-episode psychosis via statistical tests and attribute-ranking techniques.

Student's t-test

In this subsection, we comprehend the predictive power of the cannabis-related attributes over first-episode psychosis via statistical tests by re-fitting our performing models but with the cannabis-related attributes, represented in Table 1, removed from the dataset. Then, we compared the performances with and without the cannabis-related attributes using Student's t-test. We thereby demonstrated the predictive value of cannabis-related attributes with respect to first-episode psychosis by showing that there is a statistically significant difference between the performances of the predictive models built with and without the cannabis variables.

Our analysis also showed that the accuracy decreased by 5% for single-layer neural networks and by 6% for the multi-neural networks, if the cannabis-related attributes were removed from the process of building the predictive models. Then, we compared the accuracies of the single-layer neural networks models built on the data sets with and without the cannabis-related attributes using one-tailed t-test. The p-value obtained for the t-test was 5.51 × 10 -195 . As for the multi-layer neural networks models built on the data sets with and without the cannabis use attributes, the p-value obtained for the one-tailed t-test was and 2 × 10 -16 . This means that the predictive models with cannabis attributes have higher predictive accuracy than the models that were built without the cannabis attributes. In other words, the additional cannabis variables jointly account for predictive information over first-episode psychosis. These results are consistent with findings from [START_REF] Alghamdi | A Prediction Modelling and Pattern Detection Approach for the First-Episode Psychosis Associated to Cannabis Use[END_REF].

Ranking attributes' importance with the ROC curve approach

This subsection proposes the use of the ROC curves to determine the relevant variables affecting first-episode psychosis as introduced in [START_REF] Kuhn | Applied Predictive Modelling[END_REF]. We measure the individual importance of every attributes in the dataset to discover the attributes that yield significant improvements in the model predictivity power. To do so, the ROC curve is considered on each attribute. Then, a series of cutoffs is applied to the data to predict the class. The sensitivity and specificity are calculated for each cutoff, and the ROC curve is computed. Finally, the area under the curve is used as a measure of variable importance. Table 3 shows the top 10 attributes ranked by the ROC curve approach. The results in Table 3 support prior evidence that cannabis attributes, such as the type of the cannabis used and the frequency of usage, have significant power in predicting first-episode psychosis. For example, the results in Table 3 support findings from [START_REF] Di Forti | Proportion of Patients in South London with First-Episode Psychosis Attributable to Use of High Potency Cannabis: a Case-Control Study[END_REF] by associating the type of cannabis, especially high-potency cannabis, with the onset of psychosis. In addition, bindur.3 in Table 3, which represents the duration of cannabis use, is consistent with findings from [START_REF] Alghamdi | A Prediction Modelling and Pattern Detection Approach for the First-Episode Psychosis Associated to Cannabis Use[END_REF].

Conclusion and future work

This paper proposes a novel machine learning approach to developing predictive models for the onset of first-episode psychosis using artificial neural networks. We explored two types of artificial neural networks, each of which was able to recognise patterns differentiating patients from controls at an acceptable level. We based our approach on a novel methodology for optimising and post-processing predictive models. We also proposed several sampling methods and several methods for choosing the optimal cutting point on the ROC curve to improve the prediction models' performances. The models were then further tested using Monte Carlo experiments, and they consistently yielded adequate predictive power and stability. The best-performing models were multi-layer neural networks, which achieved accuracies as high as 88% in some cases and an average accuracy of 81% in Monte Carlo simulations with 1000 repetitions. The scored performances were above all performances achieved in previous studies such as [START_REF] Di Forti | Proportion of Patients in South London with First-Episode Psychosis Attributable to Use of High Potency Cannabis: a Case-Control Study[END_REF]. This paper extends on previous work as [START_REF] Di Forti | Proportion of Patients in South London with First-Episode Psychosis Attributable to Use of High Potency Cannabis: a Case-Control Study[END_REF] by proposing a new machine learning framework based on a novel methodology in which models are post-processed based on optimized cut-off point selection on the ROC curve and evaluated with the recent method of k-fold cross testing which we adapt after [START_REF] Katrinecz | Predicting Psychosis Using the Experience Sampling Method with Mobile Apps[END_REF]. Moreover, in this new framework, we developed optimized models with other powerful techniques such as artificial neural networks not addressed in [START_REF] Di Forti | Proportion of Patients in South London with First-Episode Psychosis Attributable to Use of High Potency Cannabis: a Case-Control Study[END_REF]. Also, the predictive power of cannabis-use attributes was tested via statistical tests and ranking methods to demonstrate statistically that their presence in the dataset enhances the prediction performance of the neural network models. The proposed approach proves the high potential applicability of machine learning and, particularly, artificial neural networks in psychiatry and enables researchers and doctors to predict and evaluate risks for first-episode psychosis.

One possible direction for future work is to further investigate how this prediction performance variation evolves by limiting the uncertainty in the data, represented by the high proportion of missing values. The second possible work direction is redefining the predictive modelling approach by considering more high-dimensionality data, such as genotype data. A third future work direction, which we are currently investigating, involves enhancing the power to predict first-episode psychosis using deep learning approaches.

Fig. 1 .

 1 Fig. 1. Left: The ROC curves for 2 of our optimised neural networks (NN) models: single-layer NN and multi-layer NN. Right: ROC optimisation post-processing of the multi-layer NN model, with 3 optimal cutting points: maximum accuracy, Youden and top-left methods.Numerous methods exist for finding a new cut-off. First, one can find the point on the ROC curve that is closest to the perfect model (100% sensitivity and 100% specificity), which is the point with shortest distance value from the point (0,1) as shown in the left image in Fig.1. To find the shortest distance, [(1 -sensitivity) 2 + (1specificity) 2 ] was calculated and minimised[START_REF] Pepe | The Statistical Evaluation of Medical Tests for Classification and Prediction[END_REF]. Another approach for finding an optimal cut-off point on the ROC curve is finding the largest distance from the diagonal

Fig. 2 .

 2 Fig. 2. 1000 Monte Carlo simulation for artificial neural networks. Where "SY" is single-layer NN with Youden cutting point, "ST" is single-layer NN with top-left cutting point, SM" is single-layer NN with maximum accuracy cutting point, "MY" is multi-layer NN with Youden cutting point, "MT" is multi-layer NN with top-left cutting point, and "MM" is multi-layer NN with maximum accuracy cutting point.

Table 1 .

 1 Cannabis use attributes in the analysed dataset[START_REF] Alghamdi | A Prediction Modelling and Pattern Detection Approach for the First-Episode Psychosis Associated to Cannabis Use[END_REF].

	Attribute	Description
	lifetime_cannabis_user	Ever used cannabis: yes or no
	age_first_cannabis	Age upon the first use of cannabis: 7 to 50
	age_first_cannabis_under15	Age less than 15 when first used cannabis: yes, no or
		never used
	age_first_cannabis_under14	Age less than 14 when first used cannabis: yes, no or
		never used
	current_cannabis_user	Current cannabis user: yes or no
	cannabis_fqcy	Pattern of cannabis use: never used, only on weekends
		or daily
	cannabis_measure	Cannabis usage measure: none, hash less than once per
		week, hash on weekends, hash daily, skunk less than
		once per week, skunk on weekends, skunk daily
	cannabis_type	Cannabis type: never used, hash or skunk
	duration	Cannabis use duration: 0 to 41 (months)

Table 2 .

 2 Estimations of the predictive neural networks' performances.

	Model	Accuracy mean SD mean SD mean SD mean SD Kappa Sensitivity Specificity
	Single-layers NN (Youden)	0.76 0.03 0.52	0.05 0.76 0.05 0.75 0.05
	Single-layers NN (top-left)	0.76 0.02 0.52 0.05 0.77 0.04 0.76 0.04
	Single-layers NN (max accuracy)	0.79 0.02 0.56 0.05 0.83 0.04 0.73 0.05
	Multi-layers NN (Youden)	0.8 0.02 0.59 0.05 0.84 0.02 0.75 0.04
	Multi-layers NN (top-left)	0.81 0.02 0.61 0.04 0.84 0.02 0.76 0.03
	Multi-layers NN (max accuracy)	0.80 0.08 0.59 0.13 0.87 0.04 0.74 0.08

Table 3 .

 3 ROC curve attribute importance

	Attribute	Importance
	typefreq.1	100.00
	cantype2.1	90.38
	type_use.hash	88.53
	totfreq.1	87.74
	bindur.3	85.84
	education.university.professiolqualifications	85.55
	frequenc.2	84.24
	bullying.no	78.55
	white	72.89
	homeless.1	72.86