Hi everyone!
There are already dozens of blogs on linear recurrences, why not make another one? In this article, my main goal is to highlight the possible approaches to solving linear recurrence relations, their applications and implications. I will try to derive the results with different approaches independently from each other, but will also highlight similarities between them after they're derived.
Definitions
Def. 1. An order $$$d$$$ homogeneous linear recurrence with constant coefficients (or linear recurrence) is an equation of the form
Def. 2. In the equation above, the coefficients $$$a_1, \dots, a_d \in R$$$ are called the recurrence parameters,
Def. 3. and a sequence $$$F_0, F_1, \dots \in R$$$ is called an order $$$d$$$ linear recurrence sequence.
The most common task with linear recurrences is, given initial coefficients $$$F_0, F_1, \dots, F_{d-1}$$$, to find the value of $$$F_n$$$.
Example 1. A famous Fibonacci sequence $$$F_n = F_{n-1} + F_{n-2}$$$ is an order 2 linear recurrence sequence.
Example 2. Let $$$F_n = n^2$$$. One can prove that $$$F_n = 3 F_{n-1} - 3 F_{n-2} + F_{n-3}$$$.
Example 3. Moreover, for $$$F_n = P(n)$$$, where $$$P(n)$$$ is a degree $$$d$$$ polynomial, it holds that
If this fact is not obvious to you, do not worry as it will be explained further below.
Finally, before proceeding to next sections, we'll need one more definition.
Def. 4. A polynomial
is called the characteristic polynomial of the linear recurrence defined by $$$a_1, \dots, a_d$$$.
Example 4. For Fibonacci sequence, the characteristic polynomial is $$$A(x) = x^2-x-1$$$.
Matrix approach
For linear recurrence sequences it holds that
\left(\begin{bmatrix} a_1 & a_2 & \dots & a_{d-1} & a_d \\ 1 & 0 & \dots & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \dots & 1 & 0 \end{bmatrix} \right)^{n} \begin{bmatrix}F_{d-1} \\ F_{d-2} \\ \dots \\ F_{0}\end{bmatrix} = A^n \mathbf{F}_0. $$$
This representation allows to find $$$F_n$$$ in $$$O(d^3 \log n)$$$ with binary matrix exponentiation.
Improving the complexity of this solution is possible with the following theorem:
Theorem 1 (Cayley–Hamilton theorem). Every square matrix over a commutative ring satisfies its own characteristic polynomial.
For the matrix $$$A$$$ used above, it can be proven that its characteristic polynomial is $$$(-1)^dA(x)$$$, thus
For the purpose of calculating $$$F_n$$$ it means that we may find $$$b_0, b_1, \dots, b_{d-1}$$$ such that
and this equivalence would also hold if we substitute $$$x$$$ with $$$A$$$. This also implies that
thus allowing to compute $$$F_n$$$ in $$$O(d \log d \log n)$$$ using the binary lifting modulo $$$A(x)$$$ (see on CP-Algorithms).
Approach above is intuitive and good if you want just to find the $$$n$$$-th element of the recurrence. However it might be hard to analyze sequences with it (e.g. estimate the asymptotic behavior or find a recurrence for a sum of sequences).
Problem example: CodeChef — Random Number Generator. Just find $$$F_n$$$, where $$$F_k$$$ is a given linear recurrence.
Generating function approach
Consider the generating function $$$G(x) = F_0 + F_1 x + F_2 x^2 + \dots$$$ of the sequence. For $$$G(x)$$$ it holds that
where $$$P(x)$$$ is some polynomial of degree less than $$$d$$$ used to calibrate the fact that $$$F_0, \dots, F_{d-1}$$$ are pretty much arbitrary and might not adhere to the recurrence with lower terms. From this equation we see that
where $$$Q(x)$$$ is a polynomial obtained by reversing the coefficients of $$$A(x)$$$, that is
Side note: If you know $$$F_0, \dots, F_{d-1}$$$ and $$$a_1, \dots, a_d$$$, you can find $$$P(x)$$$ from
in other words,
Now, if we allow negative powers near $$$x$$$, the following equation will hold:
where $$$[x^k] P(x)$$$ is the coefficient near $$$x^k$$$ in the series $$$P(x)$$$. We should notice that
has only positive coefficients, so is $$$x^k A(x) G(x^{-1})$$$ and so is all their linear combinations. In other words,
whenever $$$X(x)-Y(x)$$$ is divisible by $$$A(x)$$$. From this, we conclude that
which is essentially the same result as in the matrix approach.
Graeffe's method
This approach was described by Golovanov399 in this comment.
Alternatively to finding $$$x^n \bmod A(x)$$$, you may consider the following transform:
For even $$$n$$$, it is equal to $$$[x^{n/2}]\frac{S_0(x)}{T(x)}$$$ and for odd $$$k$$$ it is equal to $$$[x^{(n-1)/2}]\frac{S_1(x)}{T(x)}$$$.
It works with the same complexity, but has a better constant.
Side note: Same approach could be used to find the polynomial inverse in $$$O(n \log n)$$$:
Closed-form solution
While it could be hard to analyze the asymptotic behavior of $$$F_n$$$ with matrices (you probably can, but you'd need to work with eigenvalues of $$$A$$$), it is way simpler with the generating function representation, thanks to the existence of the partial fraction decomposition:
where $$$\lambda_1, \dots, \lambda_m$$$ are the roots of $$$A(x)$$$ with multiplicities $$$d_1,\dots,d_m$$$, and $$$C_i(x)$$$ is a polynomial of degree less than $$$d_i$$$. In other words,
It is generally known that
from which we may conclude that the general solution to the linear recurrence is given by
where $$$p_i(n)$$$ is a polynomial of degree less than $$$d_i$$$, determined from $$$F_0, \dots, F_{d-1}$$$.
Example 3 (follow-up): in particular, the formula above means that generating function for $$$p(0), p(1), \dots$$$, where $$$p(x)$$$ is a polynomial of degree $$$d$$$, has a denominator $$$(1-x)^{d+1}$$$. This proves the formula from the example 3.
Problem example: Timus — 2125. Continue the Sequence. You're given a sequence $$$a_1,\dots,a_n$$$ and a number $$$m$$$. Let $$$p(x)$$$ be a minimum-degree polynomial such that $$$p(1)=a_1,\dots, p(n)=a_n$$$. Find $$$p(n+1),\dots,p(n+m)$$$ in $$$O(n \log n)$$$.
Umbral calculus approach
Thanks to ftiasch for letting me know about umbral calculus! This approach makes so much more sense to me now.
Let's do something esoteric. We're given the recurrence
Consider formal variable $$$b$$$ such that
Using the formula above we may express any $$$b^n$$$ as a linear combination of $$$b^0, \dots, b^{d-1}$$$ by calculating it modulo $$$A(b)$$$. If, after that, we do a back-substitution of $$$F_k$$$ instead of $$$b^k$$$, we will have a formula that is true for any sequence $$$F$$$ adhering to the recurrence.
What we did is known as the umbral calculus. Formal theory behind it is as follows.
Consider a linear functional $$$L : R[b] \to R$$$ such that
By the definition of the functional,
As it is a linear functional, all possible linear combinations of $$$b^k A(b)$$$ would also give zero. In other words,
when $$$X(b)-Y(b)$$$ is divisible by $$$A(b)$$$. In other words, $$$\langle A(b) \rangle \subset \ker L$$$, where $$$\langle A(b) \rangle$$$ is the ideal generated by $$$A(b)$$$.
In particular, it is true for $$$X(b) = b^n$$$ and $$$Y(b) = b^n \bmod A(b)$$$, leading us to the $$$b^n \bmod A(b)$$$ solution again.
Side note: Comparing this approach with the genfunc one, you may notice that $$$L(P)=[x^0]P(x)G(x^{-1})$$$.
This is the approach you could've seen earlier in TLE's blog. Special thanks to islingr for a very clear explanation.
I wonder if it can be used to derive closed-form solution...
Shift operator approach
It is of little use in competitive programming, but I like it because of its genericity.
Let $$$\Delta : R^\mathbb{N} \to R^\mathbb{N}$$$ be a linear operator defined on sequences over $$$R$$$, such that $$$(\Delta a)_n = a_{n+1}$$$.
With this operator, the recurrence rewrites as
or, which is equivalent, as
Finding all possible $$$F$$$ is equivalent to parametrizing the elements of $$$\ker A(\Delta)$$$. The following result holds:
Theorem 2. Let $$$A: X \to X$$$ be a linear operator defined on a complex vector space $$$X$$$. Let $$$\lambda_1, \dots, \lambda_m$$$ be different complex numbers and $$$d_1, \dots, d_m$$$ be positive integer numbers. Then the following holds:
Here $$$\ker A = \{x : Ax = 0\}$$$ is the kernel of the linear operator $$$A$$$ and $$$\Sigma$$$ denotes the Minkowski sum of corresponding subspaces.
Proof sketch: If $$$X$$$ is finite-dimensional space, one can prove this by considering the Jordan normal form of the corresponding matrix. For infinite-dimension case, it is clear that right-hand side is a subset of left-hand side. For any vector $$$x$$$ from the left-hand side, the span of $$$x, Ax, A^2 x, \dots$$$ (known as a Krylov subspace) is finite-dimensional and you can prove that $$$x$$$ is also in the right-hand side by reducing all operators to this finite-dimensional subspace (it's invariant for all operators in the equation).
Thanks to tranquility for kindly proving this result to me!
What this means for us is that we may solve the problem independently for all $$$(\Delta-\lambda)^{d}$$$, where $$$\lambda$$$ are roots of $$$A(x)$$$ and $$$d$$$ are their multiplicities, and the final solution will be a linear combination of these basic solutions.
For $$$d = 1$$$, the equation here is of form
or, equivalently,
The solution to this is simply $$$F_n = C \cdot \lambda^n$$$, where $$$C$$$ is arbitrary constant.
As for higher-dimensional case, it can be reduced to the lower-dimensional as
In other words, let's denote by $$$F^{(k)}$$$ the generic solution to $$$(\Delta-\lambda)^k F^{(k)} = 0$$$. Then the equation above rewrites as
with the base case $$$F^{(1)}_n = C_1 \cdot \lambda^n$$$. For example,
which has a generic solution of form $$$F^{(2)}_n = (a+bn)\lambda^n$$$.
For this series of equations it can be proven that $$$F_n^{(k)} = C_k(n) \cdot \lambda^n$$$, where $$$C_k(n)$$$ is a polynomial of degree less than $$$k$$$. This essentially proves the closed-form solution to linear recurrences.
If you have some intuition about Jordan normal form and generalized eigenvectors, you can compare them to the construction above!
As you may see, $$$F^{(k)}$$$ is the generalized eigenvector of $$$\Delta$$$ of rank $$$k$$$, and the construction itself is very similar with Jordan chains.
In a very similar manner and using the same theorem you may prove that the solution to the linear differential equation
is given as
where $$$\lambda_i$$$ are the roots of differential equation's characteristic polynomial and $$$p_i(x)$$$ are polynomials of degree less than $$$d_i$$$.
To get it, one should consider operator $$$A = \frac{\partial}{\partial x}$$$, hence the equation is
Comments (0)
