  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 constrained setting remains
relatively underexplored. We present first-order linearly constrained
optimization methods with finite-time hypergradient stationarity guarantees.
For linear equality constraints, we attain $\epsilon$-stationarity in
$\widetilde{O}(\epsilon^{-2})$ gradient oracle calls, which is nearly-optimal.
For linear inequality constraints, we attain $(\delta,\epsilon)$-Goldstein
stationarity in $\widetilde{O}(d{\delta^{-1} \epsilon^{-3}})$ gradient oracle
calls, where $d$ is the upper-level dimension. Finally, we obtain for the
linear inequality setting dimension-free rates of $\widetilde{O}({\delta^{-1}
\epsilon^{-4}})$ oracle complexity under the additional assumption of oracle
access to the optimal dual variable. Along the way, we develop new nonsmooth
nonconvex optimization methods with inexact oracles. We verify these guarantees
with preliminary numerical experiments.
