# Getting Started with Lazylinop

First of all, thank you for your interest in the Lazylinop project.

This notebook intends to guide your first steps into using this Python package. It starts with the core concept before presenting the different ways to create a `LazyLinOp`

, and ends on more advanced examples. Along the way it also shows the computational interest of Lazylinop. We hope all these examples will help you to create your own sophisticated `LazyLinOp`

-s!

**Prerequisites:**

`Faust`

objects (multi-layer sparse linear operators) and SciPy that provides sparse linear operators. However no specific knowledge is expected from the reader, in particular about Pyfaust. It only serves to explain the interoperability of Lazylinop.You obviously have to install the `lazylinop`

package by a simple:

```
[1]:
```

```
#pip install lazylinop # uncomment this line if not already installed
```

More information about installing Lazylinop, including (Ana)conda install or virtual environments, are available on the web site. Note however that `numpy`

as well as `scipy`

packages, are dependencies of Lazylinop, so they are automatically installed. For `pyfaust`

you have to install it as follows:

```
[2]:
```

```
#pip install pyfaust # uncomment this line if not already installed
```

## Contents

## 1. Core concept of Lazylinop: what is a LazyLinOp?

You probably already read the project presentation but for a good start let us recall and detail the core concept of Lazylinop: the lazy encapsulation.

### The encapsulation

First, at the core of Lazylinop, is the `LazyLinOp`

class. It defines a high-level linear operator. High-level in the sense that it encapsulates any compatible implementation. The implementation can be:

a fully defined linear operator such as a NumPy array, a SciPy matrix, a Faust or any Python compatible linear operator.

functions that simply define the multiplication by a given vector/matrix and its adjoint.

The main advantage of this encapsulation is to transparently use building blocks of linear operators including efficient third-party libraries. Indeed, once defined, a `LazyLinOp`

object is then used through a convenient NumPy-like API.

### The laziness

`LazyLinOp`

until the operator is finally applied to a vector/matrix. It means that whatever is made on the object before the multiplication, potentially no heavy computation but also no substantial memory allocation happen. It allows a high parsimony all the way from the instantiation, to the transformation and aggregation of the `LazyLinOp`

with other operators to get more complex ones. When and only when the multiplication is performed the necessary resources are consumed to proceed.Let us observe this concept at work on a very simple example. Assume for some computations you need to multiply vectors by an array of ones. Using the NumPy syntax `numpy.ones`

builds an array full of ones, which can be memory costly if its `shape = (m, n)`

is of large dimension. This can be revisited with the lazy paradigm: when building an operator using `lazylinop.ones`

, not a single `1`

is ever written in memory. Yet, this operator can be efficiently applied to any vector/matrix, as
the multiplication of an array of ones by a matrix is equivalent to summing the matrix rows. This is exactly what is encapsulated in the lazy implementation of ones, offering a great opportunity for better efficiency while preserving concise code. As we will see, LazyLinop also somehow provides a language to code complex linear operators by aggregating simpler ones with NumPy-like expressions. A convenient side effect is that when an operator L is coded, its adjoint comes handy as L.H, that’s
all.

Not clear enough yet? The code below is self-explanatory:

```
[3]:
```

```
import numpy as np
import lazylinop as lz
from sys import getsizeof
m = n = 4096
np_ones = np.ones((m, n))
lz_ones = lz.ones((m, n))
print("lz_ones takes almost no memory:", getsizeof(lz_ones), "bytes")
print("while np_ones is massive with:", np_ones.nbytes // 2**20, "Mbytes")
```

```
lz_ones takes almost no memory: 56 bytes
while np_ones is massive with: 128 Mbytes
```

`ones`

.```
[4]:
```

```
ra = np.random.rand(n, n)
%timeit np_ones @ ra
```

```
4.83 s ± 161 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
```

```
[5]:
```

```
%timeit -n 10 lz_ones @ ra
```

```
601 ms ± 70.6 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
```

Feel free to increase the values of `m`

and `n`

to experiment the lazylinop performance.

## 2. Creating a LazyLinOp from a Python linear operator

Now that the key concepts have been introduced as well as the computational interest, let us discuss how to create your own `LazyLinOp`

.

There are two ways to proceed in the Lazylinop API:

Using the aslazylinop function,

Using the LazyLinOp constructor.

`aslazylinop`

gives more details about what exactly a Python compatible linear operator is.`matmat`

, `matvec`

and `rmatmat`

, `rmatvec`

functions – only two of them are required). This constructor is presented further down in this notebook with the FFT use case (see section 5).In the example below, the function `aslazylinop`

instantiates a `LazyLinOp`

that encapsulates a `Faust`

(only considered here as a kind of linear operator).

```
[6]:
```

```
from lazylinop import aslazylinop
import pyfaust as pf
import numpy as np
n = 3000
# create a random Faust
F = pf.rand(n, n, density=.001, num_factors=10)
# create a LazyLinOp upon it
lF = aslazylinop(F)
print(lF)
```

```
<3000x3000 LazyLinOp with dtype=float64>
```

Note that our operator has a `dtype`

(the type of its scalars), it is inherited from `F`

(a `Faust`

).

As already mentioned, it is also possible to create `LazyLinOp`

operators based on NumPy arrays or SciPy matrices (because they are compatible Python linear operators).

```
[7]:
```

```
from scipy.sparse import random
from numpy.random import rand
M = rand(n, n) + rand(n,n)*1j # numpy complex array
lM = aslazylinop(M)
S = random(n, n, .2, format='csc') # scipy matrix
lS = aslazylinop(S)
```

It’s worth noting that a `LazyLinOp`

must have two dimensions, not more, not less. Trying to instantiate a `LazyLinOp`

from a vector (defined with one dimension) would raise an exception, as the example below shows.

```
[8]:
```

```
try:
# rand(n) is one-dimensional
aslazylinop(np.random.rand(n))
except:
print("A LazyLinOp must have two dimensions")
```

```
A LazyLinOp must have two dimensions
```

As a matter of fact, vectors can be interpreted as matrices (hence as linear operators) but there is always an ambiguity whether this matrix has a single row or a single column. Arrays with three dimensions or more cannot be unambiguously interpreted as linear operators.

## 3. Transformation & Aggregation: from simple to complex LazyLinOp-s

Now that we have several `LazyLinOp`

objects we can start to build more complex operators. For example, let’s multiply `lF`

by a scalar:

```
[9]:
```

```
lF = 2 * lF
print(lF)
```

```
<3000x3000 LazyLinOp with dtype=float64>
```

`LazyLinOp`

. As wanted by the lazy evaluation, no operation is really computed. Only the track of the operations is kept in a new `LazyLinOp`

object.```
[10]:
```

```
import lazylinop as lz
lFs = lF @ lF
print("lF shape before concatenation:", lFs.shape)
lFc = lz.vstack((lFs, lFs))
print("lF shape after concatenation:", lFc.shape)
```

```
lF shape before concatenation: (3000, 3000)
lF shape after concatenation: (6000, 3000)
```

Note that we know the `shape`

of the resulting `LazyLinOp`

without the need to evaluate it.

If we want to increase even more the complexity of our operator we can try other operations with `lM`

and `lS`

. Both are `LazyLinOp`

objects and their dimensions are compatible for a multiplication.

```
[11]:
```

```
lMSF = lFc[:n, :] @ (2 * lM.conj().T @ lS)
# then get back the imaginary part of the LazyLinOp
lMSF_imag = lMSF.imag
```

`lMSF`

contains multiplications (and conjugation / transposition / scalar multiplication / slicing) no evaluation/computation happen. Indeed, the computation of the multiplication happens only when the `LazyLinOp`

is applied to a matrix or a vector but here we multiply `LazyLinOp`

objects so, as a result, we get again a`LazyLinOp`

(yet it is more complex regarding all the operators previously aggregated).As you have seen the `LazyLinOp`

available operations are very similar if not exactly the same to what we are used to with NumPy. So make yourself at home and take a tour to see the other supported operations on `LazyLinOp`

objects in the LazyLinOp API reference.

Let us mention most importantly:

lazy scalar multiplication,

lazy addition and subtraction,

lazy operator multiplication,

lazy operator concatenation,

lazy slicing,

lazy real/imaginary part,

lazy operator tranpose / conjugate / adjoint.

## 4. Vector/matrix multiplication & Computation times

Now that we have created and aggregated some `LazyLinOp`

objects, it is time to apply one to a vector or a matrix, represented by a NumPy array.

```
[12]:
```

```
import numpy as np
v = np.arange(n)
lMSF_imag @ v
```

```
[12]:
```

```
array([-1.46025618e+13, -1.24211699e+13, -2.93889310e+13, ...,
-2.05780932e+13, -1.37805553e+13, -1.19268034e+13])
```

Note also that the vector doesn’t have to be dense, a sparse one can of course be used.

```
[13]:
```

```
from scipy.sparse import random as srand
s = srand(n, 1, density=0.25)
lMSF_imag @ s
```

```
[13]:
```

```
array([[-1.18629901e+09],
[-1.00908424e+09],
[-2.38753070e+09],
...,
[-1.67174562e+09],
[-1.11951918e+09],
[-9.68922342e+08]])
```

If we had multiplied the `LazyLinOp`

by the identity matrix (yes, we can of course apply a `LazyLinop`

to matrices, not only to vectors), we would have obtained the array version of the `LazyLinOp`

. As this can be useful sometimes for testing or debugging purpose, this can also be done very easily by simply calling `LazyLinOp.toarray`

(to be used with care with `LazyLinOp`

of large `shape`

as this can use a lot of memory). An example comes soon in which we compare the resulting
computation times.

Note the difference between multiplication by a vector and *lazy* multiplication by another random vector taken as a `LazyLinOp`

.

```
[14]:
```

```
lMSF_imag @ aslazylinop(np.random.rand(n,1))
```

```
[14]:
```

```
<3000x1 LazyLinOp with unspecified dtype>
```

Instead of computing the resulting vector, this lazy multiplication gives yet another `LazyLinOp`

.

`LazyLinOp`

multiplication: it also supports the multiplication by an N-dimensional array (for any \(N > 2\)). It is computed exactly as made in NumPy (see numpy.matmul).```
[15]:
```

```
r = np.random.rand(1, 2, 3, 4) # random 4-dimensional array
np.allclose(lz.ones((2, 3)) @ r,
np.ones((2, 3)) @ r)
```

```
[15]:
```

```
True
```

**Comparing computation times**

As a next step in this notebook, we shall verify how much computation time it takes to use a `LazyLinOp`

compared to a NumPy array. Of course it depends on the underlying objects and implementations used behind (in the operations encoded in the `LazyLinOp`

). Here we make the measurement on `lF`

which was initialized before upon a `Faust`

object.

```
[16]:
```

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

```
3.95 s ± 669 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
1.34 s ± 173 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
8.48 s ± 372 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
consistent results: True
```

Great! As expected `lFs`

is faster to apply than its numpy array counterpart `FD`

. Speedups are even more important when the `LazyLinOp`

is associated to a fast transform, such as the FFT below.

## 5. Creating a LazyLinOp from custom matmat/vec: the FFT example

`Faust`

or a NumPy array to create a `LazyLinOp`

, it is also possible provide functions that defines how to apply the operator (and its adjoint) on a vector or on a matrix. This kind of functions are named `matvec`

, `matmat`

etc. They are not all mandatory, for example if only `matvec`

is defined then `matmat`

is implicitly defined using `matvec`

(and vice-versa), except if you provide your own definition (usually for more computational efficiency).`LazyLinOp`

class extends the LinearOperator class. Here is an example of usage of the `LazyLinOp`

constructor to build an
operator that represents the DFT/FFT (using again SciPy).```
[17]:
```

```
from lazylinop import LazyLinOp, aslazylinop
import numpy as np
from scipy.fft import fft, ifft
n = 1024
lfft = LazyLinOp(matvec=lambda x: fft(x, axis=0), rmatvec=lambda x: n * ifft(x, axis=0), shape=(n, n))
```

`lfft`

is a `LazyLinOp`

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 be found in the SciPy documentation fft doc. (look at the `norm`

argument, there are several ways to do).`LazyLinOp`

constructor, please look at the API documentation.We can compare this `lfft`

operator to the equivalent operator based this time on the DFT Faust.

```
[18]:
```

```
from pyfaust import dft
F = dft(n, normed=False)
lF = aslazylinop(F)
```

```
[19]:
```

```
x = np.random.rand(n)
np.allclose(lfft @ x, lF @ x)
```

```
[19]:
```

```
True
```

The two operators give the same results when applying them to a vector.

When creating your own operator, if you doubt it is properly defined, don’t hesitate to rely on the method LazyLinOp.check. This method is helpful to verify your operator is consistent. For example, its shape must obviously be consistent with the shape of the array that results from the multiplication. Otherwise it means that something is wrong in your underlying implementation (in the `matmat`

function or another). The documentation of LazyLinOp.check details which verifications are performed using basic linear algebra specifications.

## 6. Compatibility with SciPy: the SVD example

Because, again, a `LazyLinOp`

is a kind of SciPy `LinearOperator`

, it is straightforward to use many operations provided by SciPy on this type of object. Note that while we did not yet extensively test that all `LinearOperator`

functions work with a `LazyLinOp`

we are working toward this direction (we are, by the way, very interested in knowing about any compatibility bug).

For example, we try here the SVD decomposition of a `LazyLinOp`

on the one hand, and of the corresponding NumPy array on the other hand. We use the SciPy function svds function for this purpose.

```
[20]:
```

```
from scipy.sparse.linalg import svds
import warnings
warnings.filterwarnings("ignore")
F = pf.rand(25, 25) # a random Faust
lF = aslazylinop(F)
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)
```

```
[20]:
```

```
True
```

It works the same!

Note that, on the contrary, the `Faust`

`F`

encapsulated in the `LazyLinOp`

`lF`

in the previous code is not eligible by itself to the SciPy `svds`

. We verify it below.

```
[21]:
```

```
try:
svds(F)
except:
print("svds didn't work on F")
```

```
svds didn't work on F
```

## 7. Prebuilt Lazylinop operators

Lazylinop provide many prebuilt linear operators all thought in the lazy way. We present one of them here and we secondly invite you to look at the Lazylinop documentation for others.

### 7.1 The Kronecker product

Another operation we can try is the Kronecker product. This time we use the `numpy.kron`

function on `A`

and `B`

NumPy arrays and we compare this function to `lazylinop.kron`

which computes the Kronecker product too but as a `LazyLinOp`

. Precisely, we compare these functions on the multiplication of the Kronecker product by a vector.

```
[22]:
```

```
import lazylinop as lz
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 = lz.kron(A, B)
x = np.random.rand(AxB.shape[1], 1)
np.allclose(AxB@x, lAxB@x)
```

```
[22]:
```

```
True
```

The two versions of `kron`

give the same result. Now let’s compare the computation times.

```
[23]:
```

```
%timeit AxB @ x
```

```
99.8 ms ± 22 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
```

```
[24]:
```

```
%timeit lAxB @ x
```

```
1.91 ms ± 549 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
```

`lz.kron`

function is much faster! Indeed, it is optimized for `LazyLinOp`

-s. The array for the Kronecker product \(A \otimes B\) never exists in memory for `lz.kron`

, only the memory space needed to compute `lAxB @ x`

is used. This computation is made on the fly using the mixed-product property \((A \otimes B) x = vec(A ~\text{reshape}(x) B^t\)).### 7.2 Plenty of other prebuilt operators

Many pre-defined `LazyLinOp`

-s are available in the `lazylinop`

API, as well as functions to aggregate them. One can notably build:

but also merge existing `LazyLinOp`

using:

For even more advanced `LinearOperator`

-s you can take a look at the documentation and specialized `LazyLinOp`

-s written for signal processing use cases or built from polynomials. Some of them are still in Work-in-Progress status but feel free to try them and provide us with any feedback. https://faustgrp.gitlabpages.inria.fr/lazylinop

This notebook comes to its end. We have seen quickly what is the concept of a `LazyLinOp`

and how it can be created, transformed, aggregated and evaluated. We’ve tried a bunch of operations we can call on this kind of objects and how they can be faster or memory more efficient. For more information about `LazyLinOp`

objects you can take a look at the API documentation here. The website provides also other tutorials.

**NOTE**: this notebook was executed using the pyfaust version:

```
[25]:
```

```
pf.version()
```

```
[25]:
```

```
'3.39.19'
```