My First Kaggle Competition -- Plant Seedlings Classification

Kaggle.com is a website famous for hosting machine learning competitions. Recently I participated one competition named plant seedings classification (https://www.kaggle.com/c/plant-seedlings-classification/) and would like to share what I have learned here. The goal of this competition is for image classification. For training set, the host provides over four thousand images including 12 different plant species. The goal is to categorize the test image into one of these 12 species. This problem is a typical computer vision challenge which can be solved by deep learning. The neural network used is convolutional neural network (CNN). After CNN network trained by training images, it will be applied to test images for classification.

In order to achieve good performance, I tried various approaches and found the following two tips to be particularly useful:
1. Start from pre-trained network, and retrain it using the new data.
2. Ensemble by combining results of multiple models.

Before diving into more details, let me list the results of some approaches that I have tried.

Tip 1 is not my first choice. In fact, I started from designing a 10-layer network (conv/conv/pool/conv/conv/pool/conv/conv/pool/dense) by myself (row 'My own CNN network' in the table with 0.95717 accuracy). But then I found that if building our classifier based on some well-designed deep networks like VGG, Xception and Densenet, better accuracy can be achieved. Explaining what is VGG or Densenet is not the scope of this article. However, it is fair to say that these networks have more layers than the network I designed and also some of them such as densenet or InceptionV3 apply more exotic topology. Training these networks from scratch can be challenging, but we can start from pre-trained weight and after that it is much easier to train. An example of Keras code is shown below for building a model based on pre-trained Xception network. The base model is initialized from a Xception network with pre-trained weight based on Imagenet but not including top part. Then a pooling layer and a dense layer are added upon the base model. The output of the dense layer predicts the likelihood of which plant the image belongs to. This example can be easily extended to other network such as VGG. Weight setting corresponding to the best result are saved and used for testing. Using the pre-trained network boosts accuracy by at least two percents compared with self-designed networks in the first row. The best accuracy of single network is achieved by densenet201 and densenet169 as 0.98236.

base_model = Xception(weights='imagenet', input_shape=(img_size, img_size, 3), include_top=False)
x = base_model.output
x = GlobalAveragePooling2D()(x)
x = Dropout(0.5)(x)
x = Dense(1024, activation='relu')(x)
x = Dropout(0.5)(x)
predictions = Dense(12, activation='softmax')(x)
model = Model(inputs=base_model.input, outputs=predictions)

model.compile(optimizer='Adadelta',
              loss='categorical_crossentropy',
              metrics=['accuracy'])

model.fit_generator(datagen.flow(train_img, train_label, batch_size=batch_size), steps_per_epoch=len(train_img)//batch_size, epochs=100, verbose=1)
model.save_weights('Xception.h5')

To further improve performance, we need tip 2. If you ask machine learning practitioners the best way to improve performance, many of them probably tell you "ensemble". Ensemble method is to combine together the prediction results of multiple models. Since the prediction results generated by these models are more or less independent, combining them together reduce the noise and enhance the results. However, the non-trivial thing about ensemble is that combining more models does not mean better results. For example, in our experiments Densenet201 + Densenet169 + Densenet121 + InceptV3 + Xcep + InceptRosV2 + VGG16 + VGG19 has the most combining and not the best accuracy. Without straightforward rule, finding the best ensemble model takes lots of trials. Another problem of ensemble approach is how much weight should be assigned to each base model. Not all base models are equal; some base models behave better than others. Most likely base model with better performance should be assigned larger weight. In communication, there are schemes such as Maximum-Ratio Combining (MRC), which combines multiple models and determines combining weight using information of signal and noise level. But in deep learning, it is not easy to obtain metrics such as signal and noise level. Thus in my work, I assign each weight to each base model. But please put in your mind that other ways of weight assignment in deep learning do exist. After many trials, eventually we found that Densenet201 + Densenet169 + InceptV3 + Xcep  +InceptRosV2 gives the best performance (0.98992). I should say that by ensemble method alone, this ensemble model gives the accuracy of 0.9874. But after analyzing multiple submissions and correlated with their results, we are able to correct two entries which boosts the final accuracy to 0.98992. So we game the system a little bit but not too much. This score of 0.98992 ranks top 3% accuracy. The top player of leader board has accuracy of 1.0!

In summary, participating the competition gives lots of fun. Although we participated after the competition closed and our score does not count, working on this problem still teach us quite a few things. We also learned from many discussions and kernels in Kaggle website. Here we want to pay back by sharing our experience and hopefully it can benefit our audience.

[Update on 11/16/2018] Code examples have been added to https://github.com/legendzhangn/blog/tree/master/plant_seedlings

What Keras Source Code Tells Us about LSTM Processing

Recurrent Neural Network (RNN) is a type of neural network which can save its states and process input sequences. With this recurrent architecture, RNN can be used for applications such as language analysis, time series prediction and speech recognition. Sometimes it can be challenging to understand how RNN processes it inputs. I met with this problem when I started to learn RNN and reading Keras source code helps to clarify this puzzle.

Under Keras, the input to RNN has three dimensions (batch_size, timesteps, input_dim). batch_size controlled how often the network is updated. The network is updated once after batch_size sequences are processed. For example, if there are 256 sequences and batch_size is set as 32, that means after processing the whole set, the network will be updated 8 times. The parameter of timesteps is how many times a RNN cell runs. In each run, it starts from saved states and generates new states and outputs. input_dim is the number of features contained in an input. Another important parameter is the number of RNN units. For example, the instruction below creates 4 units of LSTM while LSTM is a popular type of RNN. Here, timesteps for input is set as 1 and input_dim is set as 2. This instruction does not define batch_size and it means that batch_size will be defined later. From now on, we name num_units as the number of units.

model.add(LSTM(4, input_shape=(1, 2)))

Keras RNN source code can be found in here. Reading the source code reveals how RNN processes its input. How LSTM operates can be found in call() function of class LSTMCell. This part of code (shown below) are only related to input_dim and num_unit but not other parameters. For instance, the dimension of self.kernel_i/self.kernel_f/self.kernel_c/self.kernel_o is (input_dim, num_units). inputs is a vector of input_dim elements. Therefore, after dot product between kernel and inputs, input_dim does not exist any more and x_i/x_f/x_c/x_o become vectors with the length of num_units, which match with the dimensionality of bias. Another observation is that in this stage, each LSTM unit works independently without interference from other units.

            if 0 < self.dropout < 1.:
                inputs_i = inputs * dp_mask[0]
                inputs_f = inputs * dp_mask[1]
                inputs_c = inputs * dp_mask[2]
                inputs_o = inputs * dp_mask[3]
            else:
                inputs_i = inputs
                inputs_f = inputs
                inputs_c = inputs
                inputs_o = inputs
            x_i = K.dot(inputs_i, self.kernel_i)
            x_f = K.dot(inputs_f, self.kernel_f)
            x_c = K.dot(inputs_c, self.kernel_c)
            x_o = K.dot(inputs_o, self.kernel_o)
            if self.use_bias:
                x_i = K.bias_add(x_i, self.bias_i)
                x_f = K.bias_add(x_f, self.bias_f)
                x_c = K.bias_add(x_c, self.bias_c)
                x_o = K.bias_add(x_o, self.bias_o)

            if 0 < self.recurrent_dropout < 1.:
                h_tm1_i = h_tm1 * rec_dp_mask[0]
                h_tm1_f = h_tm1 * rec_dp_mask[1]
                h_tm1_c = h_tm1 * rec_dp_mask[2]
                h_tm1_o = h_tm1 * rec_dp_mask[3]
            else:
                h_tm1_i = h_tm1
                h_tm1_f = h_tm1
                h_tm1_c = h_tm1
                h_tm1_o = h_tm1
            i = self.recurrent_activation(x_i + K.dot(h_tm1_i,
                                                      self.recurrent_kernel_i))
            f = self.recurrent_activation(x_f + K.dot(h_tm1_f,
                                                      self.recurrent_kernel_f))
            c = f * c_tm1 + i * self.activation(x_c + K.dot(h_tm1_c,
                                                            self.recurrent_kernel_c))
            o = self.recurrent_activation(x_o + K.dot(h_tm1_o,
                                                      self.recurrent_kernel_o))

Code above has nothing to do with time concept, which is signature of RNN. So how is time getting introduced into RNN implementation? The class of RNN(layer) defined in the same script of recurrent.py is the base class for RNN. In call() function of this class, we can see a call to backend function K.rnn where timesteps is defined as input_length to this function. We call it backend function since Keras is based on tensorflow or Theano, and rnn() function is provided not in Keras but in tensorflow or Theano library depending on which one Keras replies. Although having not yet checked the source code of rnn function, I believe what happened is that in the backend function, the above logic of LSTMCell is called timesteps times to generate the final output.

        last_output, outputs, states = K.rnn(step,
                                             inputs,
                                             initial_state,
                                             constants=constants,
                                             go_backwards=self.go_backwards,
                                             mask=mask,
                                             unroll=self.unroll,
                                             input_length=timesteps)

For reference purpose, equations for LSTM are provided here:
Input gate: \(i_{t}=\sigma(W^{(i)}x_{t}+U^{(i)}h_{t-1})\)
Forget: \(f_{t}=\sigma(W^{(f)}x_{t}+U^{(f)}h_{t-1})\)
Output: \(o_{t}=\sigma(W^{(o)}x_{t}+U^{(o)}h_{t-1})\)
New memory cell: \(\tilde{c}_{t}=tanh(W^{(c)}x_{t}+U^{(c)}h_{t-1})\)
Final memory cell: \(c_{t}=f_{t}\circ c_{t-1}+i_{t}\circ \tilde{c}_{t}\)
Final hidden state: \(h_{t}=o_{t}\circ tanh(c_{t})\)
whereas \(\circ\) means pointwise operation


In this blog, I summarized what I have found in Keras source code. Hopefully it can be useful to you.

The largest difference between ordinary RNN versus LSTM is that LSTM has an extra hidden state.

How To Load/Save Data in Python

Loading/Saving data is a common task in numerical analysis. In python, the most convenient way is to use the save/load function provided in numpy package. An example is as the follows

import numpy as np
rand_array = np.random.rand(1,10);
np.save('rand.npy',rand_array);
b=np.load('rand.npy');

This example shows how to save an array in a .npy file and then load. This save/load function is quite similar to the save/load in Matlab except for that Matlab allows saving multiple arrays in single file while in Python one file can hold only one array. For comparison purpose, Matlab save/load code is shown below where two arrays with random numbers are save in the same file

rand_array=rand(1,10);
rand_array2 = rand(1,10);
save('rand_file','rand_array','rand_array2');
myFile = load('rand_file');

3D LIDAR DIY Using PulsedLight LIDAR-Lite Device and Servo

LIDAR stands for LIght Detection And Ranging. LIDAR uses laser to measure the distance. LIDAR is capable to do things no many other sensors can do. It can accurately measure distance (tens to hundreds of meters) of objects around, which generates a good picture of surrounding environment. Because of this functionality, LIDAR is a must-have for many self-driving car designs. Many tech hobbyists want to put their fingers on a LIDAR device. However, commercial LIDAR devices are quite expensive and they are usually well above one thousand dollars. Good news is that with availability of cheap laser component like PulsedLight LIDAR-Lite and open source platform like Arduino, now we can make by ourselves a 3D LIDAR scanner with cost well below 1000$. In this blog, I will introduce such a device made by myself.

PulsedLight LIDAR-Lite is used in this device to measure the distance using laser. The maximum distance LIDAR-Lite measures is 40 meters. Two servos, one for vertical direction and another for horizontal direction, are used for scanning the 3D space. Both horizontal and vertical servos can rotate up to 180 degrees. This defines the coverage of this scanner. Arduino platform controls all of them. By connecting Arduino to PC, 3D scan data can be downloaded to PC to control point cloud file. At the same time, PC can power the device through USB. Thus the device does not need to have its own power. Figure below shows how the device looks like.

























When the device scans, it first fixes vertical angle and sweeps horizontally. Then it moves one vertical step and sweeps horizontally again. This Youtube video shows the scanning process

After scanning finished, a point cloud will be generated. Often several point clouds need to be stitched together to generate the whole picture for a building. The residential house model below is generated in that way.

3D model can be produced based on point cloud (Thanks Igor!). Below is the 3D model for the same residential house.

After 3D printing, we have a plastic model of the house. We can view it together with Google street map side-by-side.

Signal Processing Magic (3) -- Farrow Structure

Farrow structure was invented by C. W Farrow. It minimizes the number of multiplication used in fractional interpolation. Farrow structure is often deployed when there is fractional sampling rate change in the signal path. For example, ADC rate is set as 120MHz but baseband signal processing only supports 100MHz. Therefore, a filter needs to be designed for supporting rate conversion from 120MHz to 100MHz. Farrow structure is one candidate for this kind of job.

The setting of Farrow depends on the polynomial it supports. Let us take an example of using cubic Lagrange for fractional interpolation. The output y[n] can be written as combination of input x[n] as:
\[y[n]=C_{-2}x[n-2]+C_{-1}x[n-1]+C_{0}x[n]+C_{1}x[n+1]\]

Assuming fractional delay is \(\tau\), using Lagrange formula, \(C_{-2}\) equals to the division of \((-\tau-(-1))(-\tau)(-\tau-1)\) by \((-2-(-1))(-2)(-2-1)\). Thus,
\[C_{-2}=\frac{\tau^{3}}{6}-\frac{\tau}{6}\]
Extend the same formula to other coefficients, we have
\[C_{-1}=-\frac{\tau^{3}}{2}+\frac{\tau^{2}}{2}+\tau\]
\[C_{0}=\frac{\tau^{3}}{2}-\tau^{2}-\frac{\tau}{2}+1\]
\[C_{1}=-\frac{\tau^{3}}{6}+\frac{\tau^{2}}{2}-\frac{\tau}{3}\]
Rewriting the y[n] equation, we have
\[y[n]=\frac{1}{6}[(x[n-2]-3x[n-1]+3x[n]-x[n+1])\tau^{3}+ \\
(3x[n-1]-6x[n]+3x[n+1])\tau^2+ \\
(-x[n-2]+6x[n-1]-3x[n]-2x[n+1])\tau+6x[n]] \]

Corresponding Farrow structure is shown below. It only needs three multipliers and \(\tau\) is input to all multipliers. This simplifies the hardware design. When polynomial is different, the coefficients of FIR filter change accordingly. We should also note that due to the difference between input and output rates, one input sample not necessarily generates one output sample. A control logic is needed to dictate when to generate outputs.

Signal Processing Magic (2) -- Sigma Delta ADC

Sigma Delta ADC is a commonly used architecture for converting analog signal to digital. The conception of ADC is quite simple: let us say we have 4 bits with bit3 representing 1volt, bit2 for 0.5 v, bit1 for 0.25v and bit0 for 0.125v, I can quantize an analog input in volt to these 4 bits. For example, if the input is 1.2 volts, then it becomes [bit3 bit2 bit1 bit0] = [1010] = 1.25v. The quantization error is equal to 1.2-1.25 = -0.05v. Quantization noise comes from quantization error. A straightforward thought is to use binary search to convert analog input to digital. This is the basic idea behind SAR (Successive Approximation Register) ADC. 

Sigma Delta ADC, our main topic in this article, is constructed differently. Assuming x(n) is the input and y(n) is the output of sigma delta modulator, y(n) has a much higher sampling rate than x(n) and y(n) is also only 1 bit. So in term of y(n), sigma delta can be seen as a trade off between sampling rate and bitwidth. y(n) will then pass through lowpass filter and decimator to produce the multiple-bit ADC output we commonly work with.

Now let us do some math to figure out what really happens in sigma delta. First we calculate the transfer function between v(n) and u(n)

\[z^{-1}(U(z)+V(z))=V(z)\]
which gives
\[H(z)=\frac{V(z)}{U(z)}=\frac{z^{-1}}{1-z^{-1}}\]

The transfer function from input, x(n), to output, y(n), is equal to
\[\frac{Y(z)}{X(z)}=\frac{H(z)}{1+H(z)}=z^{-1}\]
which tell us there is signal delay but no signal distortion.

e(n) added to 1-bit ADC in the diagram represents quantization noise. The transfer function from e(n) to y(n) is
\[\frac{Y(z)}{E(z)}=\frac{1}{1+H(z)}=1-z^{-1}\]
Frequency response of Y(z)/E(z) is shown below. Its amplitude increases with frequency which means sigma delta modulator pushes out noise from in-band to out-of-band. Reducing in-band noise level is the major benefit of sigma delta. What we discussed here is first-order sigma delta. If we goes to higher order sigma delta, the in-band noise level will be further reduced but the cost paid is more sophisticated hardware and higher out-of-band noise level.

Signal Processing Magic (1) -- CIC Filter

In the world of signal processing, there are many elegant schemes. They are beautiful and pragmatic at the same time. Sometime you are wondering how the inventors create them. I feel obliged to introduce them to my audience. Today, let me start from a filter named CIC (Cascaded Integrator–Comb). CIC filter is mainly used for integer up-sampling and down-sampling. The main benefit of this scheme is that no multiplication is needed, which is welcomed in term of hardware complexity.

Figure below shows CIC filter architecture with N stages and decimation rate of R. M is an extra parameter and it can be 1 or 2. Normally M is set as 1. CIC filter can also be used as interpolation and the flow is opposite when used as decimation.

Accordingly, we can write the transfer function of decimation CIC filter as:

\[(\frac{1-z^{-RM}}{1-z^{-1}})^{N}\]


By increasing N with fixed R, the filter can have more rejection.

By increasing R with fixed N, the filter's cutoff frequency decreases proportionally.