Using the polynomial module
In this notebook we propose a simple tutorial on how to use the Lazylinop polynomial module.
First of all, you need to install Lazylinop if not done already, see the installation guide. If you’re not yet introduced to the basic use of Lazylinop, we invite you to jump into the library by following this quick start guide.
0. Needed imports
The notebook relies of course on Lazylinop and additionally on NumPy and SciPy. So let’s make the imports.
[1]:
import lazylinop as lz
from lazylinop.polynomial import Polynomial, Chebyshev
import numpy as np
import scipy as sp
1. Introduction
The polynomial module provides a set of classes and functions dedicated to the manipulation of polynomials able to act on linear operators.
Mathematically, given a polynomial \(P(x) = c_0 + c_1x+ \cdots + c_n x^n\) and any square matrix \(M\) one can define \(P(M) = c_0 Id + c_1 M + \cdots c_n M^n\). Convenient polynomial classes are implemented in NumPy, yet they cannot be readily applied to NumPy arrays or SciPy LinearOperators. The polynomial module fixes this by implementing the needed classes, which allow to define and manipulate \(K = P\left(L\right)\), where \(L\) is any square LazyLinOp, and the
resulting \(K\) is again a LazyLinOp.
Relationship with NumPy polynomials
polynomial module actually relies a lot on NumPy representation of polynomials but offers the support of LazyLinOp as indeterminates. As in NumPy, it allows to perform standard operations on polynomials (e.g. p1 + p2) but also to lazily evaluate a polynomial on an indeterminate (naturally written P(L) in the code). Pretty much as simple as what we are used to write basically in maths. Note however that contrary to NumPy’s the Lazylinop polynomials are defined according the
matrix power, that is, L**n = L @ L @ ... @ L (with L a LazyLinOp) while NumPy’s are defined on array elementwise power A**n = A * A * ... * A (with A a NumPy array).LazyLinOp objects of dimensions \(1 \times 1\). Indeed, it is possible to do so but to the best of our knowledge NumPy is faster for that use case.Lazy polynomials use cases & Quick API tour
Of course many use cases exist. An example is the approximation of an operator exponential \(e^{L}\approx c_0Id + c_1L + c_2L^2 + \cdots + c_nL^n\) where the \(c_i\) are the coefficients of the polynomial. The main interest of this approximation is to spare compute time. In section 5. we say more about this kind of approximation. More generally, the polynomial module follows the lazy and sparse computing way of Lazylinop. You can compute the result of \(P\left(L\right)\) applied to a batch \(X\) of column-vectors without ever explicitly computing the matrix associated to \(P(L)\). This saves memory, and delays all computations to when \(P\left(L\right) @ X\) is actually evaluated.
A bunch of functions and objects are hence available to compute lazily polynomials, including the most important listed here with their API documentations:
For general polynomials:
Polynomial inherited from NumPy Polynomial,
xpoly is the function for creating all the Lazylinop polynomials listed above or available in the module. The notebook does not show use case of
xpolyfunction. Therefore, please refer to the documentation for more details.
For the evaluation of polynomials in monomial form:
Polynomial(coef)(L) @ xorxpoly(coef, kind='monomial')(L) @ x
For the evaluation of polynomials in monomial form from roots:
Polynomial([roots, leading_coef], fromRoots=True)(L) @ xorxpoly([roots, leading_coef], kind='roots')(L) @ x
The notebook does not show use case of form='roots'. Therefore, please refer to the documentation for more details.
For the evaluation of polynomials in Chebyshev form:
Chebyshev(coef)(L) @ xorxpoly(coef, kind='chebyshev')(L) @ x
For the exponentiation:
power or LazyLinOp.__pow__ (aka
**).
Those APIs are demonstrated in the next of this notebook. In particular, the polynomial available operations are shown more smoothly than in the API doc.
Many examples will be given in this notebook. From the simplest to more advanced ones.
2. N-th power of an operator
The most simple syntax to create a polynomial is p = Polynomial([c_0,..., c_n]) which creates a polynomial instance such that \(p(L) = c_0 Id + c_1 L + \cdots + c_n L^n\) for any LazyLinOp \(L\). Hence, the polynomial instances p1 and p2 below are such that \(p_1(L) = L\) and \(p_2(L) = L^2\). Let us check that \((p_1^2)(L)\), \(p_2(L)\) and \(L@L\) all implement the same LazyLinOp: applying them to any vector \(x\) yields the same result. For
simplicity we illustrate this for a simple LazyLinOp L associated to a square matrix full of ones, and to a random vector x.
[2]:
p1 = Polynomial([0.0, 1.0])
print("p1 is", p1)
p2 = Polynomial([0.0, 0.0, 1.0])
print("p2 is", p2)
print("p1 and p2 are instances of:", p2.__class__)
N = 5000
L = lz.ones((N, N)) # L is a LazyLinOp built by lazylinop.ones
M = np.ones((N, N)) # NumPy array of ones
x = np.random.randn(N)
print("(A) lazylinop time for (p1 ** 2)(L) @ x:")
%timeit -r 100 -n 10 (p1 ** 2)(L) @ x
print("(B) lazylinop time for p2(L) @ x:")
%timeit -r 100 -n 10 p2(L) @ x
print("(C) numpy time for M @ (M @ x):")
%timeit -r 100 -n 10 M @ (M @ x)
p1 is 0.0 + 1.0·x
p2 is 0.0 + 0.0·x + 1.0·x²
p1 and p2 are instances of: <class 'lazylinop.polynomial.Polynomial'>
(A) lazylinop time for (p1 ** 2)(L) @ x:
4.85 ms ± 501 µs per loop (mean ± std. dev. of 100 runs, 10 loops each)
(B) lazylinop time for p2(L) @ x:
4.77 ms ± 138 µs per loop (mean ± std. dev. of 100 runs, 10 loops each)
(C) numpy time for M @ (M @ x):
10.4 ms ± 934 µs per loop (mean ± std. dev. of 100 runs, 10 loops each)
The result of the following quick benchmark shows that in this specific example the polynomial of a LazyLinOp ones is faster than NumPy’s. The reason why this is faster is actually not because of the instance of Poly class but rather the Lazylinop ones operator. However the two instances p1 and p2 of Polynomial class allow us to use this operator as a variable, which is needed to
perform the achieved speedup.
3. Basics operations
We have already seen that given \(p_1\) it is possible to compute the n-th power \(p_1^n\) via the usual ** syntax. Other simple operations are the addition, subtraction, or composition of polynomials. Please have a look to the following lines of code that build a polynomial p5, apply it to a linear operator L and apply the resulting operator p5(L) to a vector x.
[3]:
p1 = Polynomial([1.0, -2.0]) # First polynomial is 1 - 2 * x
p2 = Polynomial([0.0, 0.0, 1.0]) # Second polynomial is x^2
p3 = Polynomial([-1.0, 1.0]) # Third polynomial is x - 1
p4 = Polynomial([0.0, 0.0, 0.0, 0.0, 2.0]) # Fourth polynomial is 2 * x^4
print("p1 instance is", p1)
print("p2 instance is", p2)
print("p3 instance is", p3)
print("p4 instance is", p4)
p5 = 2 * (p1 + p2) ** 2 - p4(p3)
print("p5 instance is", p5)
N = 1000
L = lz.ones((N, N))
x = np.random.randn(N)
print("p5(L) @ x is all zero:", np.allclose(0.0, p5(L) @ x), "(as expected).")
p1 instance is 1.0 - 2.0·x
p2 instance is 0.0 + 0.0·x + 1.0·x²
p3 instance is -1.0 + 1.0·x
p4 instance is 0.0 + 0.0·x + 0.0·x² + 0.0·x³ + 2.0·x⁴
p5 instance is 0.0
p5(L) @ x is all zero: True (as expected).
Indeed p5(L) == 0 for any L because p4(p3)(L) == 2 * (L - 1)**4 and 2 * ((p1 + p2)**2)(L) == 2 * (L**2 - 2 * L - 1)**2 == 2 * (L - 1)**4.
Now let’s experiment the composition p = p1(t1) of two polynomials p1 and t1. The evaluation p(L) @ x can be written either as (p1(t1))(L) @ x or p1(t1(L)) @ x. It might be possible that the two evaluations are not equivalent in term of computation time even if the results are the same as you can see in this example:
[4]:
p1 = Polynomial(np.random.randn(10)) # First polynomial
t1 = Chebyshev(np.random.randn(100)) # Second polynomial in Chebyshev form
N = 10
L = lz.kron(
(1.0 / N ** 2) * np.random.randn(N, N),
(1.0 / N ** 2) * np.random.randn(N, N)
)
x = np.random.randn(L.shape[1])
print("duration of (t1(p1))(L) @ x:")
%timeit -r 40 -n 10 (t1(p1))(L) @ x
print("duration of t1(p1(L)) @ x:")
%timeit -r 40 -n 10 t1(p1(L)) @ x
print("Are the two results the same ?" +
" {0:s}".format(
str(np.allclose(
(t1(p1))(L) @ x,
t1(p1(L)) @ x))))
duration of (t1(p1))(L) @ x:
27.1 ms ± 915 µs per loop (mean ± std. dev. of 40 runs, 10 loops each)
duration of t1(p1(L)) @ x:
21.8 ms ± 573 µs per loop (mean ± std. dev. of 40 runs, 10 loops each)
Are the two results the same ? True
Even though the calculations give the same results in both cases, the computation times are not the same. The evaluation (t1(p1))(L) @ x is slower compared to the evaluation of t1(p1(L)) @ x. That is something to know when you are using Lazylinop polynomials.
4. Changing a polynomial to Chebyshev basis
As you may have noticed in our last example, one can create polynomials with coefficients expressed in different polynomial basis. We have already seen the syntax p = Polynomial([c0,..., c_n]) in the standard representation as a sum of monomials, but one can also use chebyshev, hermite (physicists polynomials), hermite_e (probabilists polynomials), laguerre, legendre and roots. For example, q = Chebyshev([c0,...,c_n]) will create a polynomial such that
\(q(L) = c_0 I + c_1 T_1(L) + \cdots + c_n T_n(L)\) where \(T_i\) is the \(i\)-th Chebyshev polynomial.
Following our previous section we now show how p1 and p2, which are defined according to a monomial basis, can be converted to a representation in the Chebyshev basis using the NumPy syntax convert(kind=Chebyshev). We then apply as before the resulting polynomial of L to x.
[5]:
print(isinstance(p1 + p2, Polynomial))
print(isinstance((p1 + p2).convert(kind=Chebyshev), Chebyshev))
print(isinstance(p1 + p2, Polynomial))
L = lz.ones((x.shape[0], x.shape[0]))
ypoly = (p1 + p2)(L) @ x
ycheb = ((p1 + p2).convert(kind=Chebyshev))(L) @ x
np.linalg.norm(ypoly - ycheb) / np.linalg.norm(ypoly)
True
True
True
[5]:
1.2341363840556378e-16
We maintain the same accuracy using Chebyshev polynomials, all is fine!
5. Approaching the exponential of a Kronecker product using Polynomial
In the above examples, we mainly manipulated p(L) with low-degree polynomials p and a simple LazyLinOp L corresponding to a squared matrix filled with ones. If you found these examples too trivial, you may be interested in the next example: we now use polynomials to efficiently approximate the exponential \(e^{\tau K}\) of a (lazy) Kronecker product \(K=A\otimes B\) using the polynomial module. We will test it and then compute an approximation of exponential action of
K on x (that is \(e^{\tau K} x\)). First we construct the NumPy arrays A,B, the NumPy array K corresponding to their Kronecker product, and the LazyLinOp lK obtained via a lazy Kronecker product.
[8]:
A = np.random.randn(2, 2)
B = np.random.randn(4, 4)
K = np.kron(A, B)
lK = lz.kron(A, B)
Then we apply SciPy’s expm to compute the matrix exponential of K and apply it to a random vector. This will serve as our reference. Of note, it matters to first normalize our operator \(K\) to avoid large eigenvalues associated to numerical issues in the computation of the exponential.
[9]:
norm = np.linalg.norm(K)
lK = lK / norm
K = K / norm
x = np.random.rand(lK.shape[1])
y = sp.linalg.expm(K) @ x
Now, we construct an instance p1 of the Polynomial class corresponding the Taylor approximation \(1+x+\frac{x^2}{2}\) of degree 2 of the exponential function \(e^x\). This allows to approximate the action of the exponential of lK on \(x\) using the simple instruction p1(lK) @ x.
[10]:
p1 = Polynomial([1.0, 1.0, 0.5])
approx = p1(lK) @ x
np.linalg.norm(y - approx) / np.linalg.norm(y)
[10]:
0.0017207308680879968
Well, the approximation is not so good, but we can increase the precision by increasing the degree of the Taylor approximation. To do so, we create a new polynomial gathering the third, fourth and fifth degree terms of the Taylor approximation \(\frac{x^3}{3!}+\frac{x^4}{4!}+\frac{x^5}{5!}+\frac{x^6}{6!}\).
[11]:
p2 = Polynomial([0.0, 0.0, 0.0, 1.0 / 6.0, 1.0 / 24.0, 1.0 / 120.0, 1.0 / 720.0])
approx = (p1 + p2)(lK) @ x
np.linalg.norm(y - approx) / np.linalg.norm(y)
[11]:
2.8228073277124983e-08
This is not perfect but alreay quite better.
6. Benchmarking the exponential of a Kronecker product
We know thanks to the previous section that the polynomial approximation works, now we want to know if it works faster. So, as before we create our Kronecker product \(K=A\otimes B\) of matrices \(A\) and \(B\) as a LazyLinOp lK and as a NumPy array K.
Again, we normalize our operator \(K\) in order to avoid a diverging approximation. The norm should be chosen carefully to reduce the cost. You might also divide by the largest eigenvalue which comes at a lower cost if you use SciPy scipy.sparse.linalg.eigsh(A, k=1).
[12]:
M, N = 16, 16
A = np.random.randn(M, N)
B = np.random.randn(M, N)
K = np.kron(A, B)
norm = np.linalg.norm(K)
lK = 1 / norm * lz.kron(A, B)
K = 1 / norm * np.kron(A, B)
Then we create a polynomial p9 in monomial form with the first ninth coefficients of the Taylor series of the exponential \(e^x\approx\sum_{i=0}^8\frac{x^i}{i!}\).
We check that our result and the one given by SciPy are close.
[13]:
x = np.random.randn(lK.shape[1])
p9 = Polynomial(1.0 / sp.special.factorial(np.arange(9), exact=True))
y = p9(lK) @ x
z1 = sp.linalg.expm(K) @ x
z2 = sp.sparse.linalg.expm_multiply(K, x)
print(np.linalg.norm(y - z1) / np.linalg.norm(z1))
print(np.linalg.norm(y - z2) / np.linalg.norm(z2))
3.9053842693899996e-16
1.7580767299374295e-16
Perfect! Now let’s measure the computation time. We confront our results to scipy.linalg.expm and scipy.sparse.linalg.expm_multiply functions.
[14]:
%timeit -r 10 -n 2 p9(lK) @ x
The slowest run took 4.43 times longer than the fastest. This could mean that an intermediate result is being cached.
979 µs ± 732 µs per loop (mean ± std. dev. of 10 runs, 2 loops each)
[15]:
%timeit -r 10 -n 2 sp.linalg.expm(K) @ x
The slowest run took 5.61 times longer than the fastest. This could mean that an intermediate result is being cached.
24.3 ms ± 22 ms per loop (mean ± std. dev. of 10 runs, 2 loops each)
[16]:
%timeit -r 10 -n 2 sp.sparse.linalg.expm_multiply(K, x)
The slowest run took 9.63 times longer than the fastest. This could mean that an intermediate result is being cached.
1.98 ms ± 1.81 ms per loop (mean ± std. dev. of 10 runs, 2 loops each)
On this very example, it’s clear that Lazylinop approximation gives better performance than SciPy!
It’s noteworthy that of course larger are the matrices better is the performance of the polynomial approximation relatively to SciPy.
You can follow the same idea to compute a Pascal matrix that involves a nilpotent matrix exponential.
Lazylinop provides many similar polynomial approximation functions with arbitrary order of approximation:
Feel free to test and explore and write yours.
Thanks for reading this notebook! Find others on the Lazylinop website.