pyfaust.poly module

The module for pyfaust Chebyshev polynomials.

pyfaust.poly.Chebyshev(L, K, dev='cpu', T0=None, impl='native')

Builds the Faust of the Chebyshev polynomial basis defined on the sparse matrix L.

Args:

L: (scipy.sparse.csr_matrix or Faust)

the sparse square matrix or square Faust (in the last case impl must be ‘py’)

K: (int)

the degree of the last polynomial, i.e. the K+1 first polynomials are built.

dev: (str)

the device to instantiate the returned Faust (‘cpu’ or ‘gpu’).

T0: (Faust or scipy.sparse.csr_matrix)

to define the 0-degree polynomial as something else than the identity. If a Faust only the first factor is take into account.

impl: (str)

‘native’ (by default) for the C++ impl., “py” for the Python impl.

Returns:

The Faust of the K+1 Chebyshev polynomials.

See pyfaust.poly.basis which is pretty the same (e.g.: calling Chebyshev(L, K) is equivalent to basis(L, K, ‘chebyshev’)

class pyfaust.poly.FaustPoly(*args, **kwargs)

Subclass of Faust specialized for orthogonal polynomial basis.

This class is used only for the native implementation of the poly functions.

Note

it is not advisable to use this class directly.

pyfaust.poly.basis(L, K, basis_name, dev='cpu', T0=None, **kwargs)

Builds the Faust of the polynomial basis defined on the sparse matrix L.

Args:

L: (scipy.sparse.csr_matrix)

the sparse square matrix. The dtype must be float32, float64 or complex128 (the dtype might have a large impact on performance).

K: (int)

the degree of the last polynomial, i.e. the K+1 first polynomials are built.

basis_name:

‘chebyshev’, and others yet to come.

dev: (str)

the device to instantiate the returned Faust (‘cpu’ or ‘gpu’).

T0: (Faust or scipy.sparse.csr_matrix)

to define the 0-degree polynomial as something else than the identity. If a Faust only the first factor is take into account.

Returns:

The Faust G of the basis composed of the K+1 orthogonal polynomials. Note that the Faust returned is also a generator: calling next(G) will return the basis of dimension K+1.

Example:

>>> from pyfaust.poly import basis
>>> from scipy.sparse import random
>>> from numpy.random import seed
>>> seed(42) # for reproducibility
>>> L = random(50, 50, .02, format='csr')
>>> L = L@L.T
>>> K = 3
>>> F = basis(L, K, 'chebyshev')
>>> F
Faust size 200x50, density 0.0699, nnz_sum 699, 4 factor(s):
- FACTOR 0 (double) SPARSE, size 200x150, density 0.00943333, nnz 283
- FACTOR 1 (double) SPARSE, size 150x100, density 0.0155333, nnz 233
- FACTOR 2 (double) SPARSE, size 100x50, density 0.0266, nnz 133
- FACTOR 3 (double) SPARSE, size 50x50, density 0.02, nnz 50
 identity matrix flag

Generate the next basis (the one with one additional dimension, whose the polynomial greatest degree is K+1 = 4)

>>> G = next(F)
>>> G
Faust size 250x50, density 0.08256, nnz_sum 1032, 5 factor(s):
- FACTOR 0 (double) SPARSE, size 250x200, density 0.00666, nnz 333
- FACTOR 1 (double) SPARSE, size 200x150, density 0.00943333, nnz 283
- FACTOR 2 (double) SPARSE, size 150x100, density 0.0155333, nnz 233
- FACTOR 3 (double) SPARSE, size 100x50, density 0.0266, nnz 133
- FACTOR 4 (double) SPARSE, size 50x50, density 0.02, nnz 50
 identity matrix flag

The factors 0 to 3 of G are views of the same factors of F.

By default, the 0-degree polynomial is the identity. However it is possible to replace the corresponding matrix by any csr sparse matrix T0 of your choice (with the only constraint that T0.shape[0] == L.shape[0]. In that purpose, do as follows:

>>> F2 = basis(L, K, 'chebyshev', T0=random(50,2, .3, format='csr'))
>>> F2
Faust size 200x2, density 1.7475, nnz_sum 699, 4 factor(s):
- FACTOR 0 (double) SPARSE, size 200x150, density 0.00943333, nnz 283
- FACTOR 1 (double) SPARSE, size 150x100, density 0.0155333, nnz 233
- FACTOR 2 (double) SPARSE, size 100x50, density 0.0266, nnz 133
- FACTOR 3 (double) SPARSE, size 50x2, density 0.5, nnz 50

pyfaust.poly.expm_multiply(A, B, t, K=10, tradeoff='time', dev='cpu', **kwargs)

Computes an approximate of the action of the matrix exponential of A on B using series of Chebyshev polynomials.

Note

This function is very similar to scipy.sparse.linalg.expm_multiply

with three major differences though:

1. A must be symmetric positive definite.

2. The time points are directly passed to the function rather to be defined as a numpy.linspace.

3. The time values must be negative.

Args:

A: (scipy.sparse.csr_matrix)

the operator whose exponential is of interest (must be a symmetric positive definite matrix).

B: (np.ndarray)

the matrix or vector to be multiplied by the matrix exponential of A.

t: (list[float])

the time points.

dev: (str)

the device (‘cpu’ or ‘gpu’) on which to compute (currently only cpu is supported).

K: (int)

the greatest polynomial degree of the Chebyshev polynomial basis. The greater it is, the better is the approximate accuracy but note that a larger K increases the computational cost.

tradeoff: (str)

‘memory’ or ‘time’ to specify what matters the most: a small memory footprint or a small time of execution. It changes the implementation of pyfaust.poly.poly used behind. It can help when the memory size is limited relatively to the value of K or the size of A and B.

Returns:

expm_A_B the result of \(e^{t_k A} B\).

Example:

>>> import numpy as np
>>> from scipy.sparse import random
>>> from pyfaust.poly import expm_multiply as fexpm_multiply
>>> np.random.seed(42) # for reproducibility
>>> L = random(5, 5, .2, format='csr')
>>> L = L@L.T
>>> x = np.random.rand(L.shape[1])
>>> t = np.linspace(start=-.5, stop=-0.1, num=3, endpoint=True)
>>> y = fexpm_multiply(L, x, t)
>>> y
array([[0.12706927, 0.21111065, 0.36636184, 0.41736344, 0.78517596],
       [0.13190047, 0.24040598, 0.36636184, 0.4324354 , 0.78517596],
       [0.13691536, 0.27376655, 0.36636184, 0.44805164, 0.78517596]])

pyfaust.poly.poly(coeffs, basis='chebyshev', L=None, X=None, dev='cpu', out=None, **kwargs)

Computes the linear combination of the polynomials defined by basis.

Args:

coeffs: (np.ndarray)

the linear combination coefficients (vector). Must be in the same dtype as L and X if they are set.

basis: (FaustPoly or np.ndarray)

either the name of the polynomial basis to build on L or the basis if already built externally (as a FaustPoly or an equivalent np.ndarray – if X is not None, basis can only be a FaustPoly).

L: (scipy.sparse.csr_matrix)

The square matrix on which the polynomial basis is built if basis is not already a Faust or a numpy.ndarray. The dtype must be float32, float64 or complex128 (the dtype might have a large impact on performance). It can’t be None if basis is not a FaustPoly or a numpy.ndarray.

X: (np.darray)

if X is not None, the linear combination of basis@X is computed (note that the memory space is optimized compared to the manual way of doing first B = basis@X and then calling poly on B with X at None).

dev: (str)

the computing device (‘cpu’ or ‘gpu’).

out: (np.ndarray)

if not None the function result is put into this np.ndarray. Note that out.flags[‘F_CONTINUOUS’] must be True. Note that this can’t work if the function returns a Faust.

Returns:

The linear combination Faust or np.ndarray depending on if basis is itself a Faust or a np.ndarray. If X is set the result is always a np.ndarray whatever is the basis type.

Example:

>>> import numpy as np
>>> from pyfaust.poly import basis, poly
>>> from scipy.sparse import random
>>> np.random.seed(42) # for reproducibility
>>> L = random(50, 50, .02, format='csr')
>>> L = L@L.T
>>> coeffs = np.array([.5, 1, 2, 3])
>>> G = poly(coeffs, 'chebyshev', L)
>>> G
Faust size 50x50, density 0.3596, nnz_sum 899, 5 factor(s):
- FACTOR 0 (double) SPARSE, size 50x200, density 0.02, nnz 200
- FACTOR 1 (double) SPARSE, size 200x150, density 0.00943333, nnz 283
- FACTOR 2 (double) SPARSE, size 150x100, density 0.0155333, nnz 233
- FACTOR 3 (double) SPARSE, size 100x50, density 0.0266, nnz 133
- FACTOR 4 (double) SPARSE, size 50x50, density 0.02, nnz 50
 identity matrix flag

Which is equivalent to do as below (in two times):

>>> K = 3
>>> F = basis(L, K, 'chebyshev')
>>> coeffs = np.array([.5, 1, 2, 3])
>>> G = poly(coeffs, F)
>>> G
Faust size 50x50, density 0.3596, nnz_sum 899, 5 factor(s):
- FACTOR 0 (double) SPARSE, size 50x200, density 0.02, nnz 200
- FACTOR 1 (double) SPARSE, size 200x150, density 0.00943333, nnz 283
- FACTOR 2 (double) SPARSE, size 150x100, density 0.0155333, nnz 233
- FACTOR 3 (double) SPARSE, size 100x50, density 0.0266, nnz 133
- FACTOR 4 (double) SPARSE, size 50x50, density 0.02, nnz 50
 identity matrix flag

Above G is a Faust because F is too. Below the full array of the Faust F is passed, so an array is returned into GA. >>> GA = poly(coeffs, F.toarray()) >>> type(GA) <class ‘numpy.ndarray’>

But of course they are equal:

>>> np.allclose(GA, G.toarray())
True