Butterfly module

Butterfly module#

Preamble#

This module gathers utilities to manipulate and create a particular family of LazyLinOp associated to matrices with the so-called butterfly structure, i.e., which are expressed as a product of a few compatible Kronecker-Sparse factors. Classical examples of such operators are the Hadamard transform and the Discrete-Fourier-Transform.

Construction of a butterfly LazyLinOp from a given array.

The main function of the module is ksd() (Kronecker-sparse decomposition), which computes a near best approximation of a given array or of a LazyLinOp by an operator with a prescribed butterfly structure [1].

Example: Consider an \(N\times N\) Hadamard matrix (with \(N\) a power of two):

>>> import scipy as sp
>>> H = sp.linalg.hadamard(N)

Turning this NumPy array into an efficient butterfly implementation requires a simple call L = ksd(H, chain) where chain is a Chain object specifying the nature of the Kronecker-sparse factors involved in the decomposition. The classical chain associated to the Hadamard matrix is square-dyadic, and can be retrieved as follows:

>>> chain = Chain.square_dyadic(H.shape)
>>> L = ksd(H, chain)

It is also possible to use other chains (see below), leading to implementations with different tradeoffs between memory consumption, energy consumption and computational speed.

>>> chain2 = Chain.monarch(H.shape, ...)
>>> L2 = ksd(H, chain2)

Construction of individual or multiple Kronecker-sparse factors.

Under the hood, the operator L resulting from a call to ksd() is a (lazy) product of Kronecker-sparse factors, each of them benefitting from an optimized implementation on various backends via the function ksm() (Kronecker-sparse multiplication operator). A single Kronecker-sparse factor is a structure determined by a 4D array \(T\) of shape (a, b, c, d), and is fully determined as L = ksm(T) given a 4-dimensional array of values \(T_{i,j,k,l}\) of size \(a\times b\times c\times d\). The shape of L is then \(abd\times acd\), and converting it to an array would yield a structured sparse matrix with support (the set of possible nonzero entries) \(I_a\otimes 1_{b\times c}\otimes I_d\). See [1] and the documentation of ksm() for further details and illustrations.

Implementation parameters.

Options of ksm() (and of ksd()) notably allow to use either the CPU or the GPU with various implementation parameters, or to directly build a lazy product of Kronecker-sparse factors L = ksm([T1, ..., Tn], ...) to minimize memory movements between CPU and GPU when computing L @ x.

Saving and loading butterfly operators.

The LazyLinOp resulting from a call to ksm() or ksd() can be saved to disk (as a .json file) and load using:

>>> lazylinop.butterfly.save(L, 'filename')
>>> L = lazylinop.butterfly.load('filename')

Plot butterfly operators.

The LazyLinOp resulting from a call to ksm() or ksd() can be plotted or saved to disk (as a .png file and a .svg file) using:

>>> N = 16
>>> L = lazylinop.butterfly.dft(N)
>>> lazylinop.butterfly.plot(L)
>>> lazylinop.butterfly.plot(L, "dft_16x16")

It is also possible to save to disk as a MATLAB-style file and load using:

>>> lazylinop.butterfly.savemat(L, 'filename')
>>> L = lazylinop.butterfly.loadmat('filename')

Pre-built Discrete-Fourier-Transform.

We provide dft() a pre-built LazyLinOp L = B @ P with the Butterfly structure B multiplied by a bit-reversal permutation P corresponding to the Discrete-Fourier-Transform [2]. The L corresponding to a DFT of a signal of size \(256\) with \(8\) factors (square-dyadic decomposition) is given by:

>>> L = lazylinop.butterfly.dft(256)

Pre-built Fast-Walsh-Hadamard-Transform.

We provide fwht() a pre-built LazyLinOp L with the Butterfly structure corresponding to the Fast-Walsh-Hadamard-Transform (FWHT). The L corresponding to a FWHT of a signal of size \(256\) with \(8\) factors (square-dyadic decomposition) is given by:

>>> L = lazylinop.butterfly.fwht(256)

Chains and their manipulation.

A chain encodes the structure of the Butterfly operator L and corresponds to a sequence of tuples

\[\begin{equation} \left(\left(a_l,~b_l,~c_l,~d_l\right)\right)_{l=1}^n \end{equation}\]

such that Kronecker-sparse operators \(L_l\) with the corresponding structure have compatible sizes, i.e. such that the product \(L = L_1\cdots L_n\) is well-defined. The shape of \(L\) is given by the attribute chain.shape. The concatenation of compatible chains is implemented via chain = chain1 @ chain2. Specifying it explicitly is simple:

>>> chain = Chain([(a1, b1, c1, d1), ..., (an, bn, cn, dn)])

The most standard chain is probably the so-called square-dyadic chain, which is behind usual fast implementations of the Hadamard transform and of DFTs. It is only defined when \(\mathtt{shape}=\left(N,~N\right)\) with \(N\) a power of two as:

>>> sd_chain = Chain.square_dyadic(shape)

Another common chain with arbitrary non-prime shape is Monarch chain with two factors.

>>> chain = Chain.monarch(shape, ...)

You can visualize the structure of all Kronecker-sparse factors in these chains as follows:

>>> sd_chain.plot()

Output file of ``sd_chain.plot("square_dyadic")``.

>>> chain.plot()

Output file of ``chain.plot("monarch")``.

and anticipate the memory footprint of the corresponding operators:

>>> sd_chain = Chain.square_dyadic((32, 32))
>>> sd_chain.mem(np.float32)
1280
>>> chain = Chain.monarch((32, 32))
>>> chain.mem(np.float32)
1536

As you can see, the memory footprint of a monarch chain is higher than that of a square-dyadic one, but when it fits in the memory of a GPUs it can lead to a faster implementation.

Chainability.

This above chains satisfy a “chainability” property, which is crucial to ensure that ksd() can be run, with approximation guarantees [1]. Chainability can be checked using the attribute chain.chainable, and ksd() will generate an error if the provided chain is not chainable.

References#

[1] Butterfly Factorization with Error Guarantees. Leon Zheng, Quoc-Tung Le, Elisa Riccietti, and Remi Gribonval https://hal.science/hal-04763712v1/document

[2] Learning Fast Algorithms for Linear Transforms Using Butterfly Factorizations. Dao T, Gu A, Eichhorn M, Rudra A, Re C. Proc Mach Learn Res. 2019 Jun;97:1517-1527. PMID: 31777847; PMCID: PMC6879380.

lazylinop.butterfly.ksd(A, chain, ortho=True, order='l2r', svd_backend=None, **kwargs)#

Returns a LazyLinOp corresponding to the (often called “butterfly”) factorization of A into Kronecker-sparse factors with sparsity patterns determined by a chainable instance chain of Chain.

L = ksd(...) returns a LazyLinOp corresponding to the factorization of A where L = ksm(...) @ ksm(...) @ ... @ ksm(...).

Note

L.ks_values contains the ks_values of the factorization (see ksm() for more details).

As an example, the DFT matrix is factorized as follows:

M = 16
# Build DFT matrix using lazylinop.signal.fft function.
from lazylinop.signal import fft
V = fft(M).toarray()
# Use square-dyadic decomposition.
sd_chain = Chain.square dyadic(V.shape)
# Multiply the DFT matrix with the bit-reversal permutation matrix.
from lazylinop.basicops import bitrev
P = bitrev(M)
L = ksd(V @ P.T, chain=sd_chain)

Args:

A: np.ndarray, cp.array, torch.Tensor or LazyLinOp

Matrix or LazyLinOp to factorize. when A is a small LazyLinOp it is often faster to perform ksd(A.toarray(), ...) than ksd(A, ...) if the dense array A.toarray() fits in memory.

chain: Chain

Chainable instance of the Chain class. See Chain documentation for more details.

ortho: bool, optional

Whether to use orthonormalisation or not during the algorithm, see [1] for more details. Default is True.

order: str, optional

Determines in which order partial factorizations are performed, see [1] for more details.

'l2r' Left-to-right decomposition (default).
'balanced'

svd_backend: str, optional

See documentation of linalg.svds() for more details.

Kwargs:

Additional arguments ksm_backend, ksm_params to pass to ksm() function. See ksm() for more details.

Returns:

L is a LazyLinOp that corresponds to the product of chain.n_patterns LazyLinOp Kronecker-sparse factors.

The namespace and device of the ks_values of all factors are determined as follows:

If A is an array (or aslazylinop(array)) then its namespace and device are used
otherwize, svd_backend determines the namespace and device

See also

Chain,
ksm().

References:

[1] Butterfly Factorization with Error Guarantees. Léon Zheng, Quoc-Tung Le, Elisa Riccietti, and Rémi Gribonval https://hal.science/hal-04763712v1/document

Examples:

>>> import numpy as np
>>> from lazylinop.butterfly import Chain, ksd
>>> from lazylinop.basicops import bitrev
>>> from lazylinop.signal import fft
>>> N = 256
>>> M = N
>>> V = fft(M).toarray()
>>> chain = Chain.square_dyadic(V.shape)
>>> # Use bit reversal permutations matrix.
>>> P = bitrev(N)
>>> approx = (ksd(V @ P.T, chain) @ P).toarray()
>>> error = np.linalg.norm(V - approx) / np.linalg.norm(V)
>>> np.allclose(error, 0.0)
True

lazylinop.butterfly.ksm(ks_values, backend='xp', params=None)#

Returns a specialization L of LazyLinOp for Kronecker Sparse Matrix Multiplication (KSMM see [1]). The sparsity pattern (or support) of a Kronecker-Sparse factor is defined as \(I_a\otimes 1_{b,c}\otimes I_d\) while its values are given by either a 4D NumPy, CuPy array or torch tensor of shape (a, b, c, d).

The shape of L is \(\left(abd,~acd\right)\).

To fill a ks_values and its Kronecker-Sparse factor M:

M = np.zeros((a * b * d, a * c * d), dtype=np.float32)
ks_values = np.empty((a, b, c, d), dtype=M.dtype)
for i in range(a):
    for j in range(b):
        for k in range(c):
            for l in range(d):
                tmp = np.random.randn()
                ks_values[i, j, k, l] = tmp
                M[i * b * d + j * d + l,
                  i * c * d + k * d + l] = tmp

For a = 3, b = 3, c = 5, d = 5 we have the following pattern.

Note

You can access the ks_values of L = ksm(...) using L.ks_values.

With OpenCL and CUDA backend, L @ X will implicitly cast X to:

match the dtype of L.ks_values
be of contiguous type

This can incur a loss of performance, as-well-as a loss of precision if the dtype of X was initially of higher precision than that of L.ks_values.

Args:

ks_values: CuPy/NumPy arrays, torch tensors or list of arrays

It could be:

A 4D array of values of the Kronecker-Sparse factor.
List of values of the Kronecker-Sparse factors. The length of the list corresponds to the number of Kronecker-Sparse factors.

The dtype of each ks_values is either torch.bfloat16 (torch only), 'float16', 'float32', 'float64', 'complex64' or 'complex128'. See code above for details on the expected indexing of ks_values.

The chain infered by the shape of ks_values must be valid and chainable.

backend: optional

The available backends depend on the namespace (see array-api-compat for more details) of the ks_values. By default, use a namespace-based implementation.

For torch namespace:
- backend='ksmm' to run the algorithm of [1]. It uses a CUDA device determined by ks_values and relies on torch.utils.cpp_extension.load_inline.
For cupy namespace:
- backend='ksmm' to run the algorithm of [1]. It uses a CUDA device determined by ks_values (must be on GPU) and relies on cp.RawModule.
- backend='cupyx' uses the cupyx.scipy.sparse.csr_matrix function.
For numpy namespace:
- backend='scipy' uses the SciPy sparse functions scipy.sparse.block_diag and scipy.sparse.csr_matrix.
- backend=(platform, device) to use OpenCL to run the algorithm of [1].
  (None, None) uses the first platform and device.
  
  (None, 'cpu') use the first platform and CPU device.
  
  (None, 'gpu') use the first platform and GPU device.
  Please consider the following piece of code for advanced choices:
  import pyopencl as cl # Get your favorite platform. platform = cl.get_platforms()[pl_id] # To get your favorite CPU device. device = platform.get_devices(device_type=cl.device_type.CPU)[dev_id] # To get your favorite GPU device. device = platform.get_devices(device_type=cl.device_type.GPU)[dev_id]
  To check platforms and devices of your system you can also run the command line clinfo -a. See PyOpenCL documentation for more details.
- backend=pycuda.driver.Device(id) uses a CUDA device to run the algorithm of [1]. See PyCUDA documentation for more details.

params: tuple or list of tuple, optional

Argument params only works for OpenCL and CUDA backends. It could be:

A tuple params of tuples where params[0] and params[1] expect a tuple of ten elements (TILEX, TILEK, TILEY, TX, TY, VSIZE, WARP_TILEX, WARP_TILEY, WARP_TX, WARP_TY) (see [1] for more details). Do not use warp-tiling if WARP_TILEX = WARP_TX = WARP_TILEY = WARP_TY = None. If (None, None) (default), the choice of hyper-parameters for multiplication L @ X and the multiplication L.H @ X is automatic. Because we did not run a fine-tuning for all the possible \(\left(a,~b,~c,~d\right)\) and \(\left(a,~c,~b,~d\right)\) tuples, automatic does not always correspond to the best choice.
List of tuple of length the number of factors. params[i][0] and params[i][1] expect a tuple of ten elements (TILEX, TILEK, TILEY, TX, TY, VSIZE, WARP_TILEX, WARP_TILEY, WARP_TX, WARP_TY) (see [1] for more details). Do not use warp-tiling if WARP_TILEX = WARP_TX = WARP_TILEY = WARP_TY = None. See How to Optimize a CUDA Matmul Kernel for cuBLAS-like Performance: a Worklog for more details about warp-tiling. If None (default), the choice of hyper-parameters for multiplication L @ X and the multiplication L.H @ X is automatic. Because we did not run a fine-tuning for all the possible \(\left(a,~b,~c,~d\right)\) and \(\left(a,~c,~b,~d\right)\) tuples, automatic does not always correspond to the best choice.

List of assertions the tuple (TILEX, TILEK, TILEY, TX, TY, VSIZE) must satisfy:

TILEX = X * TX
TILEY = Y * TY
batch size % TILEX == 0 for performance reason. Consider zero-padding of the batch.
TILEX < batch size
TILEK <= c and c % TILEK == 0 for performance reason.
TILEX > TILEK and TILEY > TILEK
(VSIZE * X * Y) % TILEX == 0
TILEK % strideInput == 0
(VSIZE * X * Y) % TILEK == 0
TILEY % strideValues == 0
TILEY <= b
(b * d) % (d * TILEY) == 0
ks_values.dtype.itemsize * 2 * (TILEY * TILEK + TILEK * TILEX) < smem

where smem is the shared memory of the hardware used to compute, VSIZE ranges from \(1\) to \(4\), strideValues = VSIZE * X * Y / TILEK and strideInput = VSIZE * X * Y / TILEX.

Returns:

An instance L of lazylinop.butterfly.KsmLazyLinOp class that is a specialization of LazyLinOp for Kronecker Sparse Matrix Multiplication (KSMM). You can access the ks-values (resp. the chain) of L using L.ks_values (resp. L.chain).

Examples:

>>> from lazylinop.butterfly.ksm import ksm
>>> import numpy as np
>>> a, b, c, d = 2, 4, 3, 2
>>> ks_values = np.full((a, b, c, d), 1.0, dtype=np.float32)
>>> L = ksm(ks_values)
>>> L.toarray(dtype='float')
array([[1., 0., 1., 0., 1., 0., 0., 0., 0., 0., 0., 0.],
       [0., 1., 0., 1., 0., 1., 0., 0., 0., 0., 0., 0.],
       [1., 0., 1., 0., 1., 0., 0., 0., 0., 0., 0., 0.],
       [0., 1., 0., 1., 0., 1., 0., 0., 0., 0., 0., 0.],
       [1., 0., 1., 0., 1., 0., 0., 0., 0., 0., 0., 0.],
       [0., 1., 0., 1., 0., 1., 0., 0., 0., 0., 0., 0.],
       [1., 0., 1., 0., 1., 0., 0., 0., 0., 0., 0., 0.],
       [0., 1., 0., 1., 0., 1., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 1., 0., 1., 0., 1., 0.],
       [0., 0., 0., 0., 0., 0., 0., 1., 0., 1., 0., 1.],
       [0., 0., 0., 0., 0., 0., 1., 0., 1., 0., 1., 0.],
       [0., 0., 0., 0., 0., 0., 0., 1., 0., 1., 0., 1.],
       [0., 0., 0., 0., 0., 0., 1., 0., 1., 0., 1., 0.],
       [0., 0., 0., 0., 0., 0., 0., 1., 0., 1., 0., 1.],
       [0., 0., 0., 0., 0., 0., 1., 0., 1., 0., 1., 0.],
       [0., 0., 0., 0., 0., 0., 0., 1., 0., 1., 0., 1.]])
>>> # List of Kronecker-Sparse factors.
>>> a1, b1, c1, d1 = 2, 4, 3, 3
>>> ks_values1 = np.full((a1, b1, c1, d1), 1.0, dtype=np.float32)
>>> a2, b2, c2, d2 = 3, 3, 5, 2
>>> ks_values2 = np.full((a2, b2, c2, d2), 1.0, dtype=np.float32)
>>> L = ksm(ks_values1) @ ksm(ks_values2)
>>> M = ksm([ks_values1, ks_values2])
>>> np.allclose(L.toarray(dtype='float'), M.toarray(dtype='float'))
True

References:

[1] Fast inference with Kronecker-sparse matrices. Antoine Gonon and Léon Zheng and Pascal Carrivain and Quoc-Tung Le https://arxiv.org/abs/2405.15013

See also

lazylinop.butterfly.dft(N, optimize=False, **kwargs)#

Return a LazyLinOp L with the Butterfly structure corresponding to the Discrete-Fourier-Transform (DFT).

Shape of L is \(\left(N,~N\right)\) where \(N=2^n\) must be a power of two.

Args:

N: int: Size of the DFT. \(N\) must be a power of two.
optimize: bool, optional: Optimization usually results in the most optimized implementation, but the optimization process may be time consuming, depending notably on the size of the factors, on the backend (see lazylinop.butterfly.ksm() for more details), on the dtype and on the device of the ks_values. The default value is False.

Kwargs:

Additional arguments ksm_backend, ksm_params, dtype and device of the ks_values to pass to ksm() function. The default values are 'xp', None, complex128 and 'cpu'. See ksm() for more details. Infer the namespace (see array-api-compat for more details) of L.ks_values from dtype and device arguments. By default, namespace of L.ks_values is numpy and dtype is 'complex128'.

Returns:

LazyLinOp L corresponding to the DFT.

Examples:

>>> import numpy as np
>>> from lazylinop.butterfly import dft as bdft
>>> from lazylinop.signal import dft as sdft
>>> N = 2 ** 5
>>> x = np.random.randn(N)
>>> y = bdft(N) @ x
>>> z = sdft(N) @ x
>>> np.allclose(y, z)
True

lazylinop.butterfly.fwht(N, optimize=False, **kwargs)#

Return a LazyLinOp L with the Butterfly structure corresponding to the Fast-Walsh-Hadamard-Transform (FWHT).

Shape of L is \(\left(N,~N\right)\) where \(N=2^n\) must be a power of two.

L is orthogonal and symmetric and the inverse WHT operator is L.T = L.

Args:

N: int: Size of the FWHT. \(N\) must be a power of two.
optimize: bool, optional: Optimization usually results in the most optimized implementation, but the optimization process may be time consuming, depending notably on the size of the factors, on the backend (see lazylinop.butterfly.ksm() for more details), on the dtype and on the device of the ks_values. The default value is False.

Kwargs:

Additional arguments ksm_backend, ksm_params, dtype and device of the ks_values to pass to ksm() function. See ksm() for more details. Infer the namespace (see array-api-compat for more details) of L.ks_values from dtype and device arguments. By default, namespace of L.ks_values is numpy and dtype is 'float64'.

Returns:

LazyLinOp L corresponding to the FWHT. L is equivalent to hadamard(N, backend, dtype) / sqrt(N).

Examples:

>>> import numpy as np
>>> from lazylinop.butterfly import fwht as bfwht
>>> from lazylinop.signal import fwht as sfwht
>>> N = 2 ** 5
>>> x = np.random.randn(N).astype('float64')
>>> y = bfwht(N) @ x
>>> z = sfwht(N) @ x
>>> np.allclose(y, z)
True

See also

hadamard()

lazylinop.butterfly.hadamard(N, optimize=False, **kwargs)#

Return a LazyLinOp L with the Butterfly structure corresponding to the Hadamard matrix.

Shape of L is \(\left(N,~N\right)\) where \(N=2^n\) must be a power of two.

L is not orthogonal and its inverse is L.T / N. L is symmetric.

The number of factors \(n\) of the square-dyadic decomposition is given by \(n=\log_2\left(N\right)\).

Infer the namespace (see array-api-compat for more details) of L.ks_values from dtype and device arguments. By default, namespace of L.ks_values is numpy and dtype is 'float64'.

Args:

N: int: Size of the Hadamard matrix. \(N\) must be a power of two.
optimize: bool, optional: Optimization usually results in the most optimized implementation, but the optimization process may be time consuming, depending notably on the size of the factors, on the backend (see lazylinop.butterfly.ksm() for more details), on the dtype and on the device of the ks_values. The default value is False.

Kwargs:

Additional arguments ksm_backend, ksm_params, dtype and device of the ks_values to pass to ksm() function. The default values are 'xp', None, float64 and 'cpu'. See ksm() for more details.

Returns:

LazyLinOp L corresponding to the Hadamard matrix.

Examples:

>>> import numpy as np
>>> import scipy as sp
>>> from lazylinop.butterfly import hadamard
>>> N = 2 ** 5
>>> x = np.random.randn(N).astype('float64')
>>> y = hadamard(N) @ x
>>> z = sp.linalg.hadamard(N) @ x
>>> np.allclose(y, z)
True

See also

Wikipedia,
fwht().

lazylinop.butterfly.suksd(A, chain, rtol=1e-06, max_iter=100, verbose=True)#

To approximate a 2d array A by a sum of Kronecker-Sparse decomposition.

Args:

A: 2d array: Matrix or LazyLinOp to factorize. when A is a small LazyLinOp it is often faster to perform suksd(A.toarray(), ...) than suksd(A, ...) if the dense array A.toarray() fits in memory.
chain: Chain: Chainable instance of the Chain class. See lazylinop.butterfly.Chain documentation for more details.
rtol: float, optional: Stop if the residual Frobenius norm of A - approx divided by the Frobenius norm of A is smaller than rtol.
max_iter: int, optional: Stop if the number of iterations is greater than max_iter.
verbose: bool, optional: Print number of iterations and relative error.

Returns:

A tuple (L, residual, rerr) where L is a LazyLinOp resulting from the approximation of A, residual is an array given by A - L.toarray(...) and rerr is the list of relative errors. The array namespace of residual is equal to:

the array namespace of A if A is an array.
NumPy if L is a LazyLinOp.

The length of rerr is the number of iterations of the algorithm.

Examples:

>>> from lazylinop.butterfly.suksd import suksd
>>> from lazylinop.butterfly import Chain
>>> import numpy as np
>>> N = 16
>>> A = np.random.randn(N, N)
>>> chain = Chain.monarch((N, N))
>>> L, res, rerr = suksd(A, chain, 1e-6, 100, False)
>>> bool(rerr[-1] < 1e-6 or len(rerr) == 100)
True

See also

lazylinop.butterfly.Chain(),
lazylinop.butterfly.ksd().

lazylinop.butterfly.fnt(N, optimize=False, **kwargs)#

Return a LazyLinOp L with the Butterfly structure corresponding to the Fast-Noiselet-Transform (FNT) (see [1] for more details).

Shape of L is \(\left(N,~N\right)\) where \(N=2^n\) must be a power of two.

L is orthonormal, and the LazyLinOp for the inverse FNT is L.H.

Args:

N: int: Size of the FNT. \(N\) must be a power of two.
optimize: bool, optional: Optimization usually results in the most optimized implementation, but the optimization process may be time consuming, depending notably on the size of the factors, on the backend (see lazylinop.butterfly.ksm() for more details), on the dtype and on the device of the ks_values. The default value is False.

Kwargs:

Returns:

LazyLinOp L corresponding to the FNT.

Examples:

>>> import numpy as np
>>> from lazylinop.butterfly import fnt as bfnt
>>> from lazylinop.signal import fnt as sfnt
>>> N = 2 ** 5
>>> x = np.random.randn(N)
>>> L = bfnt(N)
>>> y = L @ x
>>> z = sfnt(N) @ x
>>> np.allclose(y, z)
True
>>> # L is orthonormal.
>>> np.allclose(L.H @ y, x)
True

References:

[1] Noiselets. R. Coifman, F. Geshwind, Y. Meyer https://www.sciencedirect.com/science/article/pii/S1063520300903130

See also

lazylinop.butterfly.dct(N, type=2, optimize=False, real_valued=False, **kwargs)#

Returns a LazyLinOp `L for the Direct Cosine Transform (DCT).

Shape of L is \(\left(N,~N\right)\) where \(N=2^n\) must be a power of two except for DCT I (see below).

L is orthonormal, and the LazyLinOp for the inverse DCT is L.T.

Args:

N: int

Size of the input (N > 0). \(N\) must be:

a power of two for DCT II, III and IV.
a power of two plus one for DCT I.

type: int, optional

1, 2, 3, 4 (I, II, III, IV). Defaut is 2. See SciPy DCT and CuPy DCT for more details.

optimize: bool, optional

Optimization usually results in the most optimized implementation, but the optimization process may be time consuming, depending notably on the size of the factors, on the backend (see lazylinop.butterfly.ksm() for more details), on the dtype and on the device of the ks_values. The default value is False.

real_valued: bool, optional

According the a seminal paper of Makhoul the DCT matrix of size N, \(\mathbf{C}_N\), can be expressed in terms of the DFT matrix of size M, \(\mathbf{F}_{M}\) where \(M\) depends on the DCT-type. The most well-known expression is as \(\mathbf{C}_N=\mathtt{Re}(\mathbf{L}\mathbf{F}_MR)\) where the left and right matrices \(\mathbf{L}\) (of size \(N\times M\)) and \(\mathbf{R}\) (of size \(M\times N\)) are real-valued, but the DFT matrix \(\mathbf{F}_M\) is complex-valued.

This shows the existence of a “butterfly implementation” of the DCT \(\mathbf{C}_N=\mathtt{Re}(\mathbf{L}\mathbf{B}_M\mathbf{R})\) however with two caveats:

the Kronecker-sparse factors in \(\mathbf{B}_M\) are complex-valued real_valued = False
because of the \(\mathtt{Re()}\) operator, this actually expresses the DCT via a sum of two butterfly matrices \(\mathbf{C}_N=\mathbf{L}(\mathbf{B}_M+\overline{\mathbf{B}_M})\mathbf{R}\).

Since the DCT is real-valued, it would be desirable to have instead an expression in terms of a single butterfly matrix \(\mathbf{B}\) with real-valued factors real_valued = True, \(\mathbf{C}_N=\mathbf{LB'R}\) where \(\mathbf{B}':=\mathtt{real}(\mathbf{B}_M)=(\mathbf{B}_M+\overline{\mathbf{B}_M})/2\).

By the complementary low-rank characterization of butterfly matrices it is possible to show that \(\mathbf{B}'\) admits a Kronecker-sparse factorization with a different chain where each Kronecker-sparse pattern \((a,b,c,d)\) is replaced by \((a,2b,2c,d)\) except the lefmost (which is replaced by \((a,b,2c,d)\)) and the rightmost (replaced by \((a,2b,c,d)\)) ones.

The creation of the operator L = dct(N, ..., real_valued=True) can take more than one minute when \(N\ge 8192\) due to a call to the decomposition function ksd().

Kwargs:

Additional arguments ksm_backend, ksm_params, dtype and device of the ks_values to pass to ksm() function. The default values are 'xp', None, complex128 and 'cpu'. Infer the namespace (see array-api-compat for more details) of L.ks_values from dtype and device arguments. By default, namespace of L.ks_values is numpy and dtype is 'complex128'.

Returns:

LazyLinOp

Example:

>>> from lazylinop.butterfly import dct
>>> from scipy.fft import dct as sp_dct
>>> import numpy as np
>>> N = 32
>>> x = np.random.randn(N)
>>> F = dct(N, 2)
>>> y = F @ x
>>> np.allclose(y, sp_dct(x, norm='ortho'))
True
>>> # compute the inverse DCT
>>> z = F.T @ y
>>> np.allclose(z, x)
True
>>> # To mimick SciPy DCT II norm='ortho' and orthogonalize=False
>>> from lazylinop.basicops import diag
>>> v = np.full(N, 1.0)
>>> v[0] = np.sqrt(2.0)
>>> y = diag(v) @ F @ x
>>> z = sp_dct(x, 2, N, 0, 'ortho', False, 1, orthogonalize=False)
>>> np.allclose(y, z)
True

References:

[1] A Fast Cosine Transform in One and Two Dimensions,: by J. Makhoul, IEEE Transactions on acoustics, speech and signal processing vol. 28(1), pp. 27-34, :doi:`10.1109/TASSP.1980.1163351` (1980).

See also

lazylinop.butterfly.dst(N, type=2, optimize=False, real_valued=False, **kwargs)#

Returns a LazyLinOp `L for the Direct Sine Transform (DST).

Shape of L is \(\left(N,~N\right)\) where \(N=2^n\) must be a power of two except for DST I (see below).

L is orthonormal, and the LazyLinOp for the inverse DST is L.T.

Args:

N: int

Size of the input (N > 0).

\(N\) must be:

a power of two for DCT II, III and IV.
a power of two minus one for DCT I.

type: int, optional

1, 2, 3, 4 (I, II, III, IV). Defaut is 2. See SciPy DST and CuPy DST for more details.

optimize: bool, optional

real_valued: bool, optional

According the a seminal paper of Makhoul the DST matrix of size N, \(\mathbf{S}_N\), can be expressed in terms of the DFT matrix of size M, \(\mathbf{F}_{M}\) where \(M\) depends on the DST-type. The most well-known expression is as \(\mathbf{S}_N=-\mathtt{Im}(\mathbf{L}\mathbf{F}_MR)\) where the left and right matrices \(\mathbf{L}\) (of size \(N\times M\)) and \(\mathbf{R}\) (of size \(M\times N\)) are real-valued, but the DFT matrix \(\mathbf{F}_M\) is complex-valued.

This shows the existence of a “butterfly implementation” of the DST \(\mathbf{S}_N=-\mathtt{Im}(\mathbf{L}\mathbf{B}_M\mathbf{R})\) however with two caveats:

the Kronecker-sparse factors in \(\mathbf{B}_M\) are complex-valued real_valued = False
because of the \(\mathtt{Im()}\) operator, this actually expresses the DST via a sum of two butterfly matrices \(\mathbf{S}_N=\mathbf{L}(\mathbf{B}_M-\overline{\mathbf{B}_M})\mathbf{R}\).

Since the DST is real-valued, it would be desirable to have instead an expression in terms of a single butterfly matrix \(\mathbf{B}\) with real-valued factors real_valued = True, \(\mathbf{S}_N=\mathbf{LB'R}\) where \(\mathbf{B}':=\mathtt{real}(\mathbf{B}_M)=(\mathbf{B}_M-\overline{\mathbf{B}_M})/2\).

The creation of the operator L = dst(N, ..., real_valued=True) can take more than one minute when \(N\ge 8192\) due to a call to the decomposition function ksd().

Kwargs:

Additional arguments ksm_backend, ksm_params, dtype and device of the ks_values to pass to ksm() function. The default values are 'xp', None, complex128 and 'cpu'. Infer the namespace (see array-api-compat for more details) of L.ks_values from dtype and device arguments. By default, namespace of L.ks_values is numpy and dtype is 'complex128'.

Returns:

LazyLinOp

Example:

>>> from lazylinop.butterfly import dst
>>> from scipy.fft import dst as sp_dst
>>> import numpy as np
>>> N = 32
>>> x = np.random.randn(N)
>>> F = dst(N, 2)
>>> y = F @ x
>>> np.allclose(y, sp_dst(x, norm='ortho'))
True
>>> # compute the inverse DST
>>> z = F.T @ y
>>> np.allclose(z, x)
True

See also

Chain class#

class lazylinop.butterfly.Chain(ks_patterns)#

A Chain instance is built by calling the constructor Chain(ks_patterns).

Args:

ks_patterns: list: List of tuples \(((a_l,~b_l,~c_l,~d_l))_{l=1}^n\) each being called a pattern. The tuples \(((a_l,~b_l,~c_l,~d_l))_{l=1}^n\) must satisfy \(a_lc_ld_l=a_{l+1}b_{l+1}d_{l+1}\).

Attributes:

ks_patterns: list of tuple

List of patterns \((a_i,~b_i,~c_i,~d_i)\).

n_patterns: int

Equal to len(ks_patterns).

shape: tuple

shape is equal to \((a_1b_1d_1,~a_nc_nd_n)\) with n = n_patterns.

chainable: bool

True if for each \(l\):

and \(a_l\) divides \(a_{l+1}\)
and \(d_{l+1}\) divides \(d_l\)

See [1] for more details.

Return:

chain with ks_patterns, n_patterns, shape and chainable attributes.

Examples:

>>> from lazylinop.butterfly import Chain
>>> chain = Chain([(2, 1, 1, 2), (2, 1, 1, 2)])
>>> chain.ks_patterns
((2, 1, 1, 2), (2, 1, 1, 2))
>>> # Concatenation of two chains.
>>> chain1 = Chain([(1, 4, 4, 2)])
>>> chain2 = Chain([(4, 2, 2, 1)])
>>> chain = chain1 @ chain2
>>> chain.shape
(8, 8)
>>> chain.n_patterns
2
>>> chain.ks_patterns
((1, 4, 4, 2), (4, 2, 2, 1))

References:

[1] Butterfly Factorization with Error Guarantees. Lu00E9on Zheng, Quoc-Tung Le, Elisa Riccietti, and Ru00E9mi Gribonval https://hal.science/hal-04763712v1/document

lazylinop.butterfly.Chain.square_dyadic(shape)#

Build a square-dyadic chain from shape.

shape must satisfy shape[0] = shape[1] = N with \(N=2^n\) a power of two. Number of ks_patterns is equal to \(n\). The l-th pattern is given by (2 ** (l - 1), 2, 2, shape[0] // 2 ** l) where 1 <= l <= n.

We can draw the square-dyadic decomposition for \(N=16\):

using:

sq_chain = Chain.square_dyadic((16, 16))
sq_chain.plot()

Args:

shape: tuple: Shape of the input matrix must be \((N,~N)\) with \(N=2^n\).

lazylinop.butterfly.Chain.monarch(shape, p=None, q=None)#

Build a Monarch chain \(((1,~p,~q,~\frac{M}{p}),~(q,~\frac{M}{p},~\frac{N}{q},~1))\) from shape. \(p\) must divide \(M\) and \(q\) must divide \(N\). See [1] for more details.

Args:

shape: tuple: Shape of the input matrix \((M,~N)\). \(M\) and \(N\) must not be prime numbers.
p, q: int, optional: \(p\) must divide \(M\) and \(q\) must divide \(N\). If p (resp. q) is None, p (resp. q) is chosen such \(p\simeq\sqrt{m}\) (resp. \(q\simeq\sqrt{n}\)).

Examples:

>>> from lazylinop.butterfly.chain import Chain
>>> M, N = 21, 16
>>> chain = Chain.monarch((M, N))
>>> chain.ks_patterns
((1, 3, 4, 7), (4, 7, 4, 1))
>>> M, N = 12, 25
>>> chain = Chain.monarch((M, N))
>>> chain.ks_patterns
((1, 4, 5, 3), (5, 3, 5, 1))
>>> M, N = 12, 16
>>> chain = Chain.monarch((M, N))
>>> chain.ks_patterns
((1, 4, 4, 3), (4, 3, 4, 1))
>>> chain = Chain.monarch((M, N), p=2, q=4)
>>> chain.ks_patterns
((1, 2, 4, 6), (4, 6, 4, 1))

References:

[1] Monarch: Expressive structured matrices for efficient and accurate training. In International Conference on Machine Learning, pages 4690-4721. PMLR, 2022. T. Dao, B. Chen, N. S. Sohomi, A. D. Desai, M. Poli, J. Grogan, A. Liu, A. Rao, A. Rudra, and C. Ru00E9.

lazylinop.butterfly.Chain.wip_non_redundant(shape)#

This class method is still work-in-progress.

Build a non redundant chain from shape using Lemma 4.25 from [1]. The function uses a prime factorization \((q_l)_{l=1}^L\) of M = shape[0] and a prime factorization \((p_l)_{l=1}^L\) of N = shape[0] where \(L\) is the length of the prime factorization. If the two factorization do not have the same length, merge the smallest elements of the larger factorization. Raise an Exception if shape[0] > 2 and shape[1] > 2 are prime numbers.

Args:

shape: tuple: Shape of the input matrix, expect a tuple \((M,~N)\). M is equal to the number of rows \(a_1b_1d_1\) of the first factor while N is equal to the number of columns \(a_nc_nd_n\) of the last factor.

Examples:

>>> from lazylinop.butterfly import Chain
>>> sq_chain = Chain.square_dyadic((8, 8))
>>> sq_chain.ks_patterns
((1, 2, 2, 4), (2, 2, 2, 2), (4, 2, 2, 1))
>>> nr_chain = Chain.wip_non_redundant((8, 8))
>>> nr_chain.ks_patterns
((1, 2, 2, 4), (2, 2, 2, 2), (4, 2, 2, 1))

Chain methods#

lazylinop.butterfly.Chain.mem(self, dtype)#

Return the memory in bytes of the ks_values of dtype dtype to be created corresponding to self.

Args:

dtype: str: dtype of the ks_values.

Examples:

>>> import numpy as np
>>> from lazylinop.butterfly import Chain
>>> sd_chain = Chain.square_dyadic((32, 32))
>>> sd_chain.ks_patterns
((1, 2, 2, 16), (2, 2, 2, 8), (4, 2, 2, 4), (8, 2, 2, 2), (16, 2, 2, 1))
>>> sd_chain.mem(np.float32)
1280
>>> chain = Chain.monarch((32, 32))
>>> chain.ks_patterns
((1, 4, 4, 8), (4, 8, 8, 1))
>>> chain.mem(np.float32)
1536

lazylinop.butterfly.Chain.plot(self, name=None)#

Plot self.ks_patterns. The colors are randomly chosen. Matplotlib package must be installed.

Args:

name: str: Save the plot in both PNG file name + '.png' and SVG file name + '.svg'. Default value is None (it only draws the self.ks_patterns).

Examples:

>>> from lazylinop.butterfly import Chain
>>> sq_chain = Chain.square_dyadic((32, 32))
>>> sq_chain.plot("square_dyadic")
>>> chain = Chain.monarch((32, 32))
>>> chain.plot("monarch")

Chain attributes#

Chain.T#

The lazylinop.butterfly.Chain transpose.

Examples:

>>> from lazylinop.butterfly import Chain
>>> chain = Chain.monarch((64, 16))
>>> chain.ks_patterns
((1, 8, 4, 8), (4, 8, 4, 1))
>>> chain.T.ks_patterns
((4, 4, 8, 1), (1, 4, 8, 8))

Utilities#

lazylinop.butterfly.fuse()
lazylinop.butterfly.optimize()

lazylinop.butterfly.fuse(ks_values)#

Fuse a list ks_values of chainable 4D arrays (see lazylinop.butterfly.Chain for more details) into a single 4D array res. The shape of the ith ks_values[i] is \(\left(a_i,~b_i,~c_i,~d_i\right)\). The resulting 4D array of res is of shape \(\left(a_1,~\frac{b_1d_1}{d_n},~\frac{a_nc_n}{a_1},~d_n\right)\) and satisfies ksm(ks_values).toarray() == ksm(res).toarray().

This function is a utility to optimize the implementation of butterfly matrices: while being frugal in memory, butterfly matrices implemented via many (chainable) Kronecker sparse factors are sometimes less computationally efficient than their (exact) fused implementation using fewer factors.

Args:

ks_values: a list of 4D arrays: Their dimensions must correspond to a chainable chain. See lazylinop.butterfly.Chain for more details.

Returns:

A single 4D array res of shape \(\left(a_1,~\frac{b_1d_1}{d_n},~\frac{a_nc_n}{a_1},~d_n\right)\).

See also

ksm().

Examples:

>>> import numpy as np
>>> from lazylinop.butterfly import ksm, fuse
>>> a1, b1, c1, d1 = 2, 2, 2, 4
>>> a2, b2, c2, d2 = 4, 2, 2, 2
>>> a3, b3, c3, d3 = 2, 4, 4, 2
>>> v1 = np.random.randn(a1, b1, c1, d1)
>>> v2 = np.random.randn(a2, b2, c2, d2)
>>> v3 = np.random.randn(a3, b3, c3, d3)
>>> v = fuse([v1, v2, v3])
>>> v.shape
(2, 4, 4, 2)
>>> L = ksm(v)
>>> L1 = ksm(v1)
>>> L2 = ksm(v2)
>>> L3 = ksm(v3)
>>> x = np.random.randn(L.shape[1])
>>> np.allclose(L @ x, L1 @ L2 @ L3 @ x)
True

lazylinop.butterfly.optimize(input, strategy='benchmark', params={'backend': 'xp'}, verbose=True)#

Adaptively fuse:

a list input of chainable 4D arrays (see lazylinop.butterfly.Chain() for more details)
an instance of lazylinop.butterfly.KsmLazyLinOp
an instance of lazylinop.butterfly.KsmPermLazyLinOp
an instance of lazylinop.butterfly.Chain according to strategy argument and return:
a new list res of ks_values. The resulting list res of 4D arrays is of length n_factors <= len(input) and satisfies ksm(input).toarray() == ksm(res).toarray().
a new chain.

The input can be either a lazylinop.butterfly.Chain describing the structure of the Kronecker sparse factors, or an explicit list of (chainable) 4D arrays of ks_values fully specifying these factors. In the first case the output is an optimized chain, in the second one it is an optimized list of ks_values resulting from the needed hierarchical fusions, see lazylinop.butterfly.fuse().

Optimization using strategy='benchmark' usually results in the most optimized implementation, but the optimization process may be time consuming, depending notably on the size of the factors, on the backend (see lazylinop.butterfly.ksm() for more details), on the dtype and on the device of the ks_values.

Args:

input: lazylinop.butterfly.Chain or list of 4D arrays

An instance of lazylinop.butterfly.Chain or a list of chainable 4D arrays or.

strategy: str

To be chosen among the following:

'benchmark' (default): recursively fuses two consecutive elements if the resulting fused factor gives better matrix-vector multiplication performance than the succession of the two, non-fuses ones. Do not fuse if the resulting factor is dense.
'memory' iteratively fuse the two consecutive Kronecker-sparse factors that minimize the memory of the resulting fused one, until a target number of factors has been reached.
'sparsity': iteratively fuse the two consecutive Kronecker-sparse factors that minimize the ratio \(\frac{1}{ad}\) of the resulting fused one, until a target number of factors has been reached.
'speed': iteratively fuse the two consecutive Kronecker-sparse factors that minimize the ratio \(\frac{bc}{b+c}\) of the resulting fused one, until a target number of factors has been reached.
l2r: iteratively fuse factors from left to right until a target number of factors has been reached.
r2l: iteratively fuse factors from right to left until a target number of factors has been reached.
balanced: iteratively and alternatively fuse the two leftmost factors / the two rightmost ones, until a target number of factors has been reached.

Examples:

strategy='balanced':
- Case n = 6 and n_factors = 2:
  - step 0: 0 1 2 3 4 5
  - step 1: 01 2 3 45
  - step 2: 012 345
- Case n = 7 and n_factors = 2:
  - step 0: 0 1 2 3 4 5 6
  - step 1: 01 2 3 4 56
  - step 2: 012 3 456
  - step 3: 0123 456
- Case n = 7 and n_factors = 3:
  - step 0: 0 1 2 3 4 5 6
  - step 1: 01 2 3 4 56
  - step 2: 012 3 456
strategy='l2r':
- Case n = 6 and n_factors = 3:
  - step 0: 0 1 2 3 4 5
  - step 1: 01 2 3 4 5
  - step 2: 01 23 4 5
  - step 3: 01 23 45
- Case n = 7 and n_factors = 2:
  - step 0: 0 1 2 3 4 5 6
  - step 1: 01 2 3 4 5 6
  - step 2: 01 23 4 5 6
  - step 3: 01 23 45 6
  - step 4: 0123 45 6
  - step 5: 0123 456

params: dict or int

When strategy='memory', 'sparsity', 'speed' or 'balanced', 'l2r', 'r2l', params = n_factors where n_factors is the target number of factors.
When strategy='benchmark' the value of params depends on the type of input:
- If input is a chain, the four arguments params = {"backend": backend, "array_namespace": array_namespace, "dtype": dtype, "device": device} are mandatory.
- If input is a list of ks_values, params = {"backend": backend} argument is mandatory while array_namespace, dtype, device are determined using the ks_values.
Possible values of array_namespace are
- cp using import cupy as cp
- np using import numpy as np
- torch using import torch
Possible values of dtype and device depend on the array namespace. Please refer to the documentation of CuPy, NumPy and PyTorch.

verbose: bool, optional

Print fuse steps. If strategy='benchmark' writes the duration time of all possible combinations of successive fuses of the input in the output file optimize.csv. Default value is True.

Returns:

If input is a list of 4d arrays return a list res of chainable 4D arrays that satisfies ksm(input).toarray() == ksm(res).toarray().
If input is lazylinop.butterfly.Chain() return a new chain res that satisfies res.n_patterns = n_factors.

See also

lazylinop.butterfly.fuse(),
lazylinop.butterfly.ksm(),
lazylinop.butterfly.Chain.

Examples:

>>> import numpy as np
>>> from lazylinop.butterfly import ksm, optimize
>>> a1, b1, c1, d1 = 2, 2, 2, 4
>>> a2, b2, c2, d2 = 4, 2, 2, 2
>>> v1 = np.random.randn(a1, b1, c1, d1)
>>> v2 = np.random.randn(a2, b2, c2, d2)
>>> v = optimize([v1, v2], strategy='memory', params=1)
       ['0', '1']
step=0 ['01']
>>> v[0].shape
(2, 4, 4, 2)
>>> L = ksm(v)
>>> L1 = ksm(v1)
>>> L2 = ksm(v2)
>>> x = np.random.randn(L.shape[1])
>>> np.allclose(L @ x, L1 @ L2 @ x)
True
>>> # Left-to-right strategy.
>>> n = 5
>>> ksv = [np.random.randn(2, 2, 2, 2)] * n
>>> v = optimize(ksv, strategy='l2r', params=2)
       ['0', '1', '2', '3', '4']
step=0 ['01', '2', '3', '4']
step=1 ['01', '23', '4']
step=2 ['0123', '4']
>>> # Balanced strategy.
>>> n = 3
>>> ksv = [np.random.randn(2, 2, 2, 2)] * n
>>> v = optimize(ksv, strategy='balanced', params=2)
       ['0', '1', '2']
step=0 ['01', '2']
>>> n = 5
>>> ksv = [np.random.randn(2, 2, 2, 2)] * n
>>> v = optimize(ksv, strategy='balanced', params=2)
       ['0', '1', '2', '3', '4']
step=0 ['01', '2', '3', '4']
step=1 ['01', '2', '34']
step=2 ['012', '34']

I/O (load/save)#

lazylinop.butterfly.load()
lazylinop.butterfly.save()
lazylinop.butterfly.loadmat()
lazylinop.butterfly.savemat()

lazylinop.butterfly.load(name)#

Load L from file name.json. The file name.json has been created by lazylinop.butterfly.save().

Args:

name: str: Name of the .json file where to load L.

Returns:

L a lazy linear operator resulting from a call to ksm(...) or ksd(...) function. If file does not exist, return None.

See also

lazylinop.butterfly.save().

Examples:

>>> import scipy as sp
>>> import numpy as np
>>> from lazylinop.butterfly import Chain, ksd, load, save
>>> H = sp.linalg.hadamard(8)
>>> x = np.random.randn(8)
>>> chain = Chain.square_dyadic(H.shape)
>>> A = ksd(H, chain)
>>> save(A, "hadamard_8x8")
>>> A_ = load("hadamard_8x8")
>>> y = A @ x
>>> y_ = A_ @ x
>>> np.allclose(y, y_)
True

lazylinop.butterfly.save(L, name)#

Save L = ksm(...) or L = ksd(...) in a json file name + '.json'.

Args:

L:: A lazy linear operator returned by L = ksm(...) or L = ksd(...).
name: str: Name of the file.

See also

lazylinop.butterfly.load().

Examples:

>>> import scipy as sp
>>> import numpy as np
>>> from lazylinop.butterfly import Chain, ksd, load, save
>>> H = sp.linalg.hadamard(8)
>>> x = np.random.randn(8)
>>> chain = Chain.square_dyadic(H.shape)
>>> L = ksd(H, chain)
>>> save(L, "hadamard_8x8")
>>> L_ = load("hadamard_8x8")
>>> y = L @ x
>>> y_ = L_ @ x
>>> np.allclose(y, y_)
True

lazylinop.butterfly.loadmat(name)#

Load an instance L of lazylinop.butterfly.KsmLazyLinOp from file name.mat. The file name.mat has been created by savemat().

Args:

name: str: Name of the .mat file where to load L.

Returns:

An instance L of lazylinop.butterfly.KsmLazyLinOp. If file does not exist, return None.

See also

Examples:

>>> import numpy as np
>>> from lazylinop.butterfly import KsmLazyLinOp
>>> from lazylinop.butterfly import loadmat, savemat
>>> ksv = [np.random.randn(2, 2, 2, 2) for _ in range(3)]
>>> L = KsmLazyLinOp(ksv)
>>> x = np.random.randn(8)
>>> savemat(L, "savemat_L")
>>> L_load = loadmat("savemat_L")
>>> y = L @ x
>>> y_load = L_load @ x
>>> np.allclose(y, y_load)
True

lazylinop.butterfly.savemat(L, name)#

Save an instance L of lazylinop.butterfly.KsmLazyLinOp into a MATLAB file name + '.mat'.

Args:

L: lazylinop.butterfly.KsmLazyLinOp: An instance L of lazylinop.butterfly.KsmLazyLinOp to save.
name: str: Name of the file.

See also

Examples:

>>> import numpy as np
>>> from lazylinop.butterfly import KsmLazyLinOp
>>> from lazylinop.butterfly import loadmat, savemat
>>> ksv = [np.random.randn(2, 2, 2, 2) for _ in range(3)]
>>> L = KsmLazyLinOp(ksv)
>>> x = np.random.randn(8)
>>> savemat(L, "savemat_L")
>>> L_load = loadmat("savemat_L")
>>> y = L @ x
>>> y_load = L_load @ x
>>> np.allclose(y, y_load)
True

Plot#

lazylinop.butterfly.plot()

lazylinop.butterfly.plot(L, name=None, colormap='rainbow', log_scale=True)#

Plot L.ks_values on a logarithmic or linear scale. If ks_values is complex plot xp.sqrt(ks_values * xp.conj(ks_values)). Matplotlib package must be installed.

Args:

L: lazylinop.butterfly.KsmLazyLinOp

Plot L.ks_values on a logarithmic scale.

name: str

Save the plot in both PNG file name + '.png' and SVG file name + '.svg' and show the figure. Default value (None) only returns a tuple (fig, ax) for further manipulations of the plot.

colormap: str, optional

Use Matplotlib colormap. See colormaps documentation for more details.

log_scale: bool, optional

True: plot ks_values on a logarithmic scale (default value).
False: plot ks_values on a linear scale.

Returns:

If name is a str, save, show and return (None, None), otherwize return (fig, ax) corresponding to the plot of L.ks_values where fig is the matplotlib.figure.Figure instance and ax is the matplotlib.axes.Axes instance corresponding to the plot of the ks-values.

Examples:

>>> from lazylinop.butterfly import dft
>>> from lazylinop.butterfly import plot
>>> import matplotlib.image as mpl_img
>>> import matplotlib.pyplot as plt
>>> L = dft(16)
>>> # Get fig and ax of the plot.
>>> fig, ax = plot(L)
>>> axs = fig.get_axes()
>>> axs[0].set_ylim(4, 12)
(4.0, 12.0)
>>> fig.savefig("my_dft.png")
>>> fig.clf()
>>> plt.close("my_dft")
>>> # Save the plot option.
>>> fig, ax = plot(L, name="my_dft")
>>> fig is None and ax is None
True
>>> img = mpl_img.imread("my_dft.png", format='png')
>>> _ = plt.imshow(img)
>>> plt.show(block=True)
>>> plt.close()

Butterfly module

Contents

Butterfly module#

Preamble#

References#

Butterfly construction#

Chain class#

Chain methods#

Chain attributes#

Utilities#

I/O (load/save)#

Plot#