Update the lazylinop notebook with explanations on creation with matmat/matvec...

Update the lazylinop notebook with explanations on creation with matmat/matvec functions and an example based on scipy fft.

Update the lazylinop notebook with explanations on creation with matmat/matvec...
92e4dec0 · hhakim · 35cc5c42 · 92e4dec0
Commit 92e4dec0 authored 2 years ago by hhakim
--- a/gen_doc/lazylinop.ipynb
+++ b/gen_doc/lazylinop.ipynb
@@ -12,7 +12,7 @@
    "\n",
    "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).\n",
    "\n",
-    "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",
+    "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",
    "## 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"
@@ -22,9 +22,11 @@
   "cell_type": "markdown",
   "metadata": {},
   "source": [
-    "### 1. Creating a LazyLinearOp\n",
+    "### 1. How to create and use a LazyLinearOp\n",
    "\n",
-    "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 ``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.\n",
+    "\n",
+    "In the example below, the function ``aslazylinearoperator`` allows us to instantiate a ``LazyLinearOp`` that encapsulates a ``Faust``."
   ]
  },
  {
@@ -347,7 +349,62 @@
   "cell_type": "markdown",
   "metadata": {},
   "source": [
-    "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``.\n",
+    "\n",
+    "#### 4.3 The FFT and its inverse as a LazyLinearOp\n",
+    "\n",
+    "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",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from pyfaust import dft\n",
+    "from pyfaust.lazylinop import LazyLinearOperator, aslazylinearoperator\n",
+    "import numpy as np\n",
+    "from scipy.fft import fft, ifft\n",
+    "n = 1024\n",
+    "lfft = LazyLinearOperator(matvec=lambda x: fft(x, axis=0), rmatvec=lambda x: n * ifft(x, axis=0), shape=(n, n))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "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.\n",
+    "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).\n",
+    "\n",
+    "We can compare this ``lfft`` operator to the equivalent operator based this time on the DFT Faust."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "F = dft(n, normed=False)\n",
+    "lF = aslazylinearoperator(F)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The two operators give the same results when applying them to a vector."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "x = np.random.rand(n)\n",
+    "np.allclose(lfft @ x, lF @ x)"
   ]
  },
  {
@@ -355,7 +412,7 @@
   "metadata": {},
   "source": [
    "\n",
-    "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).\n",
+    "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).\n",
    "\n",
    "**NOTE**: this notebook was executed using the pyfaust version:"
   ]
@@ -386,7 +443,7 @@
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
-   "version": "3.10.7"
+   "version": "3.10.8"
  }
 },
 "nbformat": 4,

 %% 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.
+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
+### 1. How to create and use 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 ``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).
 Let us mention most importantly:
 * lazy scalar multiplication
 * lazy addition
 * lazy operator multiplication
 * lazy operator concatenation
 * lazy slicing
 * lazy real/imaginary part
 * lazy operator tranpose / conjugate / transconjugate

 %% Cell type:markdown id: tags:

 ### 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@aslazylinearoperator(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.

 %% Cell type:code id: tags:

 ``` python
 %timeit lFs@M
 %timeit lFs.toarray()
 FD = lFs.toarray() # FD is a numpy array
 %timeit FD@M
 print("consistent results:", np.allclose(lFs@M, FD@M))
 ```

 %% Cell type:markdown id: tags:

 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.

 %% Cell type:code id: tags:

 ``` python
 from scipy.sparse.linalg import svds
 lF = aslazylinearoperator(pf.rand(25, 25))
 U1, S1, Vh1 = svds(lF)
 U2, S2, Vh2 = svds(lF.toarray())
 np.allclose(U1 @ np.diag(S1) @ Vh1, U2 @ np.diag(S2) @ Vh2, atol=1e-8)
 ```

 %% Cell type:markdown id: tags:

 It works the same!

 #### 4.2 The Kronecker product

 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
+import numpy as np
+from scipy.fft import fft, ifft
+n = 1024
+lfft = LazyLinearOperator(matvec=lambda x: fft(x, axis=0), rmatvec=lambda x: n * ifft(x, axis=0), shape=(n, n))
+```
+
+%% Cell type:markdown id: tags:
+
+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 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 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:

 %% Cell type:code id: tags:

 ``` python
 pf.version()
 ```