Showing posts with label algorithms. Show all posts
Showing posts with label algorithms. Show all posts

Thursday, April 30, 2020

Rolling neighbourhood histograms: Computing histograms of sliding windows quickly

In an earlier blog post, I described local binary patterns, which are vectors that describe the texture of image regions. In segmentation tasks where you need to group pixels together by texture, you need to compute the local binary pattern of the neighbourhoods of pixels around each pixel. This means that if you want to find the texture descriptor for a particular pixel 'P', first you need to extract a square of pixels centered on 'P' called the neighbourhood of 'P', compute the local binary codes of all the pixels in the neighbourhood, compute the histogram of the local binary codes, and that histogram is the texture descriptor of pixel 'P'. This is illustrated in the figure below.



Now in order to compute the texture descriptor of every pixel in the image you would need to repeat this process of extracting neighbourhoods and computing histograms for every pixel. But this is inefficient because it ends up recomputing the same values multiple times. First of all, the local binary codes do not need to be computed for each neighbourhood but can instead be computed once for the whole image and then the neighbourhoods be extracted from the local binary codes image instead of from the original image. But we can also make the histogram computations faster as well. This is because adjacent pixels have a lot of common pixels in their neighbourhoods as shown in the figure below.



We can take advantage of this commonality by computing the histogram of the new neighbourhood based on the histogram of the previous neighbourhood. Assuming the the histograms are simple frequency histograms, the new histogram is the previous histogram plus the frequencies of the values in the blue vertical column (see last row in the above figure) minus the frequencies of the values in the green vertical column. This means that instead of having to compute the frequencies of the whole neighbourhood you would only need to compute the frequencies for two columns and then add or subtract from the previous histogram. This reduces the time complexity from quadratic time to linear time with respect to the side of the neighbourhood square. Assuming that a neighbourhood is described by the (x,y) coordinate of the central pixel and the (w,h) width and height of the neighbourhood, and that histograms can be added together by vector addition, here is the relationship of the new histogram to the previous histogram:

histogram(x,y,w,h) = histogram(x-1,y,w,h) + histogram(x-w/2,y,1,h) - histogram(x+w/2,y,1,h)

This lets us compute the histogram of each neighbourhood by iteratively using the previous histogram in the same row, requiring only the first pixel in the row to be computed fully. But we can also avoid computing the full histogram of the first pixel in a row after the first row by using the histogram of the above pixel in the previous row like this:

histogram(x,y,w,h) = histogram(x,y-1,w,h) + histogram(x,y-h/2,w,1) - histogram(x,y+h/2,w,1)

This is similar to the way rolling hashes are computed, which is where the name rolling histogram comes from.

Monday, February 26, 2018

Finding the Greatest Common Divisor (GCD) / Highest Common Factor (HCF) of two numbers

The greatest common divisor or highest common factor of two numbers can be found using a very simple algorithm described since the time of Euclid, known as the Euclidean algorithm. The nice thing about it is that it can be computed visually by drawing squares in a rectangle as follows.


Let's say that you have a line segment of length 6 units and a smaller line segment of length 2 units.




We can check if the smaller line's length is a divisor of the larger line's length if we can place several of the smaller lines next to each other and reach the exact length of the larger line.




Likewise we can check if a number is a common divisor of two other numbers by checking if a square can fill exactly a rectangle.


In the above diagram, the big blue rectangle has sides of 6 units by 4 units and the small red squares have sides 2 units by 2 units. Since the red square can exactly fill the blue rectangle by placing copies of it next to each other, 2 is a common divisor of 4 and 6.

The greatest common divisor is the largest square that can exactly fill a rectangle. The largest square that can do this has its side equal to the smallest side of the rectangle (otherwise not even one square will fit in the rectangle) whilst the smallest square has its side equal to 1 unit (assuming that sides are always whole number lengths). Let's take the above rectangle as an example and use the Euclidean algorithm to find the largest square that fits a 6 by 4 rectangle.

The Euclidean algorithm takes advantage of the fact that the greatest common divisor of two numbers is also a divisor of the difference of the two numbers. If two numbers "a" and "b" have "c" as a common divisor, then they can be rewritten as "c×p" and "c×q". The difference "a-b" is then equal to "c×p - c×q" which is equal to "c×(p - q)" which is therefore a divisor of "c" since one of its factors is "c". Notice that "c" can be any divisor of "a" and "b", which means that "a-b" should also have the greatest common divisor of "a" and "b" as a divisor.

This means that the largest square to exactly fill the original 4 by 6 square will also be the largest square to exactly fill a 2 by 4 rectangle, since 2 is the difference between 4 and 6.

This piece of information allows us to substitute the difficult problem of finding the largest square that exactly fills a rectangle with an easier problem of finding the same square but on a smaller rectangle. The smaller the rectangle, the easier it is to find the largest square to exactly fill it and the Euclidean algorithm allows us to progressively chop off parts of the rectangle until all that remains is a small square which would be the square that exactly fills the original rectangle.

So the algorithm goes as follows:

  1. If the rectangle has both sides equal then it itself is the square that is the largest square that exactly fills the rectangle and we are done.
  2. Otherwise put a square inside the rectangle of length equal to the shortest side of the rectangle.
  3. Chop off the part of the rectangle that is covered by the square and repeat.

The 6 by 4 rectangle is trivial so let's try this algorithm on a slightly harder problem of 217 by 124. This time we'll be using exact pixels to represent the rectangle lengths.

This is a 217 by 124 rectangle.



This is the overlap it has with a square of length equal to its shortest side.




We now chop off that part of the rectangle and find the largest square which fits the remaining 93 by 124 rectangle.




This is the new square's overlap.




This is the remaining 93 by 31 rectangle.




This is the new square's overlap.




This is the remaining 62 by 31 rectangle.




This is the new square's overlap.




This is the remaining 31 by 31 rectangle.




We have finally reached a square of side 31, which means that 31 is the greatest common divisor of 217 and 124.


Of course we don't really need to visualise anything. We can do this all numerically by just repeatedly subtracting the smallest number from the biggest and replacing the biggest with the difference. Here's a Python function that does that:

def gcd(a, b):
    while a != b:
        if a < b:
            b = b - a
        else:
            a = a - b
    return a

We can make this shorter by using pattern matching:

def gcd(a, b):
    while a != b:
        (a, b) = (min(a,b), max(a,b)-min(a,b))
    return a

If one of the numbers is much bigger than the other, this program will waste a lot of time repeatedly subtracting the same number until the bigger number becomes smaller than the other. Fortunately this process can be done in a single operation by dividing and getting the remainder. The remainder is what's left after you can't subtract the second number from the first any more. This speeds the algorithm up considerably. The operation is called the modulo operation. With modulo the subtractions will not stop when both numbers are equal as modulo assumes that you can do one more subtraction and make the first number zero. Instead we will now keep on looping until the remainder is zero, in which case the last number you used to get the remainder was the greatest common divisor:

def gcd(a, b):
    while b > 0:
        (a, b) = (min(a,b), max(a,b)%min(a,b))
    return a

Monday, April 3, 2017

Getting the top n most probable sentences using beam search

This is a continuation from a previous blog post on single sentence beam search.

Sometimes it is not enough to just generate the most probable sentence using a language model. Sometimes you want to generate the top 3 most probable sentences instead. In that case we need to modify our beam search a bit. We will make the function a generator that returns a sequence of sentences in order of probability instead of just returning a single most probable sentence. Here are the changes we need to make:

In the single sentence version, we were getting the most probable prefix in the current beam and checking if it is complete. If it is, then we return it and stop there. Instead, we will now not stop until the current beam is empty (or until the caller stops requesting for more sentences). After returning the most probable prefix we will check the second most probable prefix and keep on returning complete prefixes until we either find one which is not complete or we return all the beam. In the case that we return the whole beam then the algorithm stops there as there is nothing left with which to generate new prefixes. This means that the beam width gives a limit on the number of sentences that can be returned. If we do not return all the beam then we continue generating prefixes with the remainder.

In the case that some complete sentences were returned, they need to also be removed from the beam before we continue generating. Since the beam is implemented as a min-first heap queue (min-first because we want to pop the least probable prefix quickly when the beam becomes bigger than the beam width) then we cannot remove the highest probable complete sentence quickly as well. In order to do this, we first turn the heap into a list which is sorted by probability and then start popping out the items at the end if they are complete sentences. Following this, we will then heapify the list back into a min-first heap queue and continue as normal. This sorting and reheapifying should not impact on the performance too much if the beam width is relatively small.

If the clip length is reached then the whole beam is immediately returned in order of probability. This is because as soon as one prefix is equal to the allowed maximum then that means that the entire beam consists of
  1. incomplete sentences that are also as long as the allowed maximum (since all the prefixes grow together)
  2. complete sentences that were found before but which do not have a maximum probability
After returning all the incomplete sentences (they have to be returned since they cannot be extended further) then the complete sentences will also be returned as they will become the most probable sentences of what is left in the beam. The assumption is that if the complete sentence was a nonsensical one, then it wouldn't have remained in the beam so might as well return it as well rather than lose it.

Here is the modified Python 3 code:

import heapq

class Beam(object):

    def __init__(self, beam_width, init_beam=None):
        if init_beam is None:
            self.heap = list()
        else:
            self.heap = init_beam
            heapq.heapify(self.heap) #make the list a heap
        self.beam_width = beam_width

    def add(self, prob, complete, prefix):
        heapq.heappush(self.heap, (prob, complete, prefix))
        if len(self.heap) > self.beam_width:
            heapq.heappop(self.heap)

    def __iter__(self):
        return iter(self.heap)


def beamsearch(probabilities_function, beam_width=10, clip_len=-1):
    prev_beam = Beam(beam_width)
    prev_beam.add(1.0, False, [ '<start>' ])

    while True:
        curr_beam = Beam(beam_width)

        #Add complete sentences that do not yet have the best
probability to the current beam, the rest prepare to add more words to
them.
        for (prefix_prob, complete, prefix) in prev_beam:
            if complete == True:
                curr_beam.add(prefix_prob, True, prefix)
            else:
                #Get probability of each possible next word for the
incomplete prefix.
                for (next_prob, next_word) in probabilities_function(prefix):
                    if next_word == '<end>': #if next word is the end
token then mark prefix as complete and leave out the end token
                        curr_beam.add(prefix_prob*next_prob, True, prefix)
                    else: #if next word is a non-end token then mark
prefix as incomplete
                        curr_beam.add(prefix_prob*next_prob, False,
prefix+[next_word])

        sorted_beam = sorted(curr_beam) #get all prefixes in current beam sorted by probability
        any_removals = False
        while True:
            (best_prob, best_complete, best_prefix) = sorted_beam[-1] #get highest probability prefix
            if best_complete == True or len(best_prefix)-1 ==
clip_len: #if most probable prefix is a complete sentence or has a length that exceeds the clip length (ignoring the start token) then yield it
                yield (best_prefix[1:], best_prob) #yield best sentence without the start token and together with its probability
                sorted_beam.pop() #remove the yielded sentence and check the next highest probability prefix
                any_removals = True
                if len(sorted_beam) == 0: #if there are no more sentences in the beam then stop checking
                    break
            else:
                break

        if any_removals == True: #if there were any changes in the current beam then...
            if len(sorted_beam) == 0: #if there are no more prefixes in the current beam (due to clip length being reached) then end the beam search
                break
            else: #otherwise set the previous beam to the modified current beam
                prev_beam = Beam(beam_width, sorted_beam)
        else: #if the current beam was left unmodified then assign it to the previous beam as is
            prev_beam = curr_beam

Sunday, October 30, 2016

Using beam search to generate the most probable sentence

This blog post continues in a second blog post about how to generate the top n most probable sentences.

In my last blog post I talked about how to generate random text using a language model that gives the probability of a particular word following a prefix of a sentence. For example, given the prefix "The dog", a language model might tell you that "barked" has a 5% chance of being the next word whilst "meowed" has a 0.3%. It's one thing generating random text in a way that is guided by the probabilities of the words but it is an entirely different thing to generate the most probable text according to the probabilities. By most probable text we mean that if you multiply the probabilities of all the words in the sentence together you get the maximum product. This is useful for conditioned language models which give you different probabilities depending on some extra input, such as an image description generator which accepts an image apart from the prefix and returns probabilities for the next word depending on what's in the image. In this case we'd like to find the most probable description for a particular image.

You might think that the solution is to always pick the most probable word and add it to the prefix. This is called a greedy search and is known to not give the optimal solution. The reason is because it might be the case that every combination of words following the best first word might not be as good as those following the second best word. We need to use a more exploratory search than greedy search. We can do this by thinking of the problem as searching a probability tree like this:



The tree shows a probability tree of which words can follow a sequence of words together with their probabilities. To find the probability of a sentence you multiply every probability in the sentence's path from <start> to <end>. For example, the sentence "the dog barked" has a probability of 75% × 73% × 25% × 100% = 13.7%. The problem we want to solve is how to find the sentence that has the highest probability.

One way to do this is to use a breadth first search. Starting from the <start> node, go through every node connected to it, then to every node connected to those nodes and so on. Each node represents a prefix of a sentence. For each prefix compute its probability, which is the product of all the probabilities on its path from the <start> node. As soon as the most probable prefix found is a complete sentence, that would be the most probable sentence. The reason why no other less probable prefixes could ever result in more probable sentences is because as a prefix grows, its probability shrinks. This is because any additional multiplications with probabilities made to any prefix probability will only make it smaller. For example, if a prefix has a probability of 20% and another word is added to it which has a probability of 99%, then the new probability will be 20% × 99% which is the smaller probability of 19.8%.

Of course a breadth first search is impractical on any language model that has a realistic vocabulary and sentence length since it would take too long to check all the prefixes in the tree. We can instead opt to take a more approximate approach where we only check a subset of the prefixes. The subset would be the top 10 most probable prefixes found at that point. We do a breadth first search as explained before but this time only the top 10 most probable prefixes are kept and we stop when the most probable prefix in these top 10 prefixes is a complete sentence.

This is practical but it's important that the way we find the top 10 prefixes is fast. We can't sort all the prefixes and choose the first 10 as there would be too many. We can instead use a heap data structure. This data structure is designed to quickly take in a bunch of numbers and quickly pop out the smallest number. With this you can insert the prefix probabilities one by one until there are 10 prefixes in it. After that start popping out the smallest probability immediately after inserting a new one in order to only keep the 10 largest ones.

Python provides a library called "heapq" (heap queue) just for this situation. Here is a class that makes use of heapq in order to create a beam of prefix probabilities:

import heapq

class Beam(object):
#For comparison of prefixes, the tuple (prefix_probability, complete_sentence) is used.
#This is so that if two prefixes have equal probabilities then a complete sentence is preferred over an incomplete one since (0.5, False) < (0.5, True)

    def __init__(self, beam_width):
        self.heap = list()
        self.beam_width = beam_width

    def add(self, prob, complete, prefix):
        heapq.heappush(self.heap, (prob, complete, prefix))
        if len(self.heap) > self.beam_width:
            heapq.heappop(self.heap)
    
    def __iter__(self):
        return iter(self.heap)

The code to perform the actual beam search is this:

def beamsearch(probabilities_function, beam_width=10, clip_len=-1):
    prev_beam = Beam(beam_width)
    prev_beam.add(1.0, False, [ '<start>' ])
    while True:
        curr_beam = Beam(beam_width)
        
        #Add complete sentences that do not yet have the best probability to the current beam, the rest prepare to add more words to them.
        for (prefix_prob, complete, prefix) in prev_beam:
            if complete == True:
                curr_beam.add(prefix_prob, True, prefix)
            else:
                #Get probability of each possible next word for the incomplete prefix.
                for (next_prob, next_word) in probabilities_function(prefix):
                    if next_word == '<end>': #if next word is the end token then mark prefix as complete and leave out the end token
                        curr_beam.add(prefix_prob*next_prob, True, prefix)
                    else: #if next word is a non-end token then mark prefix as incomplete
                        curr_beam.add(prefix_prob*next_prob, False, prefix+[next_word])
        
        (best_prob, best_complete, best_prefix) = max(curr_beam)
        if best_complete == True or len(best_prefix)-1 == clip_len: #if most probable prefix is a complete sentence or has a length that exceeds the clip length (ignoring the start token) then return it
            return (best_prefix[1:], best_prob) #return best sentence without the start token and together with its probability
            
        prev_beam = curr_beam

"probabilities_function" returns a list of word/probability pairs given a prefix. "beam_width" is the number of prefixes to keep (so that instead of keeping the top 10 prefixes you can keep the top 100 for example). By making the beam search bigger you can get closer to the actual most probable sentence but it would also take longer to process. "clip_len" is a maximum length to tolerate, beyond which the most probable prefix is returned as an incomplete sentence. Without a maximum length, a faulty probabilities function which does not return a highly probable end token will lead to an infinite loop or excessively long garbage sentences.

Tuesday, June 21, 2016

The Backpropagation Algorithm for Artificial Neural Networks (ANNs)

In a previous post, I talked about how to use the gradient descent algorithm to optimize the weights and biases of an artificial neural network in order to give selected outputs for selected inputs. However plain gradient descent is inefficient on large neural nets since it recomputes a lot of values for each neuron, plus it has to be rewritten for every change in architecture (number of neurons and layers) and requires a lot of code. The standard solution is to use an optimized version of the gradient descent algorithm for neural nets called the backpropagation algorithm.

I will assume that you have read the previous post. Whereas the previous post was very specific using a specified architecture and a specified activation and cost function, here I will keep things as general as possible such that they can be applied on any feed forward neural network. Here are the definitions of symbols we shall be using, similar to last post's definitions:



In order to help explain the algorithm and equations, I shall be applying it to the last post's architecture, just to help you understand. So here is last post's neural net:



In the example, we have a 3 layer neural net with 2 neurons in each layer and 2 input neurons. It uses the sigmoid function as an activation function and the mean square error as the cost function.

The main point of interest in the backpropagation algorithm is the way you find the derivative of the cost function with respect to a weight or bias in any layer. Let's look at how to find these derivatives for each layer in the example neural net.

Weights and bias in layer L (output layer)

We begin by finding the derivatives for the output layer, which are the simplest.

d Cost / d Weight L
The derivative of the cost with respect to any weight in layer L can be found using



Here's how we can arrive to this equation based on a weight in layer 3 in the example:



d Cost / d Bias L

The derivative of the cost with respect to any bias in layer L can be found using



Here's how we can arrive to this equation based on a bias in layer 3 in the example:



Using delta

A part of the weights and biases equations is repeated. It will be convenient when finding derivatives of other layers to use an intermediate letter, called lower case delta, to represent this part of these equations.



Weights and bias in layer L-1

Now we find the derivative with respect to the weights and biases in layer L-1. The derivations will be continuations of the previous derivations and we won't be repeating the first bit of the derivations again.

d Cost / d Weight L-1
The derivative of the cost with respect to any weight in layer L-1 can be found using



Here's how we can arrive to this equation based on a weight in layer 2 in the example:



d Cost / d Bias L-1

The derivative of the cost with respect to any bias in layer L-1 can be found using



Here's how we can arrive to this equation based on a bias in layer 2 in the example:



Using delta

These equations can be shortened by using delta again and now we can use the previous delta to shorten this one even more.



Notice how there is a recursive definition of delta. This is the basis for the backpropagation algorithm.

Weights and bias in layer L-2

Now we find the derivative with respect to the weights and biases in layer L-2. Like before, we won't be repeating the first bit of the derivations again.

d Cost / d Weight L-2
The derivative of the cost with respect to any weight in layer L-2 can be found using



Here's how we can arrive to this equation based on a weight in layer 1 in the example (notice that there is a trick we'll be using with the summations that is explained below):



At the end of the derivation we did a trick with the sigmas (summations) in order to turn the equation into a recursive one.

First, we moved the delta to the inside of the sigmas. This can be done because it's just a common term that is factored out and we're factoring it back in. For example if the equation is "𝛿(a + b + c)" then we can turn that into "𝛿a + 𝛿b + 𝛿c".

Second, we swapped the "p" summation with the "q" summation. This does not change anything if you think about it. All it does is change the order of the summations, which does not change anything.

Finally we replaced the "p" summation with the previous delta, allowing us to shorten the equation.

d Cost / d Bias L-2

The derivative of the cost with respect to any bias in layer L-2 can be found using



Here's how we can arrive to this equation based on a bias in layer 1 in the example:



Using delta

Again, we can shorten the equations using delta.



Notice that this is exactly the same as the previous delta. This is because beyond L-1, all the deltas will be identical and can be defined using a single recursive relation.

In general: using indices

Up to now we saw what the derivatives for layer L, L-1, and L-2 are. We can now generalize these equations for any layer. Given the recursive pattern in the deltas, we can create an equation that gives the deltas of any layer provided you have the deltas of the previous layers. The base case of the recursion would be for layer L, which is defined on its own.



Now we have a way to find the derivatives for any layer, no matter how deep the neural network is. Here is how you'd use these equations on the example:



In general: using matrices

As things stand, there are a lot of index letters littering our equations (i,j,p). We can get rid of them however if we use matrix operations that work on all the indices at once. If we do this we can even get shorter code when programming the backpropagation algorithm by making use of a fast matrix library such as numpy.

Using matrices in neural nets

Let's look at how we can use matrices to compute the output of a neural net. A whole layer of activations can be calculated as a whole by treating it as a vector. We'll treat vectors as being horizontal by default and which need to be transposed in order to be made vertical.

Here's an example of what we want to calculate:


Of course we're still using indices just because we're organising our activations in a vector. In order to get rid of the indices we need to treat the weights as a matrix, where the first index is the row number and the second index is the column number. That way we can have a single letter "w" which represents all of the weights of a layer. When multiplied by a vector of the previous layer's activations, the matrix multiplication will result in another vector of weighted sums which can be added to a bias vector and passed through an activation function to compute the next vector of activations. Here is an example of the calculation based on the above example:



Final clean generalized deltas and cost gradients

And now we can finally see what the clean matrix-based general derivatives are for any layer:



The backpropagation algorithm

Finding the gradients is one step (the most complex one) of the backpropagation algorithm. The full backpropagation algorithm goes as follows:

  1. For each input-target pair in the training set,
    1. Compute the activations and the "z" of each layer when passing the input through the network. This is called a forward pass.
    2. Use the values collected from the forward pass to calculate the gradient of the cost with respect to each weight and bias, starting from the output layer and using the delta equations to compute the gradients for previous layers. This is called the backward pass.
  2. Following the backward pass, you have the gradient of the cost with respect to each weight and bias for each input in the training set. Add up all the corresponding gradients of each input in the training set (all the gradients with respect to the same weight or bias).
  3. Now use the gradient to perform gradient descent by multiplying the gradient by a step size, or learning rate as it's called in the backpropagation algorithm, and subtracting it from the corresponding weights and biases. In algebra,
    weight = weight - learningrate*dCost/dWeight

In code
And this is what the full program looks like in Python 3 using numpy to perform matrix operations. Again, we shall be training the neural net to perform the function of a half adder. A half adder adds together two binary digits and returns the sum and carry. So 0 + 1 in binary gives 1 carry 0 whilst 1 + 1 in binary gives 0 carry 1.

import math, random
import numpy as np

class ANN:
    '''General artificial neural network class.'''
    
    def __init__(self, num_in, num_hiddens, num_out, trainingset):
        '''Create a neural network with a set number of layers and neurons and with initialised weights and biases.'''
        if not all(len(x)==num_in and len(y)==num_out for (x,y) in trainingset.items()):
            raise ValueError()

        self._num_in = num_in
        self._num_hiddens = num_hiddens
        self._num_out = num_out
        self._trainingset = trainingset

        self._num_layers = len(num_hiddens) + 1
        
        self._weights = dict([
                (1, np.array([
                    [self.weight_init(1, i, j) for j in range(num_hiddens[0])]
                    for i in range(num_in)
                ]))
            ]
            + [
                (l+1, np.array([
                    [self.weight_init(l+1, i, j) for j in range(num_hiddens[l])]
                    for i in range(num_hiddens[l-1])
                ]))
                for l in range(1, len(num_hiddens))
            ]
            + [
                (self._num_layers, np.array([
                    [self.weight_init(self._num_layers, i, j) for j in range(num_out)]
                    for i in range(num_hiddens[-1])
                ]))
            ])
        
        self._biases = dict([
                (l+1, np.array([
                    self.bias_init(l+1, j) for j in range(num_hiddens[l])
                ]))
                for l in range(len(num_hiddens))
            ]
            + [
                (self._num_layers, np.array([
                    self.bias_init(self._num_layers, j) for j in range(num_out)
                ]))
            ])

    def weight_init(self, layer, i, j):
        '''How to initialise weights.'''
        return 2*random.random()-1

    def bias_init(self, layer, j):
        '''How to initialise biases.'''
        return 2*random.random()-1
        
    def a(self, z):
        '''The activation function.'''
        return 1/(1 + np.vectorize(np.exp)(-z))

    def da_dz(self, z):
        '''The derivative of the activation function.'''
        return self.a(z)*(1-self.a(z))

    def out(self, x):
        '''Compute the output of the neural network given an input vector.'''
        n = x
        for l in range(1, self._num_layers+1):
            n = self.a(np.dot(n, self._weights[l]) + self._biases[l])
        return n

    def cost(self):
        '''The cost function based on the training set.'''
        return sum(sum((self.out(x) - t)**2) for (x,t) in self._trainingset.items())/(2*len(self._trainingset))

    def dCxt_dnL(self, x, t):
        '''The derivative of the cost function for a particular input-target pair with respect to the output of the network.'''
        return self.out(x) - t

    def get_cost_gradients(self):
        '''Calculate and return the gradient of the cost function with respect to each weight and bias in the network.'''
        dC_dws = dict()
        dC_dbs = dict()
        for l in range(1, self._num_layers+1):
            dC_dws[l] = np.zeros_like(self._weights[l])
            dC_dbs[l] = np.zeros_like(self._biases[l])
        
        for (x,t) in self._trainingset.items():
            #forward pass
            zs = dict()
            ns = dict()
            ns_T = dict()
            ns[0] = np.array(x)
            ns_T[0] = np.array([np.array(x)]).T
            for l in range(1, self._num_layers+1):
                z = np.dot(ns[l-1], self._weights[l]) + self._biases[l]
                zs[l] = z
                n = self.a(z)
                ns[l] = n
                ns_T[l] = np.array([n]).T

            #backward pass
            d = self.dCxt_dnL(x, t)*self.da_dz(zs[self._num_layers])
            dC_dws[self._num_layers] += np.dot(ns_T[self._num_layers-1], np.array([d]))
            dC_dbs[self._num_layers] += d
            for l in range(self._num_layers-1, 1-1, -1):
                d = np.dot(d, self._weights[l+1].T)*self.da_dz(zs[l])
                dC_dws[l] += np.dot(ns_T[l-1], np.array([d]))
                dC_dbs[l] += d

        for l in range(1, self._num_layers+1):
            dC_dws[l] /= len(self._trainingset)
            dC_dbs[l] /= len(self._trainingset)

        return (dC_dws, dC_dbs)

    def epoch_train(self, learning_rate):
        '''Train the neural network for one epoch.'''
        (dC_dws, dC_dbs) = self.get_cost_gradients()
        for l in range(1, self._num_layers+1):
            self._weights[l] -= learning_rate*dC_dws[l]
            self._biases[l] -= learning_rate*dC_dbs[l]
    
    def check_gradient(self, epsilon=1e-5):
        '''Check if the gradients are being calculated correctly. This is done according to http://ufldl.stanford.edu/tutorial/supervised/DebuggingGradientChecking/ . This method calculates the difference between each calculated gradient (according to get_cost_gradients()) and an estimated gradient using a small number epsilon. If the gradients are calculated correctly then the returned numbers should all be very small (smaller than 1e-10.'''
        (predicted_dC_dws, predicted_dC_dbs) = self.get_cost_gradients()

        approx_dC_dws = dict()
        approx_dC_dbs = dict()
        for l in range(1, self._num_layers+1):
            approx_dC_dws[l] = np.zeros_like(self._weights[l])
            approx_dC_dbs[l] = np.zeros_like(self._biases[l])
            
        for l in range(1, self._num_layers+1):
            (rows, cols) = self._weights[l].shape
            for r in range(rows):
                for c in range(cols):
                    orig = self._weights[l][r][c]
                    self._weights[l][r][c] = orig + epsilon
                    cost_plus = self.cost()
                    self._weights[l][r][c] = orig - epsilon
                    cost_minus = self.cost()
                    self._weights[l][r][c] = orig
                    approx_dC_dws[l][r][c] = (cost_plus - cost_minus)/(2*epsilon)
                    
            (cols,) = self._biases[l].shape
            for c in range(cols):
                orig = self._biases[l][c]
                self._biases[l][c] = orig + epsilon
                cost_plus = self.cost()
                self._biases[l][c] = orig - epsilon
                cost_minus = self.cost()
                self._biases[l][c] = orig
                approx_dC_dbs[l][c] = (cost_plus - cost_minus)/(2*epsilon)

        errors_w = dict()
        errors_b = dict()
        for l in range(1, self._num_layers+1):
            errors_w[l] = abs(predicted_dC_dws[l] - approx_dC_dws[l])
            errors_b[l] = abs(predicted_dC_dbs[l] - approx_dC_dbs[l])
        return (errors_w, errors_b)

################################################################
nn = ANN(
    2, #number of inputs
    [ 2, 2 ], #number of hidden neurons (each number is a hidden layer)
    2, #number of outputs
    { #training set
        (0,0):(0,0),
        (0,1):(1,0),
        (1,0):(1,0),
        (1,1):(0,1)
    }
)

print('Initial cost:', nn.cost())
print('Starting training')
for epoch in range(1, 20000+1):
    nn.epoch_train(0.5)
    if epoch%1000 == 0:
        print(' Epoch:', epoch, ', Current cost:', nn.cost())
print('Training finished')

print('Neural net output:')
for (x,t) in sorted(nn._trainingset.items()):
    print(' ', x, ':', nn.out(x))

After 20,000 iterations I got this output:
Initial cost: 0.232581149941
Starting training
 Epoch: 1000 , Current cost: 0.218450291043
 Epoch: 2000 , Current cost: 0.215777328824
 Epoch: 3000 , Current cost: 0.178593050179
 Epoch: 4000 , Current cost: 0.102704613785
 Epoch: 5000 , Current cost: 0.0268721186425
 Epoch: 6000 , Current cost: 0.00819827730488
 Epoch: 7000 , Current cost: 0.00444660869981
 Epoch: 8000 , Current cost: 0.00298786209459
 Epoch: 9000 , Current cost: 0.00223137023455
 Epoch: 10000 , Current cost: 0.00177306322414
 Epoch: 11000 , Current cost: 0.00146724847349
 Epoch: 12000 , Current cost: 0.00124934223594
 Epoch: 13000 , Current cost: 0.00108652952015
 Epoch: 14000 , Current cost: 0.000960438140866
 Epoch: 15000 , Current cost: 0.000860007442299
 Epoch: 16000 , Current cost: 0.000778192177283
 Epoch: 17000 , Current cost: 0.000710297773182
 Epoch: 18000 , Current cost: 0.000653078479503
 Epoch: 19000 , Current cost: 0.000604220195232
 Epoch: 20000 , Current cost: 0.0005620295804
Training finished
Neural net output:
  (0, 0) : [ 0.02792593  0.00058427]
  (0, 1) : [ 0.97284142  0.01665468]
  (1, 0) : [ 0.97352071  0.01661805]
  (1, 1) : [ 0.03061993  0.97196112]

Monday, March 14, 2016

Displaying the digits of pi in your web browser

Happy pi day. Here's how you can display as many digits of pi as you want (subject to your computer's resources).

We're going to be using Gibbon's unbounded spigot algorithm which spits out consecutive digits of pi in decimal without iterative approximations. It's really cool when you see it in action.

In the above linked paper, the algorithm is given in Haskell. I won't explain the mathematical derivation in this post, but here is a translation into Python 3:
def gibbons():
    (q, r, t, k, n, l) = (1, 0, 1, 1, 3, 3)
    while True:
        if 4*q + r - t < n*t:
            yield n
            (q, r, n) = (
                10*q,
                10*(r - n*t),
                10*(3*q + r)//t - 10*n
            )
        else:
            (q, r, t, k, n, l) = (
                q*k,
                (2*q + r)*l,
                t*l,
                k+1,
                (q*(7*k + 2) + r*l)//(t*l),
                l + 2
            )
To use this code run this:
digit_seq = gibbons()
print(next(digit_seq))
print(next(digit_seq))
print(next(digit_seq))
print(next(digit_seq))
This will output
3
1
4
1
If you want it to print out a thousand digits, do this:
digit_seq = gibbons()
for (position, digit) in enumerate(digit_seq):
    if position == 1000:
        break
    print(digit, end='')
Now this code will let you print as many digits as you want; however being in Python it's not something you can just tell your non-geek friend to try out. So for pi day, here's a javascript implementation that anyone can just paste in their browser's web console and get a nice stream of pi digits. One thing you should keep in mind is that in Python you can have numbers which are as big as you like. Javascript is not like that, so we'll need to use an external library called bignumber.js. This is the full code we'll be using to generate the digits:
var script = document.createElement('script');
script.src = "https://cdnjs.cloudflare.com/ajax/libs/bignumber.js/2.3.0/bignumber.min.js";
document.getElementsByTagName('head')[0].appendChild(script);

function gibbons() {
    var q=new BigNumber(1), r=new BigNumber(0), t=new BigNumber(1), k=new BigNumber(1), n=new BigNumber(3), l=new BigNumber(3);
    function f() {
        while (true) {
            if (q.mul(4).add(r).sub(t).lt(n.mul(t))) {
                var digit = n.toString();
                var q_ = q.mul(10);
                var r_ = r.sub(n.mul(t)).mul(10);
                var n_ = q.mul(3).add(r).mul(10).divToInt(t).sub(n.mul(10));
                q=q_, r=r_, n=n_;
                return digit;
            }
            else {
                var q_ = q.mul(k);
                var r_ = q.mul(2).add(r).mul(l);
                var t_ = t.mul(l);
                var k_ = k.add(1);
                var n_ = (q.mul(k.mul(7).add(2)).add(r.mul(l))).divToInt(t.mul(l));
                var l_ = l.add(2);
                q=q_, r=r_, t=t_, k=k_, n=n_, l=l_;
            }
        }
    }
    return f;
}

function show_digits(max_num_digits) {
    var get_digit = gibbons();
    var num_digits = 0;
    function f() {
        var digits = '';
        for (var i = 1; i <= 50; i++) {
            digits += get_digit();
            num_digits++;
            if (num_digits >= max_num_digits) {
                break;
            }
        }
        console.log(digits);
        if (num_digits < max_num_digits) {
            setTimeout(f, 100);
        }
        else {
            console.log('\n'+max_num_digits+' digits of pi displayed!');
        }
    }
    f();
}

show_digits(1000);
The first line is to import the bignumber library in a way that works in the web console. After that there are two functions. The first is for the Gibbons algorithm which lets us get the digits of pi. The way it works is by returning a function which returns the next digit each time it is called. The second function outputs 50 digits of pi every 100 milliseconds.
The code to copy-paste
But this is not something you'll send you friend. Instead send them this minified code (of course you shouldn't trust obscure code sent by strangers as it might be evil code, but you can minify the above code youself if you don't trust me).
var script=document.createElement("script");script.src="https://cdnjs.cloudflare.com/ajax/libs/bignumber.js/2.3.0/bignumber.min.js",document.getElementsByTagName("head")[0].appendChild(script);

function gibbons(){function u(){for(;;){if(n.mul(4).add(i).sub(m).lt(d.mul(m))){var u=d.toString(),e=n.mul(10),r=i.sub(d.mul(m)).mul(10),t=n.mul(3).add(i).mul(10).divToInt(m).sub(d.mul(10));return n=e,i=r,d=t,u}var e=n.mul(l),r=n.mul(2).add(i).mul(o),a=m.mul(o),b=l.add(1),t=n.mul(l.mul(7).add(2)).add(i.mul(o)).divToInt(m.mul(o)),g=o.add(2);n=e,i=r,m=a,l=b,d=t,o=g}}var n=new BigNumber(1),i=new BigNumber(0),m=new BigNumber(1),l=new BigNumber(1),d=new BigNumber(3),o=new BigNumber(3);return u}function show_digits(u){function n(){for(var l="",d=1;50>=d&&(l+=i(),m++,!(m>=u));d++);console.log(l),u>m?setTimeout(n,100):console.log("\n"+u+" digits of pi displayed!")}var i=gibbons(),m=0;n()}

show_digits(1000);

What you have to do is open your browser's web console on any web page by pressing F12 (even a blank tab works, although Facebook will block it), then paste the first line of code at the bottom. After pressing enter, paste the second line of code. Finally paste the "show_digits(1000);" line and press enter. You'll see a stream of digits displayed in the console panel. I think this is a nice way to celebrate pi day. Maybe someone else can make the digits somehow interact with the web page you are on instead of just display them in the console. Until then, enjoy your pie today!

Saturday, January 23, 2016

Gradient Descent

Algebraic method: 1 variable

Say you have an equation like
y = 2x^2 - 3x + 2
that you needed to find the minimum of:


The quadratic equation is minimum (the turning point at the bottom) somewhere between x = 0.5 and x = 1. How can we find where exactly it is minimum? The standard solution is by algebraically finding where its gradient is 0. The gradient is 0 only at a turning point, which is where the minimum is. You can find the gradient at any point by taking the equation's derivative. The derivative with respect to "x" of y = 2x^2 - 3x + 2 is
dy/dx = 4x - 3


The derivative equation is 0 exactly at the same point where the first equation has a minimum turning point. So we algebraically find at which "x" "4x - 3" is 0:
0 = 4x - 3
x = 3/4 = 0.75
Therefore "x" has to be 0.75.

Algebraic method: 2 variables

What if we have 2 variables instead of just "x"? This is what the equation
z = 2x^2 + 2y^2 + 2xy - 6x
looks like:


It clearly has a minimum point somewhere close to (x,y) = (0,0). The way to find the minimum is to take it's partial derivative with respect to "x" and with respect to "y". A partial derivative is when you consider only one of the variables as a variable whilst considering the other variables as constant. This allows you to find how the gradient of the curve changes with respect to one variable. Here is the partial derivative with respect to "x":
∂z/∂x = 4x + 2y - 6
and here is the partial derivative with respect to "y":
∂z/∂y = 4y + 2x

The idea is now to find where the gradient with respect to both variables is 0. The minimum has a gradient of zero from all directions ("x" and "y") so we find where this is so:
0 = 4x + 2y - 6 [A]
0 = 4y + 2x [B]

From [B]
x = -2y [C]

Substituting [C] in [A]
0 = 4(-2y) + 2y - 6
y = -1

Therefore from [C]
x = 2

Therefore both derivatives are equal to 0 when (x,y) = (2,-1) which means that that is where the equation is minimum.

Gradient Descent

Notice how finding the the minimum of a 2 variable equation required the use of solving simultaneous equations. Imagine having a 100 variables and solving 100 simultaneous equations. Plus there is no guarantee that the equations are solvable algebraically. So sometimes we use approximate methods to find the minimum. In this case we can use a method called Gradient Descent. You start by finding the partial derivatives with respect to each variable, then you pick a point on the equation at random and follow the gradient downwards. Let's see an example with the previous equation.

∂z/∂x = 4x + 2y - 6
∂z/∂y = 4y + 2x

Let's pick (x,y) = (0,0) as a point.

∂z/∂x @ (0,0) = 4(0) + 2(0) - 6 = -6
∂z/∂y @ (0,0) = 4(0) + 2(0) = 0

The gradients suggest that at (0,0) the curve is flat at the "y" axis and sliding downwards in the positive direction of the "x" axis. Here is what that means.

If we fix "y" to be 0, "z" would vary with respect to "x" like this:

Notice how at "x" = 0 the slope steeps down towards the positive direction, which means that "x" = 0 is too small and needs to be bigger for us to find the minimum. The gradient at "x" = 0 is -6, which is negative therefore the minimum is at a greater "x" than 0.

If we fix "x" to be 0, "z" would vary with respect to "y" like this:

Notice how at "y" = 0 the slope is flat, which means that "y" is not important to find the minimum at this point and can be left as is. The gradient at "y" = 0 is 0, therefore we can leave "y" as is.

In that case, we stay where we are in the "y" axis direction and increase the "x". By how much should we increase? Is the point (1,0) closer to the minimum than (0,0)? Is it too small an increase? Would (10,0) be a better new point? Or is it too big and should be (0.1,0) instead? How can we get a new point which is closer to the minimum than (0,0) is?

As we approach the minimum, the gradients will start flattening out in both the "x" and "y" directions and start getting closer to 0. This means that we can, as a rule of thumb, increase or decrease the "x" and "y" in the point in proportion to the gradient. The larger the gradient, the steeper the slope, so the more promising moving in that direction would be to find the minimum. On the other hand, the smaller the gradient, the flatter the slope, so moving in that direction would not be very profitable as we're already close to the minimum there and might overshoot. If the gradient is zero, then we are right where we want to be and shouldn't move in that axis at all. So the greater the gradient, the greater the distance we should move.

It's important to understand that the step size should be in proportion to the gradient, not exactly the gradient. We should move at a fraction of the gradient since the gradient will usually be really big compared to the distance required to move and will probably overshoot all the time if used as is. Let's use 0.1 of the gradient.

new x = 0 - 0.1(-6) = 0.6
new y = 0 - 0.1(0) = 0

Notice that we subtracted the gradient so that if the gradient is negative we move towards the positive direction whilst if the gradient is positive we move towards the negative direction since that is how the slope will be oriented.

We check the gradients again to see if we're there yet:

∂z/∂x @ (0.6,0) = 4(0.6) + 2(0) - 6 = -3.6
∂z/∂y @ (0.6,0) = 4(0) + 2(0.6) = 1.2

At our new point we're at a slope sliding downwards towards the positive direction in the "x" and towards the negative direction in the "y" direction. Here are the graphs to confirm.

If we fix "y" to be 0, "z" would vary with respect to "x" like this:

Notice how at "x" = 0.6 the slope steeps down towards the positive direction. The gradient at "x" = 0 is -3.6.

If we fix "x" to be 0.6, "z" would vary with respect to "y" like this:

Notice how at "y" = 0 the slope is steeps down towards the negative directions. The gradient at "y" = 0 is 1.2.

So we need to make "x" bigger and "y" smaller.

new x = 0.6 - 0.1(-3.6) = 0.96
new y = 0 - 0.1(1.2) = -0.12

And we keep on reiterating until the gradient is sufficiently close to 0 in both the "x" and "y" directions or until the change in both "x" and "y" is small enough.

Here's a list of the first 20 points we'll get:
(0, 0)
(0.6, 0.0)
(0.96, -0.12)
(1.2, -0.26)
(1.37, -0.4)
(1.5, -0.51)
(1.6, -0.61)
(1.68, -0.69)
(1.75, -0.75)
(1.8, -0.8)
(1.84, -0.84)
(1.87, -0.87)
(1.9, -0.9)
(1.92, -0.92)
(1.93, -0.93)
(1.95, -0.95)
(1.96, -0.96)
(1.97, -0.97)
(1.97, -0.97)
(1.98, -0.98)

The point gets closer and closer to the actual minimum which is (2,-1). Here are the illustrated points:


See how as the point approaches the minimum it starts moving less and less?

This is called the Gradient Descent algorithm. We can fine-tune its performance by varying the starting point and the step size fraction by which we're multiplying the gradient.

Here's a Python 3 program that performs Gradient Descent given a list of gradient functions (for each variable):
def grad_desc(gradients, initial_values, step_size, threshold):
    old_values = initial_values
    while True:
        new_values = [ value - step_size*gradient(*old_values) for (value, gradient) in zip(old_values, gradients) ]
        if all(abs(new - old) < threshold for (new, old) in zip(new_values, old_values)):
            return new_values
        old_values = new_values
Here's how you would use the previous example:
grad_desc([ lambda x,y:4*x+2*y-6, lambda x,y:4*y+2*x ], [0, 0], 0.1, 0.001)

Sunday, December 13, 2015

(k-)Nearest Neighbour(s) Classification

You want a computer to learn to assign objects into categories, such as the genre of a book. You happen to have a bunch of books with a known category. One of the simplest ways to make the computer assign an unknown book's category is to find the most similar book in the bunch to the unknown book and assume that the two books share the same category. For example, you want to find what genre "Harry Potter" is and find that it is most similar to a book you have called "The Hobbit" which is tagged as fantasy, so you conclude that "Happy Potter" is also fantasy. Of course this only makes sense if you have a big collection of reference books since there might not be any books which are similar otherwise, and the most similar book would be of a genre which is significantly different.

This is called the nearest neighbours classification algorithm, in particular the 1-nearest neighbour, because you only take into consideration the most similar book. Alternatively you can take the top 10 most similar books and use the most frequent genre among the 10 books. This would be called 10-nearest neighbours classification. In general it's called k-nearest neighbours classification.

This is a simple algorithm but its advantage is in its simplicity since it makes no assumptions about the data you give it. Whereas other machine learning algorithms assume that there is some simple pattern to decide which genre a book belongs to, the nearest neighbour classifier can discriminate between very complex patterns and will adapt to any data you train it with, provided that there is enough variety of data. The more complex the relationship between the books and the genre, the more variety of books you need to train it with.

The way it works is by first converting each book in your bunch into a bunch of lists of numbers called a vectors. Each vector would be a point in space (a vector of 2 numbers is a 2D point, of 3 numbers is a 3D point, and the rest are of high dimensional space). For example, in order to convert a book into a point, each number could be the number of times a particular word occurs. Create a vocabulary of words that matter, such as "wizard" and "gun", and then create a point consisting of the number of times each word occurs in the book. So if "Happy Potter" had "wizard" appearing 100 times and "gun" appearing 0 times, then it's 2D point would be (100, 0).

Next, compare the point version of the book in question to the point versions of every book in the bunch. Use some similarity measure to quantify how similar the points are. Similarity measures include Euclidean distance (normal distance between points) and Cosine similarity (difference in the angle of the points from the origin).



Which is the most similar point to the purple one in the above diagram (the purple point is (100, 0) which represents "Harry Potter")? The purple point will be of the same colour as the closest point.

Of course comparing to every point is slow, which is a problem given that nearest neighbour classification requires a lot of points to compare to. There are nearest neighbour search algorithms but they are not very efficient when the points have a lot of dimensions (many numbers in the vector). In some cases it is enough to use approximate search algorithms that do not give exact nearest point, but will find a reasonably close point quickly. The paper "Scaling Distributional Similarity to Large Corpora" gives an overview of such algorithms for finding words that have similar meanings.

If you do not have the genres of the books but still want to categorize similar books together you can use a clustering algorithm such as k-means clustering into order to group books by similarity and then use nearest neighbour classification to associate the new book with the group of the nearest book.

Saturday, January 31, 2015

The Lempel Ziv Welch (LZW) compression algorithm

The Lempel Ziv Welch algorithm (LZW) is a classic compression algorithm published in 1984. It's a simple but practical algorithm that should be under every geek's belt and is often used in combination with other techniques.

The basic idea
Let's start with a plain English description of how this algorithm works. Let's say that we want to compress the following input string:

xxxxyyyyxxxxxxxxxxxx



Compression

If we're lazy, all we have to do to produce a valid output is to represent each letter using 2 characters.

0x0x0x0x0y0y0y0y0x0x0x0x0x0x0x0x0x0x0x0x

According to the compressed language of LZW, when a character pair starts with a "0", that means that the second character is the original letter. This has doubled the size of the sequence, but hopefully this is not the actual output. When the first character is something other than "0", the character pair becomes a reference to some prior long sequence. These references are what will compress the sequence.



In order to compress, we need to use a special table called a "dictionary" which maps 2 character values to the string they represent. Right off the bat, the dictionary will start with the following values:

Actual string2 character value
a0a
b0b
......
x0x
y0y
z0z

We start scanning through the input string from the first letter and find the longest sequence of letters which are in the dictionary. At this point, that would obviously be the first letter, since the dictionary only have single letters.

xxxxyyyyxxxxxxxxxxxx

After finding the longest string from the start which is in the dictionary, "x", including also the next letter in the input will make the shortest string which is not in the dictionary. This unregistered string is "xx".

We replace the longest string found in the dictionary with the corresponding 2 character value:

0xxxxyyyyxxxxxxxxxxxx

The next shortest string not in the dictionary is then added to the dictionary.

Actual string2 character value
a0a
b0b
......
z0z
xx1a

Now we continue scanning the input string from after the last replacement.



Again, we look for the longest string that is in the dictionary. That would be "xx", which we added in the previous step.

0xxxxyyyyxxxxxxxxxxxx

We replace this string with the corresponding 2 character value and add into the dictionary this string plus the next letter ("xxx").

0x1axyyyyxxxxxxxxxxxx

Actual string2 character value
a0a
......
z0z
xx1a
xxx1b

Repeat the process from right after the last replacement.



0x1axyyyyxxxxxxxxxxxx

This time it was "x" on its own that was the longest string in the dictionary since "xy" was not in the dictionary.

0x1a0xyyyyxxxxxxxxxxxx

Actual string2 character value
......
z0z
xx1a
xxx1b
xy1c

Repeat.



0x1a0xyyyyxxxxxxxxxxxx

Longest string in dictionary was "y", shortest string not in dictionary was "yy".

0x1a0x0yyyyxxxxxxxxxxxx

Actual string2 character value
......
z0z
xx1a
xxx1b
xy1c
yy1d

Repeat.



0x1a0x0yyyyxxxxxxxxxxxx

Longest string in dictionary was "yy", shortest string not in dictionary was "yyy".

0x1a0x0y1dyxxxxxxxxxxxx

Actual string2 character value
......
z0z
xx1a
xxx1b
xy1c
yy1d
yyy1e

Repeat.



0x1a0x0y1dyxxxxxxxxxxxx

Longest string in dictionary was "y", shortest string not in dictionary was "yx".

0x1a0x0y1d0yxxxxxxxxxxxx

Actual string2 character value
......
xx1a
xxx1b
xy1c
yy1d
yyy1e
yx1e

Repeat.



0x1a0x0y1d0yxxxxxxxxxxxx

Longest string in dictionary was "xxx", shortest string not in dictionary was "xxxx".

0x1a0x0y1d0y1bxxxxxxxxx

Actual string2 character value
......
xx1a
xxx1b
xy1c
yy1d
yyy1e
yx1e
xxxx1f

Repeat.



0x1a0x0y1d0y1bxxxxxxxxx

Longest string in dictionary was "xxxx", shortest string not in dictionary was "xxxxx".

0x1a0x0y1d0y1b1fxxxxx

Actual string2 character value
......
xx1a
xxx1b
xy1c
yy1d
yyy1e
yx1e
xxxx1f
xxxxx1g

Repeat.



0x1a0x0y1d0y1b1fxxxxx

Longest string in dictionary was "xxxxx" which resulted in the whole input string being consumed, hence there is nothing left to add to the dictionary.

0x1a0x0y1d0y1b1f1g



0x1a0x0y1d0y1b1f1g

This has resulted in an output string which is 2 characters shorter. Obviously the string was engineered to be compressed. Had it been a longer string then the dictionary would have contained longer strings which would lead to more compression.



Decompression
Compression is quite straightforward and so should decompression. Except that it's a little less straightforward because of a special case that can sneak up on you if you don't know about it (some online sources don't mention it).

If we knew what the dictionary contains then we can simply replace each 2 character value with its corresponding string. But in order to know from the start what the dictionary is we would have to include it with the output string, which would be a considerable amount of extra bytes. Instead, decompression basically consists of guessing what the dictionary contained one 2 character value at a time.



We start with the obvious. The dictionary surely contained all the single letters.

Actual string2 character value
0aa
0bb
......
0xx
0yy
0zz

The 2 character value is always one of the above initial dictionary (since that is the first shortest string not in the dictionary), so we go ahead and take care of that.

0x1a0x0y1d0y1b1f1g

x1a0x0y1d0y1b1f1g

From here on, if a 2 character value is in the dictionary then we just replace it with its corresponding string. After each replacement we use the replaced string to update the dictionary (as will be shown in an examples further down). The problem is if the 2 character value is not in the dictionary, as is the case now. This is the special case.



The next 2 character value is not in the dictionary. But notice that it is the very next 2 character value that will enter the dictionary, "1a". When this is the case, the following scenario must have taken place:

The input string has been determined to start with an "x", but the following letters are unknown.
x _ _ _

We know that the following letters were replaced with the next available 2 character value (from the compressed string), which means that it must be the shortest string that was not in the dictionary after "x" was replaced with "0x".

This shortest string must have been the last string that was replaced with a 2 character value ("x"), plus the letter after it. So the dictionary must have looked something like this:

Actual string2 character value
......
0xx
0yy
0zz
1ax_

(Notice that the blank is an unknown letter)

So then the first letter of this unknown string being referred to by the 2 character value "1a" is "x". In that case then we know what the second letter is in the input string: an "x", according to the dictionary we're constructing.

x x _ _

But wait, since the 2 character value "1a" is referring to the first "x" followed by the next letter, and since we have determined that the next letter was "x", then "1a" must be referring to "xx".

Actual string2 character value
......
0xx
0yy
0zz
1axx

In general, every time we encounter this situation, where 2 character value is not in the dictionary and is the next value to be added to the dictionary, the string being referred to is the previous replaced string followed by the same string's first letter. In this case, the previous replaced string was "x" whose first letter is "x", so the string referred to by "1a" is "xx".

x1a0x0y1d0y1b1f1g

xxx0x0y1d0y1b1f1g

After every replacement after the very first (this is the second), we need to update the dictionary with a new 2 character value. The update needs to reflect what was added during compression when the previous replacement (the first in this case) took place. This is because you need to know what letter follows the replacement in order to know which string was added to the dictionary.

Remember that during compression we were adding to the dictionary the shortest string that was not in the dictionary, which consisted of the longest string found in the dictionary (the replaced string) followed by the next letter. After the very first replacement we made, the "x", the next letter in the string is the first letter of the second replacement we made, the "xx". So we add to the dictionary the previous replacement followed by the first letter of the current replacement.

Actual string2 character value
......
0xx
0yy
0zz
1axx

Finally we've covered all the steps needed to get repeatin'. Let's continue.



xxx0x0y1d0y1b1f1g

This one is easy as it is already in the dictionary.

xxxx0y1d0y1b1f1g

The previous replacement was "xx", the current replacement was "x". So to the dictionary we add "xxx".

Actual string2 character value
......
0xx
0yy
0zz
1axx
1bxxx

Repeat.



xxxx0y1d0y1b1f1g

In the dictionary.

xxxxy1d0y1b1f1g

The previous replacement was "x", the current replacement was "y". So to the dictionary we add "xy".

Actual string2 character value
......
0xx
0yy
0zz
1axx
1bxxx
1cxy

Repeat.



xxxxy1d0y1b1f1g

Next value is the next 2 character value to be added to the dictionary. Add previous replacement, "y", followed by its own first letter. So add "yy".

xxxxyyy0y1b1f1g

The previous replacement was "y", the current replacement was "yy". So to the dictionary we add "yy".

Actual string2 character value
......
0xx
0yy
0zz
1axx
1bxxx
1cxy
1dyy

Repeat.



xxxxyyy0y1b1f1g

In the dictionary.

xxxxyyyy1b1f1g

The previous replacement was "yy", the current replacement was "y". So to the dictionary we add "yyy".

Actual string2 character value
......
0xx
0yy
0zz
1axx
1bxxx
1cxy
1dyy
1eyyy

Repeat.



xxxxyyyy1b1f1g

In the dictionary.

xxxxyyyyxxx1f1g

The previous replacement was "y", the current replacement was "xxx". So to the dictionary we add "yx".

Actual string2 character value
......
0xx
0yy
0zz
1axx
1bxxx
1cxy
1dyyy
1eyx

Repeat.



xxxxyyyyxxx1f1g

Next value is the next 2 character value to be added to the dictionary. Add previous replacement, "xxx", followed by its own first letter. So add "xxxx".

xxxxyyyyxxxxxxx1g

The previous replacement was "xxx", the current replacement was "xxxx". So to the dictionary we add "xxxx".

Actual string2 character value
......
0xx
0yy
0zz
1axx
1bxxx
1cxy
1dyyy
1eyx
1fxxxx

Repeat.



xxxxyyyyxxxxxxx1g

Next value is the next 2 character value to be added to the dictionary. Add previous replacement, "xxxx", followed by its own first letter. So add "xxxxx".

xxxxyyyyxxxxxxxxxxxx

Complete.



xxxxyyyyxxxxxxxxxxxx

In practice

In practice, we do not use 2 character values as that would have a large amount of overhead which reduces the amount of compression possible. Instead we work at the bit level and use 12 bit values, a little over 1 byte. The more bits are used, the bigger the dictionary can be and the longer the strings that are added will be, but this needs to be compromised with the overhead.

You can find code for different programming languages here.