Kinematic finite-fault source inversions aim at resolving the spatio-temporal evolution of slip on a fault given ground motion recorded on the Earth’s surface. This type of inverse problem is inherently ill-posed due to two main factors. First, the number of model parameters is typically greater than the number of independent observed data. Second, small singular-values are generated by the discretization of the physical rupture process and amplify the effect of noise in the inversion. As a result, one can find different slip distributions that fit the data equally well. This ill-posedness can be mitigated by decreasing the number of model parameters, hence improving their relationship to the observed data. In this article, we propose a fuzzy function approximation approach to describe the spatial slip function. In particular, we use an Adaptive Network-based Fuzzy Inference System (ANFIS) to find the most adequate discretization for the spatial variation of slip on the fault. The fuzzy basis functions and their respective amplitudes are optimized through hybrid learning. We solve this earthquake source problem in the frequency domain, searching for independent optimal spatial slip distributions for each frequency. The approximated frequency-dependent spatial slip functions are then used to compute the forward relationship between slip on the fault and ground motion. The method is constrained through Tikhonov regularization, requiring a smooth spatial slip variation. We discuss how the number of model parameters can be decreased while keeping the inversion stable and achieving an adequate resolution. The proposed inversion method is tested using the SIV1-benchmark exercise.