Linear algebra module

Module for linear algebra related LazyLinOps.

lazylinop.wip.linalg.coshm(L, scale=1.0, nmax=8, backend='scipy', use_numba=False)

Constructs an hyperbolic cosine of linear operator L as a lazy linear operator C(L). It uses the equation 2 * cosh(z) = exp(scale * L) + exp(-scale * L). Of note, it is only an approximation (see nmax argument).

Args:

L: 2d array Linear operator. scale: float, optional Scale factor cosh(scale * L) (default is 1). nmax: int, optional Stop the serie expansion after nmax (default is 8). backend: str, optional It can be ‘scipy’ (default) to use scipy.linalg.coshm function. nmax parameter is useless if backend is ‘scipy’. It can be ‘serie’ to use a serie expansion of coshm(scale * L).

Returns:

LazyLinOp

Raises:
ValueError

L @ x does not work because # of columns of L is not equal to the # of rows of x.

ValueError

backend value is either ‘scipy’ or ‘serie’.

ValueError

nmax must be >= 1.

Examples:
>>> import numpy as np
>>> import scipy as sp
>>> from lazylinop import eye
>>> from lazylinop.wip.linear_algebra import coshm
>>> scale = 0.01
>>> N = 10
>>> L = np.eye(N, M=N, k=0)
>>> E1 = coshm(L, scale=scale, nmax=32, backend='serie')
>>> E2 = sp.linalg.coshm(scale * L)
>>> E3 = coshm(eye(N, n=N, k=0), scale=scale, nmax=32, backend='serie')
>>> X = np.random.rand(N, 2 * N)
>>> np.allclose(E1 @ X, E2 @ X)
True
>>> np.allclose(E2 @ X, E3 @ X)
True
References:

See also scipy.linalg.coshm function.

lazylinop.wip.linalg.cosm(L, scale=1.0, nmax=8, backend='scipy')

Constructs a cosinus of linear operator L as a lazy linear operator C(L). It uses the equation expm(i * scale * L) = cos(scale * L) + i * sin(scale * L) where i^2 = -1 and returns real part. Of note, it is only an approximation (see nmax argument).

Args:

L: 2d array Linear operator. scale: float, optional Scale factor cosm(scale * L) (default is 1). nmax: int, optional Stop the serie expansion after nmax (default is 8). backend: str, optional It can be ‘scipy’ (default) to use scipy.linalg.cosm function. nmax parameter is useless if backend is ‘scipy’. It can be ‘serie’ to use a serie expansion of cosm(scale * L).

Returns:

LazyLinOp

Raises:
ValueError

L @ x does not work because # of columns of L is not equal to the # of rows of x.

ValueError

backend value is either ‘scipy’ or ‘serie’.

ValueError

If L is a 2d array, backend must be ‘scipy’.

Examples:
>>> import numpy as np
>>> import scipy as sp
>>> from lazylinop import eye
>>> from lazylinop.wip.linear_algebra import cosm
>>> scale = 0.01
>>> coefficients = np.array([1.0, scale, 0.5 * scale ** 2])
>>> N = 10
>>> L = np.eye(N, n=N, k=0)
>>> E1 = cosm(L, scale=scale, nmax=4)
>>> E2 = sp.linalg.cosm(scale * L)
>>> np.allclose(E1.toarray(), E2)
>>> E3 = eye(N, n=N, k=0)
>>> X = np.random.rand(N, 2 * N)
>>> np.allclose(E2 @ X, E3 @ X)
References:

See also scipy.linalg.cosm function. See also expm().

lazylinop.wip.linalg.expm(L, scale=1.0, nmax=8, backend='scipy', use_numba=False)

Constructs exponentiation of linear operator L as a lazy linear operator E(L). Of note, it is only an approximation E(L) @ X ~= sum((scale * L)^i / factorial(i), i=0 to nmax) @ X.

Args:
L: 2d array

Linear operator.

scale: float, optional

Scale factor expm(scale * L) (default is 1).

nmax: int, optional

Stop the serie expansion after nmax (default is 8).

backend: str, optional

It can be ‘scipy’ (default) to use scipy.linalg.expm function. nmax parameter is useless if backend is ‘scipy’. It can be ‘serie’ to use a serie expansion of expm(scale * L).

use_numba: bool, optional

If True, use prange from Numba (default is False) to compute batch in parallel. It is useful only and only if the batch size and the length of a are big enough.

Returns:

LazyLinOp

Raises:
ValueError

L @ x does not work because # of columns of L is not equal to the # of rows of x.

ValueError

backend value is either ‘scipy’ or ‘serie’.

ValueError

If L is a 2d array, backend must be ‘scipy’.

Examples:
>>> import numpy as np
>>> import scipy as sp
>>> from lazylinop import eye
>>> from lazylinop.wip.linear_algebra import expm
>>> scale = 0.01
>>> coefficients = np.array([1.0, scale, 0.5 * scale ** 2])
>>> N = 10
>>> L = np.eye(N, n=N, k=0)
>>> E1 = expm(L, scale=scale, nmax=4)
>>> E2 = sp.linalg.expm(scale * L)
>>> np.allclose(E1.toarray(), E2)
>>> E3 = eye(N, n=N, k=0)
>>> X = np.random.rand(N, 2 * N)
>>> np.allclose(E2 @ X, E3 @ X)
References:

See also scipy.linalg.expm function.

lazylinop.wip.linalg.logm(L, scale=1.0, nmax=8, backend='scipy', use_numba=False)

Constructs logarithm of linear operator L as a lazy linear operator Log(L). Of note, it is only an approximation Log(L) @ X ~= sum((-1)^(n + 1) * (L - Id)^n / n, n=1 to nmax) @ X.

Args:

L: 2d array Linear operator. scale: float, optional Scale factor logm(scale * L) (default is 1). nmax: int, optional Stop the serie expansion after nmax (default is 8). backend: str, optional It can be ‘scipy’ (default) to use scipy.linalg.logm function. nmax parameter is useless if backend is ‘scipy’. It can be ‘serie’ to use a serie expansion of logm(scale * L).

Returns:

LazyLinOp

Raises:
ValueError

L @ x does not work because # of columns of L is not equal to the # of rows of x.

ValueError

backend value is either ‘scipy’ or ‘serie’.

ValueError

nmax must be >= 1.

Examples:
>>> import numpy as np
>>> import scipy as sp
>>> from lazylinop import eye
>>> from lazylinop.wip.linear_algebra import logm
>>> scale = 0.01
>>> N = 10
>>> E1 = logm(eye(N, n=N, k=0), scale=scale, nmax=4, backend='scipy')
>>> E2 = sp.linalg.logm(scale * np.eye(N, M=N, k=0))
>>> X = np.random.rand(N, 2 * N)
>>> np.allclose(E1 @ X, E2 @ X)
References:

See also scipy.linalg.logm function. See also logarithm of a matrix.

lazylinop.wip.linalg.sinhm(L, scale=1.0, nmax=8, backend='scipy', use_numba=False)

Constructs an hyperbolic cosine of linear operator L as a lazy linear operator S(L). It uses the equation 2 * sinh(z) = exp(scale * L) - exp(-scale * L). Of note, it is only an approximation (see nmax argument).

Args:

L: 2d array Linear operator. scale: float, optional Scale factor cosh(scale * L) (default is 1). nmax: int, optional Stop the serie expansion after nmax (default is 8). backend: str, optional It can be ‘scipy’ (default) to use scipy.linalg.sinhm function. nmax parameter is useless if backend is ‘scipy’. It can be ‘serie’ to use a serie expansion of sinhm(scale * L).

Returns:

LazyLinOp

Raises:
ValueError

L @ x does not work because # of columns of L is not equal to the # of rows of x.

ValueError

backend value is either ‘scipy’ or ‘serie’.

Examples:
>>> import numpy as np
>>> import scipy as sp
>>> from lazylinop import eye
>>> from lazylinop.wip.linear_algebra import sinhm
>>> scale = 0.01
>>> N = 10
>>> L = np.eye(N, M=N, k=0)
>>> E1 = sinhm(L, scale=scale, nmax=32, backend='serie')
>>> E2 = sp.linalg.sinhm(scale * L)
>>> E3 = sinhm(eye(N, n=N, k=0), scale=scale, nmax=32, backend='serie')
>>> X = np.random.rand(N, 2 * N)
>>> np.allclose(E1 @ X, E2 @ X)
True
>>> np.allclose(E2 @ X, E3 @ X)
True
References:

See also scipy.linalg.sinhm function.

lazylinop.wip.linalg.sinm(L, scale=1.0, nmax=8, backend='scipy')

Constructs a cosinus of linear operator L as a lazy linear operator C(L). It uses the equation expm(i * scale * L) = cos(scale * L) + i * sin(scale * L) where i^2 = -1 and returns imaginary part. Of note, it is only an approximation (see nmax argument).

Args:

L: 2d array Linear operator. scale: float, optional Scale factor sinm(scale * L) (default is 1). nmax: int, optional Stop the serie expansion after nmax (default is 8). backend: str, optional It can be ‘scipy’ (default) to use scipy.linalg.sinm function. nmax parameter is useless if backend is ‘scipy’. It can be ‘serie’ to use a serie expansion of sinm(scale * L).

Returns:

LazyLinOp

Raises:
ValueError

L @ x does not work because # of columns of L is not equal to the # of rows of x.

ValueError

backend value is either ‘scipy’ or ‘serie’.

ValueError

If L is a 2d array, backend must be ‘scipy’.

Examples:
>>> import numpy as np
>>> import scipy as sp
>>> from lazylinop import eye
>>> from lazylinop.wip.linear_algebra import sinm
>>> scale = 0.01
>>> coefficients = np.array([1.0, scale, 0.5 * scale ** 2])
>>> N = 10
>>> L = np.eye(N, n=N, k=0)
>>> E1 = sinm(L, scale=scale, nmax=4)
>>> E2 = sp.linalg.sinm(scale * L)
>>> np.allclose(E1.toarray(), E2)
>>> E3 = eye(N, n=N, k=0)
>>> X = np.random.rand(N, 2 * N)
>>> np.allclose(E2 @ X, E3 @ X)
References:

See also scipy.linalg.sinm function. See also expm().

lazylinop.wip.linalg.sqrtm(L, scale=1.0, nmax=8, backend='scipy')

Constructs square root of linear operator L as a lazy linear operator R(L). It uses the equation L^1/2 = sum((-1)^n * (1/2 n) * (Id - L)^n, n=0 to inf). Of note, it is only an approximation (see nmax argument).

Args:

L: 2d array Linear operator. scale: float, optional Scale factor cosh(scale * L) (default is 1). nmax: int, optional Stop the serie expansion after nmax (default is 8). backend: str, optional It can be ‘scipy’ (default) to use scipy.linalg.sqrtm function. nmax parameter is useless if backend is ‘scipy’. It can be ‘serie’ to use a serie expansion of sqrtm(scale * L).

Returns:

LazyLinOp

Raises:
ValueError

L @ x does not work because # of columns of L is not equal to the # of rows of x.

ValueError

backend value is either ‘scipy’ or ‘serie’.

ValueError

If L is a 2d array, backend must be ‘scipy’.

Examples:
>>> import numpy as np
>>> import scipy as sp
>>> from lazylinop import eye
>>> from lazylinop.wip.linear_algebra import sqrtm
>>> scale = 0.1
>>> N = 4
>>> L = np.eye(N, M=N, k=0)
>>> X = np.random.rand(N, 2 * N)
>>> E1 = sqrtm(L, scale=scale, nmax=256, backend='serie')
>>> E2 = sp.linalg.sqrtm(scale * L)
>>> np.allclose(E1 @ X, E2 @ X)
>>> L = eye(N, n=N, k=0)
>>> E3 = sqrtm(L, scale=scale, nmax=256, backend='serie')
>>> np.allclose(E2 @ X, E3 @ X)
References:

See also scipy.linalg.sqrtm function.

lazylinop.wip.linalg.spectral_norm(L, nits=20, thres=1e-06)

Computes the approximate spectral norm or 2-norm of operator L.

Args:
L: (LazyLinOperator)

The operator to compute the 2-norm.

nits: (int)

The number of iterations of the power iteration algorithm.

thres: (float)

The precision of the the prower iteration algorithm.

Returns:

The 2-norm of the operator L.

Example:
>>> import numpy as np
>>> from numpy.linalg import norm
>>> from numpy.random import rand, seed
>>> from lazylinop.wip.linalg import spectral_norm
>>> from lazylinop import aslazylinop
>>> from scipy.linalg.interpolative import estimate_spectral_norm
>>> seed(42) # reproducibility
>>> M = rand(10, 12)
>>> L = aslazylinop(M)
>>> ref_norm = norm(M, 2)
>>> l_norm = spectral_norm(L)
>>> np.round(ref_norm, 3)
5.34
>>> np.round(l_norm, 3)
5.34
>>> np.round(estimate_spectral_norm(L), 3)
5.34