.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "generated/autoexamples/example05_solvers.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_generated_autoexamples_example05_solvers.py: ============================================================ Iterative Solve: CG, LSMR, and Polynomial Preconditioning ============================================================ This example demonstrates the :meth:`solve` method exposed on every PyGROG operator (``SparseFFT``, ``MaskedFFT`` and their :mod:`pygrog.gadgets` decorators), and the :class:`pygrog.PolynomialPreconditioner` accelerator. The pipeline mirrors :doc:`example01_basic_usage`: 1. Simulate non-Cartesian k-space from a BrainWeb phantom. 2. GROG-grid into Cartesian space (sparse and dense paths). 3. Solve the regularised least-squares problem .. math:: \hat x \;=\; \arg\min_x \; \| A x - b \|_2^2 + \lambda^2 \| x \|_2^2 using :func:`pygrog.cg`, :func:`pygrog.lsmr`, and a polynomial- preconditioned :func:`pygrog.cg`. .. GENERATED FROM PYTHON SOURCE LINES 25-39 .. code-block:: Python import matplotlib.pyplot as plt import numpy as np from brainweb_dl import get_mri from mrinufft import get_operator, initialize_2D_spiral from mrinufft.density import voronoi from pygrog import PolynomialPreconditioner from pygrog.calib import GrogInterpolator from pygrog.operator import MaskedFFT, SparseFFT .. GENERATED FROM PYTHON SOURCE LINES 56-58 .. code-block:: Python image = get_mri(0, "T1") .. GENERATED FROM PYTHON SOURCE LINES 69-71 .. code-block:: Python nufft = get_operator("finufft")( .. rst-class:: sphx-glr-script-out .. code-block:: none /home/docs/checkouts/readthedocs.org/user_builds/pygrog/envs/latest/lib/python3.12/site-packages/mrinufft/_utils.py:67: UserWarning: Samples will be rescaled to [-pi, pi), assuming they were in [-0.5, 0.5) warnings.warn( k-space shape: (16, 28800) .. GENERATED FROM PYTHON SOURCE LINES 83-89 GROG calibration and gridding ============================= Set up :class:`~pygrog.calib.GrogInterpolator` to grid the non-Cartesian k-space onto Cartesian samples. Prepare both a sparse path and a dense (gridded) path for comparison. .. GENERATED FROM PYTHON SOURCE LINES 89-122 .. code-block:: Python calib = nufft.adj_op(kspace_nc) # quick low-res-ish phantom estimate calib = calib.astype(np.complex64, copy=False)[None] # (1, *shape) grog = GrogInterpolator( shape=shape, coords=samples.reshape(48, 600, 2), kernel_width=2, oversamp=2.0, image_shape=shape, ) calib_full = smaps * image.astype( np.complex64 ) # (n_coils, *shape) — ground-truth coil images grog.calc_interp_table(calib_full, lamda=0.01, precision=1) kspace_nc_shaped = kspace_nc.reshape(n_coils, 48, 600) # Sparse path sparse = grog.interpolate(kspace_nc_shaped, ret_image=False) op_sp = SparseFFT(plan=grog.plan, smaps=smaps) b_sp = sparse * np.asarray(grog.plan.pre_weights) b_sp = b_sp.reshape(*b_sp.shape[:-1], *op_sp.natural_shape) # Dense/grid path kgrid, mplan = grog.interpolate(kspace_nc_shaped, grid=True) op_m = MaskedFFT(plan=mplan, smaps=smaps) b_m = kgrid .. GENERATED FROM PYTHON SOURCE LINES 123-128 Conjugate Gradient (CG) Solver ============================== Solve the regularized normal equations using CG. PyGROG operators automatically dispatch to CG when Toeplitz acceleration is available. .. GENERATED FROM PYTHON SOURCE LINES 130-142 .. code-block:: Python residuals_cg = [] img_cg = op_m.solve( b_m, method="cg", max_iter=20, damp=1e-3, callback=_cb_cg, ) print(f"CG: {len(residuals_cg)} iterations") .. rst-class:: sphx-glr-script-out .. code-block:: none CG: 20 iterations .. GENERATED FROM PYTHON SOURCE LINES 165-170 LSMR Solver =========== Use LSMR (iterative solver for regularised least-squares) as an alternative to CG for comparison. .. GENERATED FROM PYTHON SOURCE LINES 170-179 .. code-block:: Python img_lsmr = op_m.solve( b_m, method="lsmr", max_iter=20, callback=_cb_lsmr, ) print(f"LSMR: {len(residuals_lsmr)} iterations") .. rst-class:: sphx-glr-script-out .. code-block:: none LSMR: 20 iterations .. GENERATED FROM PYTHON SOURCE LINES 180-187 Polynomial Preconditioned CG (PCG) ================================== :class:`pygrog.PolynomialPreconditioner` builds a polynomial approximation to the inverse of ``A^H A`` from Chebyshev polynomials. Each application requires ``degree`` matrix-vector products but typically accelerates CG convergence significantly. .. GENERATED FROM PYTHON SOURCE LINES 187-202 .. code-block:: Python pc = PolynomialPreconditioner(op_m, degree=3, n_power_iter=10) print(f"Preconditioner spectrum estimate: {pc.spectrum}") print(f"Polynomial coefficients: {pc.coeffs}") img_pcg = op_m.solve( b_m, method="cg", max_iter=20, damp=1e-3, preconditioner=pc, callback=_cb_pcg, ) print(f"PCG: {len(residuals_pcg)} iterations") .. rst-class:: sphx-glr-script-out .. code-block:: none Preconditioner spectrum estimate: (0.0, 1.0315111130475998) Polynomial coefficients: [11.633418048712429, -39.47312118639093, 51.02303560548171, -22.25896137426571] PCG: 20 iterations .. GENERATED FROM PYTHON SOURCE LINES 203-205 Compare convergence =================== .. GENERATED FROM PYTHON SOURCE LINES 205-207 .. code-block:: Python fig, ax = plt.subplots(figsize=(6, 4)) .. image-sg:: /generated/autoexamples/images/sphx_glr_example05_solvers_001.png :alt: Iterative solver convergence (MaskedFFT) :srcset: /generated/autoexamples/images/sphx_glr_example05_solvers_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 220-222 Visualise the reconstructions ============================= .. GENERATED FROM PYTHON SOURCE LINES 222-225 .. code-block:: Python fig, axes = plt.subplots(1, 4, figsize=(14, 3.5)) .. image-sg:: /generated/autoexamples/images/sphx_glr_example05_solvers_002.png :alt: reference, CG, LSMR, PCG (deg=3) :srcset: /generated/autoexamples/images/sphx_glr_example05_solvers_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 247-249 The same `solve()` API works on the sparse path =============================================== .. GENERATED FROM PYTHON SOURCE LINES 249-254 .. code-block:: Python img_cg_sp = op_sp.solve(b_sp, method="cg", max_iter=20, damp=1e-3) print( "SparseFFT solve image shape:", tuple(img_cg_sp.shape), ) .. rst-class:: sphx-glr-script-out .. code-block:: none SparseFFT solve image shape: (217, 181) .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 22.985 seconds) .. _sphx_glr_download_generated_autoexamples_example05_solvers.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: example05_solvers.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: example05_solvers.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: example05_solvers.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_