View source on GitHub |
LinearOperator acting like a circulant matrix.
Inherits From: LinearOperator, Module
tf.linalg.LinearOperatorCirculant( spectrum: tf.Tensor, input_output_dtype=tf.dtypes.complex64, is_non_singular: bool = None, is_self_adjoint: bool = None, is_positive_definite: bool = None, is_square: bool = True, name='LinearOperatorCirculant' ) This operator acts like a circulant matrix A with shape [B1,...,Bb, N, N] for some b >= 0. The first b indices index a batch member. For every batch index (i1,...,ib), A[i1,...,ib, : :] is an N x N matrix. This matrix A is not materialized, but for purposes of broadcasting this shape will be relevant.
Description in terms of circulant matrices
Circulant means the entries of A are generated by a single vector, the convolution kernel h: A_{mn} := h_{m-n mod N}. With h = [w, x, y, z],
A = |w z y x| |x w z y| |y x w z| |z y x w| This means that the result of matrix multiplication v = Au has Lth column given circular convolution between h with the Lth column of u.
Description in terms of the frequency spectrum
There is an equivalent description in terms of the [batch] spectrum H and Fourier transforms. Here we consider A.shape = [N, N] and ignore batch dimensions. Define the discrete Fourier transform (DFT) and its inverse by
DFT[ h[n] ] = H[k] := sum_{n = 0}^{N - 1} h_n e^{-i 2pi k n / N} IDFT[ H[k] ] = h[n] = N^{-1} sum_{k = 0}^{N - 1} H_k e^{i 2pi k n / N} From these definitions, we see that
H[0] = sum_{n = 0}^{N - 1} h_n H[1] = "the first positive frequency" H[N - 1] = "the first negative frequency" Loosely speaking, with * element-wise multiplication, matrix multiplication is equal to the action of a Fourier multiplier: A u = IDFT[ H * DFT[u] ]. Precisely speaking, given [N, R] matrix u, let DFT[u] be the [N, R] matrix with rth column equal to the DFT of the rth column of u. Define the IDFT similarly. Matrix multiplication may be expressed columnwise:
(A u)_r = IDFT[ H * (DFT[u])_r ]
Operator properties deduced from the spectrum.
Letting U be the kth Euclidean basis vector, and U = IDFT[u]. The above formulas show thatA U = H_k * U. We conclude that the elements of H are the eigenvalues of this operator. Therefore
- This operator is positive definite if and only if
Real{H} > 0.
A general property of Fourier transforms is the correspondence between Hermitian functions and real valued transforms.
Suppose H.shape = [B1,...,Bb, N]. We say that H is a Hermitian spectrum if, with % meaning modulus division,
H[..., n % N] = ComplexConjugate[ H[..., (-n) % N] ]
- This operator corresponds to a real matrix if and only if
His Hermitian. - This operator is self-adjoint if and only if
His real.
See e.g. "Discrete-Time Signal Processing", Oppenheim and Schafer.
Example of a self-adjoint positive definite operator
# spectrum is real ==> operator is self-adjoint # spectrum is positive ==> operator is positive definite spectrum = [6., 4, 2] operator = LinearOperatorCirculant(spectrum) # IFFT[spectrum] operator.convolution_kernel() ==> [4 + 0j, 1 + 0.58j, 1 - 0.58j] operator.to_dense() ==> [[4 + 0.0j, 1 - 0.6j, 1 + 0.6j], [1 + 0.6j, 4 + 0.0j, 1 - 0.6j], [1 - 0.6j, 1 + 0.6j, 4 + 0.0j]] Example of defining in terms of a real convolution kernel
# convolution_kernel is real ==> spectrum is Hermitian. convolution_kernel = [1., 2., 1.]] spectrum = tf.signal.fft(tf.cast(convolution_kernel, tf.complex64)) # spectrum is Hermitian ==> operator is real. # spectrum is shape [3] ==> operator is shape [3, 3] # We force the input/output type to be real, which allows this to operate # like a real matrix. operator = LinearOperatorCirculant(spectrum, input_output_dtype=tf.float32) operator.to_dense() ==> [[ 1, 1, 2], [ 2, 1, 1], [ 1, 2, 1]] Example of Hermitian spectrum
# spectrum is shape [3] ==> operator is shape [3, 3] # spectrum is Hermitian ==> operator is real. spectrum = [1, 1j, -1j] operator = LinearOperatorCirculant(spectrum) operator.to_dense() ==> [[ 0.33 + 0j, 0.91 + 0j, -0.24 + 0j], [-0.24 + 0j, 0.33 + 0j, 0.91 + 0j], [ 0.91 + 0j, -0.24 + 0j, 0.33 + 0j] Example of forcing real dtype when spectrum is Hermitian
# spectrum is shape [4] ==> operator is shape [4, 4] # spectrum is real ==> operator is self-adjoint # spectrum is Hermitian ==> operator is real # spectrum has positive real part ==> operator is positive-definite. spectrum = [6., 4, 2, 4] # Force the input dtype to be float32. # Cast the output to float32. This is fine because the operator will be # real due to Hermitian spectrum. operator = LinearOperatorCirculant(spectrum, input_output_dtype=tf.float32) operator.shape ==> [4, 4] operator.to_dense() ==> [[4, 1, 0, 1], [1, 4, 1, 0], [0, 1, 4, 1], [1, 0, 1, 4]] # convolution_kernel = tf.signal.ifft(spectrum) operator.convolution_kernel() ==> [4, 1, 0, 1] Performance
Suppose operator is a LinearOperatorCirculant of shape [N, N], and x.shape = [N, R]. Then
operator.matmul(x)isO(R*N*Log[N])operator.solve(x)isO(R*N*Log[N])operator.determinant()involves a sizeNreduce_prod.
If instead operator and x have shape [B1,...,Bb, N, N] and [B1,...,Bb, N, R], every operation increases in complexity by B1*...*Bb.
Matrix property hints
This LinearOperator is initialized with boolean flags of the form is_X, for X = non_singular, self_adjoint, positive_definite, square. These have the following meaning:
- If
is_X == True, callers should expect the operator to have the propertyX. This is a promise that should be fulfilled, but is not a runtime assert. For example, finite floating point precision may result in these promises being violated. - If
is_X == False, callers should expect the operator to not haveX. - If
is_X == None(the default), callers should have no expectation either way.
References | |
|---|---|
| Toeplitz and Circulant Matrices - A Review: Gray, 2006 (pdf) |
Methods
add_to_tensor
add_to_tensor( x, name='add_to_tensor' ) Add matrix represented by this operator to x. Equivalent to A + x.
| Args | |
|---|---|
x | Tensor with same dtype and shape broadcastable to self.shape. |
name | A name to give this Op. |
| Returns | |
|---|---|
A Tensor with broadcast shape and same dtype as self. |
adjoint
adjoint( name: str = 'adjoint' ) -> 'LinearOperator' Returns the adjoint of the current LinearOperator.
Given A representing this LinearOperator, return A*. Note that calling self.adjoint() and self.H are equivalent.
| Args | |
|---|---|
name | A name for this Op. |
| Returns | |
|---|---|
LinearOperator which represents the adjoint of this LinearOperator. |
assert_hermitian_spectrum
assert_hermitian_spectrum( name='assert_hermitian_spectrum' ) Returns an Op that asserts this operator has Hermitian spectrum.
This operator corresponds to a real-valued matrix if and only if its spectrum is Hermitian.
| Args | |
|---|---|
name | A name to give this Op. |
| Returns | |
|---|---|
An Op that asserts this operator has Hermitian spectrum. |
assert_non_singular
assert_non_singular( name='assert_non_singular' ) Returns an Op that asserts this operator is non singular.
This operator is considered non-singular if
ConditionNumber < max{100, range_dimension, domain_dimension} * eps, eps := np.finfo(self.dtype.as_numpy_dtype).eps | Args | |
|---|---|
name | A string name to prepend to created ops. |
| Returns | |
|---|---|
An Assert Op, that, when run, will raise an InvalidArgumentError if the operator is singular. |
assert_positive_definite
assert_positive_definite( name='assert_positive_definite' ) Returns an Op that asserts this operator is positive definite.
Here, positive definite means that the quadratic form x^H A x has positive real part for all nonzero x. Note that we do not require the operator to be self-adjoint to be positive definite.
| Args | |
|---|---|
name | A name to give this Op. |
| Returns | |
|---|---|
An Assert Op, that, when run, will raise an InvalidArgumentError if the operator is not positive definite. |
assert_self_adjoint
assert_self_adjoint( name='assert_self_adjoint' ) Returns an Op that asserts this operator is self-adjoint.
Here we check that this operator is exactly equal to its hermitian transpose.
| Args | |
|---|---|
name | A string name to prepend to created ops. |
| Returns | |
|---|---|
An Assert Op, that, when run, will raise an InvalidArgumentError if the operator is not self-adjoint. |
batch_shape_tensor
batch_shape_tensor( name='batch_shape_tensor' ) Shape of batch dimensions of this operator, determined at runtime.
If this operator acts like the batch matrix A with A.shape = [B1,...,Bb, M, N], then this returns a Tensor holding [B1,...,Bb].
| Args | |
|---|---|
name | A name for this Op. |
| Returns | |
|---|---|
int32 Tensor |
block_shape_tensor
block_shape_tensor() Shape of the block dimensions of self.spectrum.
cholesky
cholesky( name: str = 'cholesky' ) -> 'LinearOperator' Returns a Cholesky factor as a LinearOperator.
Given A representing this LinearOperator, if A is positive definite self-adjoint, return L, where A = L L^T, i.e. the cholesky decomposition.
| Args | |
|---|---|
name | A name for this Op. |
| Returns | |
|---|---|
LinearOperator which represents the lower triangular matrix in the Cholesky decomposition. |
| Raises | |
|---|---|
ValueError | When the LinearOperator is not hinted to be positive definite and self adjoint. |
cond
cond( name='cond' ) Returns the condition number of this linear operator.
| Args | |
|---|---|
name | A name for this Op. |
| Returns | |
|---|---|
Shape [B1,...,Bb] Tensor of same dtype as self. |
convolution_kernel
convolution_kernel( name='convolution_kernel' ) Convolution kernel corresponding to self.spectrum.
The D dimensional DFT of this kernel is the frequency domain spectrum of this operator.
| Args | |
|---|---|
name | A name to give this Op. |
| Returns | |
|---|---|
Tensor with dtype self.dtype. |
determinant
determinant( name='det' ) Determinant for every batch member.
| Args | |
|---|---|
name | A name for this Op. |
| Returns | |
|---|---|
Tensor with shape self.batch_shape and same dtype as self. |
| Raises | |
|---|---|
NotImplementedError | If self.is_square is False. |
diag_part
diag_part( name='diag_part' ) Efficiently get the [batch] diagonal part of this operator.
If this operator has shape [B1,...,Bb, M, N], this returns a Tensor diagonal, of shape [B1,...,Bb, min(M, N)], where diagonal[b1,...,bb, i] = self.to_dense()[b1,...,bb, i, i].
my_operator = LinearOperatorDiag([1., 2.]) # Efficiently get the diagonal my_operator.diag_part() ==> [1., 2.] # Equivalent, but inefficient method tf.linalg.diag_part(my_operator.to_dense()) ==> [1., 2.] | Args | |
|---|---|
name | A name for this Op. |
| Returns | |
|---|---|
diag_part | A Tensor of same dtype as self. |
domain_dimension_tensor
domain_dimension_tensor( name='domain_dimension_tensor' ) Dimension (in the sense of vector spaces) of the domain of this operator.
Determined at runtime.
If this operator acts like the batch matrix A with A.shape = [B1,...,Bb, M, N], then this returns N.
| Args | |
|---|---|
name | A name for this Op. |
| Returns | |
|---|---|
int32 Tensor |
eigvals
eigvals( name='eigvals' ) Returns the eigenvalues of this linear operator.
If the operator is marked as self-adjoint (via is_self_adjoint) this computation can be more efficient.
| Args | |
|---|---|
name | A name for this Op. |
| Returns | |
|---|---|
Shape [B1,...,Bb, N] Tensor of same dtype as self. |
inverse
inverse( name: str = 'inverse' ) -> 'LinearOperator' Returns the Inverse of this LinearOperator.
Given A representing this LinearOperator, return a LinearOperator representing A^-1.
| Args | |
|---|---|
name | A name scope to use for ops added by this method. |
| Returns | |
|---|---|
LinearOperator representing inverse of this matrix. |
| Raises | |
|---|---|
ValueError | When the LinearOperator is not hinted to be non_singular. |
log_abs_determinant
log_abs_determinant( name='log_abs_det' ) Log absolute value of determinant for every batch member.
| Args | |
|---|---|
name | A name for this Op. |
| Returns | |
|---|---|
Tensor with shape self.batch_shape and same dtype as self. |
| Raises | |
|---|---|
NotImplementedError | If self.is_square is False. |
matmul
matmul( x, adjoint=False, adjoint_arg=False, name='matmul' ) Transform [batch] matrix x with left multiplication: x --> Ax.
# Make an operator acting like batch matrix A. Assume A.shape = [..., M, N] operator = LinearOperator(...) operator.shape = [..., M, N] X = ... # shape [..., N, R], batch matrix, R > 0. Y = operator.matmul(X) Y.shape ==> [..., M, R] Y[..., :, r] = sum_j A[..., :, j] X[j, r] | Args | |
|---|---|
x | LinearOperator or Tensor with compatible shape and same dtype as self. See class docstring for definition of compatibility. |
adjoint | Python bool. If True, left multiply by the adjoint: A^H x. |
adjoint_arg | Python bool. If True, compute A x^H where x^H is the hermitian transpose (transposition and complex conjugation). |
name | A name for this Op. |
| Returns | |
|---|---|
A LinearOperator or Tensor with shape [..., M, R] and same dtype as self. |
matvec
matvec( x, adjoint=False, name='matvec' ) Transform [batch] vector x with left multiplication: x --> Ax.
# Make an operator acting like batch matrix A. Assume A.shape = [..., M, N] operator = LinearOperator(...) X = ... # shape [..., N], batch vector Y = operator.matvec(X) Y.shape ==> [..., M] Y[..., :] = sum_j A[..., :, j] X[..., j] | Args | |
|---|---|
x | Tensor with compatible shape and same dtype as self. x is treated as a [batch] vector meaning for every set of leading dimensions, the last dimension defines a vector. See class docstring for definition of compatibility. |
adjoint | Python bool. If True, left multiply by the adjoint: A^H x. |
name | A name for this Op. |
| Returns | |
|---|---|
A Tensor with shape [..., M] and same dtype as self. |
range_dimension_tensor
range_dimension_tensor( name='range_dimension_tensor' ) Dimension (in the sense of vector spaces) of the range of this operator.
Determined at runtime.
If this operator acts like the batch matrix A with A.shape = [B1,...,Bb, M, N], then this returns M.
| Args | |
|---|---|
name | A name for this Op. |
| Returns | |
|---|---|
int32 Tensor |
shape_tensor
shape_tensor( name='shape_tensor' ) Shape of this LinearOperator, determined at runtime.
If this operator acts like the batch matrix A with A.shape = [B1,...,Bb, M, N], then this returns a Tensor holding [B1,...,Bb, M, N], equivalent to tf.shape(A).
| Args | |
|---|---|
name | A name for this Op. |
| Returns | |
|---|---|
int32 Tensor |
solve
solve( rhs, adjoint=False, adjoint_arg=False, name='solve' ) Solve (exact or approx) R (batch) systems of equations: A X = rhs.
The returned Tensor will be close to an exact solution if A is well conditioned. Otherwise closeness will vary. See class docstring for details.
Examples:
# Make an operator acting like batch matrix A. Assume A.shape = [..., M, N] operator = LinearOperator(...) operator.shape = [..., M, N] # Solve R > 0 linear systems for every member of the batch. RHS = ... # shape [..., M, R] X = operator.solve(RHS) # X[..., :, r] is the solution to the r'th linear system # sum_j A[..., :, j] X[..., j, r] = RHS[..., :, r] operator.matmul(X) ==> RHS | Args | |
|---|---|
rhs | Tensor with same dtype as this operator and compatible shape. rhs is treated like a [batch] matrix meaning for every set of leading dimensions, the last two dimensions defines a matrix. See class docstring for definition of compatibility. |
adjoint | Python bool. If True, solve the system involving the adjoint of this LinearOperator: A^H X = rhs. |
adjoint_arg | Python bool. If True, solve A X = rhs^H where rhs^H is the hermitian transpose (transposition and complex conjugation). |
name | A name scope to use for ops added by this method. |
| Returns | |
|---|---|
Tensor with shape [...,N, R] and same dtype as rhs. |
| Raises | |
|---|---|
NotImplementedError | If self.is_non_singular or is_square is False. |
solvevec
solvevec( rhs, adjoint=False, name='solve' ) Solve single equation with best effort: A X = rhs.
The returned Tensor will be close to an exact solution if A is well conditioned. Otherwise closeness will vary. See class docstring for details.
Examples:
# Make an operator acting like batch matrix A. Assume A.shape = [..., M, N] operator = LinearOperator(...) operator.shape = [..., M, N] # Solve one linear system for every member of the batch. RHS = ... # shape [..., M] X = operator.solvevec(RHS) # X is the solution to the linear system # sum_j A[..., :, j] X[..., j] = RHS[..., :] operator.matvec(X) ==> RHS | Args | |
|---|---|
rhs | Tensor with same dtype as this operator. rhs is treated like a [batch] vector meaning for every set of leading dimensions, the last dimension defines a vector. See class docstring for definition of compatibility regarding batch dimensions. |
adjoint | Python bool. If True, solve the system involving the adjoint of this LinearOperator: A^H X = rhs. |
name | A name scope to use for ops added by this method. |
| Returns | |
|---|---|
Tensor with shape [...,N] and same dtype as rhs. |
| Raises | |
|---|---|
NotImplementedError | If self.is_non_singular or is_square is False. |
tensor_rank_tensor
tensor_rank_tensor( name='tensor_rank_tensor' ) Rank (in the sense of tensors) of matrix corresponding to this operator.
If this operator acts like the batch matrix A with A.shape = [B1,...,Bb, M, N], then this returns b + 2.
| Args | |
|---|---|
name | A name for this Op. |
| Returns | |
|---|---|
int32 Tensor, determined at runtime. |
to_dense
to_dense( name='to_dense' ) Return a dense (batch) matrix representing this operator.
trace
trace( name='trace' ) Trace of the linear operator, equal to sum of self.diag_part().
If the operator is square, this is also the sum of the eigenvalues.
| Args | |
|---|---|
name | A name for this Op. |
| Returns | |
|---|---|
Shape [B1,...,Bb] Tensor of same dtype as self. |
__getitem__
__getitem__( slices ) __matmul__
__matmul__( other )
View source on GitHub