<p>Starting from a <code>numpy</code> array, a <code>scipy</code> matrix, a <code>Faust</code> object, or potentially many other compatible linear operators with efficient implementatons, this class follows the <em>lazy evaluation paradigm</em>.</p>
<p>In short, one can <em>aggregate low-level <code>LazyLinearOp</code> objects into higher-level ones</em> using classical operations (addition, concatenation, adjoint, real part, slicing, etc.), without actually building arrays. The actual effect of these operations is delayed until the resulting linear operator is actually applied to a vector (or to a collection of vectors, seen as a matrix).</p>
<p>The main interest of this paradigm is to enable the construction of processing pipelines that exploit as building blocks efficient implementations of ``low-level'' linear operators.</p>
<p>LazyLinearOperators are complementary to other “lazy” objects such as LazyTensors in <a href="https://www.kernel-operations.io/keops/_auto_tutorials/a_LazyTensors/plot_lazytensors_a.html">Kheops</a>, or the ones of <a href="https://lazyarray.readthedocs.io/en/latest/index.html">lazyarray</a>, <a href="https://www.weld.rs/weldnumpy/]">WeldNumpy</a> libraries, which, to the best of our knowledge, primarily provide compact descriptions of arrays which entries can be evaluated efficiently on the fly.</p>
<h2 id="This-notebook">This notebook<a class="anchor-link" href="#This-notebook">¶</a></h2><p>In this notebook we shall see how to create a <code>LazyLinearOp</code> instance, create more complex instances using various lazy operations, and finally how to apply the resulting instance on vectors or matrices. We assume the reader is familiar with at least <code>numpy</code> arrays and their operations.</p>
</div>
...
...
@@ -14267,8 +14268,8 @@ Let us mention most importantly:</p>
<pre>/home/hinria/faust2/build/wrapper/python/pyfaust/__init__.py:1232: UserWarning: The numpy array to multiply is not F_CONTIGUOUS, costly conversion on the fly. Please use np.asfortranarray by yourself.
w = lambda: warnings.warn("The numpy array to multiply is not "
"The main interest of this paradigm is to enable the construction of processing pipelines that exploit as building blocks efficient implementations of ``low-level'' linear operators.\n",
"\n",
"LazyLinearOperators are complementary to other “lazy” objects such as LazyTensors in [Kheops](https://www.kernel-operations.io/keops/_auto_tutorials/a_LazyTensors/plot_lazytensors_a.html), or the ones of [lazyarray](https://lazyarray.readthedocs.io/en/latest/index.html), [WeldNumpy](https://www.weld.rs/weldnumpy/]) libraries, which, to the best of our knowledge, primarily provide compact descriptions of arrays which entries can be evaluated efficiently on the fly.\n",
"\n",
"## This notebook\n",
"In this notebook we shall see how to create a ``LazyLinearOp`` instance, create more complex instances using various lazy operations, and finally how to apply the resulting instance on vectors or matrices. We assume the reader is familiar with at least ``numpy`` arrays and their operations.\n"
]
...
...
%% Cell type:markdown id: tags:
# Using LazyLinearOp-s
The ``LazyLinearOp`` class defines a kind of linear operator extending the scipy [LinearOperator](https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.LinearOperator.html) class.
Starting from a ``numpy`` array, a ``scipy`` matrix, a ``Faust`` object, or potentially many other compatible linear operators with efficient implementatons, this class follows the *lazy evaluation paradigm*.
In short, one can *aggregate low-level ``LazyLinearOp`` objects into higher-level ones* using classical operations (addition, concatenation, adjoint, real part, slicing, etc.), without actually building arrays. The actual effect of these operations is delayed until the resulting linear operator is actually applied to a vector (or to a collection of vectors, seen as a matrix).
The main interest of this paradigm is to enable the construction of processing pipelines that exploit as building blocks efficient implementations of ``low-level'' linear operators.
LazyLinearOperators are complementary to other “lazy” objects such as LazyTensors in [Kheops](https://www.kernel-operations.io/keops/_auto_tutorials/a_LazyTensors/plot_lazytensors_a.html), or the ones of [lazyarray](https://lazyarray.readthedocs.io/en/latest/index.html), [WeldNumpy](https://www.weld.rs/weldnumpy/]) libraries, which, to the best of our knowledge, primarily provide compact descriptions of arrays which entries can be evaluated efficiently on the fly.
## This notebook
In this notebook we shall see how to create a ``LazyLinearOp`` instance, create more complex instances using various lazy operations, and finally how to apply the resulting instance on vectors or matrices. We assume the reader is familiar with at least ``numpy`` arrays and their operations.
%% Cell type:markdown id: tags:
### 1. How to create and use a LazyLinearOp
In order to create this kind of object, you simply need to use the ``aslazylinearoperator`` function. This function receives an object that represents a linear operator, for instance a ``Faust`` (but it can also be a ``numpy`` array or a ``scipy`` matrix). Besides, note that there is another way to create a LazyLinearOp using the kind of functions we call matmat or matvec as explained in 4.3 with the FFT use case.
In the example below, the function ``aslazylinearoperator`` allows us to instantiate a ``LazyLinearOp`` that encapsulates a ``Faust``.
%% Cell type:code id: tags:
``` python
from pyfaust.lazylinop import aslazylinearoperator
import pyfaust as pf
import numpy as np
n = 3000
# create a random Faust
F = pf.rand(n, n, density=.001)
# create a LazyLinearOp upon it
lF = aslazylinearoperator(F)
print(lF)
```
%% Cell type:markdown id: tags:
As said earlier, it is also possible to create ``LazyLinearOp`` operators based on ``numpy`` arrays or ``scipy`` matrices.
%% Cell type:code id: tags:
``` python
from scipy.sparse import random
from numpy.random import rand
S = random(n, n, .2, format='csc') # scipy matrix
lS = aslazylinearoperator(S)
M = rand(n, n) + rand(n,n)*1j # numpy complex array
lM = aslazylinearoperator(M)
```
%% Cell type:markdown id: tags:
It's worth noting that a ``LazyLinearOp`` must have two dimensions. Trying to instantiate a ``LazyLinearOp`` from a vector (defined with one dimension) would raise an exception, as the example below shows.
%% Cell type:code id: tags:
``` python
try:
lMSF_imag@aslazylinearoperator(np.random.rand(n))
except:
print("A LazyLinearOp must have two dimensions")
```
%% Cell type:markdown id: tags:
As a matter of fact, vectors can be interpreted as matrices but there is always an ambiguity whether this matrix has a single row or a single column.
%% Cell type:markdown id: tags:
Then we can start to build more complex ``LazyLinearOp`` objects using various operations. For example, let's multiply ``lF`` by a scalar:
%% Cell type:code id: tags:
``` python
lF = 2 * lF
print(lF)
```
%% Cell type:markdown id: tags:
As you can see the result is still a ``LazyLinearOp`` after the scalar multiplication. That's the principle of the lazy evaluation we mentioned in the beginning of this notebook. No operation is really computed, only the track of the operations is kept in a new ``LazyLinearOp`` object.
Let's continue with other possible operations. For example, a matrix multiplication and then a concatenation.
%% Cell type:code id: tags:
``` python
import pyfaust.lazylinop as lp
lFs = lF @ lF
print("lF shape before concatenation:", lFs.shape)
lFc = lp.vstack((lFs, lFs))
print("lF shape after concatenation:", lFc.shape)
```
%% Cell type:markdown id: tags:
Note that we know the ``shape`` of the resulting LazyLinearOp without the need to evaluate it.
Let's try other operations with ``lM`` and ``lS``, all ``LazyLinearOp`` are compatible with each other provided their ``shape``s are compatible.
%% Cell type:code id: tags:
``` python
lMSF = lFc[:n, :] @ (2 * lM.conj().T @ lS)
# then get back the imaginary part of the LazyLinearOp
lMSF_imag = lMSF.imag
```
%% Cell type:markdown id: tags:
For a tour of all supported operations on ``LazyLinearOp`` objects please take a look at : [LazyLinearOp reference](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/classpyfaust_1_1lazylinop_1_1LazyLinearOp.html).
Instead of computing the resulting vector it gives another ``LazyLinearOp``.
The vector doesn't have to be dense, a sparse one can totally be used in the multiplication.
%% Cell type:code id: tags:
``` python
from scipy.sparse import random as srand
s = srand(n,1, density=0.25)
lMSF_imag@s
```
%% Cell type:markdown id: tags:
One can also simply *convert* a ``LazyLinearOp`` to an equivalent numpy array using ``LazyLinearOp.toarray``. An example will come next when we compare the resulting computation times.
%% Cell type:markdown id: tags:
### 3. Comparing computation times
As a last step in this notebook, we shall verify how much computation time it takes to use a ``LazyLinearOp`` compared to a numpy array. Of course it depends on the underlying objects used behind (in the operations encoded in the ``LazyLinearOp``). Here we make the measurement on ``lFs`` which was initialized before upon a Faust object.
Great! As expected ``lFs`` is faster to apply than its numpy array counterpart ``FD``.
%% Cell type:markdown id: tags:
### 4. Higher-level operations on LazyLinearOp
%% Cell type:markdown id: tags:
Because a ``LazyLinearOp`` is a kind of scipy ``LinearOperator``, it is straightforward to use many operations provided by scipy on this type of object.
#### 4.1 The SVD
For example, let's try a SVD decomposition on a ``LazyLinearOp`` in one hand and on a numpy array on the other hand.
Another operation we can try is the Kronecker product. This time we will use the ``numpy.kron`` function on A and B numpy arrays and we will compare this function to the ``pyfaust.lazylinop.kron`` which computes the Kronecker product too but as a ``LazyLinearOp``. Precisely, we compare these functions on the multiplication of the Kronecker product by a vector.
%% Cell type:code id: tags:
``` python
from pyfaust.lazylinop import kron as lkron
import numpy as np
from pyfaust import rand
A = np.random.rand(100, 100)
B = np.random.rand(100, 100)
AxB = np.kron(A,B)
lAxB = lkron(A, B)
x = np.random.rand(AxB.shape[1], 1)
np.allclose(AxB@x, lAxB@x)
```
%% Cell type:markdown id: tags:
The two versions of kron give the same result. Now let's compare the computation times.
%% Cell type:code id: tags:
``` python
%timeit AxB @ x
```
%% Cell type:code id: tags:
``` python
%timeit lAxB @ x
```
%% Cell type:markdown id: tags:
The LazyLinearOp kron function is much faster! Indeed, it is optimized for ``LazyLinearOp-s``.
#### 4.3 The FFT and its inverse as a LazyLinearOp
Now let's explain how it is possible to create a ``LazyLinearOp`` not using a pre-defined operator like a Faust or a numpy array but a function that defines how to apply the operator on a vector (the kind of function we name matvec) or on a matrix (in which case the function is named matmat). Actually, this process mimics the scipy [LinearOperator](https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.LinearOperator.html) constructor, so it works pretty the same for a ``LazyLinearOp`` as we show in the example below for an operator that represents the FFT.
%% Cell type:code id: tags:
``` python
from pyfaust import dft
from pyfaust.lazylinop import LazyLinearOperator, aslazylinearoperator
Hence, ``lfft`` is a ``LazyLinearOp`` defined upon the scipy ``fft`` and ``ifft`` functions. Here ``fft`` is the matvec function, it defines how to apply the ``lfft`` operator to a vector. Likewise, the ``ifft``, as a ``rmatvec`` function, defines how to apply the inverse of the lfft operator to a vector.
The operator can totally be applied on matrices too. So, since we choose to apply the 1D FFT on columns instead of rows we set the axis argument of scipy fft to 0. The reason of the scaling by n of the ``ifft`` is to find on the scipy documentation [fft doc.](https://docs.scipy.org/doc/scipy/reference/generated/scipy.fft.fft.html#scipy.fft.fft) (look at the norm argument). For more details about the ``LazyLinearOperator`` constructor, please look at the [API documentation](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/namespacepyfaust_1_1lazylinop.html).
We can compare this ``lfft`` operator to the equivalent operator based this time on the DFT Faust.
%% Cell type:code id: tags:
``` python
F = dft(n, normed=False)
lF = aslazylinearoperator(F)
```
%% Cell type:markdown id: tags:
The two operators give the same results when applying them to a vector.
%% Cell type:code id: tags:
``` python
x = np.random.rand(n)
np.allclose(lfft @ x, lF @ x)
```
%% Cell type:markdown id: tags:
This notebook comes to its end. We've seen quickly how to create and evaluate ``LazyLinearOp`` objects based on numpy arrays, scipy matrices, or a Faust objects. We've also seen how to define them upon ``matmat`` and ``matvec`` functions. We've tried a bunch of operations we can call on this kind of objects. For more information about ``LazyLinearOp`` objects you can take a look to the API documentation [here](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/namespacepyfaust_1_1lazylinop.html).
**NOTE**: this notebook was executed using the pyfaust version:
