We propose a novel method for solving linear systems with large dense operators, right-hand sides, and unknown matrices. We extend the idea of performing a low-rank or \(\mathcal {H}^{2}\) approximation of the operator to simultaneously approximating of all three matrices, using a shared bases assumption. We showcase the effectiveness of our method through its application to a particularly challenging problem in the field of seismic imaging, specifically the use of Multidimensional Deconvolution (MDD) for redatuming applications. While offering more accuracy than conventional correlation-based redatuming methods, MDD faces challenges due to the ill-posed nature of the underlying inverse problem and the requirement to handle large, dense, complex-valued matrices. These obstacles have long limited the adoption of MDD in the geophysical community. Recent interest in this technology has spurred the development of new strategies to enhance the robustness of the inversion process and reduce its computational overhead. Our proposed approach can greatly alleviate the data-heavy nature of MDD. Moreover, since in 3d applications the matrices do not lend themselves to global low rank approximations, we introduce a novel \(\mathcal {H}^2\) -like approximation. With this work, we aim to streamline MDD implementations, fostering efficiency and controlling accuracy in wavefield reconstruction. This innovation holds potential for broader applications in the geophysical domain and beyond.
