"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. \n",
"\n",
"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.\n",
"\n",
"**NOTE**: all the functions introduced in this notebook are available \n",
"\n",
"## 1. The basis function\n",
"\n",
"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:\n",
"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. \n",
"\n",
"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:\n"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "d077b727",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Faust size 768x128, density 0.85612, nnz_sum 84160, 6 factor(s): \n",
"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). \n",
"\n",
"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",
"execution_count": 2,
"id": "95b00e50",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[1., 0., 0., ..., 0., 0., 0.],\n",
" [0., 1., 0., ..., 0., 0., 0.],\n",
" [0., 0., 1., ..., 0., 0., 0.],\n",
" ...,\n",
" [0., 0., 0., ..., 1., 0., 0.],\n",
" [0., 0., 0., ..., 0., 1., 0.],\n",
" [0., 0., 0., ..., 0., 0., 1.]])"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"F[:d,:].toarray()"
]
},
{
"cell_type": "markdown",
"id": "dda9bfcc",
"metadata": {},
"source": [
"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...\n",
"\n",
"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)."
]
},
{
"cell_type": "markdown",
"id": "f2c4274c",
"metadata": {},
"source": [
"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.\n",
"\n",
"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",
"execution_count": 3,
"id": "efaf2c4f",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"1.96 ms ± 181 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n",
"7.46 ms ± 100 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
]
}
],
"source": [
"G = F.clone()\n",
"from numpy.random import rand\n",
"X = rand(F.shape[1],100)\n",
"%timeit F@X\n",
"%timeit G@X"
]
},
{
"cell_type": "markdown",
"id": "4a6cbefc",
"metadata": {},
"source": [
"Now let's verify the multiplication result is accurate:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "7f031472",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"err= 1.239890671949257e-16\n"
]
}
],
"source": [
"from numpy.linalg import norm\n",
"print(\"err=\", norm(F@X-G@X)/norm(F@X))"
]
},
{
"cell_type": "markdown",
"id": "ef58bba3",
"metadata": {},
"source": [
"As you see it's alright."
]
},
{
"cell_type": "markdown",
"id": "bc2683cd",
"metadata": {},
"source": [
"## 2. The poly function\n",
"\n",
"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:"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "2f50fe04",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Faust size 128x128, density 5.18359, nnz_sum 84928, 7 factor(s): \n",
"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."
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "5bf4ca1a",
"metadata": {},
"outputs": [],
"source": [
"x = rand(F.shape[1])\n",
"way1 = lambda: poly(coeffs, F)@x # first way to do as we've done above\n",
"way2 = lambda: poly(coeffs, F@x) # second way to do (that is quicker)\n",
"way3 = lambda: poly(coeffs, F, X=x) # third way to do (it's even quicker)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "1a438a30",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"127 µs ± 6.38 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n"
]
}
],
"source": [
"%timeit way1()"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "279c52d5",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"98.3 µs ± 2.46 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n"
]
}
],
"source": [
"%timeit way2()"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "3a0d5e60",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"78.3 µs ± 430 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n"
]
}
],
"source": [
"%timeit way3()"
]
},
{
"cell_type": "markdown",
"id": "d3319d60",
"metadata": {},
"source": [
"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."
"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)."
]
},
{
"cell_type": "markdown",
"id": "e66108e8",
"metadata": {},
"source": [
"## 3. Computing the action of the inverse / exponential of a matrix on a vector\n",
"\n",
"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.\n",
"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",
"execution_count": 54,
"id": "6fd97c47",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"32.1 ms ± 563 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
"154 ms ± 1.1 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
"465 µs ± 50.9 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n"
]
}
],
"source": [
"from pyfaust.poly import invm_multiply\n",
"from scipy.sparse.linalg import inv\n",
"from numpy.linalg import inv as numpy_inv\n",
"from scipy.sparse import csc_matrix\n",
"import numpy as np\n",
"S = random(d, d, .05, format='csr')\n",
"x = rand(d, 1)\n",
"A = S@S.T # invm_multiply is only able to treat a symmetric positive definite matrix, that's its limitation!\n",
"cscA = csc_matrix(A) # scipy prefers CSC to invert a matrix (for efficiency)\n",
"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).\n",
"\n",
"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).\n",
"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. \n",
"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",
"execution_count": 15,
"id": "82c7c82d",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"pyfaust error: 7.65649062956808e-11\n",
"130 ms ± 1.74 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
"1.2 s ± 28.6 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
]
}
],
"source": [
"from scipy.sparse.linalg import expm_multiply as scipy_expm_multiply\n",
"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)."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.10"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
%% 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.
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.
**NOTE**: all the functions introduced in this notebook are available
## 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:
%% 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):
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).
%% 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.
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
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:
%% 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):
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.
%% 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.
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).
%% 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)
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
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)
%% 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).
