[proofplan]
We expand the unknown denominator $q$ and compute the coefficients of the formal product $qf$. The low-order coefficients, from degree $0$ through degree $m$, uniquely prescribe the numerator $p$ once $q$ is chosen. The remaining order conditions form $n$ homogeneous linear equations in the $n+1$ denominator coefficients, so a nonzero denominator solution exists. Finally, the condition $q(0)=1$ is exactly the condition that some homogeneous solution has nonzero constant coefficient, because such a solution can be rescaled.
[/proofplan]
[step:Expand the formal product and isolate its coefficients]
Let $(q_0,\dots,q_n) \in \mathbb{C}^{n+1}$ be arbitrary, and define the polynomial $q \in \mathbb{C}[z]$ by
\begin{align*}
q(z)=\sum_{j=0}^{n} q_j z^j.
\end{align*}
For each integer $r \geq 0$, define $c_r \in \mathbb{C}$ to be the coefficient of $z^r$ in the formal product $qf$. Since $q$ has degree at most $n$, the Cauchy product formula gives
\begin{align*}
c_r=\sum_{j=0}^{\min(n,r)} q_j a_{r-j}.
\end{align*}
Thus
\begin{align*}
q(z)f(z)=\sum_{r=0}^{\infty} c_r z^r
\end{align*}
as an identity in $\mathbb{C}[[z]]$.
[/step]
[step:Choose the numerator to cancel the coefficients through degree $m$]
For the chosen denominator coefficients $(q_0,\dots,q_n)$, define $p \in \mathbb{C}[z]$ by
\begin{align*}
p(z)=\sum_{r=0}^{m} c_r z^r.
\end{align*}
Then $\deg p \leq m$, unless $p=0$, in which case the degree bound is vacuous. By construction, the coefficients of $z^0,z^1,\dots,z^m$ in $qf-p$ vanish. Therefore $qf-p$ has formal order at least $m+1$ exactly when the coefficients $c_{m+1},\dots,c_{m+n}$ also vanish.
Using the formula for $c_r$, these remaining conditions are precisely
\begin{align*}
\sum_{j=0}^{\min(n,r)} q_j a_{r-j}=0 \qquad \text{for } r=m+1,\dots,m+n.
\end{align*}
[guided]
Once $q$ is fixed, there is no freedom left in the numerator if we want cancellation through degree $m$. Indeed, the coefficient of $z^r$ in $qf$ is the scalar
\begin{align*}
c_r=\sum_{j=0}^{\min(n,r)} q_j a_{r-j}.
\end{align*}
To make the coefficient of $z^r$ in $qf-p$ vanish for $0 \leq r \leq m$, the coefficient of $z^r$ in $p$ must be exactly $c_r$. This forces the definition
\begin{align*}
p(z)=\sum_{r=0}^{m} c_r z^r.
\end{align*}
This polynomial has degree at most $m$ because it contains no powers $z^r$ with $r>m$.
After this choice, the coefficients through degree $m$ are already cancelled. The desired formal order condition
\begin{align*}
q(z)f(z)-p(z)=O(z^{m+n+1})
\end{align*}
means that the coefficients of $z^0,z^1,\dots,z^{m+n}$ in $qf-p$ all vanish. Since the first $m+1$ of these have been killed by the definition of $p$, the only remaining conditions are
\begin{align*}
c_r=0 \qquad \text{for } r=m+1,\dots,m+n.
\end{align*}
Substituting the Cauchy product expression for $c_r$ gives the homogeneous system
\begin{align*}
\sum_{j=0}^{\min(n,r)} q_j a_{r-j}=0 \qquad \text{for } r=m+1,\dots,m+n.
\end{align*}
It is homogeneous because every equation is linear in the unknowns $q_0,\dots,q_n$ and has zero right-hand side.
[/guided]
[/step]
[step:Find a nonzero denominator solution to the homogeneous system]
The displayed system has $n$ linear equations over $\mathbb{C}$ in the $n+1$ unknowns $q_0,\dots,q_n$. A homogeneous linear system with more unknowns than equations has a nonzero solution: after row reduction, at most $n$ pivot variables can occur, so at least one variable is free, and assigning that free variable the value $1$ gives a nonzero solution.
Choose such a nonzero solution $(q_0,\dots,q_n)$. Define $q \in \mathbb{C}[z]$ and $p \in \mathbb{C}[z]$ from this solution by the formulas above. Since $(q_0,\dots,q_n)\neq 0$, the polynomial $q$ is not the zero polynomial. Hence the pair $(p,q)$ is not both zero, $\deg q \leq n$, $\deg p \leq m$, and
\begin{align*}
q(z)f(z)-p(z)=O(z^{m+n+1}).
\end{align*}
[/step]
[step:Characterize when the homogeneous solution can be normalised]
Suppose first that the homogeneous system has a solution $(q_0,\dots,q_n)$ with $q_0 \neq 0$. Define $\lambda \in \mathbb{C}$ by
\begin{align*}
\lambda=q_0^{-1}.
\end{align*}
Because the system is homogeneous, $(\lambda q_0,\dots,\lambda q_n)$ is again a solution. The corresponding denominator $\tilde q \in \mathbb{C}[z]$ satisfies
\begin{align*}
\tilde q(0)=\lambda q_0=1.
\end{align*}
Using the numerator constructed from $\tilde q$, we obtain a normalised Padé approximant satisfying the required order condition.
Conversely, suppose a normalised Padé approximant exists. Then there are polynomials $p,q \in \mathbb{C}[z]$ with $\deg p \leq m$, $\deg q \leq n$, $q(0)=1$, and
\begin{align*}
q(z)f(z)-p(z)=O(z^{m+n+1}).
\end{align*}
Writing $q(z)=\sum_{j=0}^{n} q_j z^j$, the argument above shows that $(q_0,\dots,q_n)$ satisfies the homogeneous system. Since $q_0=q(0)=1$, this solution has $q_0\neq 0$.
Therefore a normalised Padé approximant with $q(0)=1$ exists exactly when the finite homogeneous system has a solution with nonzero constant denominator coefficient.
[/step]
