The Gauss–Newton method is one of the most common choices for solving nonlinear systems . The idea is to minimize the corresponding least squares problem by solving a sequence of linearized problems. The method has been extended to account for non-smooth regularizers, leading to proximal Gauss–Newton algorithms. At each iteration, the algorithm avoids computation or storage of the Hessian, but requires two in principle costly operations: inversion of the matrix and computation of the proximal point in a variable metric. To overcome this limitation, we propose an inexact version of the proximal Gauss–Newton algorithm based on an iterative approximation of the linearized sub-problems at each iteration. Numerical experiments on both convex penalized nonlinear least squares problems arising in binary classification, as well as non-convex bound-constrained nonlinear least squares problems, show promising performance of the suggested approach.