The high cost associated with the evaluation of Hartree-Fock exchange (HFX) makes hybrid functionals computationally challenging for large systems. In this work, we present an efficient way to accelerate HFX calculations with numerical atomic basis sets. Our approach is based on the recently proposed interpolative separable density fitting (ISDF) decomposition to construct a low-rank approximation of the HFX matrix, which avoids explicit calculations of the electron repulsion integrals (ERIs) and significantly reduces the computational cost. We implement the ISDF method for hybrid functional (PBE0) calculations in the HONPAS package. We take benzene and polycyclic aromatic hydrocarbon molecules as examples and demonstrate that hybrid functionals with ISDF yield quite promising results at a significantly reduced computational cost. Especially, the ISDF approach reduces the total cost of the evaluating HFX matrix by nearly 2 orders of magnitude compared to conventional approaches of direct evaluation of ERIs.