In this paper, we propose a technique to factorize any matrix into multiple sparse factors. The resulting factorization, called Flexible Approximate MUlti-layer Sparse Transform (FAμST), yields reduced multiplication costs by the matrix and its adjoint. Such a desirable property can be used to speed up iterative algorithms commonly used to solve high dimensional linear inverse problems. The proposed approach is first motivated, introduced and related to prior art. The compromise between computational efficiency and data fidelity is then investigated, and finally the relevance of the approach is demonstrated on a problem of brain source localization using simulated magnetoencephalography (MEG) signals.