pyfaust.poly Namespace Reference

The pyfaust module for polynomial basis as Faust objects. More...

Classes
class	FaustPoly
	Subclass of Faust specialized for orthogonal polynomial basis. More...

Functions
def	basis (L, K, basis_name, dev='cpu', T0=None, **kwargs)
	Builds the Faust of the polynomial basis defined on the sparse matrix L. More...
def	poly (coeffs, basis='chebyshev', L=None, X=None, dev='cpu', out=None, **kwargs)
	Computes the linear combination of the polynomials defined by basis. More...
def	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. More...

Detailed Description

The pyfaust module for polynomial basis as Faust objects.

Note

This module is still in BETA status.

Function Documentation

basis()

def 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.

Parameters

L	the sparse scipy square matrix in CSR format (scipy.sparse.csr_matrix). The dtype must be float32, float64 or complex128 (the dtype might have a large impact on performance).
K	the degree of the last polynomial, i.e. the K+1 first polynomials are built.
basis_name	'chebyshev', and others yet to come.
dev	the device to instantiate the returned Faust ('cpu' or 'gpu').
T0	a sparse matrix to replace the identity as a 0-degree polynomial of the basis.

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.

Examples

>>> 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

expm_multiply()

def 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.

Parameters

A	the operator whose exponential is of interest (must be a symmetric positive definite csr_matrix).
B	(ndarray) the matrix or vector to be multiplied by the matrix exponential of A.
t	(list) the time points.
dev	(str, optional) the device ('cpu' or 'gpu') on which to compute (currently only cpu is supported).
K	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	'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\).

Examples

>>> 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]])

poly()

def 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.

Parameters

coeffs	the linear combination coefficients (vector as a numpy.ndarray). Must be in the same dtype as L and X if they are set.
basis	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	if X is not None, basis can only be a FaustPoly).
L	the sparse scipy square matrix in CSR format (scipy.sparse.csr_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.

Parameters

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	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.

Examples

>>> 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