Abstract The peridynamic theory provides an appropriate description of the deformation of a continuous body involving discontinuities or other singularities, which cannot be described properly by the classical theory of solid mechanics. However, the operator in the peridynamic theory is nonlocal, so the resulting numerical methods generate dense or full coefficient matrices which require O(N2) of memory where N is the number of unknowns in the discretized system. Gaussian types of direct solvers, which were traditionally used to solve these problems, require O(N3) of operations. Furthermore, due to the singularity of the kernel in the peridynamic model, the evaluation and assembly of the coefficient matrix can be very expensive. Numerous numerical experiments have shown that in many practical simulations the evaluation and assembly of the coefficient matrix often constitute the main computational cost! The significantly increased computational work and memory requirement of the peridynamic model over those for the classical partial differential equation models severely limit their applications, especially in multiple space dimensions. We develop a fast and faithful collocation method for a two-dimensional nonlocal diffusion model, which can be viewed as a scalar-valued version of a peridynamic model, without using any lossy compression, but rather, by exploiting the structure of the coefficient matrix. The new method reduces the evaluation and assembly of the coefficient matrix by O(N), reduces the computational work from O(N3) required by the traditional methods to O(Nlog2N) and the memory requirement from O(N2) to O(N). Numerical results are presented to show the utility of the fast method.