Consider the vector $y_i$ for a given time $t_i$. We are looking for the linear operator $\mathcal{B}$ that maps $y_i$ to $y_{i+1}$ by noting that
\begin{align} \dot{y}&= \mathcal{A} y\ \frac{y_{i+1}-y_i}{\Delta t}&=\mathcal{A} y_i\ \Rightarrow y_{i+1}&=(\mathcal{A}\Delta t+I) y_i\ y_{i+1}&=\mathcal{B} y_i \end{align}
For a set of vectors $y_i$, we can construct the matrices $\mathcal{Y}$ and $\mathcal{Y}'$. \begin{align} \mathcal{Y}= \begin{bmatrix} \vert & \vert & & \vert \ y_1 & y_2 & \dots & y_{n-1} \ \vert & \vert & & \vert \end{bmatrix} \ \mathcal{Y}'= \begin{bmatrix} \vert & \vert & & \vert \ y_2 & y_3 & \dots & y_n \ \vert & \vert & & \vert \end{bmatrix} \end{align} We can therefore write: \begin{align} \mathcal{Y}'=\mathcal{B}\mathcal{Y} \end{align}
The optimization then involves solving:
\begin{align} \Phi(\mathcal{B})=\min_{\mathcal{B}}\vert\vert\mathcal{B}\mathcal{Y}-\mathcal{Y}'\vert\vert_F \end{align}
Here, $\vert\vert\cdot\vert\vert_F$ denotes the <PERSON> norm.
We can, as you said, minimize the objective function via an SVD, by seeing that
\begin{align} \mathcal{B}\mathcal{Y}&=\mathcal{Y}'\ \mathcal{B}\;{U_k^{\,}S_k^{\,}V_k^{\,}}^T&=\mathcal{Y}'\ \Rightarrow \mathcal{B}&=\mathcal{Y}'{V_k^{\,}S_k^{-1}U_k^{\,}}^T \end{align}
Here, $k$ denotes the number of singular vectors/values used for the construction of the inverse operator. You could e.g. use Matlabs svds function for the subset of the first $k$ singular vectors/values.
Alternatively, we can iterate the solution via gradient descent. 804
The derivative of the with respect to $\mathcal{B}$ of the objective function is:
\begin{align} \frac{\partial \Phi}{\partial \mathcal{B}}=2(\mathcal{B}\mathcal{Y}-\mathcal{Y}')\mathcal{Y}^T \end{align}
Here, we can preallocate the products $\mathcal{Y}'\mathcal{Y}^T$ and $\mathcal{Y}\mathcal{Y}^T$. This may be faster than waiting for the SVD.
We thus iterate until we hit some stopping criterion:
\begin{align} \mathcal{B}{k+1} = \mathcal{B}{k}-\alpha(\mathcal{B}_{k}\mathcal{Y}-\mathcal{Y}')\mathcal{Y}^T \end{align}
With $\alpha$ being the step size.
Finally, we can get $\mathcal{A}$ by
\begin{align} \mathcal{A} = \frac{\mathcal{B}-I}{\Delta t} \end{align}
To prevent overfitting due to noise, we can limit the size of $\mathcal{B}$ via adding a regularizer to the objective function. Often times this is implemented via a proximal mapping onto a ball in a given metric. Then, each iteration is projected back into a feasible set. This is a form of constrained optimization.
A good book to read about this is:
Dynamic Mode Decomposition: Data-Driven Modeling of Complex Systems
Author(s): <PERSON> , <PERSON> , <PERSON> and <PERSON> 677