Word2Vec Implementation in Rust: Continuous Bag of Words

I recently started to get serious about learning rust and I needed a small to medium project to train on.
One thing I explained in this blog and walked through the code is Word2Vec. A method to learn word embeddings. In this post I will explain the method again using my new rust code on github.

Word Similarity - Basic Idea 

In the following, I will describe the basic idea behind continuous bag of words, one way to estimate word vectors. Word vectors are a dense vectorised word representation capturing similarity of words in Euclidean space. Basically we have one vector for each word in our vocabulary. The euclidean distance between two vectors of words that are similar is smaller than the distance between vectors of words that are less similar. The assumption in Word2Vec is that words that appear in a similar context are similar. We shift a sliding window over the text and try to predict the center word in the window from the word vectors of the surrounding words in the window. 

pub struct Window<'a> {
    pub words:       &'a [usize], 
    pub predict_pos: usize
}

pub struct Document {
    pub words: Vec<usize>
}

impl<'a> Window<'a> {

    pub fn label(&self) -> usize {
        self.words[self.predict_pos]
    }

}

impl Document {

    pub fn window(&self, pos: usize, size: usize) -> Option<Window> {
        let start = sub(pos, size);
        let stop  = usize::min(pos + size, self.words.len() - 1); 
        if stop - start == 2 * size {
            Some(Window {words: &self.words[start as usize .. stop as usize], predict_pos: size})
        } else {
            None
        }
    }

}

The code above defines a document as a vector of word ids and allows to shift a sliding window over it. The return type is not a vector itself but a slice, which is convenient since there is nothing copied. In the following the center word in the window is called a label $l$ and the other words are called the context: $c_1 .. c_m$.

Word Vectors

Now that we have defined the window, we can start implementing a representation for the word vectors. There will be one vector of $d$ dimensions or columns per word in the dictionary.
In code we save the word vectors as a flat vector of doubles. We can then access the $i-th$ vector as a slice from $[i * cols .. (i + 1) * cols]$

#[derive(Serialize, Deserialize)]
pub struct ParameterStore {
    pub cols: usize,
    pub weights: Vec<f64>
}

impl ParameterStore {

    pub fn seeded(rows: usize, cols: usize) -> ParameterStore {
        let mut rng = rand::thread_rng();
        let mut weights = Vec::new();
        for _i in 0 .. (rows * cols) {
            let uniform = (rng.gen_range(0.0, 1.0) - 0.5) / cols as f64;
            weights.push(uniform);
        }
        ParameterStore {cols: cols, weights: weights}
    }

    pub fn at(&self, i: usize) -> &[f64] {
        let from = i * self.cols;
        let to = (i + 1) * self.cols;
        &self.weights[from .. to]
    }

    pub fn update(&mut self, i: usize, v: &Vec<f64>) {
        let from = i * self.cols;
        for i in 0 .. self.cols {
            self.weights[i + from] += v[i];
        }
    }

    pub fn avg(&self, result: &mut Vec<f64>, indices: Vec<usize>) {
        for col in 0 .. self.cols {
            result[col] = 0.0;
            for row in indices.iter() {
                let from = row * self.cols;
                result[col] += self.weights[col + from];
            }
            result[col] /= indices.len() as f64;
        }
    }

}

We initialise the vectors by a random number in the range $\frac{Uniform([-0.5;0.5])}{d}$, which is the standard way to initialise word2vec.  The update method simply adds a gradient to a word vector. There is also a method to compute the average of multiple word vectors, needed later during learning.

Putting It All Together: Computing the Gradients

So now we can compute windows from the text and average the word vectors of the context words. 

The averaged vectors are our hidden layer's activations. From there we compute the prediction and then back propagate the error. So we need two sets of parameters. One is the word vectors initialised as above, and one for the predictions which is initialised to zeros. 

#[derive(Serialize, Deserialize)]
pub struct CBOW {
    pub embed:   ParameterStore,
    pub predict: ParameterStore
}

Given the average embedding for the context as $h$ and the vector for word $i$ from the predictions $w_i$, the forward pass on the neural network is:

$a_i = \sigma(h * w_i)$.

We are predicting the center word of the window, that means we model the output layer as a multinomial where the center word is one hotcoded. So output should be $1$ for the center word
and $0$ for all others. Then the error for the center word is:

$e_i = log(a_i) $

and for all others it is

$e_i = log(1 - a_i)$.

The gradient on the loss is simply $\delta_i = (y - a_i) * \lambda$, with $\lambda$ being the learning rate.

 fn gradient(&self, label: usize, h: &Vecf&lt64>, truth: f64, rate: f64) -> (f64, f64) {        
        let w = self.predict.at(label);
        let a = sigmoid(dot(&h, &w));
        let d = (truth - a) * rate;
        let e = -f64::ln(if truth == 1.0 {a} else {1.0 - a});    
        (d, e)
    } 

The gradient applied to the output parameters is: $\delta_i * h$. The gradient for the embeddings is: $\delta_i * w_i$ which is applied to all context words. So all words in the same context are updated with the same gradient, in turn making them more similar, which was the main goal for word embeddings. There is one more detail. Normally, we would have to compute the losses for the complete dictionary, however for a million words, that would take very long. So instead we use something called negative sampling. We compute the gradient for the center word and then sample some negative words as the negative gradient.
The final computation for learning is shown below:

 pub fn negative_sampeling(&mut self, window: &Window, rate: f64, n_samples: usize, unigrams: &Sampler) -> f64 {
        let mut gradient_embed = vec![0.0; self.embed.cols];
        let h = self.embedding(self.ids(window));        
        let (pos, pos_err) = self.gradient(window.label(), &h, 1.0, rate);        
        let mut error = pos_err;
        let mut gradient_pos_predict = vec![0.0; self.predict.cols];
        for i in 0 .. self.embed.cols {        
            gradient_embed[i] += pos * self.predict.at(window.label())[i];
            gradient_pos_predict[i] += pos * h[i];
        }
        self.predict.update(window.label(), &gradient_pos_predict);
        for _sample in 0 .. n_samples {
            let token = unigrams.multinomial();
            let (neg, neg_error) = self.gradient(token, &h, 0.0, rate);
            error += neg_error;
            let mut gradient_neg_predict = vec![0.0; self.predict.cols];
            for i in 0 .. self.embed.cols {
                gradient_embed[i]       += neg * self.predict.at(token)[i];
                gradient_neg_predict[i] += neg * h[i];
            }
            self.predict.update(token, &gradient_neg_predict);
        }
        let ids = self.ids(window);
        for i in 0 .. ids.len() {
            self.embed.update(ids[i], &gradient_embed);
        }
        error
    }

Results

I trained the net on quora questions and when searching for similar vectors we will find similar words:

Searching: India
kashmir: 11.636631273392045
singapore: 12.155846288647625
kerala: 12.31677547739736
dubai: 12.597659271137639
syria: 12.630908898646624
pradesh: 12.657158343619573
pok: 12.667834993874745
malaysia: 12.708228117809442
taiwan: 12.75190112716203
afghanistan: 12.772290465531222
africa: 12.772372266369985
bjp: 12.802207718902666
karnataka: 12.817502941769622
europe: 12.825095372204384




Lindenmayer system: Implementation

I recently got interested in Lindenmayer systems or L-system for short and I also wrote some code.  These are generative 

systems that generate strings, that in turn can be interpreted as turtle graphics. For example, 

you can model koch curves and the famous dragon curve. This book describes how to generate plants which you can see below.

Implementation

For example, given a initial string $+A-Bf$, we have terminal symbols $\{+,-,f\}$ and non terminals $\{A, B\}$ that can still be replaced by rules like $A \rightarrow AfA$, which reads: each A can be replaced in the string. Furthermore, the string $+A-Bf$ can be interpreted by a turtle as $+$: turning right, $-$: turning left, all capital letters mean moving forward drawing a line, and all lower case letters mean moving silently to the next position. We start by defining a symbol class.

case class Symbol(id: Char) {

  def isSilent: Boolean = id.isLower

  def isTerminal: Boolean = !id.isLetter || id.isLower

}

Next we define a Lindenmayer system as a start string and a set of rules. Given the start string and a depth limit, we can recursively replace all non terminal symbols with the associated rules. This enables us to generate string.

case class LSystem(start: List[Symbol], rules: Map[Symbol, List[Symbol]]) {

  private def produce(symbol: Symbol): List[Symbol] =
    if(symbol.isTerminal) List(symbol)
    else rules(symbol)

  @tailrec
  private def produce(current: List[Symbol], result: List[Symbol]): List[Symbol] = current match {
    case Nil          => result
    case head :: tail => {
      produce(tail, result ::: produce(head))
    }
  }

  @tailrec
  private def produce(current: List[Symbol], stepsLeft: Int): List[Symbol] =
    if(stepsLeft == 0) current
    else produce(produce(current, Nil), stepsLeft - 1)

  def produce(depth: Int): List[Symbol] = produce(start, depth)

}

Now we are done with the definition of the grammar system. The only thing left is to define a turtle and an interpreter that converts strings from L-systems to a sequence of turtles. I define a turtle by it's posture $[x, y, \alpha]$ where $x,y$ is the turtle's position and the orientation is defined by the angle $\alpha$. For moving forward the next positions are defined as:

• $x + \delta * cos(\alpha)$
• $y + \delta * sin(\alpha)$

A rotation basically involves increasing the angle by some certain degree.

case class Turtle(x: Int, y: Int, angle: Double, sil: Boolean) {

  def move(distance: Int): Turtle =
    Turtle(
        math.round(x + distance * math.cos(math.toRadians(angle))).toInt, 
        math.round(y + distance * math.sin(math.toRadians(angle))).toInt, angle, false)

  def moveSil(distance: Int): Turtle =
    Turtle(
        math.round(x + distance * math.cos(math.toRadians(angle))).toInt, 
        math.round(y + distance * math.sin(math.toRadians(angle))).toInt, angle, true)

  def turn(by: Double): Turtle = {
    Turtle(x, y, (angle + by) % 360, false)
  }

  def silence: Turtle = Turtle(x, y, angle, true)

Then given a string, we can recursively build a sequence of positions by applying the action associated with each symbol.

class TurtleInterpreter(moveBy: Int, turnBy: Double) {

  import TurtleInterpreter._

  private def execute(cmd: List[Symbol], state: Turtle, stack: List[Turtle], drawable: List[Turtle]): List[Turtle] = cmd match {
    case Nil => drawable
    case RotateRight :: tail => execute(tail, state.turn(turnBy), stack, drawable)
    case RotateLeft  :: tail => execute(tail, state.turn(-turnBy), stack, drawable)
    case PushStack   :: tail => execute(tail, state, state :: stack, drawable)
    case PopStack    :: tail => execute(tail, stack.head.silence, stack.tail, stack.head.silence :: drawable)
    case head        :: tail => {
      val next = if (head.isSilent) state.moveSil(moveBy) else state.move(moveBy)
      execute(tail, next, stack, next :: drawable)
    }
  }

  def execute(cmd: List[Symbol], init: Turtle): List[Turtle] =
    execute(cmd, init, List.empty[Turtle], List(init))

}

object TurtleInterpreter {

  final val RotateRight = Symbol('+')
  final val RotateLeft  = Symbol('-')
  final val PushStack   = Symbol('[')
  final val PopStack    = Symbol(']')

Given that we can define grammars and paint pictures as the ones shown above. I wrote some code that plots all the above images that can be found on github.

Simple Gaussian Processes For Regression

In this post I am implementing a simple Gaussian Process Regressor. A Gaussian Process is a
probabilistic kernel model. Basically I am following the "Machine Learning: A probabilistic Perspective" book by Kevin Murphy.

The Basics

For an input $x_i$ and output $y_i = f(x_i)$, infer a distribution over functions of the data:
$p(f | X, y)$. In order to make predictions $y*$ for a new data point $x*$ given all the previous observations, we define:

$p(y*|x*, X, y) = \int p(y*|f, x*) p(f|x, y) df$

For our data $X = x_1 ... x_n$ the Gaussian process assumption is that the the distributions over functions $p(f_1 ... f_N)$ is jointly Gaussian with mean $\mu(x)$ and a covariance function $\Sigma_{ij} = k(x_i, x_j)$ and a kernel $k$.

Covariance structure

It is quite common to leaf the mean function at zero. Given a kernel, we need to compute the covariance function between all training and test instances. However, we structure the kernel matrix in the following way. 

The sub-matrix $K$ describes the affinity between all training instances and the sub-matrix $K*$ describes the affinity between training and test instances. Finally, the submatrix $K**$ describes the affinity between test instances to itself. The image above shows the final covariance structure. With that, our function $f*$ can be inferred as: $(f, f*) ~ N(0, \Sigma)$.

In order to obtain the prediction calculate: $f* = K_*  \bullet K^{-1} \bullet f$.
If we want the find the variance of the predictions we calculate: $trace(K_{**} - K_*  \bullet K^{-1} \bullet  K_*^T)$

In my opinion the method is very similar to support vector machines and to instance beased learning methods in general. The difference is that a SVM is only using a small set of training examples called support vectors, while the Gaussian process compares every example to all others in order to make a prediction.

A Naive Gaussian Process

Below I implemented a naive Gaussian process using numpy. Basically, I am building the three components of the covariance structure and then solve the prediction equations above. I call my method naive since computing the inverse has a worst case complexity of $O(n^3)$
So it is not advisable to use this method for larger datasets.

class GaussianProcess:
    
    def __init__(self, kernel):
        self.kernel = kernel
        
    def predict(self, train_x, train_y, test_x):
        (n, _) = train_x.shape
        (m, _) = test_x.shape
        covar_train      = np.zeros((n, n))
        covar_train_test = np.zeros((n, m))
        covar_test       = np.zeros((m, m))
        for i in range(0, n):
            for j in range(0, n):
                covar_train[i, j] = self.kernel.k(train_x[i], train_x[j], i == j)
        for i in range(0, m):
            for j in range(0, m):
                covar_test[i, j] = self.kernel.k(test_x[i], test_x[j], i == j)
        for i in range(0, n):
            for j in range(0, m):
                covar_train_test[i, j] = self.kernel.k(train_x[i], test_x[j], False)
        covar_inv  = np.linalg.inv(covar_train)
        prediction = np.dot(np.dot(covar_train_test.T, covar_inv), train_y)
        confidence = covar_test - np.dot(covar_train_test.T, np.dot(covar_inv, covar_train_test))
        confidence = np.sqrt(np.array([confidence[i, i] for i in range(0, m)]))
        return prediction, confidence 

Squared Exponential Kernel

For my experiments I use a squared exponential kernel and also add a variance parameter for for noisy observations:

$k(x_i, x_j) = \sigma_f^2 e^{-\frac{1}{2l^2} || x_i - x_j ||_2^2} + \sigma_y^2 \mathbb{I}(i == j)$

The parameters are as follows:

• $\sigma_f^2$: Variance giving the scale of the output
• $\sigma_y^2 \mathbb{I}(i == j)$: Noise variance applied to diagnoal
• $l$: Scale on input

again, the implementation is quite simple in numpy.

class SquaredExponential:

    def __init__(self, variance, alpha, noise):
        self.alpha = alpha
        self.variance = variance
        self.noise = noise
        self.scaler = -1.0 / (2 * alpha**2)
        
    def k(self, a, b, on_diag):
        error  = np.sum(np.square(a - b))
        f      = np.exp(self.scaler * error)
        scaled = self.variance * f
        if on_diag:
            scaled += self.noise
        return scaled   

Results on some simple example datasets

I tried my method on three datasets, 2 times, I am intepolating stock marktet data and one time a sine wave. Each time, the index is the input to the method and the output for the stock market prediction is the closing price.
Sine Wave Interpolation:

Stock Interpolation Dax:

Stock Interpolation US:

Integrating Weka's UI into your own code

Weka is a machine learning toolkit written in Java. Despite it's capability to run as a stand alone

user interface, you can also use Weka in your own Java code. That means, all algorithms you can find in the user interface are also trainable and usable in your Java application. The cool thing is, you can also integrate part of Weka's user interface into your own. So my goal was to integrate the option to choose a Weka classifier and configure it using Weka's native components.

Choosing and configuring a Weka classifier, here a Support Vector Machine.

After browsing a little through the Weka UI code I found two classes that can be used to

open the editor to choose a classifier and to configure it. Once the window closes, the classifier

with it's configuration is ready to be used by other functions. The first UI component is the GenericObjectEditor, which is responsible for choosing a classifier. The second is the PropertyDialog which is configuring the current classifier. Blow I showed some code that displays the dialog and prints out the classifier. I found the example in the Weka GenericObjectEditor class at the bottom. For the code to work you have to make sure the weka.jar is in your CLASSPATH.

package wekaui;

import java.awt.Frame;
import java.awt.event.WindowAdapter;
import java.awt.event.WindowEvent;
import java.beans.PropertyEditor;
import weka.classifiers.Classifier;
import weka.classifiers.rules.ZeroR;
import weka.gui.GenericObjectEditor;
import weka.gui.PropertyDialog;

/**
 * Choose a classifier and configure it in your own java code.
 *
 * @author Daniel Kohlsdorf
 */
public class WekaUI {

    public static void main(String[] args) {
        Object initial = new ZeroR();
        GenericObjectEditor.registerEditors();
        GenericObjectEditor editor = new GenericObjectEditor(true);
        ce.setClassType(weka.classifiers.Classifier.class);
        ce.setValue(initial);

        // Setup the property view
        PropertyDialog propertyDialog = new PropertyDialog((Frame) null, editor, 100, 100);

        // Get the classifier when the window is closed
        pd.addWindowListener(new WindowAdapter() {
            @Override
            public void windowDeactivated(WindowEvent e) {
                PropertyEditor propertyEditor = ((PropertyDialog) e.getSource()).getEditor();
                if (propertyEditor.getValue() != null) {
                    Classifier classifier = (Classifier) propertyEditor.getValue();
                    System.out.println("CLASSIFIER: " + classifier);
                }
            }
        });
        propertyDialog.setVisible(true);
    }

}