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;
}

TC SRM446 D1 500 (AntOnGraph)

The problem statement can be accessed via:
http://www.topcoder.com/stat?c=problem_statement&pm=10516

The problem depicts an ant which travels around a weighted directed graph. The objective of this problem is to determine the maximum score in which the ant can have given that the score is defined by the sum of all the edge weights it travels. The ant is given at most timeLimit seconds to travel (it can use less than this) and each second is defined by a mandatory set of stepsPerSecond adjacent vertex traversals within the graph. An additional constraint is that the ant must end at vertex 1 at the end (the ant starts at vertex 0), it is also possible that the ant can't reach the vertex at the right interval in time or that the vertex is inaccessible.

Given the constraints it's rather obvious we can't simulate it by time (as it'll easily exceed the time limit). There are several clues in this problem which leads to an elegant solution - firstly, we must move exactly stepsPerSecond traversals each second. What does this mean? It means that we can compute the cost of moving exactly this many traversals per second. Let's assume timeLimit = 1 second, then the answer is simply the greatest value (which could be negative) from vertex 0 to vertex 1 using exactly stepsPerSecond vertex traversals.

How do we compute this? We can think of it in a dynamic programming manner. Let's define M(x,i,j) to be the maximum value we can achieve using exactly x steps from vertex i to vertex j. If this is unreachable then we define M(x,i,j) to be negative infinity. So now we define the relationship between step x and step x-1:

M(x,i,j) = max { M(x-1,i,k) + M(1,k,j) }       
     iff (i,k) and (k,j) are edges (not negative infinity)

We can iteratively compute this from 0 to stepsPerSecond traversals and the maximum score we can get is simply M(stepsPerSecond, 0, 1). Now we need to somehow incorporate this into the case where timeLimit > 1. We can now try and relate timeLimit and timeLimit-1 seconds. It can be seen that the optimal/maximal moves does not change inside each individual timeLimit (the stepsPerSecond matrix). To see this, let's say we are on vertex x and we just finished using the timeLimit-1 interval. We now wish to go from vertex x to vertex y on the next interval. Then the optimal score is simply M(stepsPerSecond, x, y) + currentScore. We simply reused the matrix we computed since it forms a "basic indivisible unit" due to the constraint that each second must have exactly stepsPerSecond vertex traversals.

We now attempt to define a recurrence between seconds (and not specific vertex traversals) using the matrix M as a basic unit:

N(t,i,j) = max { M(t-1,i,k) + M(1,k,j) }
     iff (i,k) and (k,j) are edges

There is a problem though, since the maximum timeLimit is 1 billion - this recurrence is too slow. To optimise it we need to make one observation. With basic expansion:

N(t,i,j) = N(1,i,k) + N(1,k,l) + N(1,l,m) + ... + N(1,o,j)   (t additions)

Then we can begin collapsing:

N(t,i,j) = N(2,i,l) + ... + N(2,n,j)   (t/2 additions)

etc.

As you can see we can recursively define N(t,i,j) in terms of additions of smaller values. Note that we can define it much like how you do exponentation squaring:

N(t,i,j) = N(t/2,i,k) + N(t/2,k,j)            iff (t is even)
         = N(t/2,i,k) + N(t/2,k,m) + N(1,m,j) otherwise

With this we can compute the maximal score for exactly "t" seconds. We need to make a minor modification to allow us to compute the correct score in case we don't need to use all of the allotted "t" seconds. Basically we need to add a dummy edge of 0 score going from the last vertex to itself (self loop) so it's a possible decision to terminate early. Given the recurrence and the self-loop modification on the finishing vertex we can simply implement:

class AntOnGraph {
public:
 string maximumBonus(vector <string>, vector <string>, 
   vector <string>, int, int);
};

long long tab[51][51];

#define INF (1LL<<50)

vector<vector<long long> > mat(const vector<vector<long long> 
    >& lhs, const vector<vector<long long> >& rhs) {
 vector<vector<long long> > res = vector<vector<long long> 
     >(lhs.size(), vector<long long>(lhs.size(), -INF));
 for (int i = 0; i < lhs.size(); i++) {
  for (int j = 0; j < lhs[0].size(); j++) {
   for (int k = 0; k < lhs.size(); k++) {
    if (lhs[i][k] > -INF && rhs[k][j] > -INF)
     res[i][j] >?= (lhs[i][k] + rhs[k][j]);
   }
  }
 }
 return res;
}

vector<vector<long long> > pow(const vector<vector<long long> 
    >& vec, long long n) {
 vector<vector<long long> > res;
 if (n == 1) {
  return vec;
 }
 vector<vector<long long> > base = pow(vec, n/2);
 res = mat(base, base);
 if (n % 2 != 0) {
  return mat(res, vec);
 }
 return res;
}

string AntOnGraph::maximumBonus(vector <string> p0, vector <string> p1, 
    vector <string> p2, int stepsPerSecond, int timeLimit) {
 string res;
 vector<vector<long long> > vec = vector<vector<long long> >(p0.size()
   ,vector<long long>(p0.size(),0));
 for (int i = 0; i < p0.size(); i++) {
  for (int j = 0; j < p0[i].size(); j++) {
   vec[i][j] = ((p0[i][j] - '0') * 100 + (p1[i][j] - '0') * 10 + (p2[i][j] 
       - '0')) - 500;
   if (vec[i][j] == -500) vec[i][j] = -INF;
  }
 }
 vector<vector<long long> > M = pow(vec, stepsPerSecond);
 if (M[1][1] < 0) M[1][1] = 0;
 vector<vector<long long> > M2 = pow(M, timeLimit);
 ostringstream oss; oss << M2[0][1];
 res = oss.str();
 return M2[0][1] == -INF ? "IMPOSSIBLE" : res;
}