Skip to content

Polynomial Preconditioner

Polynomial preconditioning for iterative solvers.

torchlinops.alg.polynomial_preconditioner 

polynomial_preconditioner(
    T: NamedLinop,
    degree: int,
    norm: Literal["l_2", "l_inf", "ifista"] = "l_2",
    lower_eig: float = 0.0,
    upper_eig: float = 1.0,
) -> NamedLinop

Construct a polynomial preconditioner \(P(T)\) for a linear operator \(T\).

Given eigenvalue bounds, computes optimal polynomial coefficients and returns a NamedLinop that applies \(P(T) = \sum_k c_k T^k\).

PARAMETER DESCRIPTION
T

The operator to precondition.

TYPE: NamedLinop

degree

Polynomial degree. Use -1 for no preconditioning, 0+ for polynomial preconditioning.

TYPE: int

norm

Optimization criterion for computing the polynomial coefficients:

  • "l_2" -- Minimizes \(\int |1 - x p(x)|^2 dx\).
  • "l_inf" -- Minimizes \(\sup |1 - x p(x)|\) (Chebyshev).
  • "ifista" -- Coefficients from DOI:10.1137/140970537.

TYPE: ``{"l_2", "l_inf", "ifista"}`` DEFAULT: ``"l_2"``

lower_eig

Lower bound on the eigenvalue spectrum of \(T\).

TYPE: float DEFAULT: 0.0

upper_eig

Upper bound on the eigenvalue spectrum of \(T\).

TYPE: float DEFAULT: 1.0

RETURNS DESCRIPTION
NamedLinop

The polynomial preconditioner applied to T.
Source code in src/torchlinops/alg/poly.py 
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
def polynomial_preconditioner(
    T: NamedLinop,
    degree: int,
    norm: Literal["l_2", "l_inf", "ifista"] = "l_2",
    lower_eig: float = 0.0,
    upper_eig: float = 1.0,
) -> NamedLinop:
    """Construct a polynomial preconditioner $P(T)$ for a linear operator $T$.

    Given eigenvalue bounds, computes optimal polynomial coefficients and
    returns a ``NamedLinop`` that applies $P(T) = \\sum_k c_k T^k$.

    Parameters
    ----------
    T : NamedLinop
        The operator to precondition.
    degree : int
        Polynomial degree. Use ``-1`` for no preconditioning, ``0+`` for
        polynomial preconditioning.
    norm : ``{"l_2", "l_inf", "ifista"}``, default ``"l_2"``
        Optimization criterion for computing the polynomial coefficients:

        - ``"l_2"`` -- Minimizes $\\int |1 - x p(x)|^2 dx$.
        - ``"l_inf"`` -- Minimizes $\\sup |1 - x p(x)|$ (Chebyshev).
        - ``"ifista"`` -- Coefficients from DOI:10.1137/140970537.
    lower_eig : float, default 0.0
        Lower bound on the eigenvalue spectrum of $T$.
    upper_eig : float, default 1.0
        Upper bound on the eigenvalue spectrum of $T$.

    Returns
    -------
    NamedLinop
        The polynomial preconditioner applied to *T*.
    """
    Id: NamedLinop = Identity()  # Fixing shapes
    if degree < 0:
        return Id

    if norm == "l_2":
        c, _ = l_2_opt(degree, lower_eig, upper_eig)
    elif norm == "l_inf":
        c, _ = l_inf_opt(degree, lower_eig, upper_eig)
    elif norm == "ifista":
        c = ifista_coeffs(degree)
    else:
        raise ValueError(f"Unknown polynomial preconditioning norm option: {norm}")

    def phelper(c) -> NamedLinop:
        if c.size == 1:
            return c[0] * Id
        L = c[0] * Id  # ... -> ...
        R = phelper(c[1:]) @ T @ Id  # ... -> ...
        return L + R  # ... -> ...

    P = phelper(c)
    return P