This talk introduces a new class of algorithms for solving large-scale linear inverse problems based on new flexible and inexact Golub-Kahan factorizations. The proposed methods iteratively compute regularized solutions by approximating a solution to (re)weighted least squares problems via projection onto adaptively generated subspaces, where the constraint subspaces are formally equipped with iteration-dependent preconditioners or inexactness. Although this framework is particularly suited for handling general data fidelity functionals, e.g., those expressed in the p-norm, its use is also discussed in the context of problems with unmatched transposes. The new solvers offer a flexible and inexact Krylov subspace alternative to other existing Krylov-based approaches. Numerical experiments in imaging applications, such as deblurring and computed tomography, highlight the effectiveness and competitiveness of the proposed methods with respect other popular methods.