Non-linear statistical inverse problems pose major challenges both for statistical analysis and computation. Likelihood-based estimators typically lead to non-convex and possibly multimodal optimization landscapes, and Markov chain Monte Carlo (MCMC) methods may mix exponentially slowly. We propose a class of computationally tractable estimators--plug-in and PDE-penalized M-estimators--for inverse problems defined through operator equations of the form $L_f u = g$, where $f$ is the unknown parameter and $u$ is the observed solution. The key idea is to replace the exact PDE constraint by a weakly enforced relaxation, yielding conditionally convex and, in many PDE examples, nested quadratic optimization problems that avoid evaluating the forward map $G(f)$ and do not require PDE solvers. For prototypical non-linear inverse problems arising from elliptic PDEs, including the Darcy flow model $L_f u = \nabla\!\cdot(f\nabla u)$ and a steady-state Schrödinger model, we prove that these estimators attain the best currently known statistical convergence rates while being globally computable in polynomial time. In the Darcy model, we obtain an explicit sub-quadratic $o(N^2)$ arithmetic runtime bound for estimating $f$ from $N$ noisy samples. Our analysis is based on new generalized stability estimates, extending classical stability beyond the range of the forward operator, combined with tools from nonparametric M-estimation. We also derive adaptive rates for the Darcy problem, providing a blueprint for designing provably polynomial-time statistical algorithms for a broad class of non-linear inverse problems. Our estimators also provide principled warm-start initializations for polynomial-time Bayesian computation.