Generalizing Knuth's $O(\log_2 n)$ Fibonacci algorithm to linear homogenous recurrences Knuth has a neat algorithm that uses matrix exponentiation to compute the $n$th Fibonacci number in $O(\log_2 n)$-time 1. However, there doesn't seem to be a lot of resources on generalizing his idea to other linear recurrences. Is it a generalizable idea, or is it just a special technique only applicable to the Fibonacci recurrence?

(I'm not asking for a closed-form solution, but an efficient method similar to Knuth's technique for Fibonacci numbers).

As a related question, can we speed up other recurrences by rewriting them as linear transformations, and how far can we go in applying this technique?
Yes.  This generalizes to any linear recurrence.  Suppose we have the linear recurrence

$$x_{n+1} = a_0 x_n + a_1 x_{n-1} + \dots + a_k x_{n-k}.$$

Define the column vector $v_n = (x_n,x_{n-1},\dots,x_{n-k})^\top$.  Then we have the equation

$$v_{n+1} = M v_n$$

where the matrix $M$ is given by

$$M = \begin{pmatrix}
a_0 &amp;a_1 &amp;a_2 &amp;\dots &amp;a_k\\
1 &amp;0 &amp;0 &amp;\dots &amp;0\\
0 &amp;1 &amp;0 &amp;\dots &amp;0\\
&amp;&amp;&amp;\vdots\\
\end{pmatrix}.$$

It follows that

$$v_n = M^n v_0.$$

Therefore you can compute $v_n$ using $O(\log n)$ matrix multiplications, using a fast exponentiation algorithm for computing $M^n$.  This gives an efficient way to compute $x_n$: first compute $v_n$ using matrix exponentiation, then extract the first coefficient of $v_n$ to get $x_n$.

In fact, we can take this farther.  Suppose $M$ diagonalizes, say $M=PDP^{-1}$ where $D$ is a diagonal matrix.   Then

$$v_n = M^n v_0 = (PDP^{-1})^n v_0 = P D^n P^{-1} v_0,$$

so we don't need to perform matrix exponentiation; it suffices to exponentiate each of the diagonal entries of $D$.  In fact, expanding terms, you can get an algebraic closed-form expression for $x_n$.  If you apply this to the Fibonacci sequence, you get the closed-form expression in https://en.wikipedia.org/wiki/Fibonacci_number#Closed-form_expression (in terms of powers of $(1\pm \sqrt{5})/2$, because those are the elements of $D$); for a general sequence, you get a closed-form expression of the form

$$x_n = c_1 \lambda_1^n + \dots + c_k \lambda_k^n,$$

where $\lambda_1,\dots,\lambda_k$ are the diagonal entries of $D$ (i.e., the eigenvalues of $M$).  This also gives a potentially more efficient algorithm for computing $x_n$: instead of $O(\log n)$ multiplications of $k\times k$ matrices, we get $O(k \log n)$ operations on numbers.

If $M$ does not diagonalize, you may still be able to use the Jordan normal form of $M$.

See also https://en.wikipedia.org/wiki/Recurrence_relation#Solving_homogeneous_linear_recurrence_relations_with_constant_coefficients, which does something similar using the characteristic polynomial of the linear recurrence (which I think coincides with the characteristic polynomial of $M$, or something like that).