# Exact powers of an integer matrix

Assume that I have a large, non-symmetric matrix $\bf{A}$ of zeros and ones. If I want the exact diagonal elements of $\bf{A}^{k}$ (not the trace, but the elements themselves) where $k$ is some large whole number, I can

1. Multiply the matrices out using with exact integers, possibly speeding the computation by using some Addition-chain exponentiation.

2. Diagonalize the matrix as $\bf A = V \Gamma V^{-1}$ and compute $\bf A^{k} = V \Gamma^{k} V^{-1}$

The first method is slow but accurate, while the second method can introduce numerical errors (since the exact answer is required). Is there some decomposition I am missing that does what I'm looking for?

-
What is larger, $k$ or the size of the matrix? If it is $k$, you could use the Cayley-Hamilton theorem for computing the diagonal entries of every power of $A$ once you know those of the first $n$ powers of $A$ (where $n$ is the size of $A$). – darij grinberg Feb 18 '11 at 18:47
Both cases are needed - though I never thought that the Cayley-Hamilton theorem could be used like that, it's brilliant! It doesn't however help when $k$ is less then the size of the matrix. – Hooked Feb 18 '11 at 19:11

(edit : this is for large k, so it may not be useful for k comparable to n)

A faster way to compute high powers of a matrix is by computing in the $\mathbb{Z}$-algebra $\mathbb{Z}[M]$. The space $(I_n,M,M^2, \ldots)$ is a subspace of $M_n(\mathbb{Z})$, and in fact the Cayley-Hamilton theorem tells you that it has dimension at most $n$. For example, sometimes, the $n$ coefficients on the diagonal of the matrix (or the the coefficients on the first row, or column, etc) are enough to know what the rest of the coefficients are.

If $P$ is the minimal polynomial of $M$ over $\mathbb{Z}$ and has degree $d = d_P$, the basis the best suited in order to compute products, however, is the basis $(1,M,\ldots M^{d-1})$ : we want to use the isomorphism $\mathbb{Z}[M] = \mathbb{Z}[X]/P$. Doing one multiplication in there requires $O(d^2)$ operations :

Keep a table of the coordinates of $(1,M,\ldots M^{2d-2})$ in the basis $(1,M,\ldots M^{d-1})$, then compute the coefficients $c_k$ in $$(\Sigma_{i=0}^{d-1} a_i M^i)(\Sigma_{j=0}^{d-1} b_j M^j) = \Sigma_{k=0}^{2d-2} (c_k = \Sigma_{i+j=k} a_i b_j) M^k$$ Then use that table to convert this sum into a vector $\Sigma_{k=0}^{d-1} d_k M^k$.

You can further optimize this algorithm if the polynomial $P$ is reducible. If $P = QR$ with $Q$ and $R$ having no factor in common, then $\mathbb{Z}[M] = \mathbb{Z}[X]/P \simeq \mathbb{Z}[X]/Q \times \mathbb{Z}[X]/R$. The isomorphism here sends $X$ onto the pair $X,X$. you will need to compute the inverse of the isomorphism, but once you do that, you gain in complexity because $d_Q^2 + d_R^2 < (d_Q + d_R)^2 = d_P^2$.

You can also improve if you have $P = Q^k$ for some $k$, because then, $\mathbb{Z}[X]/P \simeq \mathbb{Z}[X,Y]/(Q(X),Y^k)$, where the isomorphism sends $X$ onto $X+Y$. Once again, the computations there are shorter since you have $k$ polynomials of degree $d/k$ to compute, it will take $O(k(d/k)^2 = O(d^2/k)$.

You have to compute $P$, its factorization, and the isomorphisms only once. Once you have done all this preliminary work you get a very good algorithm. If the matrix is diagonalizable (in $M_n(\mathbb{Z})$), you will get that $P$ factors completely into distincts factors of degree 1, and this will give you an algorithm equivalent to the one you gave.

Even if you don't want to factorize $P$, simply working in $\mathbb{Z}[X]/P$ gives a good optimisation. You can compute $P$ by taking some element $x$ in $\mathbb{Z}^n$, and checking when $M^k x$ becomes linearly dependant on $1, x, \ldots x^{k-1}$. If you are lucky you get some dependancy for $k<n$, which means you have found a factor of $P$ of degree $k<n$. If you know of any subspace of $\mathbb{Z}^n$ stable by $M$, be sure to make use of it, picking $x$ in it will give you such a factor of $P$.

-
this is great, I'll have to read into this quite a bit more. Can you suggest an entry level reference (that hopefully provides an example)? – Hooked Feb 18 '11 at 21:02