Polynomial module

Module for polynomial related LazyLinOp-s.

It provides “polynomial for LazyLinOp” functions for which the polynomial variable is itself a linear operator (especially a LazyLinOp). This allows to build a LazyLinOp M=p(L) when p is an instance of a polynomial class inherited from NumPy polynomial package and matrix representation of LazyLinOp L is square.

Below are the provided classes:

Below are the provided functions:

  • xpoly() that returns an instance of a polynomial class inherited from NumPy polynomial classes. It could be (specified by xpoly(coef, kind=...)()):

    • 'chebyshev' which is specialized for Chebyshev polynomials.

    • 'hermite' which is specialized for Hermite “physicists” polynomials.

    • 'hermite_e' which is specialized for Hermite “probabilists” polynomials.

    • 'laguerre' which is specialized for Laguerre polynomials.

    • 'legendre' which is specialized for Legendre polynomials.

    • 'roots' which defines the polynomial from its roots.

  • power() for the n-th power of any linear operator.

With p1 and p2 two polynomial instances return by xpoly(coef1, kind1)() and xpoly(coef2, kind2)(), one can:

  • add/substract (if kind1=kind2: (p1 + p2)(L), (p1 - p2)(L) with L the polynomial variable). Evaluating and applying the polynomials on the fly is also possible: (p1 + p2)(L) @ x.

  • The same is possible to multiply (*), divide (//) and modulo (%) two polynomials ((p1 * p2)(L), (p1 // p2)(L), (p1 % p2)(L).

  • And compose two polynomials: (p1(p2))(L) whatever kind1 and kind2.

  • Of note, matrix representation of instance L of LazyLinOp must be square.

  • If p is in monomial form evaluation of p(L) @ x uses the Horner’s method, Clenshaw algorithm otherwise.

  • Of note, duration of (p1(p2))(L) @ x and duration of p1(p2(L)) @ x might differ when the number of polynomial coefficients is large.

More details about implementation and features

The xpoly() returns an instance of internal classes that extend NumPy polynomial classes Polynomial, Chebyshev, Hermite, HermiteE, Laguerre, Legendre.

They override the method __call__() to implement the polynomial evaluation and calculate on the fly the available operations. Under the hood evaluation is called depending on the polynomial form.

Polynomial classes

class lazylinop.polynomial.Polynomial(coef, domain=[-1.0, 1.0], window=[-1.0, 1.0], symbol='x', fromRoots=False)

This class implements a polynomial class derived from numpy.polynomial.Polynomial and so relies on NumPy polynomial package to manipulate polynomials.

See lazylinop.polynomial for an introduction to implemented operations and their basic use.

Init instance of Poly.

Be aware that Poly(fromRoots=True)(L) returns \(p(L)=c_n\prod_{i=0}^n(L-r_iId)\) and not \(p(x)=\prod_{i=0}^n(x-r_i)\) like NumPy polyvalfromroots.

Args:
coef: list

List of coefficients \(\lbrack c_0,c_1,\cdots,c_n\rbrack\) if fromRoots=False. Polynomial of LazyLinOp L is:

\[\begin{equation} p(L)=\sum_{i=0}^nc_iL^i \end{equation}\]

List of roots \(\lbrack r_0,r_1,\cdots,r_n,c_n\rbrack\) (where c_n is the leading coefficients) if fromRoots=True. Polynomial of LazyLinOp L is:

\[\begin{equation} p(L)=c_n\prod_{i=0}^n(L-r_iId) \end{equation}\]
fromRoots: bool, optional

  • If False uses polynomial coefficients.

  • If True uses polynomial roots. Last element coef[-1] of coef is the leading coefficient.

Examples:
>>> import numpy as np
>>> import lazylinop as lz
>>> p1 = lz.polynomial.Polynomial([1, 2, 1], fromRoots=False)
>>> p2 = lz.polynomial.Polynomial([-1, -1, 1], fromRoots=True)
>>> np.allclose(p1.coef, p2.coef)
True

See also

numpy.polynomial package.

class lazylinop.polynomial.Chebyshev(coef, domain=[-1.0, 1.0], window=[-1.0, 1.0], symbol='x')

This class implements a Chebyshev polynomial class derived from numpy.polynomial.Chebyshev and so relies on NumPy polynomial package to manipulate polynomials.

See lazylinop.polynomial for an introduction to implemented operations and their basic use.

Init instance of Chebyshev.

Args:
coef: list

List of Chebyshev coefficients \(\lbrack c_0,c_1,\cdots,c_n\rbrack\). Polynomial of LazyLinOp L is:

\[\begin{equation} p(L)=\sum_{i=0}^nc_iT_i(L) \end{equation}\]
Examples:
>>> from lazylinop.polynomial import Chebyshev
>>> t = Chebyshev([1.0, 2.0, 3.0])
>>> (t + t).coef
array([2., 4., 6.])

See also

numpy.polynomial package.

class lazylinop.polynomial.Hermite(coef, domain=[-1.0, 1.0], window=[-1.0, 1.0], symbol='x')

This class implements a Hermite “physicist” polynomial class derived from numpy.polynomial.Hermite and so relies on NumPy polynomial package to manipulate polynomials.

See lazylinop.polynomial for an introduction to implemented operations and their basic use.

Init instance of Hermite.

Args:
coef: list

List of Hermite coefficients \(\lbrack c_0,c_1,\cdots,c_n\rbrack\). Polynomial of LazyLinOp L is:

\[\begin{equation} p(L)=\sum_{i=0}^nc_iH_i(L) \end{equation}\]
Examples:
>>> from lazylinop.polynomial import Hermite
>>> h = Hermite([1.0, 2.0, 3.0])
>>> (h + h).coef
array([2., 4., 6.])

See also

numpy.polynomial package.

class lazylinop.polynomial.HermiteE(coef, domain=[-1.0, 1.0], window=[-1.0, 1.0], symbol='x')

This class implements a Hermite “probabilist” polynomial class derived from numpy.polynomial.HermiteE and so relies on NumPy polynomial package to manipulate polynomials.

See lazylinop.polynomial for an introduction to implemented operations and their basic use.

Init instance of HermiteE.

Args:
coef: list

List of Hermite coefficients \(\lbrack c_0,c_1,\cdots,c_n\rbrack\). Polynomial of LazyLinOp L is:

\[\begin{equation} p(L)=\sum_{i=0}^nc_iH_i(L) \end{equation}\]
Examples:
>>> from lazylinop import eye, islazylinop
>>> from lazylinop.polynomial import HermiteE
>>> h = HermiteE([1.0, 2.0, 3.0])
>>> (h + h).coef
array([2., 4., 6.])

See also

numpy.polynomial package.

class lazylinop.polynomial.Laguerre(coef, domain=[-1.0, 1.0], window=[-1.0, 1.0], symbol='x')

This class implements a Laguerre polynomial class derived from numpy.polynomial.Laguerre and so relies on NumPy polynomial package to manipulate polynomials.

See lazylinop.polynomial for an introduction to implemented operations and their basic use.

Init instance of Laguerre.

Args:
coef: list

List of Laguerre coefficients \(\lbrack c_0,c_1,\cdots,c_n\rbrack\). Polynomial of LazyLinOp L is:

\[\begin{equation} p(L)=\sum_{i=0}^nc_iL_{a,i}(L) \end{equation}\]
Examples:
>>> from lazylinop import eye, islazylinop
>>> from lazylinop.polynomial import Laguerre
>>> la = Laguerre([1.0, 2.0, 3.0])
>>> (la + la).coef
array([2., 4., 6.])
>>> L = eye(3, n=3, k=0)
>>> islazylinop(la(L))
True

See also

numpy.polynomial package.

class lazylinop.polynomial.Legendre(coef, domain=[-1.0, 1.0], window=[-1.0, 1.0], symbol='x')

This class implements a Legendre polynomial class derived from numpy.polynomial.Legendre and so relies on NumPy polynomial package to manipulate polynomials.

See lazylinop.polynomial for an introduction to implemented operations and their basic use.

Init instance of Legendre.

Args:
coef: list

List of Legendre coefficients \(\lbrack c_0,c_1,\cdots,c_n\rbrack\). Polynomial of LazyLinOp L is:

\[\begin{equation} p(L)=\sum_{i=0}^nc_iL_{e,i}(L) \end{equation}\]
Examples:
>>> from lazylinop import eye, islazylinop
>>> from lazylinop.polynomial import Legendre
>>> le = Legendre([1.0, 2.0, 3.0])
>>> (le + le).coef
array([2., 4., 6.])
>>> L = eye(3, n=3, k=0)
>>> islazylinop(le(L))
True

See also

numpy.polynomial package.

Polynomial functions

lazylinop.polynomial.xpoly(coef, domain=[-1.0, 1.0], window=[-1.0, 1.0], symbol='x', kind='monomial')

Return instance amongst numpy.polynomial.Polynomial, numpy.polynomial.Chebyshev, numpy.polynomial.Hermite, numpy.polynomial.HermiteE, numpy.polynomial.Laguerre or numpy.polynomial.Legendre according to kind. xpoly is the extended function for polynomial creation of any kind without using the specialized polynomial classes directly. Under the hood, the function create instance depending on the kind 'chebyshev', 'hermite', 'hermite_e', 'laguerre', 'legendre', 'monomial' or 'roots' of the polynomial you ask for. It is pretty simple and you can construct a polynomial as you would do with a NumPy polynomial <https://numpy.org/doc/stable/ reference/routines.polynomials.classes.html>_.

Args:
coef: list or 1d array

List of coefficients \(\lbrack c_0,c_1,\cdots,c_n\rbrack\) if kind!='roots':

\[\begin{equation} p(L)=\sum_{i=0}^nc_iQ_i(L) \end{equation}\]

List of roots \(\lbrack r_0,r_1,\cdots,r_n,c_n\rbrack\) (where c_n is the leading coefficients) if kind='roots':

\[\begin{equation} p(L)=c_n\prod_{i=0}^n(L-r_iId) \end{equation}\]
kind: str, optional

Representation of the polynomial. It could be ‘monomial’ (default), ‘chebyshev’, ‘hermite’ physicists polynomials, ‘hermite_e’ probabilists polynomials, ‘laguerre’, ‘legendre’ and ‘roots’. If kind is ‘roots’, coef is considered to be the roots of the polynomial. Leading coefficient is the last element coef[:-1] of coef argument while the first values are the roots of the polynomial. Because of the expression \((L - r_0Id)\cdots (L - r_nId)\) coefficient \(c_n\) of the highest power \(c_nL^n\) is always 1.

Raises:
ValueError

coef size must be > 0.

Exception

coef must be a 1d array.

ValueError

kind must be either monomial, cheb, herm, herme, lag, leg or roots.

Examples:
>>> import numpy as np
>>> import lazylinop as lz
>>> p1 = lz.polynomial.xpoly([1.0, 2.0, 1.0], kind='monomial')
>>> p2 = lz.polynomial.xpoly([-1.0, -1.0, 1.0], kind='roots')
>>> np.allclose(p1.coef, p2.coef)
True

See also

numpy.polynomial package.

lazylinop.polynomial.power(L, n)

Constructs the n-th power \(L^n\) of linear operator \(L\). Matrix representation of \(L\) must be square. In some cases power(L,n) @ x can be least efficient than M=np.power(L.toarray()) @ x.

Note

It is equivalent to create an instance from xpoly(coef, kind)() such that only n-th coefficient is equal to one while the others are equal to zero.

Args:
L: LazyLinOp

Linear operator (e.g. a LazyLinOp). Matrix representation must be square.

n: int

Raise the linear operator to degree n. If n is zero, return identity matrix.

Returns:

LazyLinOp \(L^n\).

Raises:
ValueError

n must be > 0.

Exception

Matrix representation of L is not square.

Examples:
>>> import numpy as np
>>> import lazylinop as lz
>>> L = lz.polynomial.power(lz.eye(3, n=3, k=0), 3)
>>> x = np.full(3, 1.0)
>>> np.allclose(L @ x, x)
True
>>> L = lz.polynomial.power(lz.eye(3, n=3, k=1), 3)
>>> # Note that L is in fact zero (nilpotent matrix)
>>> x = np.full(3, 1.0)
>>> np.allclose(L @ x, np.zeros(3, dtype=np.float_))
True

See also

NumPy power function.