Research Project

Posts about my research.

,

Distributed Machine Learning Using PySpark

Introduction

Big Data is not a new thing nowadays, and numerous platforms and applications have been created to take advantage of Big Data. Apache Spark is one of the most popular platforms for handling Machine Learning with Big Data.  By implicitly distribute data into clusters, Spark enables developers to focus more on their analytical tasks without worrying about data parallelism under the hood. To learn how Apache Spark works and how to build a distributed ML pipeline, a local cluster on my desktop is implemented with PySpark. In this post, I will talk about how to make ML pipelines with PySpark and its automated parallel-computing functionality. 

Problem Definition & Data Ingestion

One of the most exciting ML applications would be predictive modeling, which utilizes the data collected in the past to train a model that is capable of predicting the future. A predictive model can be a regression model or a classification model. For the demonstration, I will use Apache Spark to build a predictive classification model. The dataset I use is collected from The U.S. Department of Transportation’s (DOT) Bureau of Transportation Statistics, which contains more than 5 million records of the flights’ information during 2015, including the carrier of the airline, the original airport, destination airport, the distance of the flight, the time spent, the departure delay time, and arrival delay time, etc. I will use this dataset to train a model to predict if a flight will be delayed or not. 

Data Preparation

In order to do some basic EDA, the dataset needs to be preprocessed. The following code sets up a Spark working environment (Spark Session) and loads the data into the Spark Session as a Spark DataFrame:

# instantiate a spark session
from pyspark.sql import SparkSession
spark = SparkSession.builder.master('local[*]')
                            .appName('flights_delay')
                            .config('spark.driver.memory', '8g')
                            .getOrCreate()
# load raw data
df = spark.read.csv('/home/shihao/Desktop/flights/flights.csv', header=True,
                    inferSchema=True, nullValue=' ')
# clean data set
to_keep = ['MONTH', 'DAY', 'DAY_OF_WEEK', 'AIRLINE', 'FLIGHT_NUMBER',
           'ORIGIN_AIRPORT', 'DESTINATION_AIRPORT', 'TAXI_OUT',
           'SCHEDULED_DEPARTURE', 'DEPARTURE_DELAY', 'DISTANCE',
           'SCHEDULED_ARRIVAL', 'ARRIVAL_DELAY']
df = df.select(to_keep)
df = df.dropna()

 

When instantiate the Spark session in PySpark, passing 'local[*]' to .master() sets Spark to use all the available devices as executor (8-core CPU hence 8 workers). The original dataset has 31 columns, here I only keep 13 of them, since some columns cannot be acquired beforehand for the prediction, such as the wheels-off time and tail number.

After selecting all the useful columns, drop all the null values in the Spark DataFrame. Since only roughly 20% of the total flights are delayed (flights are delayed if the actual arrival time is more than 15 minutes later than the scheduled arrival time), use the following code to create the label column. The following code also balances the dataset, specifically, makes the number of delayed samples equal to the number of on-time samples.

# create label
df = df.withColumn('label', (df.ARRIVAL_DELAY > 15).cast('integer'))
# balance data set
# get number of delayed flights
num_0 = df.filter(df.label==0).count()
# calculate the sub-sampling ratio for the on-time flights
ratio = (df.count() - num_0)/num_0
# under sample the redundant class
df = df.sampleBy('label', {0: ratio, 1:1})

Explanatory Data Analysis

Now the dataset is clean and balanced, we can do some EDA to examine the dataset. First, take a look at how each airline is performing in terms of delivering flights on time. As shown in Figure 1 (left), Frontier Airlines (F9) has the highest average delay time (68.58 minutes) among 14 airline companies and Hawaiian Airlines (HA) has the lowest average delay time (40.98 minutes). 

Figure 1. The sorted average arrival delay and maximal delay of each airline.
Figure 1. The sorted average arrival delay and maximal delay of each airline.

On the other hand, American Airlines (AA) has the most extended delay with 1546 mins \approx 25.77 hrs, and Virgin America (VX) has the shorted maximal delay with 440 mins \approx 7.33 hrs as shown in Figure 1 (right).

Another interesting fact is illustrated in Figure 2, seems like all the airlines are doing pretty poorly around February and June, and doing better in April, September, and October. Therefore, Month should be a useful indicator for our predictive model. 

Figure 2. Monthly average arrival difference of each airline.
Figure 2. Monthly average arrival difference of each airline.

How about other features? Such as the distance of the flight. Does it have some relationship with the arrival timing? Figure 3 shows the average arrival difference for each airline, along with different flight distances. Arrival difference is defined as the difference between the scheduled arrival time and the actual arrival time. 

Figure 3. Average arrival difference plotted along distance for each airline.
Figure 3. Average arrival difference plotted along distance for each airline.

While there are a lot of outliers in Figure 3, we can still see a small trend that longer distance tends to have longer arrival difference. Figure 4. shows the relationship between arrival difference and taxi-out time, where we can see a more apparent linear relationship. 

Figure 4. Average arrival difference plotted along taxi-out time.
Figure 4. Average arrival difference plotted along taxi-out time.

Finally, as shown in Figure 5, departure-delay is the best predictor for our model.

Figure 5. Average arrival difference plotted along departure-delay.
Figure 5. Average arrival difference plotted along departure-delay.

Data Segregation

Now that the data-preprocessing and EDA are done, we can go ahead and build the model. I created three simple ML models with the dataset and compare their performances. Since the problem is to predict whether a flight will be delayed or not, which is a binary classification problem, I choose the Random Forest Classification model, Logistic Regression model, and Neural Network model, to implement the predictive model. 

First of all, split the dataset into training and test set before building ML pipelines to avoid leakage. Leakage happens when any of the .fit() methods are applied to the test set, which may induce unneglectable bias to our model.  Therefore, extra attention should be paid to only apply .transform() method to the test set. The dataset can be split in PySpark using DataFrame.randomSplit([training_ratio, test_ratio], random_state) method as shown here:

# create training and test set
flights_train, flights_test = df.randomSplit([.8, .2], 5)

Model training

When the data is ready, we can begin to build our machine learning pipeline and train the model on the training set. PySpark provides us powerful sub-modules to create fully functional ML pipeline object with the minimal code. Before putting up a complete pipeline, we need to build each individual part in the pipeline.

First, convert all the categorical data in the dataset into numerical data since the underlying machine is more familiar with numbers than texts. Similar to pandas.get_dummies(), use StringIndexer() combined with OneHotEncoderEstimator() sub-module in PySpark to convert all the categorical data into numerical data. For the flights dataset, there are three categorical columns, including AIRLINE, ORIGINAL_AIRPORT, and DESTINATION_AIRPORT. Therefore, construct three StringIndexer() objects, like so: 

# create parts
from pyspark.ml.feature import StringIndexer, OneHotEncoderEstimator
arr_indexer = StringIndexer(inputCol='AIRLINE', outputCol='arr_idx').setHandleInvalid("keep")
org_indexer = StringIndexer(inputCol='ORIGIN_AIRPORT', outputCol='org_idx').setHandleInvalid("keep")
des_indexer = StringIndexer(inputCol='DESTINATION_AIRPORT', outputCol='des_idx').setHandleInvalid("keep")
onehot = OneHotEncoderEstimator(inputCols=[arr_indexer.getOutputCol(), org_indexer.getOutputCol(),
                                           des_indexer.getOutputCol()],
                                outputCols=['arr_dummy', 'org_dummy', 'des_dummy'])

 

The OneHotEncoderEstimator() object will convert all the categorical columns into their corresponding numerical columns with sparse vector representations. Basically, a sparse vector representation is a more efficient way to represent a vector if the majority of the vector is zero values. 

Next, merge all the sparse vectors and other numerical values in the dataset into a single sparse vector by concatenating them one by one and call this new vector “features”, which can be done using the PySpark VectorAssembler() sub-module:

from pyspark.ml.feature import VectorAssembler
assembler = VectorAssembler(inputCols=['MONTH', 'DAY', 'DAY_OF_WEEK',
                                       'arr_dummy', 'org_dummy', 'des_dummy',
                                       'SCHEDULED_DEPARTURE','TAXI_OUT',
                                       'DEPARTURE_DELAY', 'DISTANCE',
                                       'SCHEDULED_ARRIVAL'],
                            outputCol='features')

 

Everything so far should be identical for all the three ML models since they all share the same dataset. After this point, however, the workflow will be slightly different for different models in terms of hyperparameter-tuning. For Random Forest Classifier and Logistic Regression Classifier, instantiate the ML model and attach it to the end of our pipeline:

from pyspark.ml.classification import RandomForestClassifier
from pyspark.ml import Pipeline
rf = RandomForestClassifier(labelCol='label', featuresCol='features', 
                            numTrees=30, maxDepth=30)
pipeline = Pipeline(stages=[arr_indexer, org_indexer, des_indexer,
                            onehot, assembler, rf])

 

For Logistic Regression Model, switch to the following lines of code to instantiate the pipeline:

from pyspark.ml.classification import LogisticRegression
from pyspark.ml import Pipeline
lr = LogisticRegression(maxIter=30)
pipeline = Pipeline(stages=[arr_indexer, org_indexer, des_indexer, onehot, assembler, lr])

 

During training, it is common practice to implement K-fold Cross-validation and Grid Search. Fortunately, PySpark provides sub-modules for both of them as well. To build a parameter-searching grid, use the ParamGridBuilder() sub-module and pass the parameter map into the CrossValidator(). The code for the Random Forest Classifier is shown here:

from pyspark.ml.tuning import CrossValidator, ParamGridBuilder
param = ParamGridBuilder().build()
cv = CrossValidator(estimator=pipeline,
                    estimatorParamMaps=params,
                    evaluator=evaluator,
                    numFolds=5)

 

Do not pass any parameters when building the parameter map for Random Forest (RF) Classifier since its parameters are more straight forward and easy to interpret. The hyperparameters for RF are initialized to the best case when instantiating the RF object previously. Specifically, numTrees stands for the number of single decision trees make up the RF classifier and maxDepth stands for the maximal depth for each tree in the forest. RF does not require a lot of tuning since the higher the values of these two parameters, the better the performance of the model in terms of test-accuracy. On the other hand, of course, with more, deeper trees, the training process takes linearly more time.

Whereas for Logistic Regression (LR) Classifier, build the parameter grid by doing the following:

from pyspark.ml.tuning import CrossValidator, ParamGridBuilder
params = ParamGridBuilder()
params = params.addGrid(lr.regParam, [.01, .1, 1, 10]) 
               .addGrid(lr.elasticNetParam, [0, .5, 1])
params = params.build()
cv = CrossValidator(estimator=pipeline,
                    estimatorParamMaps=params,
                    evaluator=evaluator,
                    numFolds=5)

 

As shown in the code, there are two parameters to be tuned for LR classifier. lr.regParam is the regularization coefficient which specifies that, in the range of 0.01 to 10, how much the model should be regularized. The bigger the number, the more the model is regularized. Meanwhile, lr.elasticNetParam specifies the way the model is regularized, 0 means it’s an L2 penalty and 1 means it’s an L1 penalty. And 0.5 means the penalty is a combination of both L1 and L2.

For Neural Network Classifier, the story is a little bit different. Since Spark Mllib does not have predefined models for ANNs, instead, use keras (with tensorflow as backend) to build the ANN model for comparison. Also, keras Dense layer does not take sparse vectors as input, therefore just convert the dataset in a conventional way. In other words, creating dummy variables for categorical data:

import pandas as pd
from sklearn.model_selection import train_test_split
from keras.utils import to_categorical
# convert Spark DataFrame into Pandas DataFrame
df_pd = df.toPandas()
# create feature and label set
X = df_pd.drop(['label', 'ARRIVAL_DELAY'], axis=1)
y = df_pd.label
# convert categorical data into numerical data
X = pd.get_dummies(X)
# split into training and test set
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2,
                                                    random_state=5)
# reshape the dataset to fit ANNs
X_train = X_train.values[...,None]
X_test = X_test.values[...,None]
# convert the label into categorical data
y_train = to_categorical(y_train.values)
y_test = to_categorical(y_test.values)

 

Use the following code the initialize the ANN model:

from keras.models import Sequential
from keras.layers import Dense
from keras.layers import Flatten
from keras.layers import BatchNormalization
ann = Sequential()
ann.add(Dense(10, activation='relu', input_shape=(X_train.columns.size, 1,)))
ann.add(BatchNormalization())
ann.add(Flatten())
ann.add(Dense(2, activation='softmax'))
ann.compile(optimizer='adam', metrics=['accuracy'], loss='categorical_crossentropy')

 

A simple ANN structure that only has one hidden layer with ten hidden nodes is modeled for comparison. After the hidden layer, a BatchNormalization() layer normalizes the output of the hidden layer. Since the numerical scales of the features are different, the normalization layer makes it easier for the model to converge.

The last step is to train our models. For ANN, call the following line of code: 

from keras.callbacks import ModelCheckpoint
model_path = 'ann_keras_model'
model_check = ModelCheckpoint(model_path, save_best_only=True)
history = ann.fit(X_train, y_train,
                  batch_size = 200, epochs = 25,
                  validation_split = 0.2,
                  callbacks=[model_check])

 

Add a callback function for the training so that only the best model during training will be saved. For the Random Forest and Logistic Regression models, simply call the following line for training:

cv = cv.fit(flights_train)

Model Evaluation

For a binary classification problem, it is common practice to evaluate the performance of the model by plotting the ROC (Receiver Operating Characteristic) curve or calculating precision, recall, and F1-scores. For Random Forest classifier and Logistic Regression model, use the following code to calculate the AUC (area under ROC curve), precision, and recall, etc:

from pyspark.ml.evaluate import BinaryClassificationEvaluator
# get the prediction
prediction = cv.transform(flights_test)
# initialize the evaluator
evaluator = BinaryClassificationEvaluator()
# calculate AUC
auc = evaluator.evaluate(prediction, {evaluator.metricName: 'areaUnderROC'})
print('AUC: %0.3f' % auc)
# compute TN, TP, FN, and FP
prediction.groupBy('label', 'prediction').count().show()
# Calculate the elements of the confusion matrix
TN = prediction.filter('prediction = 0 AND label = prediction').count()
TP = prediction.filter('prediction = 1 AND label = prediction').count()
FN = prediction.filter('prediction = 0 AND label <> prediction').count()
FP = prediction.filter('prediction = 1 AND label <> prediction').count()
# calculate accuracy, precision, recall, and F1-score
accuracy = (TN + TP) / (TN + TP + FN + FP)
precision = TP / (TP + FP)
recall = TP / (TP + FN)
F =  2 * (precision*recall) / (precision + recall)
print('n precision: %0.3f' % precision)
print('n recall: %0.3f' % recall)
print('n accuracy: %0.3f' % accuracy)
print('n F1 score: %0.3f' % F)

 

Furthermore, take a look at the parameters of the best Logistic Regression model selected by Grid Search:

For the Logistic Regression model, the Grid Search chooses to have a combination of L1 and L2 penalty to get the best performance. And it also chooses to use the least regularization (0.01). For the Neural Network model, Figure 6 (left) shows the training history about training loss, validation loss, training accuracy, and validation accuracy. The training takes 15 seconds for a single epoch on a GPU implemented device. The right-hand figure shows the ROC curves for the three models.

Figure 6. Training history (left) for ANN and ROC curves (right).
Figure 6. Training history (left) for ANN and ROC curves (right).

Call .transform() method to the held-out test set to further evaluate these models :

from pyspark.ml.evaluation import BinaryClassificationEvaluator
# get predictions
prediction = cv.transform(flights_test)
# Calculate the elements of the confusion matrix
TN = prediction.filter('prediction = 0 AND label = prediction').count()
TP = prediction.filter('prediction = 1 AND label = prediction').count()
FN = prediction.filter('prediction = 0 AND label <> prediction').count()
FP = prediction.filter('prediction = 1 AND label <> prediction').count()
# show confusion matrix
prediction.groupBy('label', 'prediction').count().show()
# calculate metrics by the confusion matrix
accuracy = (TN + TP) / (TN + TP + FN + FP)
precision = TP / (TP + FP)
recall = TP / (TP + FN)
F =  2 * (precision*recall) / (precision + recall)
# calculate auc
evaluator = BinaryClassificationEvaluator()
auc = evaluator.evaluate(prediction, {evaluator.metricName: 'areaUnderROC'})
print('n precision: %0.3f' % precision)
print('n recall: %0.3f' % recall)
print('n accuracy: %0.3f' % accuracy)
print('n F1 score: %0.3f' % F)
print('AUC: %0.3f' % auc)

 

The code above prints the confusion matrix of the prediction and also calculates all metrics. And the results are summarized as a table shown below:

table1

 

From the table, all three models can predict the label pretty accurately. While Neural Network gets higher accuracy, its precision and recall are relatively lower. After all, the overall difference between these three models is relatively small. 

Conlusion

With all the models providing similar performances, it is fair to say that, for machine learning, sometimes it is the data that dominates the solution, instead of the model. In other words, with enough dataset, all models can achieve a pretty good score. With that being said, the focus shifted from data to approach in Big Data era. What’s important is no longer the result, but the way to approach the result.

For our experiment, the flights’ dataset contains more than 2 million records, which can take a long time for stand-alone devices to implement the machine learning application. Apache Spark provides an automatic distributed clustering-computing framework, under which a simple Logistic Regression can achieve almost identical performance as a GPU implemented Neural Network model with shorter training time. 

In this post, I talked about how to use PySpark to build machine learning pipelines which are suitable for Big Data analysis. With a smaller size of data, though, using standard machine learning library should be sufficient and more efficient. Nevertheless, Apache Spark is definitely a reliable tool for solving challenging machine learning problems using Big Data.

The author is still learning Apache Spark, and there might be error or mistakes in this post. Welcome to correct me if I’m wrong and let me know if you have any questions!

Thanks for reading 🙂 

, , , ,

Deep Learning with Phase Retrieval

Introduction

I spent the second year and half of the third year during my Ph.D. program to develop Chemical Holographical Imaging System (CHIS). It is designed to acquire not only the intensity of the light (electromagnetic field) but also the phase information, which is usually lost when using a standard camera or detector to collect the photons. I will write another blog to explain the mechanism behind the scene. Here, I will talk about the downstream work, which happens after the phase of the field is reconstructed.

To characterize and prove the concept, we implemented a forward model and an inverse model for CHIS, simulating the process of acquiring a phase-sensitive measurement and reconstructing the sample, respectively. Here I will only talk about the inverse model, where machine learning algorithms are applied.

We use the Mie scattering model to simulate a micro-sphere whose radius is close or longer than the wavelength of the incident light. Mie scattering is the light scattered when electromagnetic wave interacts with uniform isotropic particles (spheres) in a uniform isotropic dielectric infinite medium. Mie theory is a rigorous solution to Maxwell equations that describes the propagation of electromagnetic radiation by spherical scatterers. It has been widely applied in the characterization of atmospheric absorption and cellular structure detection in spectroscopic imaging of biological tissues.

According to the superposition principle, any form of propagating electromagnetic field is equivalent to a combination of multiple plane waves. To simplify the model, the incident light going through a condenser with a low numerical aperture is considered as a plane wave. While the forward model of CHIS can be described by Mie scattering, the inverse model is equivalent to solving the inverse Mie problem.

Here we first examine the forward Mie scattering. The scattered field for a single plane wave produced by a sphere placed at the origin with radius a can be written as:

\textbf{E}_s(\textbf{p}) = \sum _{l = 0} ^ \infty B_l(\lambda, n, a) h_l^{(1)} (kr)P_l(cos\theta)

where h_l^{(1)}(\cdot) is the order-l spherical Hankel function. The coefficients B that couple the internal and external fields are given by:
B(\lambda, n, a) = (2l+1)i^l\frac{j_l(ka)j_l^{‘}(kna)n-j_l(kna)j_l^{‘}(k)}{j_l(kna)h_l^{(1)}(ka)-h_l^{(1)’}(ka)-h_l^{(1)}(ka)j_l^{‘}(kna)n}

where j_l^{‘}(\cdot) is the first derivative of the spherical Bessel function of the first kind. And h_l^{(1)’}(\cdot) is the first derivative of the spherical Hankel function.

If we take a closer look at these two equations, it is not difficult to notice that the material properties of the sphere (radius a and refractive index n) are only embedded in the B coefficient. In other words, to build an inverse model of the system, our first step is to extract the B coefficient from the field \textbf{E}_s.

Since the Hankel function and its derivatives can be derived from the Bessel functions of the first kind, the Bessel functions of the second kind, and their derivatives, the B coefficient behaves similarly to Bessel functions. In Figure 1, B coefficients of different scattered fields scattered by different spheres are plotted (only real part is shown).

Figure 1. B Coefficients of the scattered field with different sphere radius
Figure 1. B Coefficients of the scattered field with different sphere radius

The animation shows how the B coefficient changes with a \in [1, 10] \mu m. You can run this Jupyter Notebook demo to examine the property of B coefficient with Bokeh interactive plot implementation. 

The point is, the B coefficient trends to converge to zero after going up and down along increasing order. In other words, it is highly sparse and fluctuating, which makes a lot of sense: because the incident light is scattered by a sphere, the scattered field around the surface of the sphere would have the highest intensity. Then the intensity goes up and down due to scattering and interference. When the scattered field travels further away from the sphere, its intensity converges to zero eventually (scattered into broader space). In fact, on top of the sparsity of the B coefficient, the rest of the equation is highly sparse and fluctuating as well. If we consider the scattered field as a system of linear equations where the B matrix is the input and the rest of the equation is the transferring matrix and the field is the output, we can get the following linear system: Ax = b where A is the scatter matrix, A = h_l^{(1)} (kr)P_l(cos \theta) and  x is the B coefficient and b is the field E_s. If we calculate the condition number of the linear system, it can go up to 10^{21}.

In numerical analysis, the condition number with respect to an argument of a function stands for the amount of change to the output with respect to a small shift in the input argument. A low condition number means the system is well-conditioned thus more robust to the disturbance in the input. In contrast, a high condition number means that the system is sensitive to disturbance. A small change in the input could cause a dramatic change to the output hence an ill-conditioned system. A linear system with condition number \kappa = 10^k means besides the arithmetic error during the calculation, another k digits of loss of accuracy would be induced.

In our case, if we apply traditional analytical methods such as Pseudo-inverse or Gaussian elimination to solve the inversed linear system to get B coefficient from E_s (which is our measurement), the ill-conditioned property makes the recovery with noise impossible since a tiny little bit of noise in the measurement would result in massive change in the calculation of B. Another way to think about it is the error is the noise multiplied by the condition number.

Deep Learning Approach

One pathway to the solution is broken does not mean the answer is unreachable. Before we jump into machine learning, let us take a look at some of the advantages and disadvantages of different strategies:

  1. Analytical solutions are mathematical grounded. We can break down the solution into pieces and prove each of part of it rigorously. We know what we are doing and what should we expect.
  2. Analytical solutions are universal and scalable. Top level applications are grounded with similar fundamental theories, the analytical solution (if we have one) of our system can be applied to any other similar imaging systems.
  3. Analytical solutions are one-time work. They can be applied directly once the strategy is finalized whereas machine learning methods are not. A different training set has to be provided if we want to use the model for another system even if there is only a little difference, for accuracy.
  4. Machine learning is a black box, which is well-known in the community. We do not have the detailed information of what our models are doing when they give the prediction, either do we know what to expect when we are trying to tweak some fine gears (not excluding the ability to judge if the model has a high bias or high variance, of course).

At this point you might want to ask, seems like there must be better choices than using machine learning algorithms to solve the ill-conditioned inverse problem, why don’t you try something else? The answer is Yes, we did try something else. And that something is called Compressive Sensing, which is specialized for sparse reconstruction and convex optimization. The short story is 1. the linear system we have is not a convex problem, 2. accuracy of sparse reconstruction algorithms seems to suffer from the synthetic noise, 3. they all require training somehow.

Look over all the pros and cons of all methods available, we stick with the machine learning (ANN specifically, more on later) algorithms despite all the disadvantages they might have: for proof of concept, we believe a better accuracy should be weighted higher than scalability and the ability to be generalized; the training process can be optimized in various ways such as feature selection and smart sampling.

Convolutional Neural Network

Even though Convolutional Neural Networks (CNNs) were invented during the 1980s, they were not so popular until the GPU implementation was done and generalized at the early 2010s. They are mighty in terms of image processing with the function of “convolving” the images to extract higher dimensional features. One would normally take CNN as the first choice to solve the problem (at least we did). The structure of the CNN we applied to our system is shown in Figure 2:

Figure 2. CNN Structure
Figure 2. CNN Structure

As shown in Figure 2, a convolution and max-pooling based set of layers are adopted as the basic structure. Three pairs of convolutional layer and max-pooling layer are cascaded to enable the extraction of higher level abstract features from the scattering images. The number of filters for each convolutional layer is 128, 64, and 32, respectively. The window sizes are set to 3\times 3 and 2\times 2 for convolutional layers and the sequential max-pooling layers. A fully connected layer with 128 units is attached to the last max pooling layer followed by the output layer. Given a complex field as the input for the trained CNN model, the final output is represented by the three units in the output layer corresponding to the radius a, the refractive index \eta, and the attenuation coefficient \kappa, respectively.

Since we have the forward model implemented for the system, one of the advantages is that we do not need to worry about the amount of data we could use for training. However, the trade-off between model structural complexity and training time cost is not neglectable. The three-layer structure is a result of leveraging accuracy and time cost. 

CNN performance

To ensure the CNN is trained with low bias and low variance, 8000 phase-reconstructed images with 640 \times 640 spatial resolution are generated by the forward model with different material properties for the simulated spheres. Each of the feature vectors has 20 entries uniformly distributed within their predefined range. Specifically, the radius a of the sphere is ranging in [1, 2]\mu m, the refractive index \eta is varying in the range [1.1, 2.0], and the attenuation coefficient \kappa \in [0.01, 0.05]. 

Before training the model, stratified sampling is applied to split the training and testing data set since the overall data set is too large to fit in memory. The sphere radius clusters the 8000 images into 20 subsets. Within each subset, the training set and test set are randomly sampled and shuffled with a 20\% ratio, and all images are cropped to 128 \times 128 resolution during this stage. As a result, there are 6400 and 1600 images in the training and test set, respectively. And the order of the features is randomized for \eta and \kappa but ascending for a.

During the training, the training set is further divided into training and validation set with the same ratio (20% for validation), which results in 5120 images for training and 1280 images for validation. The training process contains 25 epochs (i.e., the training and cross-validation are repeated for 25 times) and the batch size is set to 50. Meanwhile, label scaling is applied to the attenuation coefficient set due to scale differences.

Figure 3. Relative RMSE of (left) complex CNN (with phase) and (right) intensity CNN (without phase)
Figure 3. Relative RMSE of (left) complex CNN (with phase) and (right) intensity CNN (without phase)

Figure 3 shows the performance of the CNN trained on 6400 complex images and tested on the rest of the 1600 images. The plots only show 50 samples uniformly down-sampled from the testing and prediction set by a factor of 32. The relative root-mean-square error (RMSE) \epsilon_R is calculated by:

\epsilon_R = \frac{\sum_{m=1}^M\frac{\sqrt{(y_{pred}^{(m)} – y_{test}^{(m)})^2}}{y_{test}^{(m)}}}{M}

where y_{pred} is the prediction given by the CNN and y_{test} is the groud truth (label) of the corresponding sample. And the total number of testing samples is denoted by M = 1600.

To demonstrate the merits of phase-sensitive measurements achieved by chemical holography, another CNN model is trained with traditional intensity images with an identical structure and the same number of training and testing samples. The intensity images are acquired by calculating the square of the absolute value of the complex images. And the performance of the intensity CNN is shown on the right (Figure 3).

Everything seems beautiful if we stop here: not only the CNN is working with a pretty lovely regression accuracy; also the one trained with complex images (phase reconstructed) is outperforming the one trained with standard intensity images. To generalize our model, or, to demonstrate the robustness of CHIS, we added more experiments to back it up.

Table 1. Relative Root Mean Square Error (RMSE) for CNNs (%)
Table 1. Relative Root Mean Square Error (RMSE) for CNNs (%)

Table 1 shows the performance of complex CNNs and intensity CNNs summarized from 30 repetitive trails. It is clear that the complex CNNs perform better than intensity CNNs, which can be partially explained by the fact that the additional phase information existing in the complex training set is no longer available for the intensity training set (during the transaction from complex to intensity).

Artificial Neural Network

It seems a little bit weird seeing Artificial Neural Network (ANN) being mentioned at this point since, for most of us, CNN is far more advanced and powerful than ANN. We already implemented the model with CNN, why bother even talking about ANN? Well, one of the most useful lessons I learn from my advisor Dr. David Mayerich is that, as a scientist, or an engineer, we often need to calm down and take a step back. Essentially, when we are doing research or trying hard to get things done, it is quite natural for us to be tunnel-visioned. Once we stuck, it is relatively easy to get too extreme and anxious as well. Only if we take a step back and make sure everything basic has been done correctly, the path to move on would show up.

Back to the topic, why we choose ANN? The answer can be as simple as, we don’t need CNN because we can achieve the same performance with ANN, even faster, easier, and better. And I will explain the reasons explicitly and rigorously.

First of all, CNNs are widely applied to image processing to take care of complicated high dimensional feature extraction. The goal is to train the network so it can distinguish features or patterns from a 2-D image, or even a 3-D volume and CNN is extremely good at this. Now let’s consider the following statement:

We have a dataset and a model trained with this dataset, if the model is CNN and it is the best model we tried, then this dataset must contains two (or higher) dimensional images.

This statement is true (at least in most of the cases, welcome to correct me if I am wrong). But how about the converse of the statement:

If we have a dataset that contains two (or higher) dimensional images, then CNN will be the best model for this dataset.

Is this converse statement also true? There are numerous examples where it does not hold, including our case. Consider an image with all the pixel values in the x axis being the same number and only varies along the y axis, then we can “compress” this 2-D image into a 1-D vector with one of the vertical lines without losing any information. Furthermore, if all the pixel values are the same, we can even compress the whole image into a single pixel, and this is where the symmetricity comes into play. In highly symmetrical situations, such as center symmetry, radial symmetry, or even reflectional symmetry, the more symmetrical the original 2-D image is, the more redundant information it has thus, the less spatial diversity encoded in the local pixels. If there is not much or no higher dimensional features to be extracted, a simpler model such as ANN can take care of all the work adequately. 

In our case, assume there is no systematic noise or the noise is low, the whole image will be radially symmetrical if we place the sphere at the center of the image. Hence we can reduce the dimension of our dataset from 2-D images to 1-D vectors to train the model without losing any spatial information.

Once we switched from CNN to ANN using 1-D vectors reduced from the 2-D images, several obvious advantages appear naturally:

  1. Dramatically reduced data size, computational cost, and time cost: if the original 2-D simulation has a N\times N resolution, there are N^2 points to be calculated. Once been reduced to 1-D, the number of calculations decreased into N/2 (one vertical or horizontal line from half of the image). Therefore, the dataset for training and testing for the 1-D ANN model is 2N times smaller than the 2-D case. Previously, our training image resolution after padding is 640\times 640 with a total size of more than 100 GB of data, and now it will be shrunk into 1/1280 of the original size which is no more than 100MB, while without losing any information. This 2N relationship also applies to the amount of computational cost and time cost reduced.
  2. Simpler network model and more sophisticated training: as one of the more fundamental machine learning models, ANNs are more mathematically grounded, especially for single layer ANNs. The weights between layers can be interpreted more concretely. That is, for each pixel in the input 1-D vector, there is a combination of weights (the hidden layer) being optimized to minimize the cost function, including the bias term. The more units in the hidden layer, the more “fit” the model will be (also more chance of overfitting). This point-to-point mapping is less abstract than CNN’s higher dimensional feature extraction mechanism thus, easier to be controlled. With that being said, ANN models can be trained more thoroughly with much less work required compares to CNN models.

Figure 5 demonstrates the radial symmetry of the simulation and the effectiveness of the dimension reduction. As we mentioned before, for the simulation to be radially symmetrical, the sphere must be placed at the center of the field of view. One thing to keep in mind is that this condition is hard to meet in real data. In other words, in the experimental dataset, the error induced by non-symmetry would impact the performance of the proposed pipeline. However, this drawback also applies to other models such as CNNs. Adding more synthetic data for training can be a potential solution.

Figure 5. (a) 2-D Mie scattering simulation in far-field (top) and near-field (bottom). (b) 1-D far-field Mie scattering simulation corresponds to the red dashed line in (a). (c) 1-D near-field Mie-scattering simulation corresponds to the green dashed-line in (a).
Figure 5. (a) 2-D Mie scattering simulation in far-field (top) and near-field (bottom). (b) 1-D far-field Mie scattering simulation corresponds to the red dashed line in (a). (c) 1-D near-field Mie-scattering simulation corresponds to the green dashed-line in (a).

With one-dimensional data, the structure of the model can be simplified dramatically, as shown in Figure 6. As mentioned before, the 1-D dataset only occupies 1/1280 of the memory space compares to the 2-D dataset; hence the bottom neck of memory limit no longer exists. The dimension of each feature is increased to 30 results in 30^3=27000 total number of samples. Each sample contains a 320\times 1 vector as the input for the neural network. 

Figure 6. ANN structure
Figure 6. ANN structure

In terms of fine-tuning ANN model, there are several things to pay attention to:

  1. Unlike CNN where processes like convolution and max-pooling generalize the local differences in the data, ANN takes all the local information and arrange combinations of weights to each of them. As a result, ANNs with more complicated structures are more likely to be struggling with overfitting, which becomes more noticeable when applying deep ANNs to real-life data where lots of noise are presented. In a word, complicated ANNs with higher accuracy are more sensitive to noise than simple ANNs assuming that they are trained without noise. How to leverage on the trade-off between noise sensitivity and accuracy is yet to be determined by the problem to be solved.
  2. Random initialization. To ensure that the number of stochastic processes is minimized while developing a model, each operation involving randomness, such as weights initialization, should be assigned with a seed or random state. 
  3. Activation function. While the most popular activation function being ReLu (Rectified Linear units), to optimize the training process, other activation functions should be considered, too. One thing to keep in mind is the dying ReLu problem. Since ReLu sets everything negative to zero, the gradient can converge towards zero. Once a neuron (unit) hits this state, the corresponding weights will stop getting updated thus so-called dead neuron. In our case, there is a fair chance that the weighted sum can hit negative numbers. Therefore a tanh (hyperbolic tangent) activation function would be a better choice. In fact, the training process does converge significantly faster with tanh function than sigmoid or ReLu when the network has a small number of hidden units. When the network grows deeper or has more hidden units, ReLu would be a better choice. If a lot of units stuck at the edge of sigmoid or tanh functions, vanishing gradients appears which can make the training unstable and more difficult to converge. On the other hand, Dying ReLu would not be a problem here since the number of units is big enough.
ANN Performance

Starting with a simple model as shown in Figure 6, there is only one hidden layer and five hidden units. And the number of parameters of the intensity ANN is shown below.

For complex cases, since the 1-D vector contains both real and complex components, it is split into two channels for the training. The model has the same structure with 1 hidden layer and 5 hidden units. And the number of parameters in the complex model is shown below.

As shown in the tables, the number of parameters only increased by 5 for the complex ANN, which is equal to the number of extra channels added to the input layer. The model is trained with Keras package (GPU implemented) and 400 epochs, which takes around 2 minutes for 19440 training samples (19440 =  27000 \times 10% \times 20% , 10% for test set and 20% for validation set). And the results averaged from 30 repetitive trials are shown in table 2.

Table 2. Relative RMSE for simple ANNs (%)
Table 2. Relative RMSE for simple ANNs (%)

Before we judge the performances of ANNs and CNNs, let’s compare the number of parameters in both models first. Recall that the number of parameters for the complex ANN is 4818, which is not a small number. However, it is just the tip of an iceberg compares to complex CNN, of which the number of parameters is shown below:

Well, there are more than 1.7 \times 10^6 parameters in the CNN model even with max-pooling layers. In contrast, 5\times 10^3 is a small number of parameters. While using no more than 1/1000 of the parameters, ANNs are able to get only several percentages of accuracy loss, which is the point we want to point out here: 

Choosing the correct model is the first priority in machine learning.

If you are not convinced yet, totally understandable. One of the questions could be: the accuracy we are getting with the current ANNs is still worse than the CNN version, which is unwanted if we do not care about resources cost to get the model trained. In other words, what if we simply want to achieve the best accuracy, should we switch back to CNN?

The answer is, not yet. The ANN we currently have is just the simplest version. Instead of turning into CNN again, what we can do is to make our ANN model more complicated, such as adding more hidden layers, or more units at each layer. Such as a deeper ANN shown in Figure 7:

Figure 7. A more complicated ANN with 3 hidden layers and more hidden units in each layer.
Figure 7. A more complicated ANN with 3 hidden layers and more hidden units in each layer.

And the number of parameters of this new ANN is shown below:

Without the max-pooling layer and dropout layer, the number of parameters in the new ANN is significantly increased due to the last fully connected layer. But still, the training process costs less time compares to the training time of a deep CNN (no convolution needed for ANN). In this case, the network was trained with 100 epochs in 13 minutes with an Nvidia GeForce GTX 970 GPU. And its performance is shown in Table 3.

Table 2. Relative RMSE for deep ANNs (%)
Table 2. Relative RMSE for deep ANNs (%)

In this single experiment, ANN trained with phase-reconstructed 1-D data outperforms CNN trained with 2-D images, despite the comparable numbers for the intensity ANN accuracy, which in return demonstrates the merits of phase retrieval.

Conclusion

Machine learning is a profound and sophisticated discipline that should be carefully considered and examined before implementing it. ANN seems weaker than CNN at first glance, or at least the author was having some trouble giving up CNN until the ANN results were carried out.

In other words, model selection should be placed at the top with priority before taking any action. Only after taking every factor related to a specific problem into account (even considering changing the strategy of data generation), should the model be determined.

There might be error or confusion exist in this blog, feel free to correct me if I am wrong or let me know if you have any questions.

At last, the code for 

  • Generating the Mie scattering data can be found here.
  • Training and evaluating the CNN model can be found here. 
  • Training and evaluating the ANN model can be found here. 

Thanks for reading 🙂