# Using LazyLinearOp-s

The `LazyLinearOp`

class defines a kind of linear operator extending the scipy LinearOperator 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.

`LazyLinearOp`

-s are complementary to other “lazy” objects such as `LazyTensor`

-s in Kheops, or the ones of lazyarray, 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. Besides, in this notebook we make use of pyfaust, a sibling project of lazylinop.

## 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`

.

```
[1]:
```

```
from 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)
```

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

As said earlier, it is also possible to create `LazyLinearOp`

operators based on `numpy`

arrays or `scipy`

matrices.

```
[2]:
```

```
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)
```

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.

```
[3]:
```

```
try:
aslazylinearoperator(np.random.rand(n))
except:
print("A LazyLinearOp must have two dimensions")
```

```
A LazyLinearOp must have two dimensions
```

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.

Then we can start to build more complex `LazyLinearOp`

objects using various operations. For example, let’s multiply `lF`

by a scalar:

```
[4]:
```

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

```
<3000x3000 LazyLinearOp with unspecified dtype>
```

`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.```
[5]:
```

```
import lazylinop as ll
lFs = lF @ lF
print("lF shape before concatenation:", lFs.shape)
lFc = ll.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 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 dimensions (`shape`

) are compatible.

```
[6]:
```

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

For a tour of all supported operations on `LazyLinearOp`

objects please take a look at : LazyLinearOp reference. 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

## 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.

```
[7]:
```

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

```
[7]:
```

```
array([-4.38679094e+11, -4.60847874e+11, -3.33996056e+11, ...,
-3.25717562e+11, -4.68327333e+11, -2.27455345e+11])
```

Note the difference with the lazy multiplication by another random vector taken as a `LazyLinearOp`

.

```
[8]:
```

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

```
[8]:
```

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

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.

```
[9]:
```

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

```
[9]:
```

```
array([[-36668142.40986451],
[-38520886.21557115],
[-27917523.94755967],
...,
[-27226239.2297844 ],
[-39146762.73634297],
[-19012509.94659325]])
```

One can also simply *convert* a `LazyLinearOp`

to an equivalent numpy array using `LazyLinearOp.toarray`

. An example will come next when we’ll compare the resulting computation times.

## 3. Comparing computation times

As a next 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.

```
[10]:
```

```
%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))
```

```
3.59 s ± 733 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
972 ms ± 31.4 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
6.36 s ± 653 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`

.

## 4. Higher-level operations on LazyLinearOp

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.

```
[11]:
```

```
from scipy.sparse.linalg import svds
import warnings
warnings.filterwarnings("ignore")
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)
```

```
[11]:
```

```
True
```

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 `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.

```
[12]:
```

```
from 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)
```

```
[12]:
```

```
True
```

The two versions of `kron`

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

```
[13]:
```

```
%timeit AxB @ x
```

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

```
[14]:
```

```
%timeit lAxB @ x
```

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

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 constructor, so it works pretty the same
for a `LazyLinearOp`

as we show in the example below for an operator that represents the FFT.

```
[15]:
```

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

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 in the scipy documentation fft doc. (look at the norm argument). For more details about the `LazyLinearOp`

constructor, please look at the API documentation.

We can compare this `lfft`

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

```
[16]:
```

```
F = dft(n, normed=False)
lF = aslazylinearoperator(F)
```

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

```
[17]:
```

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

```
[17]:
```

```
True
```

### 4.4 Other LazyLinearOp-s

Many pre-defined `LazyLinearOp`

-s are available in the `lazylinop`

API. Most of the functions to build them are listed below:

For other advanced `LinearOperator`

-s you might want to look into the pylops API. The pylops library is compatible to scipy `LinearOperator`

and so also to our `LazyLinearOp`

-s. Besides, pylops is also lazyness oriented.

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. We’ve also seen how to create structured lazy operators like a Kronecker product or FFT. For more information about `LazyLinearOp`

objects you can take a look to the API documentation
here.

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

```
[18]:
```

```
pf.version()
```

```
[18]:
```

```
'3.39.17'
```