Let $A$ be an $n\times n$ integral matrix and $\xi$ an integral vector of dimension $n$. Let $$W_i=[\xi,A\xi,\ldots,A^{i-1}\xi]$$ for each positive $i$. Let $d_{i,j}$ (where $1\le j\le i$) denote the $j$-th invariant factor of $W_i$. It seems that $d_{i,j}$ only depends on $j$, that is $d_{j,j}=d_{j+1,j}=\cdots=d_{n,j}=d_{n+1,j}=\cdots$. By Cayley-Hamilton Theorem, we know that $W_i$ is equivalent (in $\mathbb{Z}$) to $[W_n,0_{n\times (i-n)}]$ for $i>n$. Thus we only need to consider the case that $i\le n$.
For example, let \begin{equation*} A=\left(\begin{matrix} 2&2&3&3\\ 1&3&0&0\\ 3&3&0&3\\ 0&0&3&2\\ \end{matrix} \right),\quad \xi=\left(\begin{matrix} 1\\ 2\\ 3\\ 0\\ \end{matrix} \right) \end{equation*} Then \begin{equation*} W_1=\left(\begin{matrix} 1\\ 2\\ 3\\ 0\\ \end{matrix} \right), W_2=\left(\begin{matrix} 1&15\\ 2&7\\ 3&9\\ 0&9\\ \end{matrix} \right),W_3=\left(\begin{matrix} 1&15&98\\ 2&7&36\\ 3&9&93\\ 0&9&45\\ \end{matrix} \right), W_4=\left(\begin{matrix} 1&15&98&682\\ 2&7&36&206\\ 3&9&93&537\\ 0&9&45&369\\ \end{matrix} \right) \end{equation*} The Smith normal forms of these matrices are \begin{equation*} \left(\begin{matrix} 1\\ 0\\ 0\\ 0\\ \end{matrix} \right), \left(\begin{matrix} 1&0\\ 0&1\\ 0&0\\ 0&0\\ \end{matrix} \right),\left(\begin{matrix} 1&0&0\\ 0&1&0\\ 0&0&3\\ 0&0&0\\ \end{matrix} \right), \left(\begin{matrix} 1&0&0&0\\ 0&1&0&0\\ 0&0&3&0\\ 0&0&0&9090\\ \end{matrix} \right) \end{equation*} For a fixed $j$, why these $j$-th invariant factors of $W_i$'s are all the same?