Complementarity problems and variational inequalities arise in a wide variety of areas, including machine learning, planning, game theory, and physical simulation. In all of these areas, to handle large-scale problem instances, we need fast approximate solution methods. One promising idea is Galerkin approximation, in which we search for the best answer within the span of a given set of basis functions. Bertsekas proposed one possible Galerkin method for variational inequalities. However, this method can exhibit two problems in practice: its approximation error is worse than might be expected based on the ability of the basis to represent the desired solution, and each iteration requires a projection step that is not always easy to implement efficiently. So, in this paper, we present a new Galerkin method with improved behavior: our new error bounds depend directly on the distance from the true solution to the subspace spanned by our basis, and the only projections we require are onto the feasible region or onto the span of our basis.