Probabilistic Seismic Hazard Analysis (PSHA) traditionally relies on two computationally intensive approaches: (a) Riemann Sum and (b) conventional Monte Carlo (MC) integration. The former requires fine slices across magnitude, distance, and ground motion, and the latter demands extensive synthetic earthquake catalogs. Both approaches become notably resource intensive for low-probability seismic hazards, where achieving a COV of 1% for a 10−4 annual hazard probability may require 108 MC samples. We introduce Adaptive Importance Sampling (AIS) PSHA, a novel framework to approximate optimal importance sampling (IS) distributions and dramatically reduce the number of MC samples to estimate hazards. We evaluate the efficiency and accuracy of our proposed framework using Pacific Earthquake Engineering Research Center (PEER) PSHA benchmarks that cover various seismic sources, including areal, vertical, and dipping faults, as well as combined types. Our approach computes seismic hazard up to 3.7×104 and 7.1×103 times faster than Riemann Sum and traditional MC methods, respectively, maintaining COVs below 1%. We also propose an enhanced approach with a “smart” AIS PSHA variant that leverages the sampling densities from similar ground motion intensities. This variant outperforms even “smart” implementations of Riemann Sum with enhanced grid discretizations by a factor of up to 130. Moreover, we demonstrate theoretically that optimal IS distributions are equivalent to hazard disaggregation distributions. Empirically, we show the approximated optimal IS and the disaggregation distributions are closely alike, e.g., with a Kolmogorov–Smirnov statistic between 0.017 and 0.113. This approach is broadly applicable, especially for PSHA cases requiring extensive logic trees and epistemic uncertainty.