The development of Approximate Message Passing (AMP) has become a precedent demonstrating the potential of Message Passing (MP) algorithms in solving large-scale linear inverse problems of the form y = A x + w. Not only AMP is provably convergent and Bayes-optimal, but it is also a first-order iterative method that can leverage Plug-and-Play (PnP) denoisers for recovering complex data like natural images. Unfortunately, all of these properties have been shown to hold assuming the measurement operator A is, roughly speaking, an i.i.d. random matrix, which highly limits the applicability of the algorithm. The promising extension of AMP, Vector AMP (VAMP), can handle a much broader range of A while preserving most advantages of AMP, but the algorithm requires inverting a large-scale matrix at each iteration, which makes it computationally intractable. As a result, a wide range of ideas has been proposed on upscaling VAMP while preserving its optimality and generality, and in this thesis we would like to share our contributions in this regard. The first contribution is related to developing a stable and accelerated version of Conjugate Gradient (CG) VAMP (CG-VAMP) -- the VAMP algorithm, where the matrix inversion is approximated by the CG algorithm. The originally proposed version of CG-VAMP exhibits unstable dynamics when even mildly large number of CG iterations is used and in those regimes where CG-VAMP is stable, the resulting fixed point of the algorithm might be much worse than that of VAMP. To allow CG-VAMP to use by an order more CG iterations and approximate VAMP with an almost arbitrary accuracy, we constructed a series of rigorous tools that have a negligible computational cost and that lead to stable performance of the algorithm. Additionally, we developed a combination of stopping criteria for CG that ensures efficient operation of CG-VAMP and faster time-wise convergence without sacrificing the estimation accuracy. Next, we considered an alternative way of pushing the performance of CG-VAMP closer to VAMP's and developed the warm-started CG (WS-CG) that reuses the information generated at the previous outer-loop iterations of the MP algorithm. We show that when the matrix inverse in VAMP is approximated by WS-CG, a fixed point of WS-CG VAMP (WS-CG-VAMP) is a fixed point of VAMP and, therefore, is conjectured to be Bayes-optimal. Importantly, this result is invariant with respect to the number of WS-CG iterations and the resulting algorithm can have the computational cost of AMP while being general and optimal as VAMP. We extend the tools developed for CG to WS-CG and numerically demonstrate the stability and efficiency of WS-CG-VAMP. The final contribution is the development of alternative methods for estimating the divergence of a PnP denoiser used within MP algorithms. This divergence plays a crucial role in stabilizing MP algorithms and ensuring its optimality and predictability. So far, the only suggested method for constructing an estimate of the divergence of a PnP denoiser has been the Black-Box Monte Carlo method. The main drawback of this method is that it requires executing the denoiser an additional time, which, effectively, doubles the cost of most MP algorithms. In this thesis we propose two rigorous divergence estimation methods that avoid such a problem and utilize only the information circulated in every MP algorithm.