Special matrices module

Module for special matrices related LazyLinOps.

lazylinop.wip.special_matrices.alternant(x, f, use_numba=False)

Constructs alternant matrix as a lazy linear operator A. The shape of the alternant lazy linear operator is (x.shape[0], len(f)).

\[\begin{split}A = \begin{pmatrix} f_0(x_0) & f_1(x_0) & ... & f_n(x_0)\\ f_0(x_1) & f_1(x_1) & ... & f_n(x_1)\\ f_0(x_2) & f_1(x_2) & ... & f_n(x_2)\\ . & . & & . \\ . & . & & . \\ . & . & & . \\ . & . & & . \\ f_0(x_m) & f_1(x_m) & ... & f_n(x_m)\\ \end{pmatrix}\end{split}\]

Args:

x: 1d array: Array of points
f: list: A list of lambda functions.
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: f expects at least one element.
ValueError: x expects 1d-array.

Examples:

>>> import numpy as np
>>> import scipy as sp
>>> from lazylinop.wip.special_matrices import alternant
>>> M = 5
>>> N = 6
>>> x = np.random.rand(M)
>>> f = [lambda x, n=n: np.power(x, n) for n in range(N)]
>>> X = np.random.rand(N, 3)
>>> np.allclose(alternant(x, f) @ X, np.vander(x, N=N, increasing=True) @ X)
True

References:

See also Alternant matrix.

lazylinop.wip.special_matrices.companion(a)

Constructs a companion matrix as a lazy linear operator C.

Args:

a: np.ndarray 1d array of polynomial coefficients (N, ).

Returns:

LazyLinOp

Raises:

ValueError: a expects a 1d array.
ValueError: # of coefficients must be at least >= 2.
ValueError: The first coefficient a[0] must be != 0.0.

References:

See also scipy.linalg.companion. See also Companion matrix.

lazylinop.wip.special_matrices.fiedler(a, use_numba=False)

Constructs a symmetric Fiedler matrix as a lazy linear operator F. A symmetric Fiedler matrix has entry F[i, j] = np.absolute(a[i] - a[j]).

Args:

a: np.ndarray: Sequence of numbers (shape is (n, )).
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: a is empty.

References:

See also scipy.linalg.fiedler. See also Wikipedia.

lazylinop.wip.special_matrices.helmert(n, full=False, use_numba=False)

Constructs a Helmert matrix n x n as a lazy linear operator H.

\[\begin{split}\begin{pmatrix} 1/\sqrt{n} & 1/\sqrt{n} & 1/\sqrt{n} & ... & 1/\sqrt{n}\\ 1/\sqrt{2} & -1/\sqrt{2} & 0 & ... & 0\\ 1/\sqrt{6} & 1/\sqrt{6} & 2/\sqrt{6} & ... & 0\\ . & . & . & ... & 0\\ 1/\sqrt{n(n-1)} & 1/\sqrt{n(n-1)} & 1/\sqrt{n(n-1)} & ... & -(n-1)/\sqrt{n(n-1)}\\ \end{pmatrix}\end{split}\]

Args:

n: int: The size of the matrix (n, n).
full: bool, optional: If False (default) do not return the first row H[1:, :] @ x.

Returns:

LazyLinOp

Raises:

ValueError: n must be >= 2.

Examples:

>>> from lazylinop.wip.special_matrices import helmert
>>> import numpy as np
>>> import scipy as sp
>>> N = 100
>>> X = np.random.rand(N, 10)
>>> H = helmert(N)
>>> np.allclose(H @ X, sp.linalg.helmert(N) @ X)
True

References:

See also scipy.linalg.helmert. See also R-project Helmert matrix.

lazylinop.wip.special_matrices.hilbert(n, use_numba=False, bf=False)

Constructs Hilbert matrix n x n as a lazy linear operator H (H[i, j] = 1 / (i + j + 1)). Of note Hilbert matrix is positive definite and symmetric H = L + D + L^T where L is a lower triangular matrix such that L[i, i] = 0.

Args:

n: int: The size of the matrix (n, n).
use_numba: bool, optional: If True, use Numba (default is False).
bf: bool, optional: If True, use brute force (default is False) and not NumPy.

Returns:

LazyLinOp

Raises:

ValueError: n must be >= 2.

References:

See also scipy.linalg.hilbert. See also Hilbert matrix.

lazylinop.wip.special_matrices.householder_matrix(v)

Constructs an Householder matrix lazy linear operator H from a non-zero unit column vector v.

Args:

v: 1d array: non-zero vector (unit column vector)

Returns:

LazyLinOp

Raises:

ValueError: The norm of vector v is zero.

Examples:

References:: See also Householder transformation.

lazylinop.wip.special_matrices.h_multiply(a)

Constructs a Hessenberg decomposition as a lazy linear operator H. It can be used to compute the product between Hessenberg matrix and a vector x. Hessenberg decomposition writes a = Q @ H @ Q^H.

Args:

a: np.ndarray or LazyLinOp Compute Hessenberg decomposition of the matrix a of shape (M, N).

Returns:

LazyLinOp

Raises:

ValueError: Argument a expects a 2d array.
ValueError: # of rows and # of columns are differents.

References:

See also scipy.linalg.hessenberg. See also Hessenberg decomposition.

lazylinop.wip.special_matrices.lehmer(n, use_numba=False, bf=False)

Constructs Lehmer matrix n x n as a lazy linear operator L (L[i, j] = min(i, j) / max(i, j)). Of note Lehmer matrix is symmetric.

Args:

n: int: The size of the matrix (n, n).
use_numba: bool, optional: If True, use Numba (default is False).
bf: bool, optional: If True, use brute force nested loops (default is False)

Returns:

LazyLinOp

Raises:

ValueError: n must be >= 2.

Examples:

>>> import numpy as np
>>> from lazylinop.wip.special_matrices import lehmer
>>> N = 2
>>> x = np.random.rand(N)
>>> np.allclose(lehmer(N) @ X, np.array([[1, 1 / 2], [1 / 2, 1]]) @ X)
True

References:

See also Lehmer matrix.

lazylinop.wip.special_matrices.leslie(f, s, use_numba=False)

Constructs a Leslie matrix as a lazy linear operator L.

Args:

f: np.ndarray: The fecundity coefficients (N, ).
s: np.ndarray: The survival coefficients (N - 1, ).

Returns:

LazyLinOp

Raises:

ValueError: # of fecundity coefficients must be N and # of survival coefficients must be N - 1.

References:

See also scipy.linalg.leslie. See also Leslie matrix.

lazylinop.wip.special_matrices.pascal(n, kind='symmetric', exact=True)

Constructs Pascal matrix as a lazy linear operator P. It uses the formula S = exp(A) @ exp(B) where B = A^T is a matrix with entries only on the first subdiagonal. The entries are the sequence arange(1, n) (NumPy notation). Of note, A and B are nilpotent matrices A^n = B^n = 0. To compute S @ X we use the Taylor expansion S @ X = sum(A^k / k!, k=0 to n) @ sum(B^k / k!, k=0 to n) @ X. Because of A and B are nilpotent matrices, we just have to compute the first n terms of the expansion.

Args:

n: int: The size of the Pascal matrix (n, n).
kind: str, optional: If ‘lower’ constructs lower Pascal matrix L. If ‘upper’ constructs upper Pascal matrix U. If ‘symmetric’ (default) constructs L @ U.
exact: bool, optional: If exact is False the matrix coefficients are not the exact ones. If exact is True (default) the matrix coefficients will be integers.

Returns:

LazyLinOp

Raises:

ValueError: kind is either ‘symmetric’, ‘lower’ or ‘upper’.

Examples:

>>> from lazylinop.wip.special_matrices import pascal
>>> import numpy as np
>>> import scipy as sp
>>> N = 100
>>> X = np.random.rand(N, 10)
>>> P = pascal(N, kind='symmetric', exact=True)
>>> np.allclose(P @ X, sp.linalg.pascal(N, kind='symmetric', exact=True) @ X)
True

References:

See also scipy.linalg.pascal. See also Pascal matrix.

lazylinop.wip.special_matrices.redheffer(n, use_numba=False)

Constructs Redheffer matrix n x n as a lazy linear operator R. Redheffer matrix entry R[i, j] is 1 if i divides j, 0 otherwize. Redheffer matrix for n=5 looks like:

\[\begin{split}R = \begin{pmatrix} 1 & 1 & 1 & 1 & 1\\ 1 & 1 & 0 & 1 & 0\\ 1 & 0 & 1 & 0 & 0\\ 1 & 0 & 0 & 1 & 0\\ 1 & 0 & 0 & 0 & 1\\ \end{pmatrix}\end{split}\]

and its transpose looks like:

\[\begin{split}R^T = \begin{pmatrix} 1 & 1 & 1 & 1 & 1\\ 1 & 1 & 0 & 0 & 0\\ 1 & 0 & 1 & 0 & 0\\ 1 & 1 & 0 & 1 & 0\\ 1 & 0 & 0 & 0 & 1\\ \end{pmatrix}\end{split}\]

Args:

n: int: The size of the matrix (n, n).
use_numba: bool, optional: If True, use Numba (default is False).

Returns:

LazyLinOp

Raises:

ValueError: n must be >= 2.

Examples:

>>> import numpy as np
>>> from lazylinop.wip.special_matrices import redheffer
>>> N = 3
>>> x = np.random.rand(N)
>>> np.allclose(redheffer(N) @ x, np.array([[1, 1, 1], [1, 1, 0], [1, 0, 1]]) @ x)
True

References:

See also Redheffer matrix.

lazylinop.wip.special_matrices.sylvester(cp, cq)

Constructs Sylvester matrix as a lazy linear operator S_p,q. If p has a degree m=2 and q has a degree n=3 Sylvester matrix looks like:

\[\begin{split}S = \begin{pmatrix} p_2 & p_1 & p_0 & 0 & 0\\ 0 & p_2 & p_1 & p_0 & 0\\ 0 & 0 & p_2 & p_1 & p_0\\ q_3 & q_2 & q_1 & q_0 & 0\\ 0 & q_3 & q_2 & q_1 & q_0\\ \end{pmatrix}\end{split}\]

Args:

cp: list: List of coefficients (m + 1) of the first polynomial p.
cq: list: List of coefficients (n + 1) of the second polynomial q.

Returns:

LazyLinOp

Raises:

ValueError: List of coefficients should have at least one element.

Examples:

>>> import numpy as np
>>> from lazylinop.wip.special_matrices import sylvester
>>> Op = sylvester(np.random.rand(3), np.random.rand(2))
>>> check_op(Op)
True

References:

See Sylvester matrix.

lazylinop.wip.special_matrices.tri(a, factor='lower', use_numba=False)

Constructs a lower (resp. upper) triangular matrix L (resp. U) from matrix a as a lazy linear operator Op.

Args:: a: np.ndarray or LazyLinOp Shape of the input is (M, M). factor: str, optional If ‘lower’ (default) consider a as lower triangular matrix. If ‘upper’ consider a as upper triangular matrix. use_numba: str, 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 size of a are big enough.
Returns:: LazyLinOp

Raises:

Examples:

References:

lazylinop.wip.special_matrices.vander(x, N=None)

Constructs Vandermonde matrix as a lazy linear operator V. The shape of the Vandermonde lazy linear operator is (x.shape[0], deg + 1).

\[\begin{split}V = \begin{pmatrix} 1 & x[0] & x[0]^2 & ... & x[0]^{deg}\\ 1 & x[1] & x[1]^2 & ... & x[1]^{deg}\\ 1 & x[2] & x[2]^2 & ... & x[2]^{deg}\\ . & . & . & ... & . \\ . & . & . & ... & . \\ . & . & . & ... & . \\ 1 & x[n] & x[n]^2 & ... & x[n]^{deg} \end{pmatrix}\end{split}\]

Args:

x: 1d array Array of points N: int, optional Number of columns in the output. Maximum degree is N - 1.

Returns:

LazyLinOp

Raises:

ValueError: x must be a 1d-array.
ValueError: N expects an integer value >= 1.

Examples:

References:: See also NumPy polyvander. See also NumPy vander function. See also Vandermonde matrix.