-  3.39.20
Classes | Functions
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
npif 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
pyfaust.seed
def seed(s)
(Re)Initializes the pyfaust pseudo-random generator.
Definition: __init__.py:5268
matfaust::poly::next
function next(F)
Gives the next Faust basis of dimension (n+1) from the Faust F polynomial basis of dimension n.
pyfaust.poly.basis
def basis(L, K, basis_name, dev='cpu', T0=None, **kwargs)
Builds the Faust of the polynomial basis defined on the sparse matrix L.
Definition: poly.py:201
pyfaust.poly
The pyfaust module for polynomial basis as Faust objects.
Definition: poly.py:1
pyfaust.poly.poly
def poly(coeffs, basis='chebyshev', L=None, X=None, dev='cpu', out=None, **kwargs)
Computes the linear combination of the polynomials defined by basis.
Definition: poly.py:302