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:

It is assumed the reader is familiar with at least NumPy arrays and their operations. Besides, this notebook makes use of Pyfaust that provides 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:

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

Note that pyfaust, as well as numpy and scipy packages, are all dependencies of Lazylinop, so they are automatically installed. More information about installing Lazylinop, including (Ana)conda install or virtual environments, are available on the web site.

Contents

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

2. Creating a LazyLinOp from a Python linear operator

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

4. Vector/matrix multiplication & Computation times

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

6. Compatibility with SciPy: the SVD example

7. Prebuilt Lazylinop operators

7.1 The Kronecker product

7.2 Plenty of other prebuilt operators

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

The second key concept of Lazylinop is of course the laziness of its high-level operators. The lazy paradigm rules all the operations that occur on a 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.

It allows to review many operations with a strong focus put on sparing resources.

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:

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

We have just seen that LazyLinop can save memory.

Now let’s see that laziness can also lead to smaller computation times. We multiply a random array by ones.

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

4.34 s ± 606 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

%timeit -n 10 lz_ones @ ra

832 ms ± 27.5 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:

1. Using the aslazylinop function,

2. Using the LazyLinOp constructor.

The first one is the simplest as it relies only on a pre-existing compatible operator. The documentation of aslazylinop gives more details about what exactly a Python compatible linear operator is.

The second one is more advanced because it needs custom Python functions that define the multiplications of the operator (the so-called 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).

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

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.

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:

lF = 2 * lF
print(lF)

<3000x3000 LazyLinOp with dtype=float64>

As you can see the result of the scalar multiplication is still a 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.

Let us continue with other possible operations. For example, a matrix multiplication and then a concatenation.

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.

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

Note that even if the expression leading to 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 aLazyLinOp (yet it is more complex regarding all the operators previously aggregated).

In the next section we show how all this complexity collapses in a true calculation by applying the operator to vector.

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.

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

array([-1.76212682e+13, -2.86111608e+13, -7.97462802e+12, ...,
       -2.17079593e+13, -7.99438957e+12, -1.46546029e+13])

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

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

array([[-1.48856160e+09],
       [-2.41693494e+09],
       [-6.73658580e+08],
       ...,
       [-1.83378544e+09],
       [-6.75327713e+08],
       [-1.23795158e+09]])

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.

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

<3000x1 LazyLinOp with unspecified dtype>

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

Before comparing computation times let us mention the last but not the least thing to know about 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).

Just look:

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

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.

%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.62 s ± 228 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
1.37 s ± 247 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
9.4 s ± 350 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

As explained in section 2, instead of using a pre-defined operator like a 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).

Actually, this process mimics the SciPy LinearOperator constructor, since the 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).

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

Hence, 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).

For more details about the 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.

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

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

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.

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)

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.

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.

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)

True

The two versions of kron give the same result. Now let’s compare the computation times.

%timeit AxB @ x

49.9 ms ± 1.52 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

%timeit lAxB @ x

858 µs ± 223 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

The 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\)).

For more information about the Kronecker product in Lazylinop look at the kron documentation.

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:

• eye

• ones

• zeros

but also merge existing LazyLinOp using:

• add

• block_diag

• diag

• hstack / vstack

• kron

• pad

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:

pf.version()

'3.39.19'