  The inverse design of optical metasurfaces is a rapidly emerging field that
has already shown great promise in miniaturizing conventional optics as well as
developing completely new optical functionalities. Such a design process relies
on many forward simulations of a device's optical response in order to optimize
its performance. We present a data-driven forward simulation framework for the
inverse design of metasurfaces that is more accurate than methods based on the
local phase approximation, a factor of $10^4$ times faster and requires $15$
times less memory than mesh based solvers, and is not constrained to spheroidal
scatterer geometries. We explore the scattered electromagnetic field
distribution from wavelength scale cylindrical pillars, obtaining
low-dimensional representations of our data via the singular value
decomposition. We create a differentiable model fiting the input geometries and
configurations of our metasurface scatterers to the low-dimensional
representation of the output field. To validate our model, we inverse design
two optical elements: a wavelength multiplexed element that focuses light for
$\lambda=633$nm and produces an annular beam at $\lambda=400$nm and an extended
depth of focus lens.
