  Calculating the spectral function of two dimensional systems is arguably one
of the most pressing challenges in modern computational condensed matter
physics. While efficient techniques are available in lower dimensions, two
dimensional systems present insurmountable hurdles, ranging from the sign
problem in quantum Monte Carlo (MC), to the entanglement area law in tensor
network based methods. We hereby present a variational approach based on a
Chebyshev expansion of the spectral function and a neural network
representation for the wave functions. The Chebyshev moments are obtained by
recursively applying the Hamiltonian and projecting on the space of variational
states using a modified natural gradient descent method. We compare this
approach with a modified approximation of the spectral function which uses a
Krylov subspace constructed from the "Chebyshev wave-functions". We present
results for the one-dimensional and two-dimensional Heisenberg model on the
square lattice, and compare to those obtained by other methods in the
literature.
