Mentions légales du service

Skip to content
Snippets Groups Projects
Commit 512d5774 authored by hhakim's avatar hhakim
Browse files

Add matfaust.poly livescript and update pyfaust.poly jupyter notebook.

parent c5de45d5
Branches
Tags
No related merge requests found
Pipeline #834004 skipped
Expand all Hide whitespace changes
Inline Side-by-side
gen_doc/CMakeLists.txt
+ 2
− 0
View file @ 512d5774
...@@ -166,6 +166,8 @@ if(BUILD_DOCUMENTATION) ...@@ -166,6 +166,8 @@ if(BUILD_DOCUMENTATION)
configure_file(${FAUST_DOC_SRC_DIR}/faust_projectors.mlx ${PROJECT_BINARY_DIR}/doc/html/faust_projectors.mlx COPYONLY) configure_file(${FAUST_DOC_SRC_DIR}/faust_projectors.mlx ${PROJECT_BINARY_DIR}/doc/html/faust_projectors.mlx COPYONLY)
configure_file(${FAUST_DOC_SRC_DIR}/faust_projectors.mlx.html ${PROJECT_BINARY_DIR}/doc/html/faust_projectors.mlx.html COPYONLY) configure_file(${FAUST_DOC_SRC_DIR}/faust_projectors.mlx.html ${PROJECT_BINARY_DIR}/doc/html/faust_projectors.mlx.html COPYONLY)
configure_file(${FAUST_DOC_SRC_DIR}/matfaust_poly.mlx ${PROJECT_BINARY_DIR}/doc/html/matfaust_poly.mlx COPYONLY)
configure_file(${FAUST_SRC_TEST_SRC_PYTHON_DIR}/test_svd_rc_vs_err.py ${PROJECT_BINARY_DIR}/doc/html/test_svd_rc_vs_err.py COPYONLY) configure_file(${FAUST_SRC_TEST_SRC_PYTHON_DIR}/test_svd_rc_vs_err.py ${PROJECT_BINARY_DIR}/doc/html/test_svd_rc_vs_err.py COPYONLY)
archive_tutos(pyfaust notebook ipynb) archive_tutos(pyfaust notebook ipynb)
archive_tutos(matfaust livescript mlx) archive_tutos(matfaust livescript mlx)
......
gen_doc/matfaust_poly.mlx 0 → 100644
+ 0
− 0
View file @ 512d5774
File added
gen_doc/pyfaust_poly.html
+ 37
− 37
View file @ 512d5774
This diff is collapsed.
gen_doc/pyfaust_poly.ipynb
+ 14
− 13
View file @ 512d5774
%% Cell type:markdown id:6f8cc8a4 tags: %% Cell type:markdown id:6f8cc8a4 tags:
# Using the poly module # Using the poly module
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. 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 one precise use case with the action of 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 on GPU, using the ``dev='gpu'`` argument. **NOTE**: all the functions introduced in this notebook are available on GPU, using the ``dev='gpu'`` argument.
## 1. The basis function ## 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: 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)``` ```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. 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 (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: 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: %% Cell type:code id:d077b727 tags:
``` python ``` python
from pyfaust.poly import basis from pyfaust.poly import basis
from scipy.sparse import random from scipy.sparse import random
d = 128 d = 128
L = random(d, d, .2, format='csr') L = random(d, d, .2, format='csr')
L = L@L.T L = L@L.T
K=5 K = 5
F = basis(L, K=K, basis_name='chebyshev') F = basis(L, K=K, basis_name='chebyshev')
F F
``` ```
%% Cell type:markdown id:7f0533a8 tags: %% 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). 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): 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: %% Cell type:code id:95b00e50 tags:
``` python ``` python
F[:d,:].toarray() F[:d,:].toarray()
``` ```
%% Cell type:markdown id:dda9bfcc tags: %% 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... 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 [wikipedia article](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: %% 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 compared to 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). 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). Note that ``Faust.clone`` function is not used because in this case it would return a polynomial basis Faust (after all that's the role of ``clone`` to preserve the Faust properties!). To get a classic Faust the only way is to copy the Faust factor by factor.
%% Cell type:code id:efaf2c4f tags: %% Cell type:code id:efaf2c4f tags:
``` python ``` python
from numpy.random import rand from numpy.random import rand
G = F.clone() from pyfaust import Faust
factors = [F.factors(i) for i in range(F.numfactors())]
G = Faust(factors)
X = rand(F.shape[1],100) X = rand(F.shape[1],100)
%timeit F@X %timeit F@X
%timeit G@X %timeit G@X
``` ```
%% Cell type:markdown id:4a6cbefc tags: %% Cell type:markdown id:4a6cbefc tags:
Now let's verify the multiplication result is accurate: Now let's verify the multiplication result is accurate:
%% Cell type:code id:7f031472 tags: %% Cell type:code id:7f031472 tags:
``` python ``` python
from numpy.linalg import norm from numpy.linalg import norm
print("err=", norm(F@X-G@X)/norm(F@X)) print("err=", norm(F@X-G@X)/norm(F@X))
``` ```
%% Cell type:markdown id:ef58bba3 tags: %% Cell type:markdown id:ef58bba3 tags:
As you see it's alright. As you see it's alright.
%% Cell type:markdown id:bc2683cd tags: %% Cell type:markdown id:bc2683cd tags:
## 2. The poly function ## 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 linear combination coefficients (one per polynomial, in ascending degree order) 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`` (it can also be viewed as way to compute a polynomial). 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: %% Cell type:code id:2f50fe04 tags:
``` python ``` python
from pyfaust.poly import poly from pyfaust.poly import poly
from numpy import array from numpy import array
coeffs = array([rand()*100 for i in range(K+1)]) coeffs = rand(K+1)*100
lc_F = poly(coeffs, F) lc_F = poly(coeffs, F)
lc_F lc_F
``` ```
%% Cell type:markdown id:e88d2b0b tags: %% 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). 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: %% Cell type:code id:fb13f8f8 tags:
``` python ``` python
from scipy.sparse import eye from scipy.sparse import eye
from scipy.sparse import hstack from scipy.sparse import hstack
from pyfaust import Faust 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,:] lc_G = coeffs[0]*G[:d,:]
for i in range(1, K+1): for i in range(1, K+1):
lc_G += coeffs[i]*G[d*i:d*(i+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()) print("error lc_G/lc_F:", (lc_F-lc_G).norm()/(lc_G).norm())
``` ```
%% Cell type:markdown id:2a321a8b tags: %% 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.]()) 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: %% Cell type:code id:5bf4ca1a tags:
``` python ``` python
x = rand(F.shape[1]) x = rand(F.shape[1])
way1 = lambda: poly(coeffs, F)@x # first way to do as we've done above 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) 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) way3 = lambda: poly(coeffs, F, X=x) # third way to do (it's even quicker)
``` ```
%% Cell type:code id:1a438a30 tags: %% Cell type:code id:1a438a30 tags:
``` python ``` python
%timeit way1() %timeit way1()
``` ```
%% Cell type:code id:279c52d5 tags: %% Cell type:code id:279c52d5 tags:
``` python ``` python
%timeit way2() %timeit way2()
``` ```
%% Cell type:code id:3a0d5e60 tags: %% Cell type:code id:3a0d5e60 tags:
``` python ``` python
%timeit way3() %timeit way3()
``` ```
%% Cell type:markdown id:d3319d60 tags: %% 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. Depending on the situation the ``way2`` or ``way3`` might be the quickest but they are always both quicker than ``way1``.
Just in case let's verify all ways give the same results.
%% Cell type:code id:9b0c7188 tags: %% Cell type:code id:9b0c7188 tags:
``` python ``` python
print("err way2 =", norm(poly(coeffs, F)@x-poly(coeffs, F@x))/norm(poly(coeffs, F)@x)) print("err way2 =", norm(poly(coeffs, F)@x-poly(coeffs, F@x))/norm(poly(coeffs, F)@x))
``` ```
%% Cell type:code id:cc3834b0 tags: %% Cell type:code id:cc3834b0 tags:
``` python ``` python
print("err way3 =", norm(poly(coeffs, F)@x-poly(coeffs, F, X=x))/norm(poly(coeffs, F)@x)) print("err way3 =", norm(poly(coeffs, F)@x-poly(coeffs, F, X=x))/norm(poly(coeffs, F)@x))
``` ```
%% Cell type:markdown id:90c1042a tags: %% Cell type:markdown id:90c1042a tags:
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). 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 like, for more details I invite you to consult the [API documentation](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/namespacepyfaust_1_1poly.html#afc6e14fb360a1650cddf8a7646b9209c).
%% Cell type:markdown id:e66108e8 tags: %% Cell type:markdown id:e66108e8 tags:
## 3. Computing the action of the exponential of a matrix on a vector ## 3. Computing the action of the exponential of a matrix on a vector
%% Cell type:markdown id:be5be836 tags: %% Cell type:markdown id:be5be836 tags:
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). 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 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). 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: %% Cell type:code id:82c7c82d tags:
``` python ``` python
from scipy.sparse.linalg import expm_multiply as scipy_expm_multiply from scipy.sparse.linalg import expm_multiply as scipy_expm_multiply
from pyfaust.poly import expm_multiply from pyfaust.poly import expm_multiply
from numpy import exp, linspace from numpy import exp, linspace
from numpy.linalg import norm from numpy.linalg import norm
from numpy.random import rand from numpy.random import rand
from scipy.sparse import random from scipy.sparse import random
S = random(1024, 1024, .002, format='csr') S = random(1024, 1024, .002, format='csr')
A = S@S.T A = S@S.T
X = rand(A.shape[1], 1) X = rand(A.shape[1], 1)
pts_args = {'start':-.5, 'stop':-.2, 'num':1000} pts_args = {'start':-.5, 'stop':-.2, 'num':1000}
pts = linspace(**pts_args) pts = linspace(**pts_args)
y1 = expm_multiply(A, X, pts) y1 = expm_multiply(A, X, pts)
y2 = scipy_expm_multiply(A, X, **pts_args) y2 = scipy_expm_multiply(A, X, **pts_args)
print("pyfaust error:", norm(y1-y2)/norm(y2)) print("pyfaust error:", norm(y1-y2)/norm(y2))
%timeit expm_multiply(A, X, pts) %timeit expm_multiply(A, X, pts)
%timeit scipy_expm_multiply(A, X, **pts_args) %timeit scipy_expm_multiply(A, X, **pts_args)
``` ```
%% Cell type:markdown id:c4312be4 tags: %% Cell type:markdown id:c4312be4 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). 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_1_1poly.html#aebbce7977632e3c85260738604f01104).
%% Cell type:markdown id:dcf619c8 tags: %% Cell type:markdown id:dcf619c8 tags:
Thanks for reading this notebook! Many other are available at [faust.inria.fr](https://faust.inria.fr). 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: **Note**: this notebook was executed using the following pyfaust version:
%% Cell type:code id:325a2040 tags: %% Cell type:code id:325a2040 tags:
``` python ``` python
import pyfaust import pyfaust
pyfaust.version() pyfaust.version()
``` ```
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment