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.
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
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
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
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>
LazyLinearOpafter 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
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
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
<3000x1 LazyLinearOp with unspecified dtype>
Instead of computing the resulting vector it gives another
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
4. Higher-level operations on LazyLinearOp
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)
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) np.allclose(AxB@x, lAxB@x)
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)
kron function is much faster! Indeed, it is optimized for
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
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))
lfft is a
LazyLinearOp defined upon the scipy
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)
4.4 Other LazyLinearOp-s
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
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
NOTE: this notebook was executed using the pyfaust version: