Optimizers ========== Introduction ------------ Pyralysis provides **gradient-based optimizers** for smooth (or smooth + handled-by-projection) objective functions built from visibility terms and differentiable regularizers. Each optimizer iteratively updates the image :math:`x_k` using gradients and a :doc:`linesearchers` instance to choose the step size along the search direction. For objectives with **non-differentiable** sparsity or TV terms, use proximal methods in :doc:`compressed_sensing` (FISTA, SDMM) instead of the classes documented here. .. figure:: diagrams/png/fig05_optimizers-1.png :alt: Optimizer hierarchy :width: 100% **Optimizer hierarchy** (``Optimizer``, gradient/proximal optimizers, and CG/LBFGS families). .. note:: Proximal optimizers ``FISTA`` and ``SDMM`` are covered in :doc:`compressed_sensing`. They use ``FISTABacktracking`` and optional :doc:`step_size_seeders` rather than the Armijo line searchers attached to CG/L-BFGS. Overview -------- .. list-table:: Gradient-based optimizers :header-rows: 1 :widths: 14 12 28 10 8 8 20 * - Optimizer - Type - Description - Robustness - Speed - Memory - Typical use case * - Hestenes-Stiefel - CG - Strong global convergence, nonlinear problems - Medium - Fair - Low - Nonlinear, large-scale * - Fletcher-Reeves - CG - Classic CG, quadratic problems - Medium - Fair - Low - Nonlinear, large-scale * - Polak-Ribiere - CG - Dynamic gradient update - Medium - Good - Low - Nonlinear, large-scale * - Liu-Storey - CG - Image restoration, unconstrained problems - Medium - Good - Low - Image restoration * - Dai-Yuan - CG - Performance-oriented CG scaling - Medium - Fair - Low - Nonlinear, large-scale * - Hager-Zhang - CG - Modern nonlinear CG, stable and fast - High - Good - Low - Preferred CG choice * - RMIL - CG - Sufficient descent, image recovery - Medium - Good - Low - Denoising, image recovery * - LBFGS - Quasi-Newton - Limited-memory BFGS - High - Best - Medium - Default for smooth imaging **Practical selection** - Start with **LBFGS** for most smooth imaging problems. - Prefer **HagerZhang** when memory for L-BFGS history is tight but a quasi-Newton step is not required. - Other CG variants are available for experimentation or when their theoretical properties match the problem. **References:** :cite:`hestenes_stiefel_1952,fletcher_reeves_1964,polak_ribiere_1969,polyak_1969,liu_storey_1991,dai_yuan_1999,hager_zhang_2005,rivaie_2012,nocedal_wright_2006` Common setup ------------ Every gradient optimizer needs: 1. A reconstructed **parameter** (usually an ``Image``). 2. An :class:`~pyralysis.optimization.objective_function.ObjectiveFunction` (see :doc:`objective_functions` and :doc:`optimization`). 3. A :doc:`linesearchers` instance (optionally with a :doc:`step_size_seeders`). .. code-block:: python from pyralysis.optimization.linesearch import BacktrackingArmijo, BarzilaiBorweinAlternating from pyralysis.optimization.optimizer import LBFGS seeder = BarzilaiBorweinAlternating(min_step=1e-10) ls = BacktrackingArmijo(objective_function=of, seeder=seeder) optim = LBFGS( parameter=image, objective_function=of, linesearch=ls, max_iter=100, io_handler=ioh, ) reconstructed = optim.optimize() Conjugate gradient methods -------------------------- All CG classes inherit from ``ConjugateGradient``. At iteration :math:`k` the search direction is .. math:: d_k = -g_k + \beta_k \, d_{k-1} where :math:`g_k = \nabla f(x_k)`. The line searcher finds :math:`\alpha_k` along :math:`d_k`; see :doc:`linesearchers`. **Gradient / line-search convention** ``ObjectiveFunction.dphi`` follows: ``search_direction = -dphi``. After ``calculate_gradient()``, ``dphi`` holds :math:`g`. CG sets ``dphi = -d`` before the line search so the 1D objective ``f1dim`` sees a consistent sign. Restart strategy ~~~~~~~~~~~~~~~~ When either condition holds, :math:`\beta_k = 0` (steepest descent step) :cite:`powell_1977`: 1. **Negative beta:** :math:`\beta_k \leq 0`. 2. **Powell restart:** :math:`g_{k+1}^T g_k < \eta \|g_{k+1}\|^2` with :math:`\eta = 0.2`. CG variants and :math:`\beta_k` formulas ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ HestenesStiefel ^^^^^^^^^^^^^^^ .. math:: \beta_k^{\mathrm{HS}} = \frac{g_{k+1}^T (g_{k+1} - g_k)}{(g_{k+1} - g_k)^T d_k} FletcherReeves ^^^^^^^^^^^^^^ .. math:: \beta_k^{\mathrm{FR}} = \frac{\|g_{k+1}\|^2}{\|g_k\|^2} PolakRibiere ^^^^^^^^^^^^ .. math:: \beta_k^{\mathrm{PRP}} = \frac{g_{k+1}^T (g_{k+1} - g_k)}{\|g_k\|^2} LiuStorey ^^^^^^^^^ .. math:: \beta_k^{\mathrm{LS}} = -\frac{(g_{k+1} - g_k)^T g_{k+1}}{g_k^T d_k} DaiYuan ^^^^^^^ .. math:: \beta_k^{\mathrm{DY}} = \frac{\|g_{k+1}\|^2}{(g_{k+1} - g_k)^T d_k} HagerZhang ^^^^^^^^^^ With :math:`y_k = g_{k+1} - g_k`: .. math:: \beta_k^{\mathrm{HZ}} = \frac{y_k^T g_{k+1}}{d_k^T y_k} - 2 \frac{\|y_k\|^2 (d_k^T g_{k+1})}{(d_k^T y_k)^2} RMIL ^^^^ .. math:: \beta_k^{\mathrm{RMIL}} = \frac{g_{k+1}^T (g_{k+1} - g_k)}{d_k^T (d_k - g_k)} Per-optimizer examples ~~~~~~~~~~~~~~~~~~~~~~ .. code-block:: python from pyralysis.optimization.optimizer import HagerZhang, FletcherReeves, RMIL optim = HagerZhang( image=image, objective_function=of, linesearch=ls, max_iter=100, io_handler=ioh, ) result = optim.optimize() LBFGS (limited-memory BFGS) --------------------------- **Class:** ``LBFGS`` Quasi-Newton method that approximates the inverse Hessian using the most recent ``max_corrections`` pairs :math:`(s_i, y_i)` without forming the full matrix. **Key parameters** +------------------------+---------------+------------------------------------------------+ | Parameter | Default | Role | +========================+===============+================================================+ | ``max_corrections`` | ``10`` | Max :math:`(s,y)` pairs (RAM + disk combined) | +------------------------+---------------+------------------------------------------------+ | ``memory_threshold_gb``| ``100.0`` | Switch vectorized vs sequential two-loop | +------------------------+---------------+------------------------------------------------+ | ``use_disk_storage`` | ``False`` | Spill pairs to Zarr when ``in_memory_limit`` | | | | is set | +------------------------+---------------+------------------------------------------------+ | ``in_memory_limit`` | ``None`` | Max pairs kept in RAM when spilling enabled | +------------------------+---------------+------------------------------------------------+ - ``max_corrections=0`` reduces to **steepest descent** with the same line search and stopping rules. - Pairs are accepted only if :math:`\Re(y^H s) > 0`; :math:`\rho = 1/\Re(y^H s)` uses a dtype-aware machine epsilon floor. **Memory optimizations** - Lazy deque allocation for correction pairs. - Hybrid vectorized/sequential two-loop recursion based on ``memory_threshold_gb``. - Optional Zarr-backed spill for very large images. - Cached :math:`\gamma` scaling and consistent Dask chunking across history. .. code-block:: python from pyralysis.optimization.optimizer import LBFGS optim = LBFGS( parameter=image, objective_function=of, linesearch=ls, max_iter=100, io_handler=ioh, max_corrections=10, ) # Large images: spill older pairs to disk optim = LBFGS( parameter=large_image, objective_function=of, linesearch=ls, max_iter=100, max_corrections=10, use_disk_storage=True, in_memory_limit=5, memory_threshold_gb=50.0, ) Saving intermediate results --------------------------- Pass an I/O handler to write intermediate images during ``optimize()`` with ``partial_parameter=True``. The optimizer passes the **live** ``Image`` to the handler; host/device conversion is the handler's job (see :ref:`partial-parameter-io` in :doc:`io_operations`). .. code-block:: python from pyralysis.io import FITS from pyralysis.io.zarr import ZarrArray # FITS or ZarrArray — both accept CuPy-backed Dask Images from GPU runs ioh = FITS(use_dask=True) # or ZarrArray() optim = LBFGS( parameter=image, objective_function=of, linesearch=ls, io_handler=ioh, ) optim.optimize( partial_parameter=True, partial_parameter_interval=10, # partial_10, partial_20, ...; final always saved ) Implementation notes -------------------- - All CG variants support **complex** parameters with stable inner products. - ``GradientNormError`` is raised when the gradient norm vanishes (convergence or degeneracy). - CG ``conjugate_gradient_parameter`` returns host scalars to avoid redundant ``.compute()`` calls. Dask iterate state (custom optimizers) ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Gradient optimizers assign a new lazy ``parameter.data`` each iteration (line search, acceleration, or quasi-Newton update). Dask then **accumulates graph lineage** on the carried iterate unless it is collapsed before the next loop. Without that step, memory and scheduler overhead grow **linearly with iteration count** and can OOM on large images. Built-in optimizers call :meth:`~pyralysis.optimization.optimizer.Optimizer._finalize_iteration_state` at the end of each iteration. It persists Dask-backed parameters in place (no-op for NumPy) via :meth:`~pyralysis.reconstruction.parameter.Parameter.persist`. **When implementing a new optimizer**, call ``_finalize_iteration_state`` on every ``Parameter`` you carry to the next iteration (primary iterate; FISTA also finalizes the accelerated point ``y_k``). For carried state stored as raw ``dask.array.Array`` objects, use :meth:`~pyralysis.optimization.optimizer.Optimizer._finalize_dask_arrays` instead (SDMM does this for ``z`` and ``u`` each outer iteration). This is **independent of** optional caches such as ``ObjectiveFunction.persist_gradient`` (aggregate gradient) or LBFGS direction / correction-pair persistence (multi-reuse within one iteration). .. code-block:: python from pyralysis.optimization.optimizer import Optimizer class MyOptimizer(Optimizer): def optimize(self): x = self.parameter for k in range(self.max_iter): # ... compute search direction, line search, assign x.data ... x = updated_param self._finalize_iteration_state(x) return x See :doc:`performance` for the full persist policy (gradient, proximal, LBFGS history, SDMM auxiliary variables). Progress vs diagnostics ----------------------- Optimizer iteration lines (``Iteration N``, function value, step size, convergence) are **progress** output, not diagnostic logs: - Use ``verbose=True/False`` on ``optimize()`` to gate them. - Internally they go through ``pyralysis.runtime.logging.progress``, which uses ``dask.distributed.print`` when a Distributed client is active (so lines appear once on the client) and built-in ``print`` otherwise. - Do **not** expect to silence iteration lines with ``logging`` levels alone. **Diagnostics** (I/O warnings, setup failures, debug traces) use ``pyralysis.runtime.logging.get_logger`` / the standard ``logging`` module. Configure the ``pyralysis`` logger from a CLI or notebook entry point with ``configure_logging(level="DEBUG")`` when you need subsystem detail without treating progress as ``INFO``. .. bibliography:: :filter: docname in docnames See also -------- - :doc:`linesearchers` — step acceptance and FISTA backtracking - :doc:`step_size_seeders` — Barzilai–Borwein and interpolation initial steps - :doc:`optimization` — objective construction and workflow - :doc:`regularization` — smooth regularizers compatible with these optimizers - :doc:`compressed_sensing` — FISTA and SDMM for non-smooth terms ---- :doc:`linesearchers` | :doc:`step_size_seeders` | :doc:`optimization`