## Recurrence Relation and Matrix Exponentiation

Recurrence relations appear many times in computer science. Using recurrence relation and dynamic programming we can calculate the *n ^{th}* term in O(n) time. But many times we need to calculate the

*n*in O(log n) time. This is where Matrix Exponentiation comes to rescue.

^{th}We will specifically look at linear recurrences. A linear recurrence is a sequence of vectors defined by the equation *X _{i+1}* = M

*X*for some constant matrix M. So our aim is to find this constant matrix M, for a given recurrence relation.

_{i}Let’s first start by looking at the common structure of our three matrices *X _{i+1}*,

*X*and M. For a recurrence relation where the next term is dependent on last k terms,

_{i}*X*and

_{i+1}*X*are matrices of size 1 x k and M is a matrix of size k x k.

_{i+1}| f(n+1) | | f(n) | | f(n) | | f(n-1) | | f(n-1) | = M x | f(n-2) | | ...... | | ...... | | f(n-k+1) | | f(n-k) |

Let’s look at different type of recurrence relations and how to find M.

**1.** Let’s start with the most common recurrence relation in computer science, The Fibonacci Sequence. *F _{i+1}* =

*F*+

_{i}*F*.

_{i-1}| f(n+1) | = M x | f(n) | | f(n) | | f(n-1) |

Now we know that M is a 2 x 2 matrix. Let it be

| f(n+1) | = | a b | x | f(n) | | f(n) | | c d | | f(n-1) |

Now a*f(n) + b*f(n-1) = f(n+1) and c*f(n) + d*f(n-1) = f(n). Solving these two equations we get a=1,b=1,c=1 and d=0. So,

| f(n+1) | = | 1 1 | x | f(n) | | f(n) | | 1 0 | | f(n-1) |

**2.** For recurrence relation f(n) = a*f(n-1) + b*f(n-2) + c*f(n-3), we get

| f(n+1) | = | a b c | x | f(n) | | f(n) | | 1 0 0 | | f(n-1) | | f(n-1) | | 0 1 0 | | f(n-2) |

**3.** What if the recurrence relation is f(n) = a*f(n-1) + b*f(n-2) + c, where c is a constant. We can also add it in the matrices as a state.

| f(n+1) | = | a b 1 | x | f(n) | | f(n) | | 1 0 0 | | f(n-1) | | c | | 0 0 1 | | c |

**4.** If a recurrence relation is given like this f(n) = f(n-1) if n is odd, f(n-2) otherwise, we can convert it to f(n) = (n&1) * f(n-1) + (!(n&1)) * f(n-2) and substitute the value accordingly in the matrix.

**5.**If there are more than one recurrence relation, g(n) = a*g(n-1) + b*g(n-2) + c*f(n), and, f(n) = d*f(n-1) + e*f(n-2). We can still define the matrix X in following way

| g(n+1) | | a b c 0 | | g(n) | | g(n) | = | 1 0 0 0 | x | g(n-1) | | f(n+2) | | 0 0 d e | | f(n+1) | | f(n+1) | | 0 0 1 0 | | f(n) |

Now that we have got our matrix M, how are we going to find the nth term.

*X _{i+1}* = M

*X*

_{i}(Multiplying M both sides)

M *

*X*= M * M

_{i+1}*X*

_{i}*X*= M^2

_{i+2}*X*

_{i}..

*X*= M^k

_{i+k}*X*

_{i}So all we need now is to find the matrix M^k to find the k-th term. For example in the case of Fibonacci Sequence,

M^2 = | 2 1 | | 1 1 |

Hence F(2) = 2.

Now we need to learn to find M^k in *O*(*n ^{3}* log

*b*) time. The brute force approach to calculate a^b takes O(b) time, but using a recursive divide-and-conquer algorithm takes only O(log b) time:

- If b = 0, then the answer is 1.
- If b = 2k is even, then
*a*= (^{b}*a*)^{k}.^{2} - If b is odd, then
*a*=^{b}*a * a*.^{b-1}

We take a similar approach for Matrix Exponentiation. The multiplication part takes the O(*n ^{3}*) time and hence the overall complexity is

*O*(

*n*log

^{3}*b*). Here’s a sample code in C++ using template class:

#include<iostream> using namespace std; template< class T > class Matrix { public: int m,n; T *data; Matrix( int m, int n ); Matrix( const Matrix< T > &matrix ); const Matrix< T > &operator=( const Matrix< T > &A ); const Matrix< T > operator*( const Matrix< T > &A ); const Matrix< T > operator^( int P ); ~Matrix(); }; template< class T > Matrix< T >::Matrix( int m, int n ) { this->m = m; this->n = n; data = new T[m*n]; } template< class T > Matrix< T >::Matrix( const Matrix< T > &A ) { this->m = A.m; this->n = A.n; data = new T[m*n]; for( int i = 0; i < m * n; i++ ) data[i] = A.data[i]; } template< class T > Matrix< T >::~Matrix() { delete [] data; } template< class T > const Matrix< T > &Matrix< T >::operator=( const Matrix< T > &A ) { if( &A != this ) { delete [] data; m = A.m; n = A.n; data = new T[m*n]; for( int i = 0; i < m * n; i++ ) data[i] = A.data[i]; } return *this; } template< class T > const Matrix< T > Matrix< T >::operator*( const Matrix< T > &A ) { Matrix C( m, A.n ); for( int i = 0; i < m; ++i ) for( int j = 0; j < A.n; ++j ) { C.data[i*C.n+j]=0; for( int k = 0; k < n; ++k ) C.data[i*C.n+j] = C.data[i*C.n+j] + data[i*n+k]*A.data[k*A.n+j]; } return C; } template< class T > const Matrix< T > Matrix< T >::operator^( int P ) { if( P == 1 ) return (*this); if( P & 1 ) return (*this) * ((*this) ^ (P-1)); Matrix B = (*this) ^ (P/2); return B*B; } int main() { Matrix<int> M(2,2); M.data[0] = 1;M.data[1] = 1; M.data[2] = 1;M.data[3] = 0; int F[2]={0,1}; int N; while (~scanf("%d",&N)) if (N>1) printf("%lld\n",(M^N).data[0]); else printf("%d\n",F[N]); }

Nice job…well explained

RanjanMay 18, 2011 at 12:53 PM

Good explanation. I have just learnt a new thing. An interesting one, indeed.

Just a question. Can not the matrix exponentiation be done with the diagonalization method? It would be faster, although I am not very sure. I am a little sleepy right now.

ArekuMay 22, 2011 at 3:53 AM

Just a question. Can not the matrix exponentiation be done with the diagonalization method? It would be faster, although I am not very sure. I am a little sleepy right now.+1

JabezJune 5, 2011 at 1:43 AM

hi how your code is high lighted in the wordpress

can you please tell me up??

i am not getting plugin option also in my word press did you paid for this??

gauravalgoJune 17, 2011 at 5:33 PM

You do not need any plugins to do that. Just have a look at the support section. More precisely http://en.support.wordpress.com/code/ for code syntax highlighter and http://en.support.wordpress.com/latex/ for latex.

fR0DDYJune 17, 2011 at 9:26 PM

[…] Using Matrix Multiplication In the last post we learned how to use Fast Matrix Multiplication to calculate functions having linear equations in […]

Combination « COME ON CODE ONJuly 31, 2011 at 8:43 PM

i want to know that can Determinant of matrix be find out without using any floating point operations.

gauravalgoAugust 16, 2011 at 3:59 PM

your code in pyton and simple c is better in understanding(also many coders understand like me) than c++ template so please make your write ups in these languages

gauravalgoAugust 27, 2011 at 10:09 PM

i want to ask at your 2 nd point of this article

“2. For recurrence relation f(n) = a*f(n-1) + b*f(n-2) + c*f(n-3)”

how do we calculate a,b,c by applying (subset sum algorithm) or by some else how?

and second thing that your code in c++ template so may be

cannot undestand by many users so i have this techical writing for those

http://algosolver.blogspot.com/2011/09/fibonacci-log-n.html

enjoy ; )

gauravalgoSeptember 1, 2011 at 6:00 PM

The first K terms of the sequence F(0), F(1), F(2), … , F(K-1) are given and the recurrence is, F(n) = F(n-1) + F(n-2) + … + F(n-K) for n >= K. Given N, find the value of F(N).

vcfvamcvmcOctober 27, 2011 at 2:29 AM

@fRODDY….i din understand how u’ve written x_i+1 matrix as g_n+1,g_n,f_n+2,f_n+1 ???why isn’t it f_n+1,f_n????please clarify???how can we deduce the no of elements in the k*1 matrix if reccurences like this are given????

LalitMay 16, 2012 at 2:02 PM

Can we calculate this recurrence equation using matrix ????

f(n)=f(n-1)+f(n-2)-2^n

if yes then how to calculate????please help

Nafis IslamSeptember 10, 2012 at 1:40 AM

Can anyone tell me how the matrix turns out to be {2,2,1,0} for f(n) = f(n-1)+4f(n-2)+2f(n-3) in the below link?

http://www.codechef.com/viewsolution/547362

unknownJuly 15, 2013 at 10:52 PM

can u give me the exact order of matrices! ie the matrice (Xi+1) and matric(Xi) ‘s order to be use while programming this approach

paralizedSeptember 27, 2014 at 2:10 PM