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.

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.

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.

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:

lF = 2 * lF
print(lF)

<3000x3000 LazyLinearOp with unspecified dtype>

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.

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.

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.

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

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.

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

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

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

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.

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

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)

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.

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)

True

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

%timeit AxB @ x

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

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

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.

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

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

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

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:

• eye

• diag

• sum

• block_diag

• hstack / vstack

• kron

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:

pf.version()

'3.39.17'