Using the poly module

A new module has been added to pyfaust version 3.1.x. Its name is poly and as expected it is dedicated to a kind of Fausts that are defined according to series of polynomials.

In this notebook we'll see how to use the main functions of this module then we'll introduce one precise use case with the action of exponential matrix on a vector / matrix.

NOTE: all the functions introduced in this notebook are available on GPU, using the dev='gpu' argument.

1. The basis function

Firstly, the poly module allows to define a polynomial basis (a FaustPoly) using the function pyfaust.poly.basis. Currently, only Chebyshev polynomials are supported but other are yet to come. Below is the prototype of the function:

basis(L, K, basis_name, dev='cpu', T0=None)

In the next, we shall see a simple example but I let you consult the documentation by typing help(pyfaust.poly.basis) to get more details.

For instance, if you want to instantiate a basis Faust of dimension K+1 (in the example below K=5) on L, which by the way must be a scipy.sparse.csr_matrix at least square (and most likely symmetric positive definite in much use cases), you'll make this call to the function:

As you can see, the last factor is followed by the mention identity matrix flag. It means that this factor is the identity matrix. This is not suprising, because the 0-degree Chebyshev polynomial is the identity. However, note that the T0 optional argument of the function is here to trick the basis by using another matrix than the identity even if eventually it might not be a proper basis it can be useful if you want to apply this basis on a vector or a matrix (hence you'll set T0 as this vector/matrix instead of multiplying the basis by this vector/matrix).

So how should we understand this Faust? You can see it as a vertical concatenation of polynomials. Indeed, the basis is composed of K+1 polynomials, the 0-degree polynomial is at the top of the stack (i.e. F[:d,:] is the identity):

This first 0-degree polynomial is followed by the next degree polynomials: hence F[d:d*2, :] is the 1-degree polynomial, F[d*2:d*3, :] is the 2-degree polynomial and so on...

For details about the Chebyshev polynomials, including their definition by a recurrence relationship (that is used here behind the scene), you can look at this wikipedia article.

One of the most important thing to note about a polynomial basis Faust is that the multiplication by a vector or a matrix is specialized in order to optimize the performance obtained compared to the generic Faust-vector/matrix multiplication. Indeed, due to the particular structure of the polynomial basis Faust, the multiplication can be optimized.

Let's verify that is true! In the code below F is cloned to a classic Faust G and the time of the multiplication by a matrix is measured in the both cases (with F, the basis Faust and G its classic Faust copy). Note that Faust.clone function is not used because in this case it would return a polynomial basis Faust (after all that's the role of clone to preserve the Faust properties!). To get a classic Faust the only way is to copy the Faust factor by factor.

Now let's verify the multiplication result is accurate:

As you see it's alright.

2. The poly function

The second function of the pyfaust.poly module is poly. This function purpose is to compute a linear combination of polynomials composing a FaustPoly (it can also be viewed as a way to compute a polynomial). So passing the FaustPoly and the linear combination coefficients (one per polynomial, in ascending degree order) you'll obtain:

To be explicit about lc_F let's show how to rebuild it manually using G (which again is a classic Faust equal to F).

Here again the FaustPoly operation is optimized compared to the Faust one. Speaking of which, there is ways to do even more optimized because the poly function is kind of matrix type agnostic, or precisely, it accepts a FaustPoly or a numpy.ndarray as the basis argument. Doing with the latter an optimized implementation is used whose the memory footprint is smaller than the one consumed with a FaustPoly. It can be particulary efficient when the use cases (as we'll see in 3.) consist to apply a linear combination of F to a vector x as it's shown below.

Depending on the situation the way2 or way3 might be the quickest but they are always both quicker than way1.

Just in case let's verify all ways give the same results.

All sounds good! We shall now introduce one use case of Chebyshev polynomial series in FAµST that allow to get interesting results compared to what we can do in the numpy/scipy ecosystem. But to note a last thing before going ahead to the part 3. is that the function poly is a little more complicated that it looks like, for more details I invite you to consult the API documentation.

3. Computing the action of the exponential of a matrix on a vector

In the next, we'll see how to compute action of the exponential matrix on a vector x. However this time we'll do the comparison with the scipy function scipy.sparse.linalg.expm_multiply. The both functions are intended to compute the action of the exponential matrix on a vector or matrix. Recalling that it consists to compute $exp(t A)x$ without computing directly the exponential let's compare the use, performance and accuracy of these functions.
The main difference between the two of them, is that in pyfaust the time points are passed directly as a list to the function, while the scipy version accepts only on np.linspace style arguments (to define the points as a regular space).

It is rather good for pyfaust but note that there is some drawbacks to its expm_multiply implementation. You'll find them among other details in the API documentation.

Thanks for reading this notebook! Many other are available at

Note: this notebook was executed using the following pyfaust version: