Spectral unmixing is an important task in the research field of hyperspectral image processing. It can be thought of as a regression problem, where the observed variable (i.e., an image pixel) is to be found as a function of the response variables (i.e., the pure materials in a scene, called endmembers). The Linear Mixing Model (LMM) has received a great deal of attention, due to its simplicity and ease of use in, e.g., optimization problems. Its biggest flaw is that it assumes that any pure material can be characterized by one unique spectrum throughout the entire scene. In many cases this is incorrect: the endmembers face a significant amount of spectral variability caused by, e.g., illumination conditions, atmospheric effects, or intrinsic variability. Researchers have suggested several generalizations of the LMM to mitigate this effect. However, most models lead to ill-posed and highly non-convex optimization problems, which are hard to solve and have hyperparameters that are difficult to tune. In this paper, we propose a two-step LMM that bridges the gap between model complexity and computational tractability. We show that this model leads to only a mildly non-convex optimization problem, which we solve with an interior-point solver. This method requires virtually no hyperparameter tuning, and can therefore be used easily and quickly in a wide range of unmixing tasks. We show that the model is competitive and in some cases superior to existing and well-established unmixing methods and algorithms. We do this through several experiments on synthetic data, real-life satellite data, and hybrid synthetic-real data.