The inversion of spectropolarimetric observations of the solar upper atmosphere is one of the most challenging goals in solar physics. If we account for all relevant ingredients of the spectral line formation process, such as the three-dimensional (3D) radiative transfer out of local thermodynamic equilibrium (NLTE), the task becomes extremely computationally expensive. Instead of generalizing 1D methods to 3D, we have developed a new approach to the inverse problem. In our meshfree method, we do not consider the requirement of 3D NLTE consistency as an obstacle, but as a natural regularization with respect to the traditional pixel-by-pixel methods. This leads to more robust and less ambiguous solutions. We solve the 3D NLTE inverse problem as an unconstrained global minimization problem that avoids repetitive evaluations of the Λ operator. Apart from the 3D NLTE consistency, the method allows us to easily include additional conditions of physical consistency such as the zero divergence of the magnetic field. Stochastic ingredients make the method less prone to ending up within the local minima of the loss function. Our method is capable of solving the inverse problem faster by several orders of magnitude than by using grid-based methods. The method can provide accurate and physically consistent results if sufficient computing time is available, along with approximate solutions in the case of very complex plasma structures or limited computing time.