The Levenshtein algorithm

The basic Levenshtein algorithm is used to count the minimum number of insertions, deletions and substitutions that are needed to convert one string into another. For instance, to convert English doubt into French doute, you need at least two operations. You could replace the b by a t and then replace the t by an e; or you could delete the b and then insert the e. As this example shows, there may be more than one way to convert one string into another using the minimum number of required operations, but this minimum number itself is unique for each pair of strings.

Implementation in Julia

I won’t cover the logic of the Levenshtein algorithm here. The following is a straightforward Julia implementation of the pseudocode found on Wikipedia, assuming a cost of 1 for all operations. The function takes two inputs (a string a that is to be converted to a string b) and outputs an array with the Levenshtein distances between all substrings of a on the one hand and all substrings of b on the other hand. The entry in the bottom right corner of this array is the Levenshtein distances between the full strings and this is output separately as well.

function levenshtein(a::String, b::String)
    # Initialise table
    distances = zeros(Int, length(a) + 1, length(b) + 1)
    distances[:, 1] = 0:length(a)
    distances[1, :] = 0:length(b)

    # Levenshtein logic
    for row in 2:(length(a) + 1)
        for col in 2:(length(b) + 1)
            distances[row, col] = min(
                distances[row - 1, col - 1] + Int(a[row - 1] != b[col - 1] ? 1 : 0)
                , distances[row, col - 1] + 1
                , distances[row - 1, col] + 1
            )
        end
    end

    return distances, distances[length(a) + 1, length(b) + 1]
end

levenshtein (generic function with 1 method)

Let’s compute the Levenshtein distance between the German word Zyklus (‘cycle’) and its Swedish counterpart cykel. Note the use of ; at the end of the line to suppress the output.

dist_matrix, lev_cost = levenshtein("zyklus", "cykel");
display(dist_matrix)

7×6 Matrix{Int64}:
 0  1  2  3  4  5
 1  1  2  3  4  5
 2  2  1  2  3  4
 3  3  2  1  2  3
 4  4  3  2  2  2
 5  5  4  3  3  3
 6  6  5  4  4  4

This checks out: you do indeed need four operations to transform Zyklus into cykel.

A vectorised function

But what if we wanted to apply our new functions to several pairs of strings? Let’s first define three Dutch-German word pairs:

dutch = ("boek", "zuster", "sneeuw");
german = ("buch", "schwester", "schnee");

We can run our levenshtein() on these three word pairs without introducing for-loops by simply appending a dot to the function name:

levenshtein.(dutch, german)

(([0 1 … 3 4; 1 0 … 2 3; … ; 3 2 … 2 3; 4 3 … 3 3], 3), ([0 1 … 8 9; 1 1 … 8 9; … ; 5 4 … 5 6; 6 5 … 6 5], 5), ([0 1 … 5 6; 1 0 … 4 5; … ; 5 4 … 4 3; 6 5 … 5 4], 4))

However, since the levenshtein() function outputs two pieces of information (both the matrix with the distances between the substrings as well as the final Levenshtein distance), this vectorised call yields a tuple of three subtuples, each subtuple containing both a matrix and the corresponding final Levenshtein distance. This is why the output above looks so messy. If we wanted to obtain just the Levenshtein distances, we could write a for-loop to extract them. But I think an easier solution is to first write a wrapper around the levenshtein() function that outputs only the final Levenshtein distance and use the vectorised version of this wrapper instead:

function lev_dist(a::String, b::String)
  return levenshtein(a, b)[2]
end

lev_dist (generic function with 1 method)

Now use the vectorised version of lev_dist():

lev_dist.(dutch, german)

(3, 5, 4)

Nice!

Obtaining the operations

We now know that we need four operations to transform Zyklus into cykel and five to transform zuster into Schwester. But which are the operations that you need for these transformations? The function lev_alignment() defined below outputs one possible set of operations that would do the job. (Unlike the minimum number of operations required to transform one string into another, the set of operations needed isn’t uniquely defined.)

function lev_alignment(a::String, b::String)
    source = Vector{Char}()
    target = Vector{Char}()
    operations = Vector{Char}()
    
    lev_matrix = levenshtein(a, b)[1]
    
    row = size(lev_matrix, 1)
    col = size(lev_matrix, 2)

    while (row > 1 && col > 1)
        if lev_matrix[row - 1, col - 1] == lev_matrix[row, col] &&
            lev_matrix[row - 1, col - 1] <= min(
              lev_matrix[row - 1, col]
              , lev_matrix[row, col - 1]
              )
            row = row - 1
            col = col - 1
            pushfirst!(source, a[row])
            pushfirst!(target, b[col])
            pushfirst!(operations, ' ')
        else 
            if lev_matrix[row - 1, col] <= min(lev_matrix[row - 1, col - 1], lev_matrix[row, col - 1])
                row = row - 1
                pushfirst!(source, a[row])
                pushfirst!(target, ' ')
                pushfirst!(operations, 'D')
            elseif lev_matrix[row, col - 1] <= lev_matrix[row - 1, col - 1]
                col = col - 1
                pushfirst!(source, ' ')
                pushfirst!(target, b[col])
                pushfirst!(operations, 'I')
            else
                row = row - 1
                col = col - 1
                pushfirst!(source, a[row])
                pushfirst!(target, b[col])
                pushfirst!(operations, 'S')
            end
        end
    end

    # If first column reached, move up.
    while (row > 1)
        row = row - 1
        pushfirst!(source, a[row])
        pushfirst!(target, ' ')
        pushfirst!(operations, 'D')
    end

    # If first row reached, move left.
    while (col > 1)
        col = col - 1
        pushfirst!(source, ' ')
        pushfirst!(target, b[col])
        pushfirst!(operations, 'I')
    end
    
    return vcat(
        reshape(source, (1, :))
        , reshape(target, (1, :))
        , reshape(operations, (1, :))
    )
end

lev_alignment (generic function with 1 method)

I won’t cover the logic behind the algorithm as this is more about learning Julia that the Levenshtein algorithm. On the Julia side, note first how empty character vectors can be initialised. Moreover, notice that the pushfirst!()
function is decorated with a ! (a ‘bang’). This communicates to whoever is reading the code that this function changes some of its input. For instance, pushfirst!(source, a[row]) means that the current character of a (i.e., a[row]) is added to the front of the source vector. That is, this command changes the source vector. Finally, the source, target and operations vectors are all column vectors. In order to display them somewhat nicely, I converted each of them to a single-row matrix using reshape(). Then, the three resulting rows are concatenated vertically using vcat() to show how the two strings can be aligned and which operations are needed to transform one into the other.

Let’s see how we can transform Zyklus into cykel:

lev_alignment("zyklus", "cykel")

3×7 Matrix{Char}:
 'z'  'y'  'k'  ' '  'l'  'u'  's'
 'c'  'y'  'k'  'e'  'l'  ' '  ' '
 'S'  ' '  ' '  'I'  ' '  'D'  'D'

So we substitute c for z, insert an e and delete the u and s. As I mentioned, this set of operations isn’t uniquely defined. Indeed, we could have also substituted c for z, e for l and l for u and then deleted the s. This also corresponds to a Levenshtein distance of four operations.

Normalised Levenshtein distances

Above, we computed raw Levenshtein distances. The problem with these is that longer string pairs will tend to have larger raw Levenshtein distances than shorter string pairs, even if they do seem more similar. To correct for this, we can computed normalised Levenshtein distances instead. There are various ways to compute these; one option is to divide the raw Levenshtein distance by the length of the alignment:

function norm_lev_dist(a::String, b::String)
  raw_dist = lev_dist(a, b)
  alignment_length = size(lev_alignment(a, b), 2)
  return raw_dist / alignment_length
end

norm_lev_dist (generic function with 1 method)

(Behind the scenes, we run the Levenshtein algorithm twice: once in lev_dist() and again in lev_alignment(). This seems wasteful – unless the Julia compiler is able to clean up the double work? I’m not sure.)

We obtain a normalised Levenshtein distance of about 0.57 for Zyklus – cykel:

norm_lev_dist("zyklus", "cykel")

0.5714285714285714

We can use a vectorised version of this function, too:

norm_lev_dist.(dutch, german)

(0.75, 0.5555555555555556, 0.5)

Of course, normalised Levenshtein distances are symmetric, so we obtain the same result when running the following command:

norm_lev_dist.(german, dutch)

(0.75, 0.5555555555555556, 0.5)

The Fibonacci sequence

The famous Fibonacci sequence is an infinite sequence of natural numbers, the first of which are 1, 1, 2, 3, 5, 8, 13, …. The sequence is defined as follows:

Let’s write some Julia functions that can generate this sequence.

Julia

You can download Julia from julialang.org. I’m currently using the Pluto.jl package that allows you to write Julia code in a reactive notebook. Check out the Pluto.jl page for more information.

First alternative: A purely recursive function

The Fibonacci sequence is defined recursively: To obtain the nth Fibonacci number, you first need to compute the n-1th and n-2th Fibonacci number and then add them. We can write a Julia function that exactly reflects the definition of the Fibonacci sequence like so:

function fibonacci(n)
  if n <= 2
    return 1
  end
  return fibonacci(n - 1) + fibonacci(n - 2)
end

fibonacci (generic function with 1 method)

This function tacitly assumes that n is a non-zero natural number. If n is equal to or lower than 2, i.e., if n is 1 or 2, it immediately returns 1, as per the definition of the sequence. If this condition isn’t met, the output is computed recursively. The function can be run as follows:

fibonacci(10)

55

Checks out! But from a computational point of view, the fibonacci() function is quite wasteful. In order to obtain fibonacci(10), we need to compute fibonacci(9) and fibonacci(8). But in order to compute fibonacci(9), we also need to compute fibonacci(8). For both fibonacci(9) and fibonacci(8), we need to compute fibonacci(7), etc. In fact, we need to compute the value of fibonacci(8) two times, that of fibonacci(7) three times, that of fibonacci(6) five times, and that of fibonacci(5) seven times. So we’d be doing lots of computations over and over again. For this reason, the fibonacci() function is hopelessly inefficient: While you can compute fibonacci(10) in a fraction of a second, it may take minutes to compute, say, fibonacci(60). Luckily, we can speed up our function considerably.

Second alternative: Recursion with memoisation

Memoisation is a programming technique where any intermediate result that you computed is stored in an array. Before computing any further intermediate results, you then first look up in the array if you haven’t in fact already computed it, saving you a lot of unnecessary computations. The following Julia function is a bit more involved that the previous one, but it’s much more efficient.

function fib_memo(n)
  known = zeros(Int64, n)
  function memoize(k)
    if known[k] != 0
      # do nothing
    elseif k == 1 || k == 2
      known[k] = 1
    else
      known[k] = memoize(k-1) + memoize(k-2)
    end
    return known[k]
  end
  return memoize(n)
end

fib_memo (generic function with 1 method)

The overall function that we’ll actually call is fib_memo(). It creates an array called known with n zeroes. Then it defines an inner function memoize(). This latter function obtains an integer k that in practice will range from 0 to n and does the following. First, it checks if the kth value in the array known is still 0. If it got changed, the function just returns the kth value in known. Otherwise, if k is equal to either 1 or 2, it sets the first or second value of known to 1. If k is greater than 2, the kth value of known is computed recursively. In all cases, the memoize() function returns the k value of the known array. The outer fib_memo() function then just returns the result of memoize(n).

Perhaps by now, your computer has finished running fibonacci(60) and you can try out the alternative implementation:

fib_memo(60)

1548008755920

Notice how much faster this new function is! Even the 200th Fibonacci number can be computed in a fraction of a second:

fib_memo(200)

-1123705814761610347

Unfortunately, we’ve ran into a different problem now: integer overflow. The result of the computations has become so large that it exceeded the range of 64-bit integers. To fix this problem, we can work with BigIntegers instead:

function fib_memo(n)
  known = zeros(BigInt, n)
  function memoize(k)
    if known[k] != 0
      # do nothing
    elseif k == 1 || k == 2
      known[k] = 1
    else
      known[k] = memoize(k-1) + memoize(k-2)
    end
    return known[k]
  end
  return memoize(n)
end

fib_memo (generic function with 1 method)

fib_memo(200)

280571172992510140037611932413038677189525

Nice!

Third alternative: Using Binet’s formula

The third alternative is more of a mathematical solution rather than a programming solution. According to Binet’s formula, the nth Fibonacci number can be computed as where , the Golden Ratio, and , its conjugate. In Julia:

function fib_binet(n)
  φ = (1 + sqrt(5))/2
  ψ = (1 - sqrt(5))/2
  fib_n = 1/sqrt(5) * (φ^n - ψ^n)
  return BigInt(round(fib_n))
end

fib_binet (generic function with 1 method)

Note that you can use mathematical symbols like and in Julia. This function runs very fast, too:

fib_binet(60)

1548008755920

fib_binet(200)

280571172992512015699912586503521287798784

Notice, however, that the result for the 200th Fibonacci number differs by 27 orders of magnitude from the one obtained using fib_memo():

fib_binet(200) - fib_memo(200)

1875662300654090482610609259

By using Binet’s formula, we’ve left the fairly neat world of integer arithmetic and entered the realm of floating point arithmetic that is rife with approximation errors. While we’re at it, we might as well compute and plot the size of these approximation errors. In the snippet below, I first use list comprehension in order to compute the first 200 Fibonacci numbers using both fib_memo() and fib_binet(). Note that I added a dot (.) to both function names. This is Julia notation for running vectorised computations. Further note that I end all lines with a semi-colon so that the results don’t get printed to the prompt. Then, I compute the absolute values of the differences between the numbers obtained by both computation methods. Note again the use of a dot in both abs.() and .- that is required to have both of these functions work on vectors. Finally, I convert these absolute differences to differences relative to the correct answers;

fib_integer = fib_memo.(1:200);
fib_math    = fib_binet.(1:200);
abs_diff = abs.(fib_math .- fib_integer);
rel_diff = abs_diff ./ fib_integer;

To wrap off this blog post, let’s now plot these absolute and relative differences using the Plots.jl package. While Figure 1 shows that the absolute error becomes huge, Figure 2 shows that these discrepancies only amount to a negligble fraction of the correct answers.

using Plots
plot(1:200, abs_diff, seriestype=:scatter,
     xlabel = "n",
     ylabel = "absolute difference",
     label = "")

Figure 1. Absolute difference between the Fibonacci numbers obtained using fib_binet() and those obtained using fib_memo().

plot(1:200, rel_diff, seriestype=:scatter, 
     xlabel = "n",
     ylabel = "relative difference",
     label = "")

Figure 2. Relative difference between the Fibonacci numbers obtained using fib_binet() and those obtained using fib_memo().