Step Size Seeders
=================

Introduction
------------

In Pyralysis, **step size seeders** (``StepSizeSeeder``) estimate an initial step size
:math:`\alpha_k` before a :doc:`linesearchers` refines it. Seeders use gradient history,
function values, or proximal residuals to provide a good starting point, which reduces
backtracking iterations and improves convergence.

.. figure:: diagrams/png/fig_step_size_trajectories.png
   :alt: Optimization paths for small, large, and seeded steps on quadratic level curves
   :width: 100%

   **Why seeders matter** — illustrative paths on an ill-conditioned quadratic: (a) step too
   small (slow zig-zag); (b) step too large (oscillation); (c) seeded / near-optimal steps.
   Full discussion: :ref:`step-size-geometry` in :doc:`linesearchers`.

.. figure:: diagrams/png/fig06_linesearch_seeders-2.png
   :alt: Step-size seeder hierarchy
   :width: 100%

   **Step-size seeder hierarchy** (``StepSizeSeeder``, Barzilai–Borwein family, and
   interpolation seeders).

.. note::
   Seeders do **not** replace line search. They only supply the initial guess passed to
   :meth:`~pyralysis.optimization.linesearch.LineSearcher._get_initial_step_size`.
   See :doc:`linesearchers` for Armijo, Brent, and related strategies.

Role in the optimization stack
------------------------------

Each gradient-based optimizer (see :doc:`optimizers`) calls a ``LineSearcher`` at every
iteration. When a seeder is attached to the line searcher, the workflow is:

1. Optimizer computes a search direction :math:`d_k`.
2. Seeder estimates :math:`\alpha_k^{(0)}` from history.
3. Line searcher accepts or backtracks from :math:`\alpha_k^{(0)}` until its conditions hold.

Barzilai–Borwein family
-----------------------

All BB seeders inherit from ``BarzilaiBorwein`` and share two secant step estimates
computed from

.. math::

   s_{k-1} = x_k - x_{k-1}, \quad y_{k-1} = g_k - g_{k-1}

where :math:`g_k = \nabla f(x_k)` :cite:`barzilai_borwein_1988`.

**Long step (BB1)**

.. math::

   \alpha^{\mathrm{LONG}}_k = \frac{s_{k-1}^T s_{k-1}}{s_{k-1}^T y_{k-1}}
   = \frac{\|s_{k-1}\|^2}{s_{k-1}^T y_{k-1}}

**Short step (BB2)**

.. math::

   \alpha^{\mathrm{SHORT}}_k = \frac{s_{k-1}^T y_{k-1}}{y_{k-1}^T y_{k-1}}
   = \frac{s_{k-1}^T y_{k-1}}{\|y_{k-1}\|^2}

On the first iteration (no history), all BB seeders return ``init_step`` (default ``1.0``).
Invalid denominators yield ``np.nan``; selection logic then falls back to the other variant
or to ``min_step``.

BarzilaiBorweinCyclic
~~~~~~~~~~~~~~~~~~~~~

**Class:** ``BarzilaiBorweinCyclic``

Cycles through BB1 and BB2 with configurable period ``cycle_length`` (default ``4``, must be
:math:`\geq 2`). At cycle position ``iteration % cycle_length``:

- **Even position:** prefer BB1, fall back to BB2.
- **Odd position:** prefer BB2, fall back to BB1.

.. math::

   \alpha_k = \begin{cases}
       \alpha^{\mathrm{LONG}}_k & \text{if cycle position is even} \\
       \alpha^{\mathrm{SHORT}}_k & \text{if cycle position is odd}
   \end{cases}

**Use when:** you want explicit control over how often long vs short BB steps are tried;
default ``cycle_length=4`` alternates two BB1 and two BB2 steps per cycle.

BarzilaiBorweinAlternating
~~~~~~~~~~~~~~~~~~~~~~~~~~

**Class:** ``BarzilaiBorweinAlternating``

Backward-compatible alias for ``BarzilaiBorweinCyclic(cycle_length=2)``: odd iterations
prefer BB2, even iterations prefer BB1. Any ``cycle_length`` passed at construction is
ignored.

BarzilaiBorweinAdaptiveMin1 (ABBmin1)
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

**Class:** ``BarzilaiBorweinAdaptiveMin1``

Adaptive rule from Dai & Fletcher :cite:`dai_fletcher_2005` using a sliding window of
recent BB2 values:

.. math::

   \alpha_k = \begin{cases}
       \min\{\alpha^{\mathrm{SHORT}}_i : i \in \text{window}\}
           & \text{if } \alpha^{\mathrm{SHORT}} / \alpha^{\mathrm{LONG}} < \tau \\
       \alpha^{\mathrm{LONG}}_k & \text{otherwise}
   \end{cases}

**Parameters:** ``tau`` (default ``0.8``), ``window`` (default ``10``).

**Use when:** ill-conditioned problems where short BB steps become overly conservative.

BarzilaiBorweinAdaptiveMin2 (ABBmin2)
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

**Class:** ``BarzilaiBorweinAdaptiveMin2``

Simpler adaptive rule:

.. math::

   \alpha_k = \begin{cases}
       \alpha_{k-1}
           & \text{if } \alpha^{\mathrm{SHORT}} / \alpha^{\mathrm{LONG}} < \tau
             \text{ and } \alpha_{k-1} \text{ exists} \\
       \alpha^{\mathrm{LONG}}_k & \text{otherwise}
   \end{cases}

**Parameters:** ``tau`` (default ``0.9``).

**Use when:** well-conditioned problems with relatively stable step sizes.

BarzilaiBorweinRetarded
~~~~~~~~~~~~~~~~~~~~~~~

**Class:** ``BarzilaiBorweinRetarded``

Uses a **delayed** gradient difference instead of the standard one-step difference:

.. math::

   y_k = g_k - g_{k-m}, \quad \alpha_k = \frac{s_k^T s_k}{s_k^T y_k}

where :math:`s_k = x_k - x_{k-1}` and ``m`` is ``delay`` (default ``3``). When
``delay=1``, this reduces to BB1 with the usual :math:`y_k = g_k - g_{k-1}`.

Falls back to ``init_step`` until the circular gradient buffer holds ``delay`` past
gradients. Only BB1 is selected (BB2 is computed but not used).

**Use when:** you want BB1 with a longer effective gradient memory; can help on noisy
or slowly varying curvature landscapes.

ProximalBarzilaiBorwein
~~~~~~~~~~~~~~~~~~~~~~~

**Class:** ``ProximalBarzilaiBorwein``

Extends ``BarzilaiBorweinCyclic`` for **proximal** objectives (FISTA). Replaces
:math:`y` with the proximal residual

.. math::

   r_k = \operatorname{prox}_{\rho}(x_k) - \operatorname{prox}_{\rho}(x_{k-1})

and uses cyclic BB1/BB2 proximal formulas:

.. math::

   \alpha^{\mathrm{LONG}}_k = \frac{s_k^T s_k}{s_k^T r_k}, \quad
   \alpha^{\mathrm{SHORT}}_k = \frac{s_k^T r_k}{r_k^T r_k}

with :math:`\rho = 1 / \alpha_{k-1}` from the previous accepted step.

**Requirements:** exactly **0 or 1** non-differentiable term in the objective. With no
non-smooth term the proximal map is the identity and gradient-based BB seeders should be
used instead.

**Not compatible with SDMM** when multiple non-smooth terms are split separately; use
gradient-based BB seeders for the inner SDMM optimizer. See :doc:`compressed_sensing`.

Interpolation seeders
---------------------

CubicInterpolationSeeder
~~~~~~~~~~~~~~~~~~~~~~~~

**Class:** ``CubicInterpolationSeeder``

Fits a cubic polynomial along the **previous** search direction using
:math:`f(x_{k-1})`, :math:`f(x_k)`, and gradients at both points, finds the minimizer
:math:`\alpha_{\mathrm{opt}}` along that direction, then scales to the current direction:

.. math::

   \alpha_{\mathrm{new}} = \alpha_{\mathrm{opt}} \cdot
   \frac{\|\nabla f(x_k)\|}{\|\nabla f(x_{k-1})\|}

**Use when:** smooth objectives with varying curvature; more accurate than quadratic but
needs reliable function and gradient history.

QuadraticInterpolationSeeder
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

**Class:** ``QuadraticInterpolationSeeder``

Same framework as the cubic seeder but with a quadratic fit (function value and
gradients at endpoints of the previous step). Closed-form minimizer when the quadratic
curvature coefficient is positive.

**Use when:** smoother, more stable alternative to cubic; good fallback when cubic
interpolation fails numerically.

Shared parameters
-----------------

All seeders support:

+------------------+---------------+------------------------------------------+
| Parameter        | Default       | Role                                     |
+==================+===============+==========================================+
| ``init_step``    | ``1.0``       | Step when history is insufficient        |
+------------------+---------------+------------------------------------------+
| ``min_step``     | ``1e-10``     | Floor on accepted estimates              |
+------------------+---------------+------------------------------------------+
| ``curvature_tol``| ``1e-12``     | Denominator guard in BB formulas         |
+------------------+---------------+------------------------------------------+

Selection guide
---------------

+---------------------------+-----------------------------------------------+
| Problem type              | Suggested seeder                              |
+===========================+===============================================+
| General smooth imaging    | ``BarzilaiBorweinAlternating`` or ``Cyclic``  |
+---------------------------+-----------------------------------------------+
| Ill-conditioned           | ``BarzilaiBorweinAdaptiveMin1``               |
+---------------------------+-----------------------------------------------+
| Well-conditioned, stable  | ``BarzilaiBorweinAdaptiveMin2``               |
+---------------------------+-----------------------------------------------+
| Longer gradient memory    | ``BarzilaiBorweinRetarded``                   |
+---------------------------+-----------------------------------------------+
| FISTA / single prox term  | ``ProximalBarzilaiBorwein``                   |
+---------------------------+-----------------------------------------------+
| Strongly quadratic local  | ``CubicInterpolationSeeder``                  |
| behavior                  |                                               |
+---------------------------+-----------------------------------------------+
| Cubic unstable            | ``QuadraticInterpolationSeeder``              |
+---------------------------+-----------------------------------------------+
| Fixed / no history        | No seeder (line searcher default step)        |
+---------------------------+-----------------------------------------------+

Implementation example
----------------------

.. code-block:: python

   from pyralysis.optimization.linesearch import (
       BacktrackingArmijo,
       BarzilaiBorweinCyclic,
       BarzilaiBorweinRetarded,
       ProximalBarzilaiBorwein,
       CubicInterpolationSeeder,
   )

   # Gradient-based BB with 4-step cycle
   seeder = BarzilaiBorweinCyclic(cycle_length=4, min_step=1e-10)
   ls = BacktrackingArmijo(objective_function=of, seeder=seeder)

   # Retarded BB1 with delay m=3
   seeder = BarzilaiBorweinRetarded(delay=3, init_step=0.5)

   # FISTA: attach ProximalBarzilaiBorwein to FISTABacktracking (see linesearchers)
   from pyralysis.optimization.linesearch import FISTABacktracking
   prox_seeder = ProximalBarzilaiBorwein(cycle_length=2)
   fista_ls = FISTABacktracking(objective_function=of, seeder=prox_seeder)

References
----------

.. bibliography::
   :filter: docname in docnames

See also
--------

- :doc:`linesearchers` — how seeders connect to line search algorithms
- :doc:`optimizers` — gradient and quasi-Newton methods that consume line searchers
- :doc:`compressed_sensing` — FISTA and proximal seeding

----

:doc:`linesearchers` | :doc:`optimizers` | :doc:`optimization`
