Matrix Exponentiation

Matrix Exponentiation is a useful technique that can be applied to a wide variety of problems. Most commonly it is used for efficiently solving linear recurrence problems and as such can be used in any problem that can be represented as a linear recurrence. The main problem lies with the fact that it is inefficient when we naively evaluate them. Let's consider an easy example that you should know about: Fibonacci numbers.

The sequence is defined by:

It's easy enough to increasingly evaluate it 1 by 1 each time until we hit our targetted n-th Fibonacci number. However, let's say n was 1 billion - we need at least 1 billion computations of the sequence to derive the answer. If we represent the sequence in terms of a matrix we discover a clever optimisation which can be applied.

The correctness can be verified by a simple matrix multiplication. What's more important to note is that the next two sequence of Fibonacci numbers can be represented like the one above.

If we expand:

This is equivalent to:

We can continue expanding and eventually we see a remarkable observation. We can always reduce it to the power of our matrix multiplied by our initial conditions. This yields the recurrence below:

 

Since exponentiation can be done in logarithmic time we have essentially reduced our Fibonacci computation speed from linear to logarithmic. A huge difference in speed for large numbers! As an exercise, write a logarithmic Fibonacci term generator - as the numbers do get quite large, feel free to modulo it with an appropriate number (like 10^8).

We now expand our knowledge to a slightly more complicated linear recurrence. This is based on the problem "Number Sequences": http://acm.tju.edu.cn/toj/showp2169.html

To solve this recurrence we simply just do the same thing as we did with Fibonacci with a slight modification. Fibonacci is really a special case of the above recurrence where x = 1 and y = 1 and a0 = a1 = 1. Note the prime observation below:

If we matrix multiply - we yield the correct recurrence for both Fn and Fn-1. Therefore we simply need to change our original Fibonacci matrix of [ 1 1, 1 0 ] to [ x y, 1 0] and the initial conditions from being always 1 and 1 (F1 and F0 respectively) to [a1 a0]. We then simply use matrix exponentiation to calculate the correct term, as always we apply modulo arithmetic to keep the number representable with integers.

The implementation below demonstrates this idea:

class Matrix {
public:
   long long a;
   long long b;
   long long c;
   long long d;
   Matrix() { }
   Matrix(long long _a, long long _b, long long _c, long long _d) : 
      a(_a), b(_b), c(_c), d(_d) { }
};

Matrix multiply(const Matrix& lhs, const Matrix& rhs) {
   long long a = lhs.a * rhs.a + lhs.b * rhs.c;
   long long b = lhs.a * rhs.b + lhs.b * rhs.d;
   long long c = lhs.c * rhs.a + lhs.d * rhs.c;
   long long d = lhs.c * rhs.b + lhs.d * rhs.d;
   return Matrix(a % 100, b % 100, c % 100, d % 100);
}

ostream& operator<<(ostream& out, const Matrix& M) {
   out << "a: " << M.a << " b: " << M.b << " c: " << M.c << " d: " << M.d;
   return out;
}

Matrix power(Matrix M, int n) {
   if (n == 1) return M;
   Matrix tmp = power(M, n/2);
   if (n % 2 != 0) {
      Matrix n = multiply(tmp,tmp);
      Matrix m = multiply(n, M);
      return m;
   }
   Matrix v = multiply(tmp,tmp);
   return v;
}

void output(int n) {
   if (n < 10) cout << "0";
   cout << n << "\n";
}

int main() {
   long long x, y, a0, a1, n;
   while (cin >> x >> y) {
      if (x == 0 && y == 0) break;
      cin >> a0 >> a1 >> n;
      if (n == 0) { output(a0%100); continue; }
      else if (n == 1) { output(a1%100); continue; }
      Matrix m = power(Matrix(x,y,1,0), n-1);
      output((m.a * a1 + m.b * a0) % 100);
   }
   return 0;
}

Now for a more advanced example, we use matrix exponentiation to determine the number of cycles with a length smaller than k in a given directed graph.

Problem: TourCounting
Source: TopCoder SRM 306 D1-1050
URL: http://www.topcoder.com/stat?c=problem_statement&pm=6386

An elegant way to calculate the number of paths from a source vertex u to a destination vertex v that takes exactly k steps is to use matrix multiplication. First let's define A to be the adjacency matrix with each entry representing the number of edges from u to v. The (u,v) entry of A^t is precisely the answer we are looking for. The way this works is similar to the way that Floyd-Warshall algorithm works. The matrix multiplication considers the connectivity between an intermediate node k and attempts to link u to k and k to v. Due to the multiplicative nature of matrix multiplication as opposed to an additive (for shortest paths) it determines the number of paths based on the multiplication principle from discrete mathematics.

To illustrate this principle, consider a snapshot of the matrix A^t and it's adjacency matrix which is represented as A^1:

If we were to calculate the number of paths that start from vertex 0 to itself (i.e. a cycle) that takes exactly t+1 steps. Note that we can access vertex 0 from either vertex 1 or vertex 2. Since at state t, we can reach vertex 1 from vertex 0 with 2 paths and we can reach vertex 2 from vertex 0 with 3 paths. It stands that we can now reach vertex 0 to itself from 2+3=5 paths. Note that this is a direct consequence of the matrix multiplication process - validate this yourself for the first row and column. If the adjacency matrix was changed such that there are two edges from vertex 1 to vertex 0 like the following:

Then we can reach vertex 0 in a total of 2*2 + 1*3 = 7 ways since we have an option of choosing one of the two edges from vertex 1 to vertex 0. Again, this is consequence of the matrix multiplication algorithm. The argument holds for every other cell in the matrix. So what's the advantage of this method? The prime factor is efficiency - we can calculate matrix powers in logarithmic number of matrix multiplications. We accomplish this through a similar algorithm to fast powering.

So now we know how to calculate the number of paths from a source vertex u to a destination vertex v that takes exactly k steps. The actual problem requires us to compute the number of cycles that are present with less than k steps. So for a given instance of the matrix A^t we want to sum the diagonal to yield the number of paths with exactly t steps. To calculate it for less than k steps we would normally have to compute each power of t and sum all the diagonals up, however this is not fast enough as the number of steps can be as much as 1 million. We can reduce this using a similar method to our matrix powering algorithm. First define a function f:

f(n) = represents the number of ways we can get from all vertices to other vertices using less than n steps (matrix form)

The main optimisation is reducing it from linear (summing all matrices less than n) to logarithmic. The huge hint in logarithmic is dividing the input space by a (usually) constant factor. Note that if we divide the space into two f(n/2) parts we are close to our answer but not quite there. The problem is that upper n/2 part is not symmetrical to the lower n/2 part. But wait! Given the simple mathematical property of powers a^(b+n) = a^b * a^n we can simply just multiply one of the f(n/2) parts by A^(n/2) to yield the correct upper half! With this our recurrence becomes:

Like our fast powering algorithm we need to distinguish between odd and even states. The easy way to make an odd number even is to subtract it by 1. We then just calculate one instance of the odd power and make all successively calls even. This maintains the correctness of the recurrence whilst also maintaining the time complexity.

Then it's simply a matter of implementation which is the easy part. Just note we need to use sensible modulo arithmetic to ensure that we don't overflow any of our computations.

class TourCounting {
public:
   int countTours(vector <string>, int, int);
};

class Matrix {
public:
   vector<vector<long long> > data;
   Matrix() { }
   Matrix(int n, int m) {
      data = vector<vector<long long> >(n, vector<long long>(m, 0));
   }
};

long long MOD;

Matrix operator*(const Matrix& lhs, const Matrix& rhs) {
   Matrix res(lhs.data.size(), rhs.data[0].size());
   for (int i = 0; i < lhs.data.size(); i++) {
      for (int j = 0; j < lhs.data[i].size(); j++) {
         for (int k = 0; k < rhs.data[0].size(); k++) {
            res.data[i][k] = (res.data[i][k] + lhs.data[i][j] * rhs.data[j][k])%MOD;
         }
      }
   }
   return res;
}

Matrix operator+(const Matrix& lhs, const Matrix& rhs) {
   Matrix res(lhs.data.size(),lhs.data[0].size());
   for (int i = 0; i < lhs.data.size(); i++) {
      for (int j = 0; j < lhs.data[0].size(); j++) {
         res.data[i][j] = (lhs.data[i][j] + rhs.data[i][j])%MOD;
      }
   }
   return res;
}

Matrix power(const Matrix& lhs, int P) {
   if (P == 1) return lhs;
   Matrix tmp = power(lhs, P/2);
   if (P % 2 != 0) {
      Matrix n = tmp * tmp;
      Matrix m = n * lhs;
      return m;
   }
   Matrix v = tmp * tmp;
   return v;
}

long long compute(const Matrix& m) {
   long long res = 0;
   for (int i = 0; i < m.data.size(); i++) res = (res + m.data[i][i]) % MOD;
   return res % MOD;
}

Matrix dp[1000001];
int visited[1000001];

Matrix f(const Matrix& ref, int n) {
   if (visited[n]) return dp[n];
   visited[n] = 1;
   if (n == 1) return dp[n] = ref;
   if (n % 2 != 0) {
      // odd
      Matrix v = f(ref, n-1) + power(ref, n);
      return dp[n] = v;
   }
   // even 
   Matrix vE = f(ref, n/2) * power(ref, n/2);
   Matrix vEE = vE + f(ref, n/2);
   return dp[n] = vEE;
}

int TourCounting::countTours(vector <string> g, int k, int m) {
   Matrix mat((int)g.size(),(int)g.size());
   for (int i = 0; i < g.size(); i++)
      for (int j = 0; j < g[i].size(); j++)
         mat.data[i][j] = g[i][j] == 'Y' ? 1 : 0;
   MOD = m;
   long long res = 0;
   Matrix rr = f(mat, k-1);
   res = compute(rr);
   return res;
}