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.
## 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. Creating a LazyLinearOp
In order to create this kind of object, you simply need to use the ``asLazyLinearOp`` 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). The function instantiates a ``LazyLinearOp`` that encapsulates the ``Faust`` you gave.
%% Cell type:code id: tags:
``` python
from pyfaust.lazylinop import asLazyLinearOp
import pyfaust as pf
import numpy as np
n = 4096
n = 3000
# create a random Faust
F = pf.rand(n, n, density=.001)
# create a LazyLinearOp upon it
lF = asLazyLinearOp(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 = asLazyLinearOp(S)
M = rand(n, n) + rand(n,n)*1j # numpy complex array
lM = asLazyLinearOp(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@asLazyLinearOp(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).
### 2. Applying a LazyLinearOp to a vector or a matrix
Now that we've seen how to create and operate a ``LazyLinearOp`` let's see how to apply it to a vector or a matrix, represented by a ``numpy`` array.
%% Cell type:code id: tags:
``` python
import numpy as np
v = np.arange(n)
lMSF_imag@v
```
%% Cell type:markdown id: tags:
Note the difference with the lazy multiplication by the another random vector taken as A ``LazyLinearOp``.
%% Cell type:code id: tags:
``` python
lMSF_imag@asLazyLinearOp(np.random.rand(n,1))
```
%% Cell type:markdown id: tags:
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``.
%% 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 seen 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:
