Lazylinop core module

Lazylinop core module#

  • LazyLinOp: the core class for lazy linear operators.

Available operations

+ (addition), - (subtraction), @ (matrix product), * (scalar multiplication), ** (matrix power for square operators), indexing, slicing and others. For a nicer introduction you might look at these tutorials.

LazyLinOp class#

class lazylinop.LazyLinOp(shape, matvec=None, matmat=None, rmatvec=None, rmatmat=None)#

The LazyLinOp class.

The lazy principle

The evaluation of any defined operation on a LazyLinOp is delayed until a multiplication by a matrix/vector or a call of LazyLinOp.toarray() is made.

Two ways to instantiate

Available operations

+ (addition), - (subtraction), @ (matrix product), * (scalar multiplication), ** (matrix power for square operators), indexing, slicing and others. For a nicer introduction you might look at these tutorials.

Recursion limit

Repeated “inplace” modifications of a LazyLinOp through any operation like a concatenation (op = vstack((op, anything))) are subject to a RecursionError if the number of recursive calls exceeds sys.getrecursionlimit(). You might change this limit if needed using sys.setrecursionlimit().

A LazyLinOp instance is defined by a shape and at least functions matvec or matmat and rmatvec or rmatmat.

Parameters#

shape: (tuple[int, int])

Operator \(L\) dimensions \((M, N)\).

matvec: (callable)

Returns \(y = L * v\) with \(v\) a vector of size \(N\). \(y\) size is \(M\) with the same number of dimension(s) as \(v\).

rmatvec: (callable)

Returns \(y = L^H * v\) with \(v\) a vector of size \(M\). \(y\) size is \(N\) with the same number of dimension(s) as \(v\).

matmat: (callable)

Returns \(L * V\). The output matrix shape is \((M, K)\).

rmatmat: (callable)

Returns \(L^H * V\). The output matrix shape is \((N, K)\).

Auto-implemented operations

  • If only matvec is defined and not matmat, an automatic naive matmat will be defined upon the given matvec but note that it might be suboptimal (in which case a matmat is useful). The same applies for rmatvec and rmatmat.

  • No need to provide the implementation of the multiplication by a LazyLinOp, or a numpy array with ndim > 2 because both of them are auto-implemented. For the latter operation, it is computed as in numpy.__matmul__.

Return:

LazyLinOp

Example:
>>> # In this example we create a LazyLinOp
>>> # for the DFT using the fft from scipy
>>> import numpy as np
>>> from scipy.fft import fft, ifft
>>> from lazylinop import LazyLinOp
>>> fft_mm = lambda x: fft(x, norm='ortho')
>>> fft_rmm = lambda x: ifft(x, norm='ortho')
>>> n = 16
>>> F = LazyLinOp((n, n), matvec=fft_mm, rmatvec=fft_rmm)
>>> x = np.random.rand(n)
>>> y = F @ x
>>> np.allclose(y, fft(x, norm='ortho'))
True
>>> np.allclose(x, F.H @ y)
True

See also

SciPy linear Operator. LazyLinOp.check() scipy fft

Attributes#

LazyLinOp.T#

The LazyLinOp transpose.

LazyLinOp.H#

The LazyLinOp adjoint/transconjugate.

LazyLinOp.ndim#

The number of dimensions of the LazyLinOp (it is always 2).

LazyLinOp.shape#
LazyLinOp.real#

Returns the LazyLinOp real part.

LazyLinOp.imag#

Returns the LazyLinOp imaginary part.

Methods#

LazyLinOp main methods: LazyLinOp.toarray(), LazyLinOp.conj(), LazyLinOp.check(), LazyLinOp.__pow__().

lazylinop.LazyLinOp.toarray(self, dtype=None, array_namespace=None, device=None)#

Returns self as an array

Internally, it computes self @ array_namespace.eye(self.shape[1], dtype=dtype), with array_namespace and dtype depending on arguments.

with smallest possible dtype and returns self as a NumPy/CuPy array or torch tensor. dtype of the output y depends on the LazyLinOp instance self.

Args:
dtype: str, NumPy/CuPy or torch dtype, optional

The dtype used eye() in eye()

Default value is None and will select the smallest possible dtype.

Note that dtype of the returned array depends on the LazyLinOp instance self.

array_namespace: namespace, optional

The type of the return array, as an array API namespace (NumPy, torch, CuPy, …) Default is None, which would use NumPy.

device: str, optional

The device where the returned array will be computed

Default value is None, the array will reside on CPU

Note: device has no effect if array_namespace is not equal to 'torch'.

Examples:
>>> import numpy as np
>>> from lazylinop import aslazylinop
>>> L = aslazylinop(np.eye(2, dtype='int'))
>>> L.toarray(array_namespace=np, dtype='float')
array([[1., 0.],
       [0., 1.]])
lazylinop.LazyLinOp.conj(self)#

Returns the LazyLinOp conjugate.

lazylinop.LazyLinOp.check(self, array_namespace=None, dtype='float64', device=None, rtol=1e-05, atol=1e-08)#

Verifies validity assertions on any LazyLinOp.

Notations:

  • Let op a LazyLinOp,

  • u, v vectors such that u.shape[0] == op.shape[1]

and v.shape[0] == op.shape[0], - X, Y 2d-arrays such that X.shape[0] == op.shape[1] and Y.shape[0] == op.shape[0].

The function verifies:

  • Consistency of operator/adjoint product shape:

      1. (op @ u).shape == (op.shape[0],),

      2. (op.H @ v).shape == (op.shape[1],),

      1. (op @ X).shape == (op.shape[0], X.shape[1]),

      2. (op.H @ Y).shape == (op.shape[1], Y.shape[1]),

  • Consistency of operator & adjoint products:

    1. (op @ u).conj().T @ v == u.conj().T @ op.H @ v

  • Consistency of operator-by-matrix & operator-by-vector products:

    1. op @ X is equal to the horizontal concatenation of all op @ X[:, j] (\(0\le j < X.shape[1]\)).

      (it implies also that (op @ X).shape[1] == X.shape[1], as previously verified in 2.a)

  • Consistency of adjoint-by-matrix & adjoint-by-vector products:

    1. op.H @ Y is equal to the horizontal concatenation of all op.H @ Y[:, j] (\(0\le j < Y.shape[1]\)).

      (it implies also that (op.H @ Y).shape[1] == Y.shape[1], as previously verified in 2.b)

  • Linearity:

    1. op @ (a1 * u1 + a2 * u2) == a1 * (op @ u1) + a2 * (op @ u2).

  • Device:

    1. array_api_compat.device(x) == array_api_compat.device(op @ x)

Raises:

  • Exception("Operator shape[0] and operator-by-vector product shape must agree") (assertion 1.a)

  • Exception("Operator shape[1] and adjoint-by-vector product shape must agree") (assertion 1.b)

  • Exception("Operator-by-matrix product shape and operator/input-matrix shape must agree") (assertion 2.a)

  • Exception("Operator-by-matrix & operator-by-vector products must agree") (assertion 2.b)

  • Exception("Operator and adjoint products do not match") (assertion 3)

  • Exception("Operator-by-matrix & operator-by-vector products must agree") (assertion 4)

  • Exception("Adjoint-by-matrix product shape and adjoint/input-matrix shape must agree") (assertion 5)

Computational cost

This function has a computational cost of several matrix products. It shouldn’t be used into an efficient implementation but only to test a LazyLinOp implementation is valid.

Necessary condition but not sufficient

This function is able to detect an inconsistent LazyLinOp according to the assertions above but it cannot ensure a particular operator computes what someone is excepted this operator to compute. In other words, the operator can be consistent but not correct at the same time. Thus, this function is not enough by itself to write unit tests for an operator, complementary tests are necessary.

Args:
self: (LazyLinOp)

Operator to test.

array_namespace: namespace, optional

Namespace of the input to test self (see Examples section).

  • None which would use NumPy (defaut value).

  • np NumPy namespace import numpy as np.

  • cp CuPy namespace import cupy as cp.

  • torch PyTorch namespace import torch.

  • list or array namespaces.

dtype: str, NumPy/CuPy or torch dtype, optional

dtype of the input that will be used to test self.

device: optional

Use device device to run self.check(...). Default value is None.

rtol: float, optional

The relative tolerance parameter used by NumPy allclose function to verify validity assertions. Default value is 1e-5. See NumPy allclose for more details.

atol: float, optional

The absolute tolerance parameter used by NumPy allclose function to verify validity assertions. Default value is 1e-8. See NumPy allclose for more details.

Example:
>>> import numpy as np
>>> from numpy.random import rand
>>> from lazylinop import aslazylinop, LazyLinOp
>>> M = rand(12, 14)
>>> # numpy array M is OK as a LazyLinOp
>>> aslazylinop(M).check(array_namespace=np)
>>> # the next LazyLinOp is not
>>> L2 = LazyLinOp((6, 7), matmat=lambda x: np.ones((6, 7)), rmatmat=lambda x: np.zeros((7,6)))
>>> L2.check(array_namespace=np) 
Traceback (most recent call last):
    ...
Exception: ...

See also

aslazylinop(), LazyLinOp.

lazylinop.LazyLinOp.__pow__(self, n)#

Returns the LazyLinOp for the n-th power of self.

  • L**n == L @ L @ ... @ L (n-1 multiplications).

Args:

n: a positive integer

Raises:

The LazyLinOp is not square.

Example:
>>> from lazylinop import aslazylinop
>>> import numpy as np
>>> M = np.random.rand(10, 10).astype('float32')
>>> lM = aslazylinop(M)
>>> lM
<10x10 ArrayBasedLazyLinOp with dtype=float32>
>>> np.allclose((lM**2).toarray(), M @ M)
True

Module utility functions#

lazylinop.islazylinop(op)#
lazylinop.aslazylinop(op)#

Basic operators#

Construction#

Functions for creating a Lazylinop from scratch (providing the related parameters).

Important

  • Note that lazylinop provides two closely named functions pad and padder:

    • pad takes a lazylinop as input and returns a padded lazylinop, and is part of a family of tools to “glue” various types of lazylinops together.

    • padder builds a lazylinop that maps a (batch of) vector(s) to their (batch of) padded versions.

  • In general, function names such as slicer, indexer, padder (see also padder2d) also refer to functions that return a LazyLinop which action is described by the verb appearing in their name.

lazylinop.anti_diag(v, k=0, extract_meth='canonical_vectors', extract_batch=1)#

Returns a LazyLinOp L that extracts an antidiagonal or builds an antidiagonal .

The shape of L is square and depends on the size of v.

Args:
v: (compatible linear operator, 1D numpy.ndarray)

  • If v is a LazyLinOp or an array-like compatible object, returns a copy of its k-th antidiagonal.

  • If v is a 1D numpy array, returns a LazyLinOp with v on the k-th antidiagonal. L @ x returns an error if v and x are not on the same device where L = anti_diag(v, k).

k: int, optional

The index of antidiagonal, 0 (default) for the main antidiagonal (the one starting from the upper right corner), k > 0 for upper antidiagonals, k < 0 for lower antidiagonals below (see anti_eye()).

extract_meth: str, optional

The method used to extract the antidiagonal vector. The interest to have several methods resides in their difference of memory and execution time costs but also on the operator capabilities (e.g. not all of them support a CSC matrix as multiplication operand).

  • 'canonical_vectors': use canonical basis vectors \(e_i\) to extract each antidiagonal element of the operator. It takes an operator-vector multiplication to extract each antidiagonal element.

  • 'canonical_vectors_csc': The same as above but using scipy CSC matrices to encode the canonical vectors. The memory cost is even smaller than that of 'canonical_vectors'. However v must be compatible with CSC matrix-vector multiplication.

  • 'slicing': extract antidiagonal elements by slicing rows and columns by blocks of shape (extract_batch, extract_batch).

  • 'toarray': use LazyLinOp.toarray() to extract the antidiagonal after a conversion to a whole numpy array.

extract_batch: int, optional

  • The size of the batch used for partial antidiagonal extraction in 'canonical_vectors', 'canonical_vectors_csc' and 'slicing ' methods.

    This argument is ignored for 'toarray' method.

Antidiagonal extraction cost

Even though the 'toarray' method is generally faster if the operator is not extremely large it has an important memory cost (\(O(v.shape[0] \times v.shape[1])\)) . Hence the default method is canonical_vectors in order to avoid memory consumption. However note that this method allows to define a memory-time trade-off with the extract_batch argument. The larger is the batch, the faster should be the execution (provided enough memory is available).

Returns:

The extracted antidiagonal numpy vector or the constructed antidiagonal LazyLinOp.

Example: (antidiagonal LazyLinOp)
>>> import lazylinop.basicops as lz
>>> import numpy as np
>>> v = np.arange(1, 6)
>>> v
array([1, 2, 3, 4, 5])
>>> ld1 = lz.anti_diag(v)
>>> ld1
<5x5 LazyLinOp with unspecified dtype>
>>> ld1.toarray('int')
array([[0, 0, 0, 0, 1],
       [0, 0, 0, 2, 0],
       [0, 0, 3, 0, 0],
       [0, 4, 0, 0, 0],
       [5, 0, 0, 0, 0]])
>>> ld2 = lz.anti_diag(v, -2)
>>> ld2
<7x7 LazyLinOp with unspecified dtype>
>>> ld2.toarray('int')
array([[0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 1],
       [0, 0, 0, 0, 0, 2, 0],
       [0, 0, 0, 0, 3, 0, 0],
       [0, 0, 0, 4, 0, 0, 0],
       [0, 0, 5, 0, 0, 0, 0]])
>>> ld3 = lz.anti_diag(v, 2)
>>> ld3
<7x7 LazyLinOp with unspecified dtype>
>>> ld3.toarray('int')
array([[0, 0, 0, 0, 1, 0, 0],
       [0, 0, 0, 2, 0, 0, 0],
       [0, 0, 3, 0, 0, 0, 0],
       [0, 4, 0, 0, 0, 0, 0],
       [5, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0]])
Example: (antidiagonal extraction)
>>> import lazylinop.basicops as lz
>>> import numpy as np
>>> lD = aslazylinop(np.random.rand(10, 12))
>>> d = lz.anti_diag(lD, -2, extract_meth='toarray', extract_batch=3)
>>> # verify d is really the antidiagonal of index -2
>>> d_ = np.diag(np.fliplr(lD.toarray()), -2)
>>> np.allclose(d, d_)
True

See also

diag() numpy.diag numpy.fliplr aslazylinop()

lazylinop.anti_eye(M, N=None, k=0)#

Constructs a LazyLinOp whose equivalent array is filled with ones on the k-th antidiagonal and zeros everywhere else.

L = anti_eye(M, N, k) is such that L.toarray() == numpy.flip(numpy.eye(M, N, k), axis=1).

Args:
M: int

Number of rows.

N: int, optional

Number of columns (default is M).

k: int, optional

Anti-diagonal to place ones on.

  • zero for the main antidiagonal (default), the one starting from the upper right corner.

  • positive integer for an upper antidiagonal.

  • negative integer for a lower antidiagonal,

if k >= N or k <= - M then anti_eye(M, N, k) is zeros((M, N)) (k-th antidiagonal is out of operator shape).

Returns:

The anti-eye LazyLinOp.

Examples:
>>> import lazylinop as lz
>>> import numpy as np
>>> x = np.arange(3)
>>> L = lz.basicops.anti_eye(3)
>>> np.allclose(L @ x, np.flip(x))
True
>>> # Check the main diagonal is the one starting
>>> # from the upper right corner.
>>> L = lz.basicops.anti_eye(3, 4, k=0)
>>> L.toarray(dtype='float')
array([[0., 0., 0., 1.],
       [0., 0., 1., 0.],
       [0., 1., 0., 0.]])
>>> L = lz.basicops.anti_eye(3, N=3, k=0)
>>> L.toarray(dtype='int')
array([[0, 0, 1],
       [0, 1, 0],
       [1, 0, 0]])
>>> L = lz.basicops.anti_eye(3, N=3, k=1)
>>> L.toarray(dtype='int')
array([[0, 1, 0],
       [1, 0, 0],
       [0, 0, 0]])
>>> L = lz.basicops.anti_eye(3, N=3, k=-1)
>>> L.toarray(dtype='int')
array([[0, 0, 0],
       [0, 0, 1],
       [0, 1, 0]])
>>> L = lz.basicops.anti_eye(3, N=4, k=0)
>>> L.toarray(dtype='int')
array([[0, 0, 0, 1],
       [0, 0, 1, 0],
       [0, 1, 0, 0]])
>>> L = lz.basicops.anti_eye(3, N=4, k=1)
>>> L.toarray(dtype='int')
array([[0, 0, 1, 0],
       [0, 1, 0, 0],
       [1, 0, 0, 0]])
>>> L = lz.basicops.anti_eye(3, N=4, k=-1)
>>> L.toarray(dtype='int')
array([[0, 0, 0, 0],
       [0, 0, 0, 1],
       [0, 0, 1, 0]])

See also

eye(), LazyLinOp.toarray(), numpy.eye, numpy.flip.

lazylinop.diag(v, k=0, extract_meth='canonical_vectors', extract_batch=1)#

Extracts a diagonal or constructs a diagonal LazyLinOp.

Args:
v: (compatible linear operator, 1D numpy.ndarray)

  • If v is a LazyLinOp or an array-like compatible object, returns a copy of its k-th diagonal.

  • If v is a 1D NumPy/CuPy array or torch tensor, returns a LazyLinOp with v on the k-th diagonal. L @ x returns an error if v and x are not on the same device where L = diag(v, k).

k: (int)

The index of diagonal, 0 for the main diagonal, k > 0 for diagonals above, k < 0 for diagonals below (see eye()).

extract_meth: (str)

The method used to extract the diagonal vector. The interest to have several methods resides in their difference of memory and execution time costs but also on the operator capabilities (e.g. not all of them support a CSC matrix as multiplication operand).

  • 'canonical_vectors': use canonical basis vectors \(e_i\) to extract each diagonal element of the operator. It takes an operator-vector multiplication to extract each diagonal element.

  • 'canonical_vectors_csc': The same as above but using scipy CSC matrices to encode the canonical vectors. The memory cost is even smaller than that of 'canonical_vectors'. However v must be compatible to CSC matrices multiplication.

  • 'slicing': extract diagonal elements by slicing rows and columns by blocks of shape (extract_batch, extract_batch).

  • 'toarray': use LazyLinOp.toarray() to extract the diagonal after a conversion to a whole numpy array.

extract_batch: (int)

  • The size of the batch used for partial diagonal extraction in 'canonical_vectors', 'canonical_vectors_csc' and 'slicing ' methods.

    This argument is ignored for 'toarray' method.

Diagonal extraction cost

Even though the 'toarray' method is generally faster if the operator is not extremely large it has an important memory cost (\(O(v.shape[0] \times v.shape[1])\)) . Hence the default method is canonical_vectors in order to avoid memory consumption. However note that this method allows to define a memory-time trade-off with the extract_batch argument. The larger is the batch, the faster should be the execution (provided enough memory is available).

Returns:

The extracted diagonal numpy vector or the constructed diagonal LazyLinOp.

Example: (diagonal LazyLinOp)
>>> import lazylinop as lz
>>> import numpy as np
>>> v = np.arange(1, 6)
>>> v
array([1, 2, 3, 4, 5])
>>> ld1 = lz.diag(v)
>>> ld1
<5x5 LazyLinOp with unspecified dtype>
>>> ld1.toarray('int')
array([[1, 0, 0, 0, 0],
       [0, 2, 0, 0, 0],
       [0, 0, 3, 0, 0],
       [0, 0, 0, 4, 0],
       [0, 0, 0, 0, 5]])
>>> ld2 = lz.diag(v, -2)
>>> ld2.toarray('int')
array([[0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0],
       [1, 0, 0, 0, 0, 0, 0],
       [0, 2, 0, 0, 0, 0, 0],
       [0, 0, 3, 0, 0, 0, 0],
       [0, 0, 0, 4, 0, 0, 0],
       [0, 0, 0, 0, 5, 0, 0]])
>>> ld3 = lz.diag(v, 2)
>>> ld3.toarray('int')
array([[0, 0, 1, 0, 0, 0, 0],
       [0, 0, 0, 2, 0, 0, 0],
       [0, 0, 0, 0, 3, 0, 0],
       [0, 0, 0, 0, 0, 4, 0],
       [0, 0, 0, 0, 0, 0, 5],
       [0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0]])
Example: (diagonal extraction)
>>> import lazylinop as lz
>>> from lazylinop import aslazylinop
>>> import numpy as np
>>> M = np.random.rand(10, 12)
>>> lD = aslazylinop(M)
>>> d = lz.diag(lD, -2)
>>> # verify d is really the diagonal of index -2
>>> d_ = np.array([M[i, i-2] for i in range(abs(-2), lD.shape[0])])
>>> np.allclose(d, d_)
True

See also

numpy.diag aslazylinop()

lazylinop.diff(N, n=1, prepend=0, append=0, backend='numpy')#

Returns a LazyLinOp L that calculates the n-th discrete difference of an input vector.

Shape of L is \((M,~N)\) where \(M = prepend + N + append - n\).

Args:
N: int

Size of the input.

nint, optional

The number of times values are differenced (default is 1). If zero, the input is returned as-is.

prepend, appendint, optional

Prepend or append input vector with a number of zeros equals to the argument, prior to performing the difference. By default it is equal to 0.

Returns:

LazyLinOp

See also

NumPy diff function.

Examples:
>>> import numpy as np
>>> import lazylinop as lz
>>> x = np.array([1, 2, 2, 4, 2, 2, 1])
>>> D = diff(x.shape[0])
>>> y = D @ x
>>> np.allclose(np.array([1,  0,  2, -2,  0, -1]), y)
True
lazylinop.eye(M, N=None, k=0)#

Returns the LazyLinOp L for eye (identity matrix and variants).

Args:
M: int

Number of rows.

N: int, optional

Number of columns. Default is M.

k: int, optional

Diagonal to place ones on.

  • zero for the main diagonal (default),

  • positive integer for an upper diagonal,

  • negative integer for a lower diagonal.

Example:
>>> import lazylinop as lz
>>> le1 = lz.eye(5)
>>> le1
<5x5 LazyLinOp with unspecified dtype>
>>> le1.toarray(dtype='float')
array([[1., 0., 0., 0., 0.],
       [0., 1., 0., 0., 0.],
       [0., 0., 1., 0., 0.],
       [0., 0., 0., 1., 0.],
       [0., 0., 0., 0., 1.]])
>>> le2 = lz.eye(5, 2)
>>> le2
<5x2 LazyLinOp with unspecified dtype>
>>> le2.toarray(dtype='int')
array([[1, 0],
       [0, 1],
       [0, 0],
       [0, 0],
       [0, 0]])
>>> le3 = lz.eye(5, 3, 1)
>>> le3
<5x3 LazyLinOp with unspecified dtype>
>>> le3.toarray(dtype='int')
array([[0, 1, 0],
       [0, 0, 1],
       [0, 0, 0],
       [0, 0, 0],
       [0, 0, 0]])
>>> le4 = lz.eye(5, 3, -1)
>>> le4
<5x3 LazyLinOp with unspecified dtype>
>>> le4.toarray(dtype='int')
array([[0, 0, 0],
       [1, 0, 0],
       [0, 1, 0],
       [0, 0, 1],
       [0, 0, 0]])

See also

scipy.sparse.eye, numpy.eye.

lazylinop.ones(shape)#

Returns a LazyLinOp ones.

Fixed memory cost

Whatever is the shape of the ones, it has the same memory cost.

Args:

shape: tuple[int, int]

Operator shape, e.g., (2, 3).

Returns:

LazyLinOp ones.

Example:
>>> import numpy as np
>>> import lazylinop as lz
>>> L = lz.ones((6, 5))
>>> v = np.arange(5)
>>> v
array([0, 1, 2, 3, 4])
>>> L @ v
array([10, 10, 10, 10, 10, 10])
>>> Oa = np.ones((6, 5)).astype('int')
>>> Oa @ v
array([10, 10, 10, 10, 10, 10])
>>> M = np.arange(5 * 4).reshape(5, 4)
>>> L @ M
array([[40, 45, 50, 55],
       [40, 45, 50, 55],
       [40, 45, 50, 55],
       [40, 45, 50, 55],
       [40, 45, 50, 55],
       [40, 45, 50, 55]])
>>> Oa @ M
array([[40, 45, 50, 55],
       [40, 45, 50, 55],
       [40, 45, 50, 55],
       [40, 45, 50, 55],
       [40, 45, 50, 55],
       [40, 45, 50, 55]])

See also

numpy.ones

lazylinop.tri(N, M=None, k=0)#

Returns a LazyLinOp L associated to a matrix such that its lower triangular part is filled with ones. Above given diagonal k, matrix is filled with zeros.

Shape of L is \((N,~M)\) where \(M = N\) by default.

LazyLinOp L = tri(N, k=0) corresponds to a lower triangular matrix of shape \((N,~N)\) and its transposition L.T corresponds to an upper triangular matrix.

Args:
N: int

Number of rows.

M: int, optional

Number of columns. By default M is equal to N.

k: int, optional

The sub-diagonal at and below which the matrix is filled with ones (the rest being filled by zeroes).

  • \(k<0\) below the main diagonal.

  • \(k=0\) corresponds to the main diagonal (default). When called with \(k=0\) (default), the resulting operator L performs a cumulative sum.

  • \(k>0\) above the main diagonal.

Returns:

LazyLinOp

Examples:
>>> import numpy as np
>>> import lazylinop as lz
>>> N = 4
>>> x = np.arange(N)
>>> x
array([0, 1, 2, 3])
>>> L = lz.basicops.tri(N)
>>> y = L @ x
>>> y_ = np.cumsum(x)
>>> np.allclose(y, y_)
True

See also

lazylinop.cumsum(N)#

Returns a LazyLinOp L associated to the cumulative sum of an input of size N.

Shape of L is \((N,~N)\).

Args:
N: int

Size of the input array.

Returns:

LazyLinOp

Examples:
>>> import numpy as np
>>> import lazylinop as lz
>>> N = 4
>>> x = np.arange(N)
>>> x
array([0, 1, 2, 3])
>>> L = lz.basicops.tri(N)
>>> y = L @ x
>>> y_ = np.cumsum(x)
>>> np.allclose(y, y_)
True
lazylinop.zeros(shape)#

Returns a zero LazyLinOp.

Fixed memory cost

Whatever is the shape of the zeros, it has the same memory cost.

Args:
shape: (tuple[int, int])

The operator shape.

Example:
>>> import numpy as np
>>> import lazylinop as lz
>>> Lz = lz.zeros((10, 12))
>>> x = np.random.rand(12)
>>> Lz @ x
array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])

See also

numpy.zeros.

lazylinop.padder(N, width=(0, 0), mode='zero')#

Returns a LazyLinOp L that extends a signal with either zero, periodic, symmetric, antisymmetric or reflect boundary conditions.

Shape of L is \((M,~N)\). By default \(M = N\).

According to mode and value of width, padder() will add width[0] // N “copies” of the signal plus width[0] % N elements before and width[1] // N “copies” of the signal plus width[1] % N elements after.

For an input x with entries \(x_0, x_1, \ldots, x_{N-1}\), the output y = L @ x depends on the mode.

In the general case, width is a multiple of input length N:

  • with a periodic boundary condition and width = (N, N), the entries of y are:

    \(x_0, x_1, ..., x_{N-1} | x_0, x_1, ..., x_{N-1} | x_0, x_1, ..., x_{N-1}\)

  • with a symmetric boundary condition and width = (N, N), y has entries:

    \(x_{N-1}, ..., x_1, x_0 | x_0, x_1, ..., x_{N-1} | x_{N-1}, ..., x_1, x_0\)

  • with an antisymmetric boundary condition and width = (N, N), y has entries:

    \(-x_{N-1}, ..., -x_1, -x_0 | x_0, x_1, ..., x_{N-1} | -x_{N-1}, ..., -x_1, -x_0\)

  • with a reflect boundary condition and width = (N - 1, N - 1), the entries of y are:

    \(x_{N-1}, ..., x_1 | x_0, x_1, ..., x_{N-1} | x_{N-2}, x_{N-3}, ..., x_0\)

In the case, width is lesser than the input length N we have:

  • with a periodic boundary condition and width = (b, a), b < N - 1 and a < N - 1, the entries of y are:

    \(x_{N-1-b+1}, ..., x_{N-1} | x_0, x_1, ..., x_{N-1} | x_0, ..., x_a\)

  • with a symmetric boundary condition and width = (b, a), b < N - 1 and a < N - 1, the entries of y are:

    \(x_b, ..., x_0 | x_0, x_1, ..., x_{N-1} | x_{N-1}, ..., x_{N-1-a+1}\)

  • with an antisymmetric boundary condition and width = (b, a), b < N - 1 and a < N - 1, the entries of y are:

    \(-x_b, ..., -x_0 | x_0, x_1, ..., x_{N-1} | -x_{N-1}, ..., -x_{N-1-a+1}\)

  • with a reflect boundary condition and width = (b, a), 1 < b < N - 2 and 1 < a < N - 2, the entries of y are:

    \(x_b, ..., x_1 | x_0, x_1, ..., x_{N-1} | x_{N-2}, ..., x_{N-1-a}\)

y = L @ X with X a 2D array (corresponding to an image) will only extend each column of X, since this is the normal behaviour of a LazyLinOp.

Args:
N: int

Size of the input.

width: tuple, optional

Number of values padded on both side of the input (before, after). By default it is equal to (0, 0). The size of the output is width[0] + N + width[1]. width[0] and width[1] must be greater or equal to zero.

mode: str, optional

zero (default), 'wrap'/'periodic', 'symm'/'symmetric', 'antisymmetric' or 'reflect' boundary condition.

Returns:

LazyLinOp of shape \((M,~N)\) where \(M = before + N + after\) (default is \((3N,~N)\)).

Examples:
>>> import lazylinop as lz
>>> import numpy as np
>>> N = 3
>>> x = np.arange(1, N + 1).astype(np.float64)
>>> x
array([1., 2., 3.])
>>> L = lz.basicops.padder(N, (N, N - 1))
>>> L @ x
array([0., 0., 0., 1., 2., 3., 0., 0.])
>>> L = lz.basicops.padder(N, (N, N), mode='periodic')
>>> L @ x
array([1., 2., 3., 1., 2., 3., 1., 2., 3.])
>>> L = lz.basicops.padder(N, (N, N), mode='symmetric')
>>> L @ x
array([3., 2., 1., 1., 2., 3., 3., 2., 1.])
>>> L = lz.basicops.padder(N, (1, N), mode='periodic')
>>> L @ x
array([3., 1., 2., 3., 1., 2., 3.])
>>> L = lz.basicops.padder(N, (2, 1), mode='symmetric')
>>> L @ x
array([2., 1., 1., 2., 3., 3.])
>>> L = lz.basicops.padder(N, (N + 2, N + 1), mode='symmetric')
>>> L @ x
array([2., 3., 3., 2., 1., 1., 2., 3., 3., 2., 1., 1.])
>>> L = lz.basicops.padder(N, (N, N), mode='antisymmetric')
>>> L @ x
array([-3., -2., -1.,  1.,  2.,  3., -3., -2., -1.])
>>> L = lz.basicops.padder(N, (N, N), mode='reflect')
>>> L @ x
array([2., 3., 2., 1., 2., 3., 2., 1., 2.])
>>> X = np.array([[0., 3.], [1., 4.], [2., 5.]])
>>> X
array([[0., 3.],
       [1., 4.],
       [2., 5.]])
>>> L = lz.basicops.padder(N, (N + 2, N + 1), mode='symmetric')
>>> L @ X
array([[1., 4.],
       [2., 5.],
       [2., 5.],
       [1., 4.],
       [0., 3.],
       [0., 3.],
       [1., 4.],
       [2., 5.],
       [2., 5.],
       [1., 4.],
       [0., 3.],
       [0., 3.]])

See also

vstack(), hstack(), padder2d()

Agglomeration & Transformation#

Functions for creating a Lazylinop from preexisting linear operators.

lazylinop.add(*ops)#

Returns a LazyLinOp L that acts as a sum of given compatible linear operators ops.

Args:
ops:

Operators (LazyLinOp-s or other compatible linear operators) to sum.

Returns:

The LazyLinOp for the sum of ops.

Examples:
>>> import numpy as np
>>> import lazylinop as lz
>>> from lazylinop import aslazylinop
>>> nt = 10
>>> d = 8
>>> v = np.random.rand(d)
>>> terms = [np.ones((d, d)) for i in range(nt)]
>>> # terms are all Fausts here
>>> ls = lz.add(*terms) # ls is the LazyLinOp add of terms
>>> np_sum = 0
>>> for i in range(nt): np_sum += terms[i]
>>> np.allclose(ls @ v, nt * np.ones((d, d)) @ v)
True

See also

aslazylinop()

lazylinop.block_diag(*ops)#

Returns a LazyLinOp L that acts as the block-diagonal concatenation of compatible linear operators ops.

Args:
ops:

Operators (LazyLinOp-s or other compatible linear operators) to concatenate block-diagonally.

Returns:

The resulting block-diagonal LazyLinOp.

Example:
>>> import numpy as np
>>> import lazylinop as lz
>>> from lazylinop import aslazylinop
>>> import scipy
>>> nt = 10
>>> d = 64
>>> v = np.random.rand(d)
>>> terms = [np.random.rand(64, 64) for _ in range(10)]
>>> ls = lz.block_diag(*terms) # ls is the block diagonal LazyLinOp
>>> np.allclose(scipy.linalg.block_diag(*terms), ls.toarray())
True

See also

scipy.linalg.block_diag, aslazylinop()

lazylinop.hstack(ops)#

Concatenates linear operators horizontally.

Args:
ops: (tuple of compatible linear operators)

For any pair i, j < len(ops), ops[i].shape[0] == ops[i].shape[0].

Returns:

A concatenation LazyLinOp.

Example:
>>> import numpy as np
>>> import lazylinop as lz
>>> from lazylinop import islazylinop
>>> A = np.ones((10, 10))
>>> B = lz.ones((10, 2))
>>> lcat = lz.hstack((A, B))
>>> islazylinop(lcat)
True
>>> np.allclose(lcat.toarray(), np.hstack((A, B.toarray())))
True

See also

vstack(), numpy.hstack, aslazylinop()

lazylinop.kron(L1, L2)#

Returns the LazyLinOp for the Kronecker product \(L_1\otimes L_2\) using as a definition the mixed Kronecker matrix-vector property

\[\begin{equation} (L_2^T\otimes L_1)x=\mathtt{vec}(L_1XL_2) \end{equation}\]

where \(X\) is a matrix of appropriate dimension such that \(\mathtt{vec}(X)=x\). Importantly here \(\mathtt{vec}(X)\) stacks the columns of \(X\).

Note

This specialization is particularly optimized for multiplying the operator by a vector.

Args:
L1: (compatible linear operator)

scaling factor,

L2: (compatible linear operator)

block factor.

Returns:

The Kronecker product LazyLinOp.

Example:
>>> import numpy as np
>>> import lazylinop as lz
>>> A = np.random.randn(100, 100)
>>> B = np.random.randn(100, 100)
>>> AxB = np.kron(A, B)
>>> lAxB = lz.kron(A, B)
>>> x = np.random.randn(AxB.shape[1], 1)
>>> print(np.allclose(AxB@x, lAxB@x))
True
>>> from timeit import timeit
>>> timeit(lambda: AxB @ x, number=10) 
0...
>>> # example: 0.4692082800902426
>>> timeit(lambda: lAxB @ x, number=10) 
0...
>>> # example 0.03464869409799576

See also

lazylinop.pad(op, pad_width, mode='constant', constant_values=0)#

Returns a LazyLinOp L that acts as a padded version of a given compatible linear operator op.

Args:
op: (scipy LinearOperator, LazyLinOperator, numpy array, torch.Tensor)

The operator/array to pad.

pad_width: (tuple, list)

Number of values padded to the edges of each axis.

  • ((B0, A0), (B1, A1)) (See Figure Padding format).

  • (B, A) is equivalent to ((B, A), (B, A)).

  • ((B0, ), (B1, )) is equivalent to ((B0, B0), (B1, B1)).

  • (B, ) is equivalent to ((B, B), (B, B)).

  • C is equivalent to ((C, C), (C, C)).

mode: (str)
  • 'constant':

    Pads with a constant value.

  • 'symmetric':

    Pads with the reflection of the vector mirrored along the edge of the array.

  • 'antisymmetric':

    Pads with the reflection of the vector mirrored and negated along the edge of the array.

  • 'reflect':

    Pads with the reflection of the vector mirrored on the first and last values of the vector along each axis.

  • 'mean':

    Pads with the mean value of all the vector along each axis.

  • 'edge':

    Pads with the edge values of LazyLinOp.

  • 'wrap':

    Pads with the wrap of the vector along the axis. The first values are used to pad the end and the end values are used to pad the beginning.

constant_values: (tuple, list, scalar)

The padded values for each axis (in mode='constant').

  • ((VB0, VA0), (VB1, VA1)): padding values before (VBi) and values after (VAi) on each dimension. In Figure Padding format value VBi (resp. VAi) goes where padding width Bi (resp. Ai) is.

  • ((VB0, VA0)) is equivalent to ((VB0, VA0), (VB0, VA0)).

  • (V,) or V is equivalent to ((V, V), (V, V)).

  • ((VB0,), (VB1,)) is equivalent to ((VB0, VB0), (VB1, VB1)).

Padding format (for an operator op)#

_images/pad_width.svg
Example mode='constant':
>>> import lazylinop as lz
>>> import numpy as np
>>> A = np.arange(18 * 2).reshape((18, 2))
>>> A
array([[ 0,  1],
       [ 2,  3],
       [ 4,  5],
       [ 6,  7],
       [ 8,  9],
       [10, 11],
       [12, 13],
       [14, 15],
       [16, 17],
       [18, 19],
       [20, 21],
       [22, 23],
       [24, 25],
       [26, 27],
       [28, 29],
       [30, 31],
       [32, 33],
       [34, 35]])
>>> lpA = lz.pad(A, (2, 3))
>>> lpA
<23x7 LazyLinOp with unspecified dtype>
>>> lpA.toarray().astype(int)
array([[ 0,  0,  0,  0,  0,  0,  0],
       [ 0,  0,  0,  0,  0,  0,  0],
       [ 0,  0,  0,  1,  0,  0,  0],
       [ 0,  0,  2,  3,  0,  0,  0],
       [ 0,  0,  4,  5,  0,  0,  0],
       [ 0,  0,  6,  7,  0,  0,  0],
       [ 0,  0,  8,  9,  0,  0,  0],
       [ 0,  0, 10, 11,  0,  0,  0],
       [ 0,  0, 12, 13,  0,  0,  0],
       [ 0,  0, 14, 15,  0,  0,  0],
       [ 0,  0, 16, 17,  0,  0,  0],
       [ 0,  0, 18, 19,  0,  0,  0],
       [ 0,  0, 20, 21,  0,  0,  0],
       [ 0,  0, 22, 23,  0,  0,  0],
       [ 0,  0, 24, 25,  0,  0,  0],
       [ 0,  0, 26, 27,  0,  0,  0],
       [ 0,  0, 28, 29,  0,  0,  0],
       [ 0,  0, 30, 31,  0,  0,  0],
       [ 0,  0, 32, 33,  0,  0,  0],
       [ 0,  0, 34, 35,  0,  0,  0],
       [ 0,  0,  0,  0,  0,  0,  0],
       [ 0,  0,  0,  0,  0,  0,  0],
       [ 0,  0,  0,  0,  0,  0,  0]])
>>> lpA2 = lz.pad(A, ((2, 3), (4, 1)))
>>> lpA2
<23x7 LazyLinOp with unspecified dtype>
>>> lpA2.toarray().astype('int')
array([[ 0,  0,  0,  0,  0,  0,  0],
       [ 0,  0,  0,  0,  0,  0,  0],
       [ 0,  0,  0,  0,  0,  1,  0],
       [ 0,  0,  0,  0,  2,  3,  0],
       [ 0,  0,  0,  0,  4,  5,  0],
       [ 0,  0,  0,  0,  6,  7,  0],
       [ 0,  0,  0,  0,  8,  9,  0],
       [ 0,  0,  0,  0, 10, 11,  0],
       [ 0,  0,  0,  0, 12, 13,  0],
       [ 0,  0,  0,  0, 14, 15,  0],
       [ 0,  0,  0,  0, 16, 17,  0],
       [ 0,  0,  0,  0, 18, 19,  0],
       [ 0,  0,  0,  0, 20, 21,  0],
       [ 0,  0,  0,  0, 22, 23,  0],
       [ 0,  0,  0,  0, 24, 25,  0],
       [ 0,  0,  0,  0, 26, 27,  0],
       [ 0,  0,  0,  0, 28, 29,  0],
       [ 0,  0,  0,  0, 30, 31,  0],
       [ 0,  0,  0,  0, 32, 33,  0],
       [ 0,  0,  0,  0, 34, 35,  0],
       [ 0,  0,  0,  0,  0,  0,  0],
       [ 0,  0,  0,  0,  0,  0,  0],
       [ 0,  0,  0,  0,  0,  0,  0]])
>>> # the same with arbitrary values
>>> pw = ((2, 3), (4, 1))
>>> cv = ((-1, -2), (-3, -4))
>>> lpA3 = lz.pad(A, pw, constant_values=cv)
>>> lpA3
<23x7 LazyLinOp with unspecified dtype>
>>> lpA3.toarray().astype('int')
array([[-3, -3, -3, -3, -1, -1, -4],
       [-3, -3, -3, -3, -1, -1, -4],
       [-3, -3, -3, -3,  0,  1, -4],
       [-3, -3, -3, -3,  2,  3, -4],
       [-3, -3, -3, -3,  4,  5, -4],
       [-3, -3, -3, -3,  6,  7, -4],
       [-3, -3, -3, -3,  8,  9, -4],
       [-3, -3, -3, -3, 10, 11, -4],
       [-3, -3, -3, -3, 12, 13, -4],
       [-3, -3, -3, -3, 14, 15, -4],
       [-3, -3, -3, -3, 16, 17, -4],
       [-3, -3, -3, -3, 18, 19, -4],
       [-3, -3, -3, -3, 20, 21, -4],
       [-3, -3, -3, -3, 22, 23, -4],
       [-3, -3, -3, -3, 24, 25, -4],
       [-3, -3, -3, -3, 26, 27, -4],
       [-3, -3, -3, -3, 28, 29, -4],
       [-3, -3, -3, -3, 30, 31, -4],
       [-3, -3, -3, -3, 32, 33, -4],
       [-3, -3, -3, -3, 34, 35, -4],
       [-3, -3, -3, -3, -2, -2, -4],
       [-3, -3, -3, -3, -2, -2, -4],
       [-3, -3, -3, -3, -2, -2, -4]])
zero-padded DFT example:
>>> import lazylinop as lz
>>> from lazylinop.signal import fft
>>> e = lz.eye(5)
>>> pe = lz.pad(e, (0, 3))
>>> pfft = fft(8) @ pe
Example mode='symmetric', mode='reflect':
>>> import lazylinop as lz
>>> a = np.arange(25).reshape(5, 5)
>>> sp_a = lz.pad(a, (2, 1), mode='symmetric')
>>> print(sp_a)
<8x8 LazyLinOp with unspecified dtype>
>>> sp_a.toarray().astype('int')
array([[ 6,  5,  5,  6,  7,  8,  9,  9],
       [ 1,  0,  0,  1,  2,  3,  4,  4],
       [ 1,  0,  0,  1,  2,  3,  4,  4],
       [ 6,  5,  5,  6,  7,  8,  9,  9],
       [11, 10, 10, 11, 12, 13, 14, 14],
       [16, 15, 15, 16, 17, 18, 19, 19],
       [21, 20, 20, 21, 22, 23, 24, 24],
       [21, 20, 20, 21, 22, 23, 24, 24]])
>>> sp_a2 = lz.pad(a, ((1, 1), (2, 1)), mode='symmetric')
>>> print(sp_a2)
<7x8 LazyLinOp with unspecified dtype>
>>> sp_a2.toarray().astype('int')
array([[ 1,  0,  0,  1,  2,  3,  4,  4],
       [ 1,  0,  0,  1,  2,  3,  4,  4],
       [ 6,  5,  5,  6,  7,  8,  9,  9],
       [11, 10, 10, 11, 12, 13, 14, 14],
       [16, 15, 15, 16, 17, 18, 19, 19],
       [21, 20, 20, 21, 22, 23, 24, 24],
       [21, 20, 20, 21, 22, 23, 24, 24]])
>>> rp_a = lz.pad(a, (2, 1), mode='reflect')
>>> print(rp_a)
<8x8 LazyLinOp with unspecified dtype>
>>> rp_a.toarray().astype('int')
array([[12, 11, 10, 11, 12, 13, 14, 13],
       [ 7,  6,  5,  6,  7,  8,  9,  8],
       [ 2,  1,  0,  1,  2,  3,  4,  3],
       [ 7,  6,  5,  6,  7,  8,  9,  8],
       [12, 11, 10, 11, 12, 13, 14, 13],
       [17, 16, 15, 16, 17, 18, 19, 18],
       [22, 21, 20, 21, 22, 23, 24, 23],
       [17, 16, 15, 16, 17, 18, 19, 18]])
>>> rp_a2 = lz.pad(a, ((1, 1), (2, 1)), mode='reflect')
>>> print(rp_a2)
<7x8 LazyLinOp with unspecified dtype>
>>> rp_a2.toarray().astype('int')
array([[ 7,  6,  5,  6,  7,  8,  9,  8],
       [ 2,  1,  0,  1,  2,  3,  4,  3],
       [ 7,  6,  5,  6,  7,  8,  9,  8],
       [12, 11, 10, 11, 12, 13, 14, 13],
       [17, 16, 15, 16, 17, 18, 19, 18],
       [22, 21, 20, 21, 22, 23, 24, 23],
       [17, 16, 15, 16, 17, 18, 19, 18]])
Example mode='mean':
>>> import lazylinop as lz
>>> a = np.arange(25).reshape(5, 5)
>>> mp_a = lz.pad(a, (2, 1), mode='mean')
>>> print(mp_a)
<8x8 LazyLinOp with unspecified dtype>
>>> mp_a.toarray()
array([[12., 12., 10., 11., 12., 13., 14., 12.],
       [12., 12., 10., 11., 12., 13., 14., 12.],
       [ 2.,  2.,  0.,  1.,  2.,  3.,  4.,  2.],
       [ 7.,  7.,  5.,  6.,  7.,  8.,  9.,  7.],
       [12., 12., 10., 11., 12., 13., 14., 12.],
       [17., 17., 15., 16., 17., 18., 19., 17.],
       [22., 22., 20., 21., 22., 23., 24., 22.],
       [12., 12., 10., 11., 12., 13., 14., 12.]])
>>> mp_a2 = lz.pad(a, ((1, 1), (2, 1)), mode='mean')
>>> print(mp_a2)
<7x8 LazyLinOp with unspecified dtype>
>>> mp_a2.toarray()
array([[12., 12., 10., 11., 12., 13., 14., 12.],
       [ 2.,  2.,  0.,  1.,  2.,  3.,  4.,  2.],
       [ 7.,  7.,  5.,  6.,  7.,  8.,  9.,  7.],
       [12., 12., 10., 11., 12., 13., 14., 12.],
       [17., 17., 15., 16., 17., 18., 19., 17.],
       [22., 22., 20., 21., 22., 23., 24., 22.],
       [12., 12., 10., 11., 12., 13., 14., 12.]])
Example mode='edge':
>>> import lazylinop as lz
>>> a = np.arange(25).reshape(5, 5)
>>> ep_a = lz.pad(a, (2, 1), mode='edge')
>>> print(ep_a)
<8x8 LazyLinOp with unspecified dtype>
>>> y = ep_a.toarray().astype('int')
>>> y
array([[ 0,  0,  0,  1,  2,  3,  4,  4],
       [ 0,  0,  0,  1,  2,  3,  4,  4],
       [ 0,  0,  0,  1,  2,  3,  4,  4],
       [ 5,  5,  5,  6,  7,  8,  9,  9],
       [10, 10, 10, 11, 12, 13, 14, 14],
       [15, 15, 15, 16, 17, 18, 19, 19],
       [20, 20, 20, 21, 22, 23, 24, 24],
       [20, 20, 20, 21, 22, 23, 24, 24]])
>>> z = np.pad(a, (2, 1), mode='edge')
>>> np.allclose(y, z)
True
>>> ep_a2 = lz.pad(a, ((1, 1), (2, 1)), mode='edge')
>>> print(ep_a2)
<7x8 LazyLinOp with unspecified dtype>
>>> ep_a2.toarray().astype('int')
array([[ 0,  0,  0,  1,  2,  3,  4,  4],
       [ 0,  0,  0,  1,  2,  3,  4,  4],
       [ 5,  5,  5,  6,  7,  8,  9,  9],
       [10, 10, 10, 11, 12, 13, 14, 14],
       [15, 15, 15, 16, 17, 18, 19, 19],
       [20, 20, 20, 21, 22, 23, 24, 24],
       [20, 20, 20, 21, 22, 23, 24, 24]])
Example mode='wrap':
>>> import lazylinop as lz
>>> a = np.arange(25).reshape(5, 5)
>>> wp_a = lz.pad(a, (2, 1), mode='wrap')
>>> print(wp_a)
<8x8 LazyLinOp with unspecified dtype>
>>> wp_a.toarray().astype('int')
array([[18, 19, 15, 16, 17, 18, 19, 15],
       [23, 24, 20, 21, 22, 23, 24, 20],
       [ 3,  4,  0,  1,  2,  3,  4,  0],
       [ 8,  9,  5,  6,  7,  8,  9,  5],
       [13, 14, 10, 11, 12, 13, 14, 10],
       [18, 19, 15, 16, 17, 18, 19, 15],
       [23, 24, 20, 21, 22, 23, 24, 20],
       [ 3,  4,  0,  1,  2,  3,  4,  0]])
>>> wp_a2 = lz.pad(a, ((1, 1), (2, 1)), mode='wrap')
>>> print(wp_a2)
<7x8 LazyLinOp with unspecified dtype>
>>> wp_a2.toarray().astype('int')
array([[23, 24, 20, 21, 22, 23, 24, 20],
       [ 3,  4,  0,  1,  2,  3,  4,  0],
       [ 8,  9,  5,  6,  7,  8,  9,  5],
       [13, 14, 10, 11, 12, 13, 14, 10],
       [18, 19, 15, 16, 17, 18, 19, 15],
       [23, 24, 20, 21, 22, 23, 24, 20],
       [ 3,  4,  0,  1,  2,  3,  4,  0]])

See also

numpy.pad, aslazylinop()

lazylinop.vstack(ops)#

Concatenates linear operators horizontally.

Args:
ops: (tuple of compatible linear operators)

For any pair i, j < len(ops), ops[i].shape[1] == ops[i].shape[1].

Returns:

A concatenation LazyLinOp.

Example:
>>> import numpy as np
>>> import lazylinop as lz
>>> from lazylinop import islazylinop
>>> A = np.ones((10, 10))
>>> B = lz.ones((2, 10))
>>> lcat = lz.vstack((A, B))
>>> islazylinop(lcat)
True
>>> np.allclose(lcat.toarray(), np.vstack((A, B.toarray())))
True

See also

hstack() numpy.vstack aslazylinop()

Slicing and indexing operators#

lazylinop.basicops.slicer(N, start=None, stop=None, step=None)#

Returns a LazyLinOp' ``L` that extracts slices from a vector of size N, such that with:

y = slicer(N, start, stop, step) @ x

y is the vector formed by concatenation of several x slices given by each element of start/stop/step parameters, i.e:

y = [ x[start[0]:stop[0]:step[0]] … x[start[i]:stop[i]:step[i]] … ]

where start[i], stop[i], step[i] are interpreted in the same way as in usual Python slice parameters (i.e. they can be a positive or negative integer or None).

If start, stop and step are integers (or None), a single slicing is performed, i.e:

slicer(N, 0, 10, 2) @ x == x[0:10:2]

A mix of arrays and integers (or None) may be provided for start/stop/step parameters. In that case, each integer parameter is expanded into an array filled with its value of the same length as the arrays provided in other parameters. (Note: all provided arrays must have the same length)

Args:
N: int

Length of the input.

start: int or np.ndarray

The slice’s first element or slices’ list of first elements. Each start element is interpreted in the same way as Python slice()’s start parameter. Default is None.

stop: int or np.ndarray

The slice’s stop index or slices’ list of stop index. Each stop element is interpreted in the same way as Python slice()’s stop parameter. Default is None.

step: int or np.ndarray

The slice’s stride or slices’s list of strides to be used. Each stop element is interpreted in the same way as Python slice()’s step parameter. Default is None.

Returns:

The slices LazyLinOp.

Examples:
>>> import numpy as np
>>> import lazylinop as lz
>>> N = 10
>>> x = np.arange(N)
>>> x
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
>>> # use integers to extract a single slice
>>> lz.basicops.slicer(N, 5, 8) @ x
array([5, 6, 7])
>>> # use list to extract multiple slices
>>> lz.basicops.slicer(N, [0, 5], [2, 8]) @ x
array([0, 1, 5, 6, 7])
>>> # pick one every two elements
>>> lz.basicops.slicer(N, step=2) @ x
array([0, 2, 4, 6, 8])
>>> # same example using a batch, in reverse order
>>> X = np.arange(3*N).reshape((N, 3))
>>> X
array([[ 0,  1,  2],
       [ 3,  4,  5],
       [ 6,  7,  8],
       [ 9, 10, 11],
       [12, 13, 14],
       [15, 16, 17],
       [18, 19, 20],
       [21, 22, 23],
       [24, 25, 26],
       [27, 28, 29]])
>>> lz.basicops.slicer(N, step=-2) @ X
array([[27, 28, 29],
       [21, 22, 23],
       [15, 16, 17],
       [ 9, 10, 11],
       [ 3,  4,  5]])
>>> # windows of 5 elements every 2 elements across the vector
>>> lz.basicops.slicer(N, np.arange(0, 5, 2), np.arange(5, N, 2)) @ x
array([0, 1, 2, 3, 4, 2, 3, 4, 5, 6, 4, 5, 6, 7, 8])
lazylinop.basicops.indexer(N, index)#

Returns a LazyLinOp' ``L` that extracts elements from a vector of size N, such that with I an array of integers or booleans:

indexer(N, I) @ x == x[I]

Args:
N: int

Length of the input.

index: An array-compatible of integers or booleans.

If index is an array of integers, it reprensents the indices of the element of the input vector to select.

If index in an array of booleans, indexer() will returns the input vector elements whose indices are True in index. In that case index size must be N or 1.

Returns:

The indexer LazyLinOp.

Examples:
>>> import numpy as np
>>> import lazylinop as lz
>>> N = 10
>>> x = np.arange(N)
>>> x
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
>>> # use integers to extract elements
>>> lz.indexer(N, [2, 1, 1]) @ x
array([2, 1, 1])
>>> # use booleans to extract elements
>>> lz.indexer(N, x > 6) @ x
array([7, 8, 9])

Permutations operators#

The following functions return LazyLinOp’s that implement permutations. These operators are all of shape \(\left(N,~N\right)\) and orthonormal.

lazylinop.basicops.bitrev(N)#

Returns a LazyLinOp L for the bit-reversal permutation. Size of the signal N must be a power of two. Bit-reversal permutation maps each item of the sequence 0 to N - 1 to the item whose bit representation has the same bits in reversed order.

Args:
N: int

Size of the signal.

Returns:

LazyLinOp

Example:
>>> import numpy as np
>>> from lazylinop.basicops import bitrev
>>> x = np.arange(4)
>>> L = bitrev(4)
>>> L @ x
array([0, 2, 1, 3])
References:
[1] Fast Bit-Reversal Algorithms, Anne Cathrine Elster.

IEEE International Conf. on Acoustics, Speech, and Signal Processing 1989 (ICASSP’89), Vol. 2, pp. 1099-1102, May 1989.

lazylinop.basicops.flip(N, start=0, end=None)#

Returns a LazyLinOp L that flips an input array or a sub-interval of the array.

For an input x with entries \(x_0, x_1, \ldots, x_{N - 1}\), the result of y = L @ x is:

  • with defaults parameters, the entries of y are:

    \(x_{N - 1}, x_{N - 2}, \ldots, x_0\)

  • with start = a and end = b (a < b), the entries of y are:

    \(x_0, x_1, \ldots, x_{b - 1}, \ldots, x_a, x_b, \ldots, x_{N - 1}\)

Shape of L is \((N,~N)\).

Args:
N: int

Size of the input.

start: int, optional

Start to flip from this value (default is 0).

end: int, optional

Stop to flip (not included, default is None).

Returns:

LazyLinOp

Examples:
>>> import numpy as np
>>> from lazylinop.basicops import flip
>>> N = 6
>>> x = np.arange(N)
>>> x
array([0, 1, 2, 3, 4, 5])
>>> y = flip(N, 0, 5) @ x
>>> y
array([4, 3, 2, 1, 0, 5])
>>> z = flip(N, 2, 4) @ x
>>> z
array([0, 1, 3, 2, 4, 5])
>>> X = np.eye(6, 5)
>>> X
array([[1., 0., 0., 0., 0.],
       [0., 1., 0., 0., 0.],
       [0., 0., 1., 0., 0.],
       [0., 0., 0., 1., 0.],
       [0., 0., 0., 0., 1.],
       [0., 0., 0., 0., 0.]])
>>> flip(N, 1, 4) @ X
array([[1., 0., 0., 0., 0.],
       [0., 0., 0., 1., 0.],
       [0., 0., 1., 0., 0.],
       [0., 1., 0., 0., 0.],
       [0., 0., 0., 0., 1.],
       [0., 0., 0., 0., 0.]])
lazylinop.basicops.roll(N, shift=0)#

Returns a LazyLinOp for rolling the elements of a vector x.

The elements that roll beyond the first position (resp. the last) re-enter at the last (resp. the first). Rolling of \(x=\left(x_0,x_1,\cdots,x_{N-1}\right)\) with shift=s is: \(x_{s}=\left(x_s,\cdots,x_{N-1},x_0,\cdots,x_{s-1}\right)\).

Args:
N: int

Size of the input.

shift: int, optional

Shift the elements by this number (to the left or to the right).

  • If negative, shift to the left.

  • If positive, shift to the right.

  • If zero (default), do nothing.

Returns:

LazyLinOp

See also

NumPy roll function.

Examples:
>>> import numpy as np
>>> import lazylinop as lz
>>> x = np.arange(4)
>>> L = roll(4, 2)
>>> y = L @ x
>>> np.allclose(np.array([2, 3, 0, 1]), y)
True
>>> L = roll(4, -1)
>>> y = L @ x
>>> np.allclose(np.array([1, 2, 3, 0]), y)
True