We introduce BiLO (Bilevel Local Operator Learning), a novel neural network-based approach for solving inverse problems in partial differential equations (PDEs). BiLO formulates the PDE inverse problem as a bilevel optimization problem: at the upper level, we optimize PDE parameters by minimizing data loss, while at the lower level, we train a neural network to locally approximate the PDE solution operator near given PDE parameters. This localized approximation enables accurate descent direction estimation for the upper-level optimization. We apply gradient descent simultaneously on both the upper and lower level optimization problems, leading to an effective and fast algorithm. Additionally, BiLO can infer unknown functions within PDEs by introducing an auxiliary variable. Extensive experiments across various PDE systems demonstrate that BiLO enforces strong PDE constraints, is robust to sparse and noisy data, and eliminates the need for manually balancing residual and data loss, a common challenge in soft PDE constraints. We also discuss how to apply the BILO for uncertainty quantification in a Bayesian framework.