Suppose one is given a (say) real sequence $(a_n)_{n\in\Bbb N}$ such that
- one can explicitly compute it
- one knows it, a priori, satisfies some (unknown) linear recurrence relation of order $\leq N$.
The problem I would like to solve is that of finding an explicit recurrence relation satisfied by $(a_n)_{n\in\Bbb N}$ : how could one find such a recurrence relation explicitly?
Such a sequence might be constructed as follows : put, for all $n\geq 0$, $a_n=\langle a\mid D^n b\rangle$ where $D\in M_N(\Bbb R)$ is some square matrix, and $a,b\in\Bbb R^N$ are column vectors. This sequence satisfies the premise of the question. Indeed, such a sequence will satisfy a recurrence relation of order $\leq N$ by virtue of the fact that the sequence $(D^n)_{n\in\Bbb N}$ satisfies one of order $\leq N$. The order $N$, however, might not be optimal, and $(a_n)$ might already satisfy a recurrence of lower order.
If one can make an educated guess as to the minimal degree of the recurrence relation, say $p\leq N$, then one can find such a recurrence relation explicitly.
Indeed, if $a=(a_n)_{n\in\Bbb N}$ satisfies a linear recurrence relation of order $p$ (which is minimal) given by $$ \forall n\in\Bbb N,\quad a_{n+p}=c_{p-1}a_{n+p-1}+\cdots{}+c_{1}a_{n+1}+c_{0}a_{n} $$ then, by minimality of $p$, the family of sequences $(a,\sigma a, \dots{},\sigma^{p-1} a)$ is free, where $\sigma$ is the shift operator that sends a sequence $(x_n)_{n\in\Bbb N}$ to $(x_{n+1})_{n\in\Bbb N}$. Since a sequence satisfying a recurrence relation of order $p$ is entirely determined by its $p$ first terms, one gets that the matrix $$ \mathbf{A}= \begin{pmatrix} a_{0} & a_{1} & a_{2} & a_{3} &\cdots{} & a_{p-2} & a_{p-1}\\ a_{1} & a_{2} & a_{3} & & &a_{p-1} &a_{p}\\ a_{2} & a_{3} & & & & & \vdots\\ a_{3}&&&&&& a_{2p-5}\\ \vdots & &&&& a_{2p-5}& a_{2p-4}\\ a_{p-2} & a_{p-1} & & & a_{2p-5}& a_{2p-4} & a_{2p-3}\\ a_{p-1} & a_{p} & \cdots{} & a_{2p-5} & a_{2p-4}& a_{2p-3} & a_{2p-2} \end{pmatrix} $$ is invertible, and for real numbers $c_0,c_1,\dots,c_{p-1}$, the sequence $(a_n)$ satisfies the recurrence relation $$ \forall n\in\Bbb N,\quad a_{n+p}=c_{p-1}a_{n+p-1}+\cdots{}+c_{1}a_{n+1}+c_{0}a_{n} $$ if and only if $$ \mathbf{A} \begin{pmatrix} c_{0}\\ c_{1}\\ c_{2}\\ c_{3}\\ \vdots\\ c_{p-2}\\ c_{p-1}\\ \end{pmatrix} = \begin{pmatrix} a_{p}\\ a_{p+1}\\ a_{p+2}\\ a_{p+3}\\ \vdots\\ a_{2p-2}\\ a_{2p-1}\\ \end{pmatrix} =A $$ Both the matrix $\mathbf{A}$ and the column vector $A$ are explicitly computable, and so is the inverse $\mathbf{A}^{-1}$, hence one can explicitly find the coefficients $c_0,\dots{},c_{p-1}$ by inverting the matrix $\mathbf{A}$.
What about the case where one only has an upper bound on the minimal order of a recurrence relation satisfied by $(a_n)$?
Two solutions follow : the first misses an important fact and is far more complicated and time consuming than need be. The last paragraph provides a much simpler and faster solution to the problem of finding the minimal recurrence relation satisfied by a sequence known to satisfy some recurrence relation of order $\leq N$. It is based on the observation that the degree of the minimal relation equals the rank of the matrix $\mathbf{A}$ introduced below.
The simple observation that follows allows us to identify sequences satisfying possibly different recurrence relations
Now suppose $a$ is known to satisfy some recurrence relation of order $p\leq N$, and let $(c_0,\dots,c_{N-1})\in\Bbb C^N$ be a solution to the matrix system (1) below $$ \mathbf{A} \begin{pmatrix} c_{0}\\ c_{1}\\ c_{2}\\ c_{3}\\ \vdots\\ c_{N-2}\\ c_{N-1}\\ \end{pmatrix} = \begin{pmatrix} a_{N}\\ a_{N+1}\\ a_{N+2}\\ a_{N+3}\\ \vdots\\ a_{2N-2}\\ a_{2N-1}\\ \end{pmatrix} =A $$ where $\mathbf{A}=(a_{i+j})_{0\leq i,j\leq N-1}$. Define $\tilde{a}=(\tilde{a}_n)_{n\in\Bbb N}$ to be the sequence satisfying the recurrence relation $$ \forall n\in\Bbb N,\quad \tilde a_{n+N}=c_{N-1}\tilde a_{n+N-1}+\cdots{}+c_{1}\tilde a_{n+1}+c_{0}\tilde a_{n} $$ with initial conditions $\tilde{a}_i=a_i$, for $i=0,\dots{},N-1$. Then by definition of the sequence $(c_0,\dots,c_{N-1})$, one has $\tilde{a}_i=a_i$, for $i=0,\dots{},2N-1$. Thus $a$ and $\tilde a$ satisfy recurrence relations of order $p$ and $N$ respectively, and coincide at the first $2N\geq N+p$ values ($N+p$ is an upper bound on the degree of the lcm of polynomials giving the recurrence relations for $a$ and $\tilde a$). By the previous observation, $\tilde a=a$.
If we furthermore assume that $p$ is the degree of the minimal polynomial of the shift operator at $a$, and write $\mu_a$ for the minimal polynomial of the shift operator at $a$ (so that $\deg(\mu_a)=p$), we get $$ \mu_a\mid X^N-\sum_{k=0}^{N-1}c_kX^k. $$ This shows that upon identifying $\Bbb C^N$ with the affine space of unital degree $N$ polynomials (via $(c_0,\dots,c_{N-1})\longmapsto X^N-\sum_{k=0}^{N-1}c_kX^k$), we get that the set of all $(c_0,\dots,c_{N-1})$ solution to the matrix system (1) is affinely equivalent to the set of all unital degree $N$ multiples of the minimal polynomial of $\sigma$ at $a$.
With enough patience one can then get to the minimal equation satisfied by $a$ (that is, the minimal polynomial of $\sigma$ at $a$) :
This (should) will yield the minimal polynomial at $a$ of $\sigma$, that is, the minimal order recurrence relation satisfied by $a$.
This does enough to convince me that the minimal recurrence relation satisfied by a concrete sequence known to satisfy a recurrence relation of order $\leq N$ is in theory effectively computable.
Note that the degree $p$ of the minimal relation can be extracted from $\mathbf{A}$ as its rank : $\mathrm{rk}(\mathbf{A})=p$ ... And one can then apply the idea exposed in the body of the question, and bypass this answer altogether : extract $\mathbf{A}_p=(a_{i+j})_{0\leq i,j\leq p-1}$, which must be invertible, and find the coefficients as was explained in the question.