We explore the connection between fractional order partial differential
equations in two or more spatial dimensions with boundary integral operators to
develop techniques that enable one to efficiently tackle the integral
fractional Laplacian. In particular, we develop techniques for the treatment of
the dense stiffness matrix including the computation of the entries, the
efficient assembly and storage of a sparse approximation and the efficient
solution of the resulting equations. The main idea consists of generalising
proven techniques for the treatment of boundary integral equations to general
fractional orders. Importantly, the approximation does not make any strong
assumptions on the shape of the underlying domain and does not rely on any
special structure of the matrix that could be exploited by fast transforms. We
demonstrate the flexibility and performance of this approach in a couple of
two-dimensional numerical examples