For some irrational $\theta\in(0,1)$, let $i_k=1$ if $\lfloor (k+1)\theta\rfloor > \lfloor k\theta\rfloor$ and $i_k=0$ otherwise. In other words, $(i_k)_{k\in\mathbb N}$ is a Sturmian sequence. For example, if $\theta=1/\sqrt 2$ we get $(i_k) = (0, 1, 1, 0, 1, 1, 0, 1, 1, 1, \dots)$.
Now, let $M_0, M_1\in\mathbb R^{k\times k}$ be matrices, and let $$M_n = \prod_{k=0}^n M_{i_k}.$$ Is there a way to compute the eigenvalues of $M_n$ symbolically?
Of course, I can do it approximate $M_n$ by doing a bunch of multiplications, but I'm looking for something similar to how the eigenvalues of $M_1^n$ are just the eigenvalues of $M_1$ to the $n$th power.
From the ergodic properties of the Sturmian sequence, I might approximate $M_n$ by $ (\theta M_1 + (1-\theta) M_0)^n$, but that approach only seems to give a lower bound on the eigenvalues.
If no symbolic approach is known, an algorithm faster than linear time would also be useful. Like the $O(\log n)$ algorithm for matrix exponentiation.
Here's a fast multiplication algorithm for (finite) Sturmian words: Consider a Sturmian word with more $1$s than $0$s. If you replace all instances of $10$ by $2$, you get a new Sturmian word made only of $1$s and $2$s, possibly prepended by a $0$. Now if you let $M_2 := M_1 M_0$, then the matrix product generated by the new word is equal to the matrix product generated by the old word. The $2$-density of the new word is easily computable from the $1$-density of the old word, and the new word is in general significantly shorter.
Instead, if the old Sturmian word had more $0$s than $1$s, then replacing all $10$s by $2$ yields a Sturmian word of $0$s and $2$s, possibly followed by a $1$. The same principle as above applies.
This way you can recursively reduce your long word to a very short word with very few iteration steps on average. This approach can still have very bad worst-case performance, but it is possible to make the algorithm perform well even worst-case by doing the following instead: If your Sturmian word has between $m$ and $m+1$ times as many $1$s as $0$s, you can replace instances of $1\ldots 10$ by $2$ (always replace exactly $m$ ones), and let $M_2 = M_1^mM_0$.
Note that this $m$ is the first entry in the continued fraction expansion of $\theta/(1-\theta)$, and the next few values of $n$ will be the subsequent entries in this expansion.
This algorithm will let you multiply $n$ matrices in $O(\log(n))$ time, assuming constant-time addition and multiplication of integers.