Mentions légales du service

Skip to content
Snippets Groups Projects
Commit 21aa7692 authored by hhakim's avatar hhakim
Browse files

Review of poly module tutorials with minor changes.

parent 05723db2
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
gen_doc/matfaust_poly.html
+ 6
− 6
View file @ 21aa7692
...@@ -55,7 +55,7 @@ ...@@ -55,7 +55,7 @@
.matrixElement .verticalEllipsis,.textElement .verticalEllipsis,.rtcDataTipElement .matrixElement .verticalEllipsis,.rtcDataTipElement .textElement .verticalEllipsis {margin-left: 35px; width: 12px; height: 30px; background-repeat: no-repeat; background-image: url("data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAAoAAAAZCAYAAAAIcL+IAAAALklEQVR42mP4//8/AzGYgWyFMECMwv8QddRS+P//KyimlmcGUOFoOI6GI/UVAgDnd8Dd4+NCwgAAAABJRU5ErkJggg==");} .matrixElement .verticalEllipsis,.textElement .verticalEllipsis,.rtcDataTipElement .matrixElement .verticalEllipsis,.rtcDataTipElement .textElement .verticalEllipsis {margin-left: 35px; width: 12px; height: 30px; background-repeat: no-repeat; background-image: url("data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAAoAAAAZCAYAAAAIcL+IAAAALklEQVR42mP4//8/AzGYgWyFMECMwv8QddRS+P//KyimlmcGUOFoOI6GI/UVAgDnd8Dd4+NCwgAAAABJRU5ErkJggg==");}
.S9 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px; padding: 6px 45px 4px 13px; line-height: 17.2339992523193px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, 'Courier New', monospace; font-size: 14px; } .S9 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px; padding: 6px 45px 4px 13px; line-height: 17.2339992523193px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, 'Courier New', monospace; font-size: 14px; }
.S10 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 0px none rgb(0, 0, 0); border-radius: 0px; padding: 6px 45px 0px 13px; line-height: 17.2339992523193px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, 'Courier New', monospace; font-size: 14px; } .S10 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 0px none rgb(0, 0, 0); border-radius: 0px; padding: 6px 45px 0px 13px; line-height: 17.2339992523193px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, 'Courier New', monospace; font-size: 14px; }
.S11 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px 0px 4px 4px; padding: 6px 45px 4px 13px; line-height: 17.2339992523193px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, 'Courier New', monospace; font-size: 14px; }</style></head><body><div class = rtcContent><h1 class = 'S0'><span>Using the poly module</span></h1><div class = 'S1'><span>A new module has been added to matfaust version 3.1.x. Its name is </span><span style=' font-family: monospace;'>poly</span><span> and as expected it is dedicated to a kind of Fausts that are defined according to series of polynomials.</span></div><div class = 'S1'><span>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.</span></div><div class = 'S1'><span style=' font-weight: bold;'>NOTE</span><span>: all the functions introduced in this notebook are available on GPU, using the '</span><span style=' font-family: monospace;'>dev', 'gpu'</span><span> arguments.</span></div><h2 class = 'S2'><span>1. The basis function</span></h2><div class = 'S1'><span></span></div><div class = 'S1'><span>Firstly, the </span><span style=' font-family: monospace;'>poly</span><span> module allows to define a polynomial basis (aka </span><span style=' font-family: monospace;'>FaustPoly</span><span>) using the function </span><span style=' font-family: monospace;'>matfaust.poly.basis</span><span>. Currently, only Chebyshev polynomials are supported but other are yet to come. Below is the prototype of the function:</span></div><div class = 'S1'><span style=' font-family: monospace;'>basis(L, K, basis_name, varargin)</span></div><div class = 'S1'><span style=' font-weight: bold; font-family: monospace;'>NOTE:</span><span style=' font-family: monospace;'> it is noteworthy that varargin contains the optional T0 argument, we'll come back to this point later.</span></div><div class = 'S1'><span>In the next, we shall see a simple example but I let you consult the documentation by typing </span><span style=' font-family: monospace;'>help matfaust.poly.basis</span><span> to get more details (alternatively, this is the </span><a href = "https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/namespacematfaust_1_1poly.html#a556d6a4eedc41529ddea692d227caaf4"><span>online doc</span></a><span>).</span></div><div class = 'S1'><span>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 sparse matrix at least square (and most likely symmetric positive definite in much use cases), you'll make this call to the function:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre;"><span>import </span><span style="color: rgb(160, 32, 240);">matfaust.poly.basis</span></span></div></div><div class="inlineWrapper"><div class = 'S4'><span style="white-space: pre;"><span>d = 128;</span></span></div></div><div class="inlineWrapper"><div class = 'S4'><span style="white-space: pre;"><span>L = sprand(d, d, .2);</span></span></div></div><div class="inlineWrapper"><div class = 'S4'><span style="white-space: pre;"><span>L = L*L';</span></span></div></div><div class="inlineWrapper"><div class = 'S4'><span style="white-space: pre;"><span>K = 5;</span></span></div></div><div class="inlineWrapper outputs"><div class = 'S5'><span style="white-space: pre;"><span>F = basis(L, K, </span><span style="color: rgb(160, 32, 240);">'chebyshev'</span><span>)</span></span></div><div class = 'S6'><div class="inlineElement eoOutputWrapper embeddedOutputsVariableStringElement" uid="F7B6A34C" data-testid="output_0" data-width="495" data-height="132" data-hashorizontaloverflow="false" style="width: 525px; max-height: 261px;"><div class="textElement"><div><span class="variableNameElement">F = </span></div><div>Faust size 768x128, density 0.847371, nnz_sum 83300, 6 factor(s): .S11 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px 0px 4px 4px; padding: 6px 45px 4px 13px; line-height: 17.2339992523193px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, 'Courier New', monospace; font-size: 14px; }</style></head><body><div class = rtcContent><h1 class = 'S0'><span>Using the poly module</span></h1><div class = 'S1'><span>A new module has been added to matfaust version 3.1.x. Its name is </span><span style=' font-family: monospace;'>poly</span><span> and as expected it is dedicated to a kind of Fausts that are defined according to series of polynomials.</span></div><div class = 'S1'><span>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.</span></div><div class = 'S1'><span style=' font-weight: bold;'>NOTE</span><span>: all the functions introduced in this notebook are available on GPU, using the '</span><span style=' font-family: monospace;'>dev', 'gpu'</span><span> arguments.</span></div><h2 class = 'S2'><span>1. The basis function</span></h2><div class = 'S1'><span></span></div><div class = 'S1'><span>Firstly, the </span><span style=' font-family: monospace;'>poly</span><span> module allows to define a polynomial basis (aka </span><span style=' font-family: monospace;'>FaustPoly</span><span>) using the function </span><span style=' font-family: monospace;'>matfaust.poly.basis</span><span>. Currently, only Chebyshev polynomials are supported but others are yet to come. Below is the prototype of the function:</span></div><div class = 'S1'><span style=' font-family: monospace;'>basis(L, K, basis_name, varargin)</span></div><div class = 'S1'><span style=' font-weight: bold; font-family: monospace;'>NOTE:</span><span style=' font-family: monospace;'> it is noteworthy that varargin contains the optional T0 argument, we'll come back to this point later.</span></div><div class = 'S1'><span>In the next, we shall see a simple example but I let you consult the documentation by typing </span><span style=' font-family: monospace;'>help matfaust.poly.basis</span><span> to get more details (alternatively, this is the </span><a href = "https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/namespacematfaust_1_1poly.html#a556d6a4eedc41529ddea692d227caaf4"><span>online doc</span></a><span>).</span></div><div class = 'S1'><span>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 sparse matrix at least square (and most likely symmetric positive definite in much use cases), you'll make this call to the function:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre;"><span>import </span><span style="color: rgb(160, 32, 240);">matfaust.poly.basis</span></span></div></div><div class="inlineWrapper"><div class = 'S4'><span style="white-space: pre;"><span>d = 128;</span></span></div></div><div class="inlineWrapper"><div class = 'S4'><span style="white-space: pre;"><span>L = sprand(d, d, .2);</span></span></div></div><div class="inlineWrapper"><div class = 'S4'><span style="white-space: pre;"><span>L = L*L';</span></span></div></div><div class="inlineWrapper"><div class = 'S4'><span style="white-space: pre;"><span>K = 5;</span></span></div></div><div class="inlineWrapper outputs"><div class = 'S5'><span style="white-space: pre;"><span>F = basis(L, K, </span><span style="color: rgb(160, 32, 240);">'chebyshev'</span><span>)</span></span></div><div class = 'S6'><div class="inlineElement eoOutputWrapper embeddedOutputsVariableStringElement" uid="F7B6A34C" data-testid="output_0" data-width="495" data-height="132" data-hashorizontaloverflow="false" style="width: 525px; max-height: 261px;"><div class="textElement"><div><span class="variableNameElement">F = </span></div><div>Faust size 768x128, density 0.847371, nnz_sum 83300, 6 factor(s):
- FACTOR 0 (real) SPARSE, size 768x640, density 0.0344157, nnz 16916 - FACTOR 0 (real) SPARSE, size 768x640, density 0.0344157, nnz 16916
- FACTOR 1 (real) SPARSE, size 640x512, density 0.0512329, nnz 16788 - FACTOR 1 (real) SPARSE, size 640x512, density 0.0512329, nnz 16788
- FACTOR 2 (real) SPARSE, size 512x384, density 0.0847371, nnz 16660 - FACTOR 2 (real) SPARSE, size 512x384, density 0.0847371, nnz 16660
...@@ -72,7 +72,7 @@ ...@@ -72,7 +72,7 @@
0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
</div><div class="horizontalEllipsis" style="white-space: nowrap; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"></div><div class="verticalEllipsis" style="white-space: nowrap; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"></div></div></div></div></div></div></div><div class = 'S7'><span>This first 0-degree polynomial is followed by the next degree polynomials: hence </span><span style=' font-family: monospace;'>F(d+1:d*2, :)</span><span> is the 1-degree polynomial, </span><span style=' font-family: monospace;'>F(d*2+1:d*3, :)</span><span> is the 2-degree polynomial and so on...</span></div><div class = 'S1'><span>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 </span><a href = "https://en.wikipedia.org/wiki/Chebyshev_polynomials"><span style=' text-decoration: underline;'>wikipedia article</span></a><span>.</span></div><div class = 'S1'><span>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.</span></div><div class = 'S1'><span>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 </span><span style=' font-family: monospace;'>Faust.clone</span><span> function is not used because in this case it would return a polynomial basis Faust (after all that's the role of </span><span style=' font-family: monospace;'>clone </span><span>to preserve the Faust properties!). To get a classic Faust the only way is to copy the Faust factor by factor.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre;"><span>facs = {};</span></span></div></div><div class="inlineWrapper"><div class = 'S4'><span style="white-space: pre;"><span style="color: rgb(0, 0, 255);">for </span><span>i=1:numfactors(F)</span></span></div></div><div class="inlineWrapper"><div class = 'S4'><span style="white-space: pre;"><span> </span><span class="warning_squiggle_rte">facs</span><span> = [facs {factors(F,i)}]; </span></span></div></div><div class="inlineWrapper"><div class = 'S4'><span style="white-space: pre;"><span style="color: rgb(0, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div class = 'S4'><span style="white-space: pre;"><span>G = matfaust.Faust(facs);</span></span></div></div><div class="inlineWrapper"><div class = 'S4'><span style="white-space: pre;"><span>X = rand(size(F, 2), 100);</span></span></div></div><div class="inlineWrapper outputs"><div class = 'S5'><span style="white-space: pre;"><span>timeit(@() F*X)</span></span></div><div class = 'S6'><div class='variableElement' style='font-family: Menlo, Monaco, Consolas, &quot;Courier New&quot;, monospace; font-size: 12px; '>ans = 0.0132</div></div></div><div class="inlineWrapper outputs"><div class = 'S9'><span style="white-space: pre;"><span>timeit(@() G*X)</span></span></div><div class = 'S6'><div class='variableElement' style='font-family: Menlo, Monaco, Consolas, &quot;Courier New&quot;, monospace; font-size: 12px; '>ans = 0.0221</div></div></div></div><div class = 'S7'><span>Now let's verify the multiplication result is accurate:</span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div class = 'S8'><span style="white-space: pre;"><span>err = norm(F*X-G*X)/norm(F*X)</span></span></div><div class = 'S6'><div class='variableElement' style='font-family: Menlo, Monaco, Consolas, &quot;Courier New&quot;, monospace; font-size: 12px; '>err = 2.7722e-17</div></div></div></div><div class = 'S7'><span>As you see it's alright.</span></div><h2 class = 'S2'><span>2. The poly function</span></h2><div class = 'S1'><span></span></div><div class = 'S1'><span>The second function of the </span><span style=' font-family: monospace;'>matfaust.poly</span><span> module is </span><span style=' font-family: monospace;'>poly</span><span>. This function purpose is to compute a linear combination of polynomials composing a polynomial basis / </span><span style=' font-family: monospace;'>FaustPoly</span><span> (it can also be viewed as a way to compute a polynomial). So passing the </span><span style=' font-family: monospace;'>FaustPoly</span><span> </span><span style=' font-family: monospace;'>F</span><span> and the linear combination coefficients (one per polynomial, in ascending degree order) you'll obtain:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre;"><span>import </span><span style="color: rgb(160, 32, 240);">matfaust.poly.poly</span></span></div></div><div class="inlineWrapper"><div class = 'S4'><span style="white-space: pre;"><span>coeffs = rand(K+1, 1)*100;</span></span></div></div><div class="inlineWrapper outputs"><div class = 'S5'><span style="white-space: pre;"><span>lc_F = poly(coeffs, F)</span></span></div><div class = 'S6'><div class="inlineElement eoOutputWrapper embeddedOutputsVariableStringElement" uid="21DA6029" data-testid="output_5" data-width="495" data-height="146" data-hashorizontaloverflow="false" style="width: 525px; max-height: 261px;"><div class="textElement"><div><span class="variableNameElement">lc_F = </span></div><div>Faust size 128x128, density 5.1311, nnz_sum 84068, 7 factor(s): </div><div class="horizontalEllipsis" style="white-space: nowrap; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"></div><div class="verticalEllipsis" style="white-space: nowrap; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"></div></div></div></div></div></div></div><div class = 'S7'><span>This first 0-degree polynomial is followed by the next degree polynomials: hence </span><span style=' font-family: monospace;'>F(d+1:d*2, :)</span><span> is the 1-degree polynomial, </span><span style=' font-family: monospace;'>F(d*2+1:d*3, :)</span><span> is the 2-degree polynomial and so on...</span></div><div class = 'S1'><span>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 </span><a href = "https://en.wikipedia.org/wiki/Chebyshev_polynomials"><span style=' text-decoration: underline;'>wikipedia article</span></a><span>.</span></div><div class = 'S1'><span>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.</span></div><div class = 'S1'><span>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 both cases (with F, the basis Faust and G its classic Faust copy). Note that </span><span style=' font-family: monospace;'>Faust.clone</span><span> function is not used because in this case it would return a polynomial basis Faust (after all that's the role of </span><span style=' font-family: monospace;'>clone </span><span>to preserve the Faust properties!). To get a classic Faust the only way is to copy the Faust factor by factor.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre;"><span>facs = {};</span></span></div></div><div class="inlineWrapper"><div class = 'S4'><span style="white-space: pre;"><span style="color: rgb(0, 0, 255);">for </span><span>i=1:numfactors(F)</span></span></div></div><div class="inlineWrapper"><div class = 'S4'><span style="white-space: pre;"><span> </span><span class="warning_squiggle_rte">facs</span><span> = [facs {factors(F,i)}]; </span></span></div></div><div class="inlineWrapper"><div class = 'S4'><span style="white-space: pre;"><span style="color: rgb(0, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div class = 'S4'><span style="white-space: pre;"><span>G = matfaust.Faust(facs);</span></span></div></div><div class="inlineWrapper"><div class = 'S4'><span style="white-space: pre;"><span>X = rand(size(F, 2), 100);</span></span></div></div><div class="inlineWrapper outputs"><div class = 'S5'><span style="white-space: pre;"><span>timeit(@() F*X)</span></span></div><div class = 'S6'><div class='variableElement' style='font-family: Menlo, Monaco, Consolas, &quot;Courier New&quot;, monospace; font-size: 12px; '>ans = 0.0132</div></div></div><div class="inlineWrapper outputs"><div class = 'S9'><span style="white-space: pre;"><span>timeit(@() G*X)</span></span></div><div class = 'S6'><div class='variableElement' style='font-family: Menlo, Monaco, Consolas, &quot;Courier New&quot;, monospace; font-size: 12px; '>ans = 0.0221</div></div></div></div><div class = 'S7'><span>Now let's verify the multiplication result is accurate:</span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div class = 'S8'><span style="white-space: pre;"><span>err = norm(F*X-G*X)/norm(F*X)</span></span></div><div class = 'S6'><div class='variableElement' style='font-family: Menlo, Monaco, Consolas, &quot;Courier New&quot;, monospace; font-size: 12px; '>err = 2.7722e-17</div></div></div></div><div class = 'S7'><span>As you see it's alright.</span></div><h2 class = 'S2'><span>2. The poly function</span></h2><div class = 'S1'><span></span></div><div class = 'S1'><span>The second function of the </span><span style=' font-family: monospace;'>matfaust.poly</span><span> module is </span><span style=' font-family: monospace;'>poly</span><span>. This function purpose is to compute a linear combination of polynomials composing a polynomial basis / </span><span style=' font-family: monospace;'>FaustPoly</span><span> (it can also be viewed as a way to compute a polynomial). So passing the </span><span style=' font-family: monospace;'>FaustPoly</span><span> </span><span style=' font-family: monospace;'>F</span><span> and the linear combination coefficients (one per polynomial, in ascending degree order) you'll obtain:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre;"><span>import </span><span style="color: rgb(160, 32, 240);">matfaust.poly.poly</span></span></div></div><div class="inlineWrapper"><div class = 'S4'><span style="white-space: pre;"><span>coeffs = rand(K+1, 1)*100;</span></span></div></div><div class="inlineWrapper outputs"><div class = 'S5'><span style="white-space: pre;"><span>lc_F = poly(coeffs, F)</span></span></div><div class = 'S6'><div class="inlineElement eoOutputWrapper embeddedOutputsVariableStringElement" uid="21DA6029" data-testid="output_5" data-width="495" data-height="146" data-hashorizontaloverflow="false" style="width: 525px; max-height: 261px;"><div class="textElement"><div><span class="variableNameElement">lc_F = </span></div><div>Faust size 128x128, density 5.1311, nnz_sum 84068, 7 factor(s):
- FACTOR 0 (real) SPARSE, size 128x768, density 0.0078125, nnz 768 - FACTOR 0 (real) SPARSE, size 128x768, density 0.0078125, nnz 768
- FACTOR 1 (real) SPARSE, size 768x640, density 0.0344157, nnz 16916 - FACTOR 1 (real) SPARSE, size 768x640, density 0.0344157, nnz 16916
- FACTOR 2 (real) SPARSE, size 640x512, density 0.0512329, nnz 16788 - FACTOR 2 (real) SPARSE, size 640x512, density 0.0512329, nnz 16788
...@@ -96,7 +96,7 @@ ...@@ -96,7 +96,7 @@
0.4435 -0.1650 -0.0126 -0.0780 0.9582 0.4745 0.3213 -0.1900 0.1681 0.8485 0.3578 0.3991 0.0178 0.3890 0.2835 -0.0715 -0.1600 0.0966 0.1111 0.1194 0.1658 0.1124 0.1484 0.7558 0.3178 0.2321 0.0195 0.0959 0.6406 -0.1509 0.0169 0.3552 0.1125 0.2538 0.2038 0.5370 0.0098 -0.2254 0.0679 -0.2711 0.4737 0.6133 0.3881 0.0179 0.7434 0.2237 -0.0931 0.0388 -0.2206 0.5995 0.4435 -0.1650 -0.0126 -0.0780 0.9582 0.4745 0.3213 -0.1900 0.1681 0.8485 0.3578 0.3991 0.0178 0.3890 0.2835 -0.0715 -0.1600 0.0966 0.1111 0.1194 0.1658 0.1124 0.1484 0.7558 0.3178 0.2321 0.0195 0.0959 0.6406 -0.1509 0.0169 0.3552 0.1125 0.2538 0.2038 0.5370 0.0098 -0.2254 0.0679 -0.2711 0.4737 0.6133 0.3881 0.0179 0.7434 0.2237 -0.0931 0.0388 -0.2206 0.5995
0.4436 -0.1650 -0.0126 -0.0779 0.9582 0.4747 0.3215 -0.1900 0.1682 0.8485 0.3580 0.3992 0.0178 0.3892 0.2835 -0.0715 -0.1599 0.0966 0.1112 0.1194 0.1659 0.1124 0.1484 0.7558 0.3178 0.2322 0.0196 0.0960 0.6407 -0.1508 0.0169 0.3553 0.1125 0.2540 0.2039 0.5371 0.0099 -0.2253 0.0679 -0.2710 0.4737 0.6134 0.3882 0.0179 0.7434 0.2239 -0.0931 0.0388 -0.2206 0.5996 0.4436 -0.1650 -0.0126 -0.0779 0.9582 0.4747 0.3215 -0.1900 0.1682 0.8485 0.3580 0.3992 0.0178 0.3892 0.2835 -0.0715 -0.1599 0.0966 0.1112 0.1194 0.1659 0.1124 0.1484 0.7558 0.3178 0.2322 0.0196 0.0960 0.6407 -0.1508 0.0169 0.3553 0.1125 0.2540 0.2039 0.5371 0.0099 -0.2253 0.0679 -0.2710 0.4737 0.6134 0.3882 0.0179 0.7434 0.2239 -0.0931 0.0388 -0.2206 0.5996
0.4438 -0.1650 -0.0126 -0.0778 0.9582 0.4748 0.3217 -0.1900 0.1684 0.8485 0.3582 0.3992 0.0178 0.3893 0.2836 -0.0715 -0.1598 0.0966 0.1114 0.1194 0.1660 0.1125 0.1484 0.7558 0.3178 0.2323 0.0196 0.0961 0.6407 -0.1508 0.0169 0.3555 0.1125 0.2542 0.2040 0.5372 0.0100 -0.2253 0.0679 -0.2709 0.4738 0.6134 0.3884 0.0179 0.7434 0.2241 -0.0931 0.0388 -0.2205 0.5997 0.4438 -0.1650 -0.0126 -0.0778 0.9582 0.4748 0.3217 -0.1900 0.1684 0.8485 0.3582 0.3992 0.0178 0.3893 0.2836 -0.0715 -0.1598 0.0966 0.1114 0.1194 0.1660 0.1125 0.1484 0.7558 0.3178 0.2323 0.0196 0.0961 0.6407 -0.1508 0.0169 0.3555 0.1125 0.2542 0.2040 0.5372 0.0100 -0.2253 0.0679 -0.2709 0.4738 0.6134 0.3884 0.0179 0.7434 0.2241 -0.0931 0.0388 -0.2205 0.5997
</div><div class="horizontalEllipsis" style="white-space: nowrap; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"></div><div class="verticalEllipsis" style="white-space: nowrap; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"></div></div></div></div></div></div><div class="inlineWrapper"><div class = 'S11'></div></div></div><div class = 'S7'><span>For more details about the use and the limitations of </span><span style=' font-family: monospace;'>expm_multiply</span><span> please see the </span><a href = "https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/namespacematfaust_1_1poly.html#a7d06aeb92187809089a9917b6670b978"><span>API documentation</span></a><span>.</span></div><div class = 'S1'><span>Thanks for reading this notebook! Many other are available at </span><a href = "https://faust.inria.fr"><span style=' text-decoration: underline;'>faust.inria.fr</span></a><span>.</span></div><div class = 'S1'><span>This live script has been executed using the following version of matfaust:</span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div class = 'S8'><span style="white-space: pre;"><span>matfaust.version()</span></span></div><div class = 'S6'><div class="inlineElement eoOutputWrapper embeddedOutputsVariableStringElement" uid="C48497BD" data-testid="output_14" data-width="495" data-height="20" data-hashorizontaloverflow="false" style="width: 525px; max-height: 261px;"><div class="textElement"><div><span class="variableNameElement">ans = </span>'3.1.2'</div></div></div></div></div></div><div class = 'S7'></div></div><br> </div><div class="horizontalEllipsis" style="white-space: nowrap; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"></div><div class="verticalEllipsis" style="white-space: nowrap; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"></div></div></div></div></div></div><div class="inlineWrapper"><div class = 'S11'></div></div></div><div class = 'S7'><span>For more details about the use and the limitations of </span><span style=' font-family: monospace;'>expm_multiply</span><span> please see the </span><a href = "https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/namespacematfaust_1_1poly.html#a7d06aeb92187809089a9917b6670b978"><span>API documentation</span></a><span>.</span></div><div class = 'S1'><span>Thanks for reading this notebook! Many others are available at </span><a href = "https://faust.inria.fr"><span style=' text-decoration: underline;'>faust.inria.fr</span></a><span>.</span></div><div class = 'S1'><span>This live script has been executed using the following version of matfaust:</span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div class = 'S8'><span style="white-space: pre;"><span>matfaust.version()</span></span></div><div class = 'S6'><div class="inlineElement eoOutputWrapper embeddedOutputsVariableStringElement" uid="C48497BD" data-testid="output_14" data-width="495" data-height="20" data-hashorizontaloverflow="false" style="width: 525px; max-height: 261px;"><div class="textElement"><div><span class="variableNameElement">ans = </span>'3.1.2'</div></div></div></div></div></div><div class = 'S7'></div></div><br>
<!-- <!--
##### SOURCE BEGIN ##### ##### SOURCE BEGIN #####
%% Using the poly module %% Using the poly module
...@@ -115,7 +115,7 @@ ...@@ -115,7 +115,7 @@
% %
% Firstly, the |poly| module allows to define a polynomial basis (aka |FaustPoly|) % Firstly, the |poly| module allows to define a polynomial basis (aka |FaustPoly|)
% using the function |matfaust.poly.basis|. Currently, only Chebyshev polynomials % using the function |matfaust.poly.basis|. Currently, only Chebyshev polynomials
% are supported but other are yet to come. Below is the prototype of the function: % are supported but others are yet to come. Below is the prototype of the function:
% %
% |basis(L, K, basis_name, varargin)| % |basis(L, K, basis_name, varargin)|
% %
...@@ -268,7 +268,7 @@ y = expm_multiply(A, X, pts) ...@@ -268,7 +268,7 @@ y = expm_multiply(A, X, pts)
% please see the <https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/namespacematfaust_1_1poly.html#a7d06aeb92187809089a9917b6670b978 % please see the <https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/namespacematfaust_1_1poly.html#a7d06aeb92187809089a9917b6670b978
% API documentation>. % API documentation>.
% %
% Thanks for reading this notebook! Many other are available at <https://faust.inria.fr % Thanks for reading this notebook! Many others are available at <https://faust.inria.fr
% faust.inria.fr>. % faust.inria.fr>.
% %
% This live script has been executed using the following version of matfaust: % This live script has been executed using the following version of matfaust:
...@@ -277,4 +277,4 @@ matfaust.version() ...@@ -277,4 +277,4 @@ matfaust.version()
%% %%
% %
##### SOURCE END ##### ##### SOURCE END #####
--></body></html> --></body></html>
\ No newline at end of file
gen_doc/pyfaust_poly.html
+ 2
− 2
View file @ 21aa7692
...@@ -14262,7 +14262,7 @@ a.anchor-link { ...@@ -14262,7 +14262,7 @@ a.anchor-link {
<h1 id="Using-the-poly-module">Using the poly module<a class="anchor-link" href="#Using-the-poly-module">&#182;</a></h1><p>A new module has been added to pyfaust version 3.1.x. Its name is <code>poly</code> and as expected it is dedicated to a kind of Fausts that are defined according to series of polynomials.</p> <h1 id="Using-the-poly-module">Using the poly module<a class="anchor-link" href="#Using-the-poly-module">&#182;</a></h1><p>A new module has been added to pyfaust version 3.1.x. Its name is <code>poly</code> and as expected it is dedicated to a kind of Fausts that are defined according to series of polynomials.</p>
<p>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.</p> <p>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.</p>
<p><strong>NOTE</strong>: all the functions introduced in this notebook are available on GPU, using the <code>dev='gpu'</code> argument.</p> <p><strong>NOTE</strong>: all the functions introduced in this notebook are available on GPU, using the <code>dev='gpu'</code> argument.</p>
<h2 id="1.-The-basis-function">1. The basis function<a class="anchor-link" href="#1.-The-basis-function">&#182;</a></h2><p>Firstly, the <code>poly</code> module allows to define a polynomial basis (a <code>FaustPoly</code>) using the function <code>pyfaust.poly.basis</code>. Currently, only Chebyshev polynomials are supported but other are yet to come. Below is the prototype of the function:</p> <h2 id="1.-The-basis-function">1. The basis function<a class="anchor-link" href="#1.-The-basis-function">&#182;</a></h2><p>Firstly, the <code>poly</code> module allows to define a polynomial basis (a <code>FaustPoly</code>) using the function <code>pyfaust.poly.basis</code>. Currently, only Chebyshev polynomials are supported but others are yet to come. Below is the prototype of the function:</p>
<p><code>basis(L, K, basis_name, dev='cpu', T0=None)</code></p> <p><code>basis(L, K, basis_name, dev='cpu', T0=None)</code></p>
<p>In the next, we shall see a simple example but I let you consult the documentation by typing <code>help(pyfaust.poly.basis)</code> to get more details.</p> <p>In the next, we shall see a simple example but I let you consult the documentation by typing <code>help(pyfaust.poly.basis)</code> to get more details.</p>
<p>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 <code>scipy.sparse.csr_matrix</code> at least square (and most likely symmetric positive definite in much use cases), you'll make this call to the function:</p> <p>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 <code>scipy.sparse.csr_matrix</code> at least square (and most likely symmetric positive definite in much use cases), you'll make this call to the function:</p>
...@@ -14380,7 +14380,7 @@ a.anchor-link { ...@@ -14380,7 +14380,7 @@ a.anchor-link {
<div class="jp-Cell-inputWrapper"><div class="jp-InputPrompt jp-InputArea-prompt"> <div class="jp-Cell-inputWrapper"><div class="jp-InputPrompt jp-InputArea-prompt">
</div><div class="jp-RenderedHTMLCommon jp-RenderedMarkdown jp-MarkdownOutput " data-mime-type="text/markdown"> </div><div class="jp-RenderedHTMLCommon jp-RenderedMarkdown jp-MarkdownOutput " data-mime-type="text/markdown">
<p>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.</p> <p>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.</p>
<p>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 <code>Faust.clone</code> function is not used because in this case it would return a polynomial basis Faust (after all that's the role of <code>clone</code> to preserve the Faust properties!). To get a classic Faust the only way is to copy the Faust factor by factor.</p> <p>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 both cases (with F, the basis Faust and G its classic Faust copy). Note that <code>Faust.clone</code> function is not used because in this case it would return a polynomial basis Faust (after all that's the role of <code>clone</code> to preserve the Faust properties!). To get a classic Faust the only way is to copy the Faust factor by factor.</p>
   
</div> </div>
</div><div class="jp-Cell jp-CodeCell jp-Notebook-cell "> </div><div class="jp-Cell jp-CodeCell jp-Notebook-cell ">
......
gen_doc/pyfaust_poly.ipynb
+ 4
− 4
View file @ 21aa7692
%% 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 others 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). 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. 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 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
from pyfaust import Faust from pyfaust import Faust
factors = [F.factors(i) for i in range(F.numfactors())] factors = [F.factors(i) for i in range(F.numfactors())]
G = Faust(factors) 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`` (it can also be viewed as a way to compute a polynomial). 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 a 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 = rand(K+1)*100 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
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:
Depending on the situation the ``way2`` or ``way3`` might be the quickest but they are always both quicker than ``way1``. 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. 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 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). 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_1_1poly.html#aebbce7977632e3c85260738604f01104). 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