An algebraic multilevel (ML) preconditioner is presented for the Helmholtz equation in heterogeneous media. It is based on a multilevel incomplete LDL(T) factorization and preserves the inherent (complex) symmetry of the Helmholtz equation. The ML preconditioner incorporates two key components for efficiency and numerical stability: symmetric maximum weight matchings and an inverse-based pivoting strategy. The former increases the block-diagonal dominance of the system, whereas the latter controls parallel to L(-1)parallel to for numerical stability. When applied recursively, their combined effect yields an algebraic coarsening strategy, similar to algebraic multigrid methods, even for highly indefinite matrices. The ML preconditioner is combined with a Krylov subspace method and applied as a "black-box" solver to a series of challenging two-and three-dimensional test problems, mainly from geophysical seismic imaging. The numerical results demonstrate the robustness and efficiency of the ML preconditioner, even at higher frequency regimes.