# Fast way of computing the action of a matrix power on a vector

For integer $$k>0$$, it is well-known that one can use binary exponentiation to evaluate the matrix power $$\mathbf A^k$$, where $$\mathbf A$$ is an $$n\times n$$ matrix.

However, it is not clear to me if that method can be readily modified to evaluate $$\mathbf A^k \mathbf v$$, where $$\mathbf v$$ is a vector of dimension $$n$$. Of course, one could just form the power as usual, and then apply it to the vector, but in cases like this (e.g. if the matrix is very large), one would wish to avoid forming a whole matrix. (At worst, one could use a few vectors of comparable size.)

Is there a good way to compute the action of a matrix power on a vector for integer exponents? My attempts to search the literature have come up empty, but perhaps I am not using the right keywords.

Solutions, or at least pointers to the literature, would be very much appreciated.

• Is something known about the structure of A (e.g. sparsity, where elements are) or could it be an entirely general matrix? Jan 15, 2022 at 16:31
• What about $A*(A*(A*...(A*v)...))$? No need to compute the matrix power, just apply A to your vector k times Jan 16, 2022 at 2:36
• I would turn @CharlieS comment on its side and suggest using repeated squaring as suggested by other comments and the start or finish with the matrix-vector product to peel off one power of $A$. Jan 16, 2022 at 15:15
• If the matrix is sparse, then you cannot compute $A^k$ efficiently because it will in general be a much denser matrix. So the repeated matrix-vector product suggested by CharlieS's comment is the only reasonable way. Jan 16, 2022 at 18:57
• If the matrix is dense, repeated squaring will take about $\log_2 k$ matrix-matrix products, which cost $n^3$ each. So the cost will be $n^3 \log_2 k$. The repeated matrix-vector product can be computed in $k$ steps, each of which takes $n^2$ operations -- so the cost of that algorithm is $n^2 k$. As a consequence, the repeated matrix-vector product wins if $k < n \log_2 k$, which will universally be the case because you should not be interested in situations where $k>n$. Jan 16, 2022 at 18:59

It's going to be hard to beat the two naive approaches, iterated squaring and repeated application associating from the right $$A(A(A(\cdots Av)))$$. Which one of the two is better depends on the value of $$k$$, sparsity, etc. as noted in the comments.
(1) treat your problem as a general $$f(A)v$$, where $$f$$ is a scalar function. This opens up a few strategies, for instance the diagonalization or Schur--Parlett method: a change of basis to reduce to diagonal or triangular form + solving explicitly the triangular/diagonal problem. This might be convenient for dense $$A$$ and huge $$k$$ since its complexity does not depend on $$k$$ (apart from computing $$k$$th powers of scalars, which can usually be done in constant time). Or the Arnoldi method for sparse $$A$$, which similarly has a very mild dependence on the exponent $$k$$. These are treated on Higham's book Functions of matrices.
(2) approximate $$A$$ with a rank-$$r$$ matrix $$UV^T$$, for instance with randomized methods, and then use $$(UV^T)^{k} = U(V^TU)^{k-1}V^T$$, which leaves you with powers of a $$r\times r$$ matrix rather than a $$n\times n$$ one. The techniques for this random approximation are sketched (pun intended) for instance in the popular review article https://arxiv.org/abs/2002.01387 . Note that the Arnoldi method suggested above is also essentially a rank-$$r$$ approximation computed in a peculiar way.