  We establish conditions for which graph Laplacians
$\Delta_{\lambda,\epsilon}$ on compact, boundaryless, smooth submanifolds
$\mathcal{M}$ of Euclidean space are semiclassical pseudodifferential operators
($\Psi$DOs): essentially, that the graph Laplacian's kernel bandwidth
($\textit{bias term}$) $\sqrt{\epsilon}$ decays faster than the semiclassical
parameter $h$, $\textit{i.e.}$, $h \gg \sqrt{\epsilon}$ and we compute the
symbol. Coupling this with Egorov's theorem and coherent states $\psi_h$
localized at $(x_0, \xi_0) \in T^*\mathcal{M}$, we show that with
$U_{\lambda,\epsilon}^t := e^{-i t \sqrt{\Delta}_{\lambda,\epsilon}}$
spectrally defined, the (co-)geodesic flow $\Gamma^t$ on $T^*\mathcal{M}$ is
approximated by $\langle U_{\lambda,\epsilon}^{-t} \operatorname{Op}_h(a)
U_{\lambda,\epsilon}^t \psi_h, \psi_h \rangle = a \circ \Gamma^t(x_0, \xi_0) +
O(h)$. Then, we turn to the discrete setting: for $\Delta_{\lambda,\epsilon,N}$
a normalized graph Laplacian defined on a set of $N$ points $x_1, \ldots, x_N$
sampled $\textit{i.i.d.}$ from a probability distribution with smooth density,
we establish Bernstein-type lower bounds on the probability that
$||U_{\lambda,\epsilon,N}^t[u] - U_{\lambda,\epsilon}^t[u]||_{L^{\infty}} \leq
\delta$ with $U_{\lambda,\epsilon,N}^t := e^{-i t
\sqrt{\Delta}_{\lambda,\epsilon,N}}$. We apply this to coherent states to show
that the geodesic flow on $\mathcal{M}$ can be approximated by matrix dynamics
on the discrete sample set, namely that $\textit{with high probability}$,
$c_{t,N}^{-1} \sum_{j=1}^N |U_{\lambda,\epsilon,N}^t[\psi_h](x_j)|^2 u(x_j) =
u(x_t) + O(h)$ for $c_{t,N} := \sum_{j=1}^N
|U_{\lambda,\epsilon,N}^t[\psi_h](x_j)|^2$ and $x_t$ the projection of
$\Gamma^t(x_0, \xi_0)$ onto $\mathcal{M}$.
