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

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.0672, nnz_sum 672, 4 factor(s):

• FACTOR 0 (double) SPARSE, size 200x150, density 0.00913333, nnz 274
• FACTOR 1 (double) SPARSE, size 150x100, density 0.0149333, nnz 224
• FACTOR 2 (double) SPARSE, size 100x50, density 0.0248, nnz 124
• 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.07968, nnz_sum 996, 5 factor(s):

• FACTOR 0 (double) SPARSE, size 250x200, density 0.00648, nnz 324
• FACTOR 1 (double) SPARSE, size 200x150, density 0.00913333, nnz 274
• FACTOR 2 (double) SPARSE, size 150x100, density 0.0149333, nnz 224
• FACTOR 3 (double) SPARSE, size 100x50, density 0.0248, nnz 124
• 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.68, nnz_sum 672, 4 factor(s):

• FACTOR 0 (double) SPARSE, size 200x150, density 0.00913333, nnz 274
• FACTOR 1 (double) SPARSE, size 150x100, density 0.0149333, nnz 224
• FACTOR 2 (double) SPARSE, size 100x50, density 0.0248, nnz 124
• 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	(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\).

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.24591891, 0.24529726, 0.35223879, 0.7409274 ],
[0.13190047, 0.26391877, 0.2890644 , 0.38977573, 0.75815737],
[0.13691536, 0.28256865, 0.33890501, 0.43252159, 0.77600955]])

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

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.3488, nnz_sum 872, 5 factor(s):

• FACTOR 0 (double) SPARSE, size 50x200, density 0.02, nnz 200
• FACTOR 1 (double) SPARSE, size 200x150, density 0.00913333, nnz 274
• FACTOR 2 (double) SPARSE, size 150x100, density 0.0149333, nnz 224
• FACTOR 3 (double) SPARSE, size 100x50, density 0.0248, nnz 124
• 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.3488, nnz_sum 872, 5 factor(s):

• FACTOR 0 (double) SPARSE, size 50x200, density 0.02, nnz 200
• FACTOR 1 (double) SPARSE, size 200x150, density 0.00913333, nnz 274
• FACTOR 2 (double) SPARSE, size 150x100, density 0.0149333, nnz 224
• FACTOR 3 (double) SPARSE, size 100x50, density 0.0248, nnz 124
• 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