Least Square fit of a system of equations

Ok, now it's too long for a comment, so I can try to answer your question right away. Let's write your matrix $DF$ as $D+F$ where $D$ is a matrix with only the last row filled with unknowns and $F$ has as only non-zero entries ones on the first subdiagonal. Then your equation becomes

$$ z^{k+1} = (D + F) z^{k} rightarrow z^{k+1} - F z^{k} = D z^{k} $$ Writing it out element-wise leads to $$ (z^{k+1} - F z^{k})_{ij} = sum_{l=1}^d D_{il} z^{k}_{lj} $$ and due to the fact that only the last row of $D$ is nonzero this becomes (using the Kronecker-Delta) $$ (z^{k+1} - F z^{k})_{ij} = delta_{id} sum_{l=1}^d D_{dl} z^{k}_{lj} $$

You see that a solution can only exists if the first $(d-1)$ rows of $(z^{k+1} - F z^{k})_{ij}$ are zero anyways. Moreover, with your choice of $D$, you have no effect in these dimensions. This means that when you want to find a least-squares solution, you are basically minimizing in the last dimension only. I'm not sure whether this is what you want.

But let's asusme it is intended. Then, for the last row, as $(Fz^{k})_{dj} = 0$, one gets $$ z^{k+1}_{dj} = sum_{l=1}^d D_{dl} z^{k}_{lj} rightarrow sum_{l=1}^d (z^{k})^T_{jl} D_{dl} = z^{k+1}_{dj} $$ which is a linear equation for the unknown coefficients $D_{dl}$ with coefficient matrix $(z^{k})^T$, where $T$ means transpose, and right-hand side $z^{k+1}_{dj}$. And this you can solve with SciPy and whatever.

And again, what you get is a form of least square solution, but quite a special one where you solve in only one dimension exactly (at least if a solution exists), and do nothing in the all the other dimensions. (Put in geometrical terms: there is only one non-vanishing singular value of the matrix $DF$)