  We present a method to construct an efficient approximation to the bare
exchange and screened direct interaction kernels of the Bethe-Salpeter
Hamiltonian for periodic solid state systems via the interpolative separable
density fitting technique. We show that the cost of constructing the
approximate Bethe-Salpeter Hamiltonian scales nearly optimally as
$\mathcal{O}(N_k)$ with respect to the number of samples in the Brillouin zone
$N_k$. In addition, we show that the cost for applying the Bethe-Salpeter
Hamiltonian to a vector scales as $\mathcal{O}(N_k \log N_k)$. Therefore the
optical absorption spectrum, as well as selected excitation energies can be
efficiently computed via iterative methods such as the Lanczos method. This is
a significant reduction from the $\mathcal{O}(N_k^2)$ and $\mathcal{O}(N_k^3)$
scaling associated with a brute force approach for constructing the Hamiltonian
and diagonalizing the Hamiltonian respectively. We demonstrate the efficiency
and accuracy of this approach with both one-dimensional model problems and
three-dimensional real materials (graphene and diamond). For the diamond system
with $N_k=2197$, it takes $6$ hours to assemble the Bethe-Salpeter Hamiltonian
and $4$ hours to fully diagonalize the Hamiltonian using $169$ cores when the
brute force approach is used. The new method takes less than $3$ minutes to set
up the Hamiltonian and $24$ minutes to compute the absorption spectrum on a
single core.
