Mentions légales du service

Skip to content
Snippets Groups Projects
Commit 51c997e8 authored by hhakim's avatar hhakim
Browse files

Update the pyfaust.poly jupyter notebook and add the HTML pre-renredered version.

parent 1fc4ddf2
Branches
Tags
No related merge requests found
Hide whitespace changes
Inline Side-by-side
gen_doc/pyfaust_poly.html 0 → 100644
+ 14914
− 0
View file @ 51c997e8
Source diff could not be displayed: it is too large. Options to address this: view the blob.
gen_doc/pyfaust_poly.ipynb
+ 59
− 233
View file @ 51c997e8
%% Cell type:markdown id:6f8cc8a4 tags:
# Using the poly module
A new module has been added to pyfaust version 3.0.x. Its name is ``poly`` and as expected it is dedicated to a kind of Fausts that are defined according to polynomials.
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 two precise use cases with the action of inverse or exponential matrix on a vector / matrix.
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
**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 (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:
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:
%% Cell type:code id:d077b727 tags:
``` python
from pyfaust.poly import basis
from scipy.sparse import random
d = 128
L = random(d, d, .2, format='csr')
L = L@L.T
K=5
F = basis(L, K=K, basis_name='chebyshev')
F
```
%% Output
Faust size 768x128, density 0.85612, nnz_sum 84160, 6 factor(s):
- FACTOR 0 (real) SPARSE, size 768x640, density 0.0347656, nnz 17088
- FACTOR 1 (real) SPARSE, size 640x512, density 0.0517578, nnz 16960
- FACTOR 2 (real) SPARSE, size 512x384, density 0.085612, nnz 16832
- FACTOR 3 (real) SPARSE, size 384x256, density 0.169922, nnz 16704
- FACTOR 4 (real) SPARSE, size 256x128, density 0.501953, nnz 16448
- FACTOR 5 (real) SPARSE, size 128x128, density 0.0078125, nnz 128
identity matrix flag
%% Cell type:markdown id:7f0533a8 tags:
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):
%% Cell type:code id:95b00e50 tags:
``` python
F[:d,:].toarray()
```
%% Output
array([[1., 0., 0., ..., 0., 0., 0.],
[0., 1., 0., ..., 0., 0., 0.],
[0., 0., 1., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 1., 0., 0.],
[0., 0., 0., ..., 0., 1., 0.],
[0., 0., 0., ..., 0., 0., 1.]])
%% Cell type:markdown id:dda9bfcc tags:
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 [page](https://en.wikipedia.org/wiki/Chebyshev_polynomials).
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](https://en.wikipedia.org/wiki/Chebyshev_polynomials).
%% Cell type:markdown id:f2c4274c tags:
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 with the generic Faust-vector/matrix multiplication. Indeed, due to the particular structure of the polynomial basis Faust, the multiplication can be optimized.
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).
%% Cell type:code id:efaf2c4f tags:
``` python
G = F.clone()
from numpy.random import rand
G = F.clone()
X = rand(F.shape[1],100)
%timeit F@X
%timeit G@X
```
%% Output
1.96 ms ± 181 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
7.46 ms ± 100 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%% Cell type:markdown id:4a6cbefc tags:
Now let's verify the multiplication result is accurate:
%% Cell type:code id:7f031472 tags:
``` python
from numpy.linalg import norm
print("err=", norm(F@X-G@X)/norm(F@X))
```
%% Output
err= 1.239890671949257e-16
%% Cell type:markdown id:ef58bba3 tags:
As you see it's alright.
%% Cell type:markdown id:bc2683cd tags:
## 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``. So passing the ``FaustPoly`` and the coefficients (one per polynomial) for the linear combination you'll obtain:
The second function of the ``pyfaust.poly`` module is ``poly``. This function purpose is to compute a linear combination of polynomials composing a ``FaustPoly``. So passing the ``FaustPoly`` and the linear combination coefficients (one per polynomial, in ascending degree order) you'll obtain:
%% Cell type:code id:2f50fe04 tags:
``` python
from pyfaust.poly import poly
from numpy import array
coeffs = array([rand()*100 for i in range(K+1)])
lc_F = poly(coeffs, F)
lc_F
```
%% Output
Faust size 128x128, density 5.18359, nnz_sum 84928, 7 factor(s):
- FACTOR 0 (real) SPARSE, size 128x768, density 0.0078125, nnz 768
- FACTOR 1 (real) SPARSE, size 768x640, density 0.0347656, nnz 17088
- FACTOR 2 (real) SPARSE, size 640x512, density 0.0517578, nnz 16960
- FACTOR 3 (real) SPARSE, size 512x384, density 0.085612, nnz 16832
- FACTOR 4 (real) SPARSE, size 384x256, density 0.169922, nnz 16704
- FACTOR 5 (real) SPARSE, size 256x128, density 0.501953, nnz 16448
- FACTOR 6 (real) SPARSE, size 128x128, density 0.0078125, nnz 128
identity matrix flag
%% Cell type:markdown id:e88d2b0b tags:
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).
%% Cell type:code id:fb13f8f8 tags:
``` python
from scipy.sparse import eye
from scipy.sparse import hstack
from pyfaust import Faust
Id = eye(d, format='csr')
#coeffs_factor = hstack(tuple(coeffs[i]*Id for i in range(K+1)), format='csr')
#lc_G = Faust(coeffs_factor)@G
lc_G = coeffs[0]*G[:d,:]
for i in range(1, K+1):
lc_G += coeffs[i]*G[d*i:d*(i+1),:]
print("error lc_G/lc_F:", (lc_F-lc_G).norm()/(lc_G).norm())
```
%% Output
error lc_G/lc_F: 4.280382463461641e-16
%% Cell type:markdown id:2a321a8b tags:
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.]()) that consist to apply a linear combination of F to a vector x as it's shown below.
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.
%% Cell type:code id:5bf4ca1a tags:
``` python
x = rand(F.shape[1])
way1 = lambda: poly(coeffs, F)@x # first way to do as we've done above
way2 = lambda: poly(coeffs, F@x) # second way to do (that is quicker)
way3 = lambda: poly(coeffs, F, X=x) # third way to do (it's even quicker)
```
%% Cell type:code id:1a438a30 tags:
``` python
%timeit way1()
```
%% Output
127 µs ± 6.38 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
%% Cell type:code id:279c52d5 tags:
``` python
%timeit way2()
```
%% Output
98.3 µs ± 2.46 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
%% Cell type:code id:3a0d5e60 tags:
``` python
%timeit way3()
```
%% Output
78.3 µs ± 430 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
%% Cell type:markdown id:d3319d60 tags:
OK! The second way to do is quicker and the third is the quickest but just in case let's verify all ways give the same results.
%% Cell type:code id:9b0c7188 tags:
``` python
print("err way2 =", norm(poly(coeffs, F)@x-poly(coeffs, F@x))/norm(poly(coeffs, F)@x))
```
%% Output
err way2 = 1.2393110086488054e-16
%% Cell type:code id:cc3834b0 tags:
``` python
print("err way3 =", norm(poly(coeffs, F)@x-poly(coeffs, F, X=x))/norm(poly(coeffs, F)@x))
```
%% Output
err way3 = 1.2393110086488054e-16
%% Cell type:markdown id:90c1042a tags:
All sounds good! We shall now introduce two use cases of Chebyshev polynomials 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, fore more details I invite you to consult the [API documentation](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/namespacepyfaust.html).
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, fore more details I invite you to consult the [API documentation](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/namespacepyfaust.html).
%% Cell type:markdown id:e66108e8 tags:
## 3. Computing the action of the inverse / exponential of a matrix on a vector
Computing the action of the inverse of a matrix A on a vector x means to compute $y = A^{-1} x$ without actually computing the inverse of A.
That's exactly what the function ``pyfaust.poly.invm_multiply`` is able to do using ``FaustPoly`` behind the scene. Let's generate a random sparse matrix (because that's the interesting use case) and compute $y$ using ``invm_multiply`` besides to computing the inverse in order to measure the relative performance and accurracy.
%% Cell type:code id:6fd97c47 tags:
``` python
from pyfaust.poly import invm_multiply
from scipy.sparse.linalg import inv
from numpy.linalg import inv as numpy_inv
from scipy.sparse import csc_matrix
import numpy as np
S = random(d, d, .05, format='csr')
x = rand(d, 1)
A = S@S.T # invm_multiply is only able to treat a symmetric positive definite matrix, that's its limitation!
cscA = csc_matrix(A) # scipy prefers CSC to invert a matrix (for efficiency)
max_K = np.inf
rel_err=1e-3
full_A = A.toarray()
%timeit inv(cscA)@x
%timeit invm_multiply(A, x, tradeoff='time', rel_err=rel_err, max_K=max_K)
%timeit numpy_inv(full_A)@x
```
%% Output
32.1 ms ± 563 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
154 ms ± 1.1 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
465 µs ± 50.9 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
%% Cell type:markdown id:101f2a3a tags:
## 3. Computing the action of the exponential of a matrix on a vector
Now that we've seen that ``invm_multiply`` works quicker, as usual we check below that it's accurate:
%% Cell type:code id:e94a7186 tags:
``` python
y1 = inv(cscA)@x
y2 = invm_multiply(A, x, tradeoff='time', rel_err=rel_err, max_K=max_K)
print("err:", norm(y2-y1)/norm(y1))
```
%% Output
err: 0.8689511053955197
%% Cell type:markdown id:be5be836 tags:
There are function arguments to govern the precision and computation time tradeoff (the most important argument being 'tradeoff' itself, which takes the 'memory' or 'time' depending on your primary criterion). You'll find them in the [API documentation](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/namespacepyfaust.html).
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``](https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.expm_multiply.html).
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).
%% Cell type:code id:82c7c82d tags:
``` python
from scipy.sparse.linalg import expm_multiply as scipy_expm_multiply
from pyfaust.poly import expm_multiply
from numpy import exp, linspace
from numpy.linalg import norm
from numpy.random import rand
from scipy.sparse import random
S = random(1024, 1024, .002, format='csr')
A = S@S.T
X = rand(A.shape[1], 1)
pts_args = {'start':-.5, 'stop':-.2, 'num':1000}
pts = linspace(**pts_args)
y1 = expm_multiply(A, X, pts)
y2 = scipy_expm_multiply(A, X, **pts_args)
print("pyfaust error:", norm(y1-y2)/norm(y2))
%timeit expm_multiply(A, X, pts)
%timeit scipy_expm_multiply(A, X, **pts_args)
```
%% Output
%% Cell type:markdown id:c4312be4 tags:
pyfaust error: 7.65649062956808e-11
130 ms ± 1.74 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
1.2 s ± 28.6 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
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](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/namespacepyfaust.html).
%% Cell type:markdown id:c4312be4 tags:
%% Cell type:markdown id:dcf619c8 tags:
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](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/namespacepyfaust.html).
Thanks for reading this notebook! Many other are available at [faust.inria.fr](https://faust.inria.fr).
**Note**: this notebook was executed using the following pyfaust version:
%% Cell type:code id:325a2040 tags:
``` python
import pyfaust
pyfaust.version()
```
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment