How to solve for discrete state space matrices given input and output

Question

I have a set of time-series data that consists of inputs $u_k$ where $ u \in R $ and $k = 1 ... T$, and outputs $ y_k $ where $ y \in R^2 $ and also $k = 1 ... T$, from a given system. I believe this system can be modeled in discrete canonical state space form as $$ x_{k+1} = Fx_k + Gu_k $$ $$ y_{k+1} = Cx_{k+1} $$ In this case, $ y_{k} = x_{k} $ so the form becomes $$ y_{k+1} = Fy_k + Gu_k $$ Given all of the $ y_k $ and $ u_k $ I am fairly sure that I should be able to solve for $ F $ and $ G $, but how do I actually do this? Some work on paper got me nowhere and I can't seem to find anything on the internet.

Answer 1

The answer to this turned out to be to frame the problem as a least-squares problem. Specifically I found & used State Space Model Identification by Least Squares. This allowed me to solve for the rows of $F$ and $G$ one-by-one.