For highly complex urban systems, we may not be able to know which variables will be present, and we may not know how to explicitly describe the underlying dynamics, or we may not even know the governing equations of city evolution. This chapter discusses how to discover the underlying governing equations of city dynamics directly from experimental or observational data. The first approach employs sparse regression methods that take advantage of the sparsity of nonlinear terms constituting the underlying governing equations to learn parsimonious models of city evolution. Sparse identification of nonlinear dynamics can be extended to model learning in a noisy environment with stochastic perturbations and controls. It can also be combined with neural networks to find the right coordinates of the dynamical systems. The second approach is the Gaussian process with mixture distributions. It is the statistical learning of hidden dynamics with complete or incomplete observations in the phase space. Neural networks constitute the third approach in which black box methods are employed to uncover nonlinear mappings between system inputs and outputs. Using urban population dynamics and spatial interaction as examples, the chapter also examines ways different approaches can be integrated for effective and efficient data-driven learning. It concludes with a discussion on intelligibility of data-driven learning.