We show how the basic idea of parabolic Jacobi relaxation can be modified to obtain a new class of hyperbolic relaxation schemes that are suitable for the solution of elliptic equations. Some of the analytic and numerical properties of hyperbolic relaxation are examined. We describe its implementation as a first order system in a pseudospectral evolution code, demonstrating that certain elliptic equations can be solved within a framework for hyperbolic evolution systems. Applications include various initial data problems in numerical general relativity. In particular we generate initial data for the evolution of a massless scalar field, a single neutron star, and binary neutron star systems.