  In this article, we propose two numerical methods, the Gaussian Process (GP)
method and the Fourier Features (FF) algorithm, to solve mean field games
(MFGs). The GP algorithm approximates the solution of a MFG with maximum a
posteriori probability estimators of GPs conditioned on the partial
differential equation (PDE) system of the MFG at a finite number of sample
points. The main bottleneck of the GP method is to compute the inverse of a
square gram matrix, whose size is proportional to the number of sample points.
To improve the performance, we introduce the FF method, whose insight comes
from the recent trend of approximating positive definite kernels with random
Fourier features. The FF algorithm seeks approximated solutions in the space
generated by sampled Fourier features. In the FF method, the size of the matrix
to be inverted depends only on the number of Fourier features selected, which
is much less than the size of sample points. Hence, the FF method reduces the
precomputation time, saves the memory, and achieves comparable accuracy to the
GP method. We give the existence and the convergence proofs for both
algorithms. The convergence argument of the GP method does not depend on any
monotonicity condition, which suggests the potential applications of the GP
method to solve MFGs with non-monotone couplings in future work. We show the
efficacy of our algorithms through experiments on a stationary MFG with a
non-local coupling and on a time-dependent planning problem. We believe that
the FF method can also serve as an alternative algorithm to solve general PDEs.
