Abstract: Algorithms for bilevel optimization often encounter Hessian computations, which are prohibitive in high dimensions. While recent works offer first-order methods for unconstrained bilevel problems, the \textit{constrained} setting remains relatively underexplored. We present fully first-order linearly constrained optimization methods with finite-time hypergradient stationarity guarantees. For linear \textit{equality} constraints, we attain $\epsilon$-stationarity in $\widetilde{O}(\epsilon^{-2})$ gradient oracle calls, which is nearly-optimal. For linear \textit{inequality} constraints, we attain $(\delta,\epsilon)$-Goldstein stationarity in $\widetilde{O}({\delta^{-1} \epsilon^{-4}})$ or $\widetilde{O}(d{\delta^{-1} \epsilon^{-3}})$ gradient oracle calls, where $d$ is the upper-level dimension. Along the way, we develop new nonsmooth nonconvex optimization methods with inexact oracles. Our numerical experiments verify these guarantees.