GitLab upgrade completed. Current version is 17.11.3.
Admin message
Exceptional GitLab maintenance is scheduled for June 19, 2025. The service will be unavailable between 12am and 2pm. Please do not work on the platform until an announcement indicates that maintenance is complete.
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.
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*.
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).
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.
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
## 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.
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:
%% Cell type:markdown id: tags:
### 1. Creating a LazyLinearOp
### 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.
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:
%% Cell type:code id: tags:
``` python
``` python
from pyfaust.lazylinop import asLazyLinearOp
from pyfaust.lazylinop import asLazyLinearOp
import pyfaust as pf
import pyfaust as pf
import numpy as np
import numpy as np
n = 4096
n = 3000
# create a random Faust
# create a random Faust
F = pf.rand(n, n, density=.001)
F = pf.rand(n, n, density=.001)
# create a LazyLinearOp upon it
# create a LazyLinearOp upon it
lF = asLazyLinearOp(F)
lF = asLazyLinearOp(F)
print(lF)
print(lF)
```
```
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
As said earlier, it is also possible to create ``LazyLinearOp`` operators based on ``numpy`` arrays or ``scipy`` matrices.
As said earlier, it is also possible to create ``LazyLinearOp`` operators based on ``numpy`` arrays or ``scipy`` matrices.
%% Cell type:code id: tags:
%% Cell type:code id: tags:
``` python
``` python
from scipy.sparse import random
from scipy.sparse import random
from numpy.random import rand
from numpy.random import rand
S = random(n, n, .2, format='csc') # scipy matrix
S = random(n, n, .2, format='csc') # scipy matrix
lS = asLazyLinearOp(S)
lS = asLazyLinearOp(S)
M = rand(n, n) + rand(n,n)*1j # numpy complex array
M = rand(n, n) + rand(n,n)*1j # numpy complex array
lM = asLazyLinearOp(M)
lM = asLazyLinearOp(M)
```
```
%% Cell type:markdown id: tags:
%% 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.
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:
%% Cell type:code id: tags:
``` python
``` python
try:
try:
lMSF_imag@asLazyLinearOp(np.random.rand(n))
lMSF_imag@asLazyLinearOp(np.random.rand(n))
except:
except:
print("A LazyLinearOp must have two dimensions")
print("A LazyLinearOp must have two dimensions")
```
```
%% Cell type:markdown id: tags:
%% 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.
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:
%% 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:
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:
%% Cell type:code id: tags:
``` python
``` python
lF = 2 * lF
lF = 2 * lF
print(lF)
print(lF)
```
```
%% Cell type:markdown id: tags:
%% 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.
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.
Let's continue with other possible operations. For example, a matrix multiplication and then a concatenation.
%% Cell type:code id: tags:
%% Cell type:code id: tags:
``` python
``` python
import pyfaust.lazylinop as lp
import pyfaust.lazylinop as lp
lFs = lF @ lF
lFs = lF @ lF
print("lF shape before concatenation:", lFs.shape)
print("lF shape before concatenation:", lFs.shape)
lFc = lp.vstack((lFs, lFs))
lFc = lp.vstack((lFs, lFs))
print("lF shape after concatenation:", lFc.shape)
print("lF shape after concatenation:", lFc.shape)
```
```
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
Note that we know the ``shape`` of the resulting LazyLinearOp without the need to evaluate it.
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.
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:
%% Cell type:code id: tags:
``` python
``` python
lMSF = lFc[:n, :] @ (2 * lM.conj().T @ lS)
lMSF = lFc[:n, :] @ (2 * lM.conj().T @ lS)
# then get back the imaginary part of the LazyLinearOp
# then get back the imaginary part of the LazyLinearOp
lMSF_imag = lMSF.imag
lMSF_imag = lMSF.imag
```
```
%% Cell type:markdown id: tags:
%% 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).
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
### 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.
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:
%% Cell type:code id: tags:
``` python
``` python
import numpy as np
import numpy as np
v = np.arange(n)
v = np.arange(n)
lMSF_imag@v
lMSF_imag@v
```
```
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
Note the difference with the lazy multiplication by the another random vector taken as A ``LazyLinearOp``.
Note the difference with the lazy multiplication by the another random vector taken as A ``LazyLinearOp``.
%% Cell type:code id: tags:
%% Cell type:code id: tags:
``` python
``` python
lMSF_imag@asLazyLinearOp(np.random.rand(n,1))
lMSF_imag@asLazyLinearOp(np.random.rand(n,1))
```
```
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
Instead of computing the resulting vector it gives another ``LazyLinearOp``.
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.
The vector doesn't have to be dense, a sparse one can totally be used in the multiplication.
%% Cell type:code id: tags:
%% Cell type:code id: tags:
``` python
``` python
from scipy.sparse import random as srand
from scipy.sparse import random as srand
s = srand(n,1, density=0.25)
s = srand(n,1, density=0.25)
lMSF_imag@s
lMSF_imag@s
```
```
%% Cell type:markdown id: tags:
%% 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.
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:
%% Cell type:markdown id: tags:
### 3. Comparing computation times
### 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.
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``.
Great! As expected ``lFs`` is faster to apply than its numpy array counterpart ``FD``.
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
### 4. Higher-level operations on LazyLinearOp
### 4. Higher-level operations on LazyLinearOp
%% Cell type:markdown id: tags:
%% 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.
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
#### 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.
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.
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:
%% Cell type:code id: tags:
``` python
``` python
from pyfaust.lazylinop import kron as lkron
from pyfaust.lazylinop import kron as lkron
import numpy as np
import numpy as np
from pyfaust import rand
from pyfaust import rand
A = np.random.rand(100, 100)
A = np.random.rand(100, 100)
B = np.random.rand(100, 100)
B = np.random.rand(100, 100)
AxB = np.kron(A,B)
AxB = np.kron(A,B)
lAxB = lkron(A, B)
lAxB = lkron(A, B)
x = np.random.rand(AxB.shape[1], 1)
x = np.random.rand(AxB.shape[1], 1)
np.allclose(AxB@x, lAxB@x)
np.allclose(AxB@x, lAxB@x)
```
```
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
The two versions of kron give the same result. Now let's compare the computation times.
The two versions of kron give the same result. Now let's compare the computation times.
%% Cell type:code id: tags:
%% Cell type:code id: tags:
``` python
``` python
%timeit AxB @ x
%timeit AxB @ x
```
```
%% Cell type:code id: tags:
%% Cell type:code id: tags:
``` python
``` python
%timeit lAxB @ x
%timeit lAxB @ x
```
```
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
The LazyLinearOp kron function is much faster! Indeed, it is optimized for ``LazyLinearOp-s``.
The LazyLinearOp kron function is much faster! Indeed, it is optimized for ``LazyLinearOp-s``.
%% Cell type:markdown id: tags:
%% 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).
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:
**NOTE**: this notebook was executed using the pyfaust version:
