Linear recurrence relations (or how to solve recursive sequences effectively)

This article will have quite a few equations but they are not as hard to follow as they look so bear with me.

Every programmer has written a program that finds the nth Fibonacci number. Some of you did it with recursion, others went a step ahead and used dynamic programming. Still, even with dynamic programming your algorithm was $O(n)$. We can go even further and write an $O(\log{n})$ algorithm for finding the nth number of any sequence of the type:

where $k$  is a positive integer, $C_i (i = 0, 1,2, \cdots, k)$ are some constants and $a_{n-i}$ is the (n - i)-th member of the sequence. A sequence defined as above is called a linear recurrence homogenous relation (LRHR) of order k.

Clearly, the Fibonacci sequence is a LRHR of order k = 2 since it is defined as:

Let me now show you how to compute the nth Fibonacci number not only in time complexity of $O(\log{n})$, but also in constant space complexity (in contrast to dynamic programming which requires linear)

To do so we must find a closed form of our relation. We'll explore how to do that by first giving the general solution to such a problem and then applying it to the Fibonacci sequence as an example.

Algorithm for solving linear recurrence homogenous relations

1. In the equation $a_n=C_1 a_{n-1}+C_2 a_{n-2}+\cdots+C_k a_{n-k}$ replace all the $a_{n-i}$ with  $x^{k-i} (i=0,1,\cdots, k)$ and you get:

Call the polynomial $x^k - (C_1 x^{k-1}+C_2 x^{k-2}+\cdots+C_{k-1} x+C_k)$ a characteristic polynomial.

Applied to our example:

Replacing the $a_{n-i}$ as described above gives us $x^2=x+1$.

Therefore, the characteristic polynomial of the Fibonacci sequence is $x^2-x-1$.

2. Find the roots of the characteristic polynomial.
Applied to our example:

Roots of the Fibonacci's polynomial are $x=\frac{1\pm\sqrt{5}}{2}$. (As a side note, the number $x=\frac{1+\sqrt{5}}{2}$ is known as the Golden Ratio)

3. Let the roots of the characteristic polynomial be $r_1,r_2,\cdots,r_k$. For now let's assume that there are no repeated roots. Then the closed form of our sequence is

where $A_i$ are some unknown (for now) constants.

We will consider the repeated roots case later.

Applied to our example:

For the Fibonacci sequence, applying step 3 gives us:

4. Given that you know the first $k$ members of the sequence, you can now calculate the unknown constants $A_i$ by substituting the values of the first members in the equation of step 3.
Applied to our example:

The first two members of the Fibonacci sequence are $a_0=0, a_1=1$.

Substitute for $n=0$ and we get:

Substitute for $n=1$ and we arrive at:

But since we just found out that $A_1=-A_2$:

Therefore:

5. Go back to the equation in step 3 and replace all the $A_i$ with the values that you just calculated to receive the final closed form.
Applied to our example:

Substitute $A_1=\frac{1}{\sqrt{5}},A_2=-\frac{1}{\sqrt{5}}$ in $a_n=A_1 \left(\frac{1+\sqrt{5}}{2}\right)^n+A_2 \left(\frac{1-\sqrt{5}}{2}\right)^n$ and:

Finally, the closed form of the Fibonacci sequence is:

Obviously the above formula can be calculated in $O(\log{n})$ so that completes our algorithm.

Here's an example implementation in C#. I tested it by computing the first 92 Fibonacci numbers a million times (the 92nd number is the last one that can be stored in a Int64). Results prove that direct computation with the formula is much faster(every number was calculated a million times in ~0.17s), while dynamic programming is slower for all numbers after the 18th. Average running time (over 10 tests):

 Direct formula Dynamic programming 15.8 s 33.575 s

Here's the same code implemented in JavaScript for you to play with.

There are a few things to note though:

1. The roots in the Fibonacci polynomial are not integers. Therefore, we cannot hope to effectively calculate every single number in the sequence due to floating – point arithmetic errors. Test the example implementations above and see for yourself that the last number to be calculated correctly is the 71st for the C# code and the 75th for the JS one. If the roots were integers it would've been a whole another story.
2. In the JS example, even the dynamic programming approach fails after 78th number due to the way JavaScript integers work. To see the error, compare the results from our example with this calculator.

Repeated Roots

Now as promised let's take a look at what happens when one of the roots of the characteristic polynomial is repeated (step 3). For simplicity's sake, let's assume that only the first root is repeated and let $m$ denote the number of times it is repeated. The last statement is equivalent to saying that the roots of our polynomial are $r_1,r_2,\cdots,r_{k-m}$. In this case, our closed form solution will look like this:

(Of course if more than one root repeats, do the same for it)

And since that last equation looks a little horrifying, let's give an example. The characteristic roots of the sequence $a_n=4a_{n-1}-4a_{n-2}$ are $r_1=r_2=2$.

Applying what we said above we get the following closed form:

Say that the first 2 members of the sequence are 0 and 1 as in the Fibonacci numbers. Substitute for $a_0=0$ and you get:

Substitute $a_1=1$:

Therefore the closed form solution of the relation is precisely:

Summary:

For every linear recurrence homogenous relation there exists a closed form solution. Hence, since you can model any recursive function like a sequence, you can find a closed form for every recursive function knowing only its first few values.

This algorithm may not be the remedy for everything, but it sure kicks your recursive function into turbo mode. Best results are achieved when the roots of the characteristic polynomial are integers and the sequence is slowly increasing (unlike the Fibonacci which quickly overflows 64-bits).

It may not always be the case you need the performance boost from using LRHR over DP but it is always a good thing to know that such an algorithm exists if a problem of this kind arises.

Things to note:

• The order $k$ must be a constant. (i.e. the factorial function cannot be solved by the above algorithm since it depends on all previous values which is a variable amount)