Radial functions are a powerful tool in many areas of multidimensional approximation, especially when dealing with scattered data. We present a fast approximate algorithm for the evaluation of linear combinations of radial functions on the sphere S2 . The approach is based on a particular rank approximation of the corresponding Gram matrix and fast algorithms for spherical Fourier transforms. The proposed method takes O(L) arithmetic operations for L arbitrarily distributed nodes on the sphere. In contrast to other methods, we do not require the nodes to be sorted or pre-processed in any way, thus the precomputation effort only depends on the particular radial function and the desired accuracy. We establish explicit error bounds for a range of radial functions and provide numerical examples covering approximation quality, speed measurements, and a comparison of our particular matrix approximation with a truncated singular value decomposition. AMS Subject Classification: 65T50, 65F30, ...