Earthquake focal mechanism data provide information about the stress state at the origin of earthquakes. The inversion methods that are commonly used to infer the stress tensor from focal mechanisms have varying complexity but always rely on a number of assumptions. We present an iterative method built upon a classic linear stress tensor inversion that allows for relaxing the assumption on shear stress magnitudes while preserving the computational simplicity of the linear problem. Every iteration of our method computes the least‐squares solution of the problem, which makes the method fast enough to estimate the inverted parameter errors with nonparametric resampling methods such as bootstrapping. Following previous studies, this method removes the fault plane ambiguity in focal mechanism data by selecting the nodal plane that best satisfies the Mohr–Coulomb failure criterion. We first test the performance and robustness to noise of the proposed method on synthetic data sets and then apply it to data from the Southern California and Geysers geothermal field data sets. We focus the study on investigating the consequences of relaxing the assumption of constant shear stress magnitudes. Our variable shear method successfully generalizes its constant shear counterpart: it is able to perform similarly when the constant shear assumption is a good approximation and provides more accurate results when it is not. We provide the Python package iterative linear stress inversion to implement the proposed method.