1# Copyright 2018 The TensorFlow Authors. All Rights Reserved. 2# 3# Licensed under the Apache License, Version 2.0 (the "License"); 4# you may not use this file except in compliance with the License. 5# You may obtain a copy of the License at 6# 7# http://www.apache.org/licenses/LICENSE-2.0 8# 9# Unless required by applicable law or agreed to in writing, software 10# distributed under the License is distributed on an "AS IS" BASIS, 11# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 12# See the License for the specific language governing permissions and 13# limitations under the License. 14# ============================================================================== 15"""`LinearOperator` coming from a [[nested] block] circulant matrix.""" 16 17from __future__ import absolute_import 18from __future__ import division 19from __future__ import print_function 20 21import numpy as np 22 23from tensorflow.python.framework import dtypes 24from tensorflow.python.framework import ops 25from tensorflow.python.framework import tensor_shape 26from tensorflow.python.ops import array_ops 27from tensorflow.python.ops import check_ops 28from tensorflow.python.ops import math_ops 29from tensorflow.python.ops.distributions import util as distribution_util 30from tensorflow.python.ops.linalg import linalg_impl as linalg 31from tensorflow.python.ops.linalg import linear_operator 32from tensorflow.python.ops.linalg import linear_operator_util 33from tensorflow.python.ops.signal import fft_ops 34from tensorflow.python.util.tf_export import tf_export 35 36__all__ = [ 37 "LinearOperatorCirculant", 38 "LinearOperatorCirculant2D", 39 "LinearOperatorCirculant3D", 40] 41 42# Different FFT Ops will be used for different block depths. 43_FFT_OP = {1: fft_ops.fft, 2: fft_ops.fft2d, 3: fft_ops.fft3d} 44_IFFT_OP = {1: fft_ops.ifft, 2: fft_ops.ifft2d, 3: fft_ops.ifft3d} 45 46 47# TODO(langmore) Add transformations that create common spectrums, e.g. 48# starting with the convolution kernel 49# start with half a spectrum, and create a Hermitian one. 50# common filters. 51# TODO(langmore) Support rectangular Toeplitz matrices. 52class _BaseLinearOperatorCirculant(linear_operator.LinearOperator): 53 """Base class for circulant operators. Not user facing. 54 55 `LinearOperator` acting like a [batch] [[nested] block] circulant matrix. 56 """ 57 58 def __init__(self, 59 spectrum, 60 block_depth, 61 input_output_dtype=dtypes.complex64, 62 is_non_singular=None, 63 is_self_adjoint=None, 64 is_positive_definite=None, 65 is_square=True, 66 parameters=None, 67 name="LinearOperatorCirculant"): 68 r"""Initialize an `_BaseLinearOperatorCirculant`. 69 70 Args: 71 spectrum: Shape `[B1,...,Bb, N]` `Tensor`. Allowed dtypes: `float16`, 72 `float32`, `float64`, `complex64`, `complex128`. Type can be different 73 than `input_output_dtype` 74 block_depth: Python integer, either 1, 2, or 3. Will be 1 for circulant, 75 2 for block circulant, and 3 for nested block circulant. 76 input_output_dtype: `dtype` for input/output. 77 is_non_singular: Expect that this operator is non-singular. 78 is_self_adjoint: Expect that this operator is equal to its hermitian 79 transpose. If `spectrum` is real, this will always be true. 80 is_positive_definite: Expect that this operator is positive definite, 81 meaning the quadratic form `x^H A x` has positive real part for all 82 nonzero `x`. Note that we do not require the operator to be 83 self-adjoint to be positive-definite. See: 84 https://en.wikipedia.org/wiki/Positive-definite_matrix\ 85 #Extension_for_non_symmetric_matrices 86 is_square: Expect that this operator acts like square [batch] matrices. 87 parameters: Python `dict` of parameters used to instantiate this 88 `LinearOperator`. 89 name: A name to prepend to all ops created by this class. 90 91 Raises: 92 ValueError: If `block_depth` is not an allowed value. 93 TypeError: If `spectrum` is not an allowed type. 94 """ 95 96 allowed_block_depths = [1, 2, 3] 97 98 self._name = name 99 100 if block_depth not in allowed_block_depths: 101 raise ValueError("Expected block_depth to be in %s. Found: %s." % 102 (allowed_block_depths, block_depth)) 103 self._block_depth = block_depth 104 105 with ops.name_scope(name, values=[spectrum]): 106 self._spectrum = self._check_spectrum_and_return_tensor(spectrum) 107 108 # Check and auto-set hints. 109 if not self.spectrum.dtype.is_complex: 110 if is_self_adjoint is False: 111 raise ValueError( 112 "A real spectrum always corresponds to a self-adjoint operator.") 113 is_self_adjoint = True 114 115 if is_square is False: 116 raise ValueError( 117 "A [[nested] block] circulant operator is always square.") 118 is_square = True 119 120 super(_BaseLinearOperatorCirculant, self).__init__( 121 dtype=dtypes.as_dtype(input_output_dtype), 122 is_non_singular=is_non_singular, 123 is_self_adjoint=is_self_adjoint, 124 is_positive_definite=is_positive_definite, 125 is_square=is_square, 126 parameters=parameters, 127 name=name) 128 # TODO(b/143910018) Remove graph_parents in V3. 129 self._set_graph_parents([self.spectrum]) 130 131 def _check_spectrum_and_return_tensor(self, spectrum): 132 """Static check of spectrum. Then return `Tensor` version.""" 133 spectrum = linear_operator_util.convert_nonref_to_tensor(spectrum, 134 name="spectrum") 135 136 if spectrum.shape.ndims is not None: 137 if spectrum.shape.ndims < self.block_depth: 138 raise ValueError( 139 "Argument spectrum must have at least %d dimensions. Found: %s" % 140 (self.block_depth, spectrum)) 141 return spectrum 142 143 @property 144 def block_depth(self): 145 """Depth of recursively defined circulant blocks defining this `Operator`. 146 147 With `A` the dense representation of this `Operator`, 148 149 `block_depth = 1` means `A` is symmetric circulant. For example, 150 151 ``` 152 A = |w z y x| 153 |x w z y| 154 |y x w z| 155 |z y x w| 156 ``` 157 158 `block_depth = 2` means `A` is block symmetric circulant with symmetric 159 circulant blocks. For example, with `W`, `X`, `Y`, `Z` symmetric circulant, 160 161 ``` 162 A = |W Z Y X| 163 |X W Z Y| 164 |Y X W Z| 165 |Z Y X W| 166 ``` 167 168 `block_depth = 3` means `A` is block symmetric circulant with block 169 symmetric circulant blocks. 170 171 Returns: 172 Python `integer`. 173 """ 174 return self._block_depth 175 176 def block_shape_tensor(self): 177 """Shape of the block dimensions of `self.spectrum`.""" 178 # If spectrum.shape = [s0, s1, s2], and block_depth = 2, 179 # block_shape = [s1, s2] 180 return self._block_shape_tensor() 181 182 def _block_shape_tensor(self, spectrum_shape=None): 183 if self.block_shape.is_fully_defined(): 184 return linear_operator_util.shape_tensor( 185 self.block_shape.as_list(), name="block_shape") 186 spectrum_shape = ( 187 array_ops.shape(self.spectrum) 188 if spectrum_shape is None else spectrum_shape) 189 return spectrum_shape[-self.block_depth:] 190 191 @property 192 def block_shape(self): 193 return self.spectrum.shape[-self.block_depth:] 194 195 @property 196 def spectrum(self): 197 return self._spectrum 198 199 def _vectorize_then_blockify(self, matrix): 200 """Shape batch matrix to batch vector, then blockify trailing dimensions.""" 201 # Suppose 202 # matrix.shape = [m0, m1, m2, m3], 203 # and matrix is a matrix because the final two dimensions are matrix dims. 204 # self.block_depth = 2, 205 # self.block_shape = [b0, b1] (note b0 * b1 = m2). 206 # We will reshape matrix to 207 # [m3, m0, m1, b0, b1]. 208 209 # Vectorize: Reshape to batch vector. 210 # [m0, m1, m2, m3] --> [m3, m0, m1, m2] 211 # This is called "vectorize" because we have taken the final two matrix dims 212 # and turned this into a size m3 batch of vectors. 213 vec = distribution_util.rotate_transpose(matrix, shift=1) 214 215 # Blockify: Blockfy trailing dimensions. 216 # [m3, m0, m1, m2] --> [m3, m0, m1, b0, b1] 217 if (vec.shape.is_fully_defined() and 218 self.block_shape.is_fully_defined()): 219 # vec_leading_shape = [m3, m0, m1], 220 # the parts of vec that will not be blockified. 221 vec_leading_shape = vec.shape[:-1] 222 final_shape = vec_leading_shape.concatenate(self.block_shape) 223 else: 224 vec_leading_shape = array_ops.shape(vec)[:-1] 225 final_shape = array_ops.concat( 226 (vec_leading_shape, self.block_shape_tensor()), 0) 227 return array_ops.reshape(vec, final_shape) 228 229 def _unblockify_then_matricize(self, vec): 230 """Flatten the block dimensions then reshape to a batch matrix.""" 231 # Suppose 232 # vec.shape = [v0, v1, v2, v3], 233 # self.block_depth = 2. 234 # Then 235 # leading shape = [v0, v1] 236 # block shape = [v2, v3]. 237 # We will reshape vec to 238 # [v1, v2*v3, v0]. 239 240 # Un-blockify: Flatten block dimensions. Reshape 241 # [v0, v1, v2, v3] --> [v0, v1, v2*v3]. 242 if vec.shape.is_fully_defined(): 243 # vec_shape = [v0, v1, v2, v3] 244 vec_shape = vec.shape.as_list() 245 # vec_leading_shape = [v0, v1] 246 vec_leading_shape = vec_shape[:-self.block_depth] 247 # vec_block_shape = [v2, v3] 248 vec_block_shape = vec_shape[-self.block_depth:] 249 # flat_shape = [v0, v1, v2*v3] 250 flat_shape = vec_leading_shape + [np.prod(vec_block_shape)] 251 else: 252 vec_shape = array_ops.shape(vec) 253 vec_leading_shape = vec_shape[:-self.block_depth] 254 vec_block_shape = vec_shape[-self.block_depth:] 255 flat_shape = array_ops.concat( 256 (vec_leading_shape, [math_ops.reduce_prod(vec_block_shape)]), 0) 257 vec_flat = array_ops.reshape(vec, flat_shape) 258 259 # Matricize: Reshape to batch matrix. 260 # [v0, v1, v2*v3] --> [v1, v2*v3, v0], 261 # representing a shape [v1] batch of [v2*v3, v0] matrices. 262 matrix = distribution_util.rotate_transpose(vec_flat, shift=-1) 263 return matrix 264 265 def _fft(self, x): 266 """FFT along the last self.block_depth dimensions of x. 267 268 Args: 269 x: `Tensor` with floating or complex `dtype`. 270 Should be in the form returned by self._vectorize_then_blockify. 271 272 Returns: 273 `Tensor` with `dtype` `complex64`. 274 """ 275 x_complex = _to_complex(x) 276 return _FFT_OP[self.block_depth](x_complex) 277 278 def _ifft(self, x): 279 """IFFT along the last self.block_depth dimensions of x. 280 281 Args: 282 x: `Tensor` with floating or complex dtype. Should be in the form 283 returned by self._vectorize_then_blockify. 284 285 Returns: 286 `Tensor` with `dtype` `complex64`. 287 """ 288 x_complex = _to_complex(x) 289 return _IFFT_OP[self.block_depth](x_complex) 290 291 def convolution_kernel(self, name="convolution_kernel"): 292 """Convolution kernel corresponding to `self.spectrum`. 293 294 The `D` dimensional DFT of this kernel is the frequency domain spectrum of 295 this operator. 296 297 Args: 298 name: A name to give this `Op`. 299 300 Returns: 301 `Tensor` with `dtype` `self.dtype`. 302 """ 303 with self._name_scope(name): # pylint: disable=not-callable 304 h = self._ifft(_to_complex(self.spectrum)) 305 return math_ops.cast(h, self.dtype) 306 307 def _shape(self): 308 s_shape = self._spectrum.shape 309 # Suppose spectrum.shape = [a, b, c, d] 310 # block_depth = 2 311 # Then: 312 # batch_shape = [a, b] 313 # N = c*d 314 # and we want to return 315 # [a, b, c*d, c*d] 316 batch_shape = s_shape[:-self.block_depth] 317 # trailing_dims = [c, d] 318 trailing_dims = s_shape[-self.block_depth:] 319 if trailing_dims.is_fully_defined(): 320 n = np.prod(trailing_dims.as_list()) 321 else: 322 n = None 323 n_x_n = tensor_shape.TensorShape([n, n]) 324 return batch_shape.concatenate(n_x_n) 325 326 def _shape_tensor(self, spectrum=None): 327 spectrum = self.spectrum if spectrum is None else spectrum 328 # See self.shape for explanation of steps 329 s_shape = array_ops.shape(spectrum) 330 batch_shape = s_shape[:-self.block_depth] 331 trailing_dims = s_shape[-self.block_depth:] 332 n = math_ops.reduce_prod(trailing_dims) 333 n_x_n = [n, n] 334 return array_ops.concat((batch_shape, n_x_n), 0) 335 336 def assert_hermitian_spectrum(self, name="assert_hermitian_spectrum"): 337 """Returns an `Op` that asserts this operator has Hermitian spectrum. 338 339 This operator corresponds to a real-valued matrix if and only if its 340 spectrum is Hermitian. 341 342 Args: 343 name: A name to give this `Op`. 344 345 Returns: 346 An `Op` that asserts this operator has Hermitian spectrum. 347 """ 348 eps = np.finfo(self.dtype.real_dtype.as_numpy_dtype).eps 349 with self._name_scope(name): # pylint: disable=not-callable 350 # Assume linear accumulation of error. 351 max_err = eps * self.domain_dimension_tensor() 352 imag_convolution_kernel = math_ops.imag(self.convolution_kernel()) 353 return check_ops.assert_less( 354 math_ops.abs(imag_convolution_kernel), 355 max_err, 356 message="Spectrum was not Hermitian") 357 358 def _assert_non_singular(self): 359 return linear_operator_util.assert_no_entries_with_modulus_zero( 360 self.spectrum, 361 message="Singular operator: Spectrum contained zero values.") 362 363 def _assert_positive_definite(self): 364 # This operator has the action Ax = F^H D F x, 365 # where D is the diagonal matrix with self.spectrum on the diag. Therefore, 366 # <x, Ax> = <Fx, DFx>, 367 # Since F is bijective, the condition for positive definite is the same as 368 # for a diagonal matrix, i.e. real part of spectrum is positive. 369 message = ( 370 "Not positive definite: Real part of spectrum was not all positive.") 371 return check_ops.assert_positive( 372 math_ops.real(self.spectrum), message=message) 373 374 def _assert_self_adjoint(self): 375 # Recall correspondence between symmetry and real transforms. See docstring 376 return linear_operator_util.assert_zero_imag_part( 377 self.spectrum, 378 message=( 379 "Not self-adjoint: The spectrum contained non-zero imaginary part." 380 )) 381 382 def _broadcast_batch_dims(self, x, spectrum): 383 """Broadcast batch dims of batch matrix `x` and spectrum.""" 384 spectrum = ops.convert_to_tensor_v2_with_dispatch(spectrum, name="spectrum") 385 # spectrum.shape = batch_shape + block_shape 386 # First make spectrum a batch matrix with 387 # spectrum.shape = batch_shape + [prod(block_shape), 1] 388 batch_shape = self._batch_shape_tensor( 389 shape=self._shape_tensor(spectrum=spectrum)) 390 spec_mat = array_ops.reshape( 391 spectrum, array_ops.concat((batch_shape, [-1, 1]), axis=0)) 392 # Second, broadcast, possibly requiring an addition of array of zeros. 393 x, spec_mat = linear_operator_util.broadcast_matrix_batch_dims((x, 394 spec_mat)) 395 # Third, put the block shape back into spectrum. 396 x_batch_shape = array_ops.shape(x)[:-2] 397 spectrum_shape = array_ops.shape(spectrum) 398 spectrum = array_ops.reshape( 399 spec_mat, 400 array_ops.concat( 401 (x_batch_shape, 402 self._block_shape_tensor(spectrum_shape=spectrum_shape)), 403 axis=0)) 404 405 return x, spectrum 406 407 def _matmul(self, x, adjoint=False, adjoint_arg=False): 408 x = linalg.adjoint(x) if adjoint_arg else x 409 # With F the matrix of a DFT, and F^{-1}, F^H the inverse and Hermitian 410 # transpose, one can show that F^{-1} = F^{H} is the IDFT matrix. Therefore 411 # matmul(x) = F^{-1} diag(spectrum) F x, 412 # = F^{H} diag(spectrum) F x, 413 # so that 414 # matmul(x, adjoint=True) = F^{H} diag(conj(spectrum)) F x. 415 spectrum = _to_complex(self.spectrum) 416 if adjoint: 417 spectrum = math_ops.conj(spectrum) 418 419 x = math_ops.cast(x, spectrum.dtype) 420 421 x, spectrum = self._broadcast_batch_dims(x, spectrum) 422 423 x_vb = self._vectorize_then_blockify(x) 424 fft_x_vb = self._fft(x_vb) 425 block_vector_result = self._ifft(spectrum * fft_x_vb) 426 y = self._unblockify_then_matricize(block_vector_result) 427 428 return math_ops.cast(y, self.dtype) 429 430 def _determinant(self): 431 axis = [-(i + 1) for i in range(self.block_depth)] 432 det = math_ops.reduce_prod(self.spectrum, axis=axis) 433 return math_ops.cast(det, self.dtype) 434 435 def _log_abs_determinant(self): 436 axis = [-(i + 1) for i in range(self.block_depth)] 437 lad = math_ops.reduce_sum( 438 math_ops.log(math_ops.abs(self.spectrum)), axis=axis) 439 return math_ops.cast(lad, self.dtype) 440 441 def _solve(self, rhs, adjoint=False, adjoint_arg=False): 442 rhs = linalg.adjoint(rhs) if adjoint_arg else rhs 443 spectrum = _to_complex(self.spectrum) 444 if adjoint: 445 spectrum = math_ops.conj(spectrum) 446 447 rhs, spectrum = self._broadcast_batch_dims(rhs, spectrum) 448 449 rhs_vb = self._vectorize_then_blockify(rhs) 450 fft_rhs_vb = self._fft(rhs_vb) 451 solution_vb = self._ifft(fft_rhs_vb / spectrum) 452 x = self._unblockify_then_matricize(solution_vb) 453 return math_ops.cast(x, self.dtype) 454 455 def _diag_part(self): 456 # Get ones in shape of diag, which is [B1,...,Bb, N] 457 # Also get the size of the diag, "N". 458 if self.shape.is_fully_defined(): 459 diag_shape = self.shape[:-1] 460 diag_size = self.domain_dimension.value 461 else: 462 diag_shape = self.shape_tensor()[:-1] 463 diag_size = self.domain_dimension_tensor() 464 ones_diag = array_ops.ones(diag_shape, dtype=self.dtype) 465 466 # As proved in comments in self._trace, the value on the diag is constant, 467 # repeated N times. This value is the trace divided by N. 468 469 # The handling of self.shape = (0, 0) is tricky, and is the reason we choose 470 # to compute trace and use that to compute diag_part, rather than computing 471 # the value on the diagonal ("diag_value") directly. Both result in a 0/0, 472 # but in different places, and the current method gives the right result in 473 # the end. 474 475 # Here, if self.shape = (0, 0), then self.trace() = 0., and then 476 # diag_value = 0. / 0. = NaN. 477 diag_value = self.trace() / math_ops.cast(diag_size, self.dtype) 478 479 # If self.shape = (0, 0), then ones_diag = [] (empty tensor), and then 480 # the following line is NaN * [] = [], as needed. 481 return diag_value[..., array_ops.newaxis] * ones_diag 482 483 def _trace(self): 484 # The diagonal of the [[nested] block] circulant operator is the mean of 485 # the spectrum. 486 # Proof: For the [0,...,0] element, this follows from the IDFT formula. 487 # Then the result follows since all diagonal elements are the same. 488 489 # Therefore, the trace is the sum of the spectrum. 490 491 # Get shape of diag along with the axis over which to reduce the spectrum. 492 # We will reduce the spectrum over all block indices. 493 if self.spectrum.shape.is_fully_defined(): 494 spec_rank = self.spectrum.shape.ndims 495 axis = np.arange(spec_rank - self.block_depth, spec_rank, dtype=np.int32) 496 else: 497 spec_rank = array_ops.rank(self.spectrum) 498 axis = math_ops.range(spec_rank - self.block_depth, spec_rank) 499 500 # Real diag part "re_d". 501 # Suppose spectrum.shape = [B1,...,Bb, N1, N2] 502 # self.shape = [B1,...,Bb, N, N], with N1 * N2 = N. 503 # re_d_value.shape = [B1,...,Bb] 504 re_d_value = math_ops.reduce_sum(math_ops.real(self.spectrum), axis=axis) 505 506 if not self.dtype.is_complex: 507 return math_ops.cast(re_d_value, self.dtype) 508 509 # Imaginary part, "im_d". 510 if self.is_self_adjoint: 511 im_d_value = array_ops.zeros_like(re_d_value) 512 else: 513 im_d_value = math_ops.reduce_sum(math_ops.imag(self.spectrum), axis=axis) 514 515 return math_ops.cast(math_ops.complex(re_d_value, im_d_value), self.dtype) 516 517 @property 518 def _composite_tensor_fields(self): 519 return ("spectrum", "input_output_dtype") 520 521 522@tf_export("linalg.LinearOperatorCirculant") 523@linear_operator.make_composite_tensor 524class LinearOperatorCirculant(_BaseLinearOperatorCirculant): 525 """`LinearOperator` acting like a circulant matrix. 526 527 This operator acts like a circulant matrix `A` with 528 shape `[B1,...,Bb, N, N]` for some `b >= 0`. The first `b` indices index a 529 batch member. For every batch index `(i1,...,ib)`, `A[i1,...,ib, : :]` is 530 an `N x N` matrix. This matrix `A` is not materialized, but for 531 purposes of broadcasting this shape will be relevant. 532 533 #### Description in terms of circulant matrices 534 535 Circulant means the entries of `A` are generated by a single vector, the 536 convolution kernel `h`: `A_{mn} := h_{m-n mod N}`. With `h = [w, x, y, z]`, 537 538 ``` 539 A = |w z y x| 540 |x w z y| 541 |y x w z| 542 |z y x w| 543 ``` 544 545 This means that the result of matrix multiplication `v = Au` has `Lth` column 546 given circular convolution between `h` with the `Lth` column of `u`. 547 548 #### Description in terms of the frequency spectrum 549 550 There is an equivalent description in terms of the [batch] spectrum `H` and 551 Fourier transforms. Here we consider `A.shape = [N, N]` and ignore batch 552 dimensions. Define the discrete Fourier transform (DFT) and its inverse by 553 554 ``` 555 DFT[ h[n] ] = H[k] := sum_{n = 0}^{N - 1} h_n e^{-i 2pi k n / N} 556 IDFT[ H[k] ] = h[n] = N^{-1} sum_{k = 0}^{N - 1} H_k e^{i 2pi k n / N} 557 ``` 558 559 From these definitions, we see that 560 561 ``` 562 H[0] = sum_{n = 0}^{N - 1} h_n 563 H[1] = "the first positive frequency" 564 H[N - 1] = "the first negative frequency" 565 ``` 566 567 Loosely speaking, with `*` element-wise multiplication, matrix multiplication 568 is equal to the action of a Fourier multiplier: `A u = IDFT[ H * DFT[u] ]`. 569 Precisely speaking, given `[N, R]` matrix `u`, let `DFT[u]` be the `[N, R]` 570 matrix with `rth` column equal to the DFT of the `rth` column of `u`. 571 Define the `IDFT` similarly. 572 Matrix multiplication may be expressed columnwise: 573 574 ```(A u)_r = IDFT[ H * (DFT[u])_r ]``` 575 576 #### Operator properties deduced from the spectrum. 577 578 Letting `U` be the `kth` Euclidean basis vector, and `U = IDFT[u]`. 579 The above formulas show that`A U = H_k * U`. We conclude that the elements 580 of `H` are the eigenvalues of this operator. Therefore 581 582 * This operator is positive definite if and only if `Real{H} > 0`. 583 584 A general property of Fourier transforms is the correspondence between 585 Hermitian functions and real valued transforms. 586 587 Suppose `H.shape = [B1,...,Bb, N]`. We say that `H` is a Hermitian spectrum 588 if, with `%` meaning modulus division, 589 590 ```H[..., n % N] = ComplexConjugate[ H[..., (-n) % N] ]``` 591 592 * This operator corresponds to a real matrix if and only if `H` is Hermitian. 593 * This operator is self-adjoint if and only if `H` is real. 594 595 See e.g. "Discrete-Time Signal Processing", Oppenheim and Schafer. 596 597 #### Example of a self-adjoint positive definite operator 598 599 ```python 600 # spectrum is real ==> operator is self-adjoint 601 # spectrum is positive ==> operator is positive definite 602 spectrum = [6., 4, 2] 603 604 operator = LinearOperatorCirculant(spectrum) 605 606 # IFFT[spectrum] 607 operator.convolution_kernel() 608 ==> [4 + 0j, 1 + 0.58j, 1 - 0.58j] 609 610 operator.to_dense() 611 ==> [[4 + 0.0j, 1 - 0.6j, 1 + 0.6j], 612 [1 + 0.6j, 4 + 0.0j, 1 - 0.6j], 613 [1 - 0.6j, 1 + 0.6j, 4 + 0.0j]] 614 ``` 615 616 #### Example of defining in terms of a real convolution kernel 617 618 ```python 619 # convolution_kernel is real ==> spectrum is Hermitian. 620 convolution_kernel = [1., 2., 1.]] 621 spectrum = tf.signal.fft(tf.cast(convolution_kernel, tf.complex64)) 622 623 # spectrum is Hermitian ==> operator is real. 624 # spectrum is shape [3] ==> operator is shape [3, 3] 625 # We force the input/output type to be real, which allows this to operate 626 # like a real matrix. 627 operator = LinearOperatorCirculant(spectrum, input_output_dtype=tf.float32) 628 629 operator.to_dense() 630 ==> [[ 1, 1, 2], 631 [ 2, 1, 1], 632 [ 1, 2, 1]] 633 ``` 634 635 #### Example of Hermitian spectrum 636 637 ```python 638 # spectrum is shape [3] ==> operator is shape [3, 3] 639 # spectrum is Hermitian ==> operator is real. 640 spectrum = [1, 1j, -1j] 641 642 operator = LinearOperatorCirculant(spectrum) 643 644 operator.to_dense() 645 ==> [[ 0.33 + 0j, 0.91 + 0j, -0.24 + 0j], 646 [-0.24 + 0j, 0.33 + 0j, 0.91 + 0j], 647 [ 0.91 + 0j, -0.24 + 0j, 0.33 + 0j] 648 ``` 649 650 #### Example of forcing real `dtype` when spectrum is Hermitian 651 652 ```python 653 # spectrum is shape [4] ==> operator is shape [4, 4] 654 # spectrum is real ==> operator is self-adjoint 655 # spectrum is Hermitian ==> operator is real 656 # spectrum has positive real part ==> operator is positive-definite. 657 spectrum = [6., 4, 2, 4] 658 659 # Force the input dtype to be float32. 660 # Cast the output to float32. This is fine because the operator will be 661 # real due to Hermitian spectrum. 662 operator = LinearOperatorCirculant(spectrum, input_output_dtype=tf.float32) 663 664 operator.shape 665 ==> [4, 4] 666 667 operator.to_dense() 668 ==> [[4, 1, 0, 1], 669 [1, 4, 1, 0], 670 [0, 1, 4, 1], 671 [1, 0, 1, 4]] 672 673 # convolution_kernel = tf.signal.ifft(spectrum) 674 operator.convolution_kernel() 675 ==> [4, 1, 0, 1] 676 ``` 677 678 #### Performance 679 680 Suppose `operator` is a `LinearOperatorCirculant` of shape `[N, N]`, 681 and `x.shape = [N, R]`. Then 682 683 * `operator.matmul(x)` is `O(R*N*Log[N])` 684 * `operator.solve(x)` is `O(R*N*Log[N])` 685 * `operator.determinant()` involves a size `N` `reduce_prod`. 686 687 If instead `operator` and `x` have shape `[B1,...,Bb, N, N]` and 688 `[B1,...,Bb, N, R]`, every operation increases in complexity by `B1*...*Bb`. 689 690 #### Matrix property hints 691 692 This `LinearOperator` is initialized with boolean flags of the form `is_X`, 693 for `X = non_singular, self_adjoint, positive_definite, square`. 694 These have the following meaning: 695 696 * If `is_X == True`, callers should expect the operator to have the 697 property `X`. This is a promise that should be fulfilled, but is *not* a 698 runtime assert. For example, finite floating point precision may result 699 in these promises being violated. 700 * If `is_X == False`, callers should expect the operator to not have `X`. 701 * If `is_X == None` (the default), callers should have no expectation either 702 way. 703 704 References: 705 Toeplitz and Circulant Matrices - A Review: 706 [Gray, 2006](https://www.nowpublishers.com/article/Details/CIT-006) 707 ([pdf](https://ee.stanford.edu/~gray/toeplitz.pdf)) 708 """ 709 710 def __init__(self, 711 spectrum, 712 input_output_dtype=dtypes.complex64, 713 is_non_singular=None, 714 is_self_adjoint=None, 715 is_positive_definite=None, 716 is_square=True, 717 name="LinearOperatorCirculant"): 718 r"""Initialize an `LinearOperatorCirculant`. 719 720 This `LinearOperator` is initialized to have shape `[B1,...,Bb, N, N]` 721 by providing `spectrum`, a `[B1,...,Bb, N]` `Tensor`. 722 723 If `input_output_dtype = DTYPE`: 724 725 * Arguments to methods such as `matmul` or `solve` must be `DTYPE`. 726 * Values returned by all methods, such as `matmul` or `determinant` will be 727 cast to `DTYPE`. 728 729 Note that if the spectrum is not Hermitian, then this operator corresponds 730 to a complex matrix with non-zero imaginary part. In this case, setting 731 `input_output_dtype` to a real type will forcibly cast the output to be 732 real, resulting in incorrect results! 733 734 If on the other hand the spectrum is Hermitian, then this operator 735 corresponds to a real-valued matrix, and setting `input_output_dtype` to 736 a real type is fine. 737 738 Args: 739 spectrum: Shape `[B1,...,Bb, N]` `Tensor`. Allowed dtypes: `float16`, 740 `float32`, `float64`, `complex64`, `complex128`. Type can be different 741 than `input_output_dtype` 742 input_output_dtype: `dtype` for input/output. 743 is_non_singular: Expect that this operator is non-singular. 744 is_self_adjoint: Expect that this operator is equal to its hermitian 745 transpose. If `spectrum` is real, this will always be true. 746 is_positive_definite: Expect that this operator is positive definite, 747 meaning the quadratic form `x^H A x` has positive real part for all 748 nonzero `x`. Note that we do not require the operator to be 749 self-adjoint to be positive-definite. See: 750 https://en.wikipedia.org/wiki/Positive-definite_matrix\ 751 #Extension_for_non_symmetric_matrices 752 is_square: Expect that this operator acts like square [batch] matrices. 753 name: A name to prepend to all ops created by this class. 754 """ 755 parameters = dict( 756 spectrum=spectrum, 757 input_output_dtype=input_output_dtype, 758 is_non_singular=is_non_singular, 759 is_self_adjoint=is_self_adjoint, 760 is_positive_definite=is_positive_definite, 761 is_square=is_square, 762 name=name 763 ) 764 super(LinearOperatorCirculant, self).__init__( 765 spectrum, 766 block_depth=1, 767 input_output_dtype=input_output_dtype, 768 is_non_singular=is_non_singular, 769 is_self_adjoint=is_self_adjoint, 770 is_positive_definite=is_positive_definite, 771 is_square=is_square, 772 parameters=parameters, 773 name=name) 774 775 def _eigvals(self): 776 return ops.convert_to_tensor_v2_with_dispatch(self.spectrum) 777 778 779@tf_export("linalg.LinearOperatorCirculant2D") 780@linear_operator.make_composite_tensor 781class LinearOperatorCirculant2D(_BaseLinearOperatorCirculant): 782 """`LinearOperator` acting like a block circulant matrix. 783 784 This operator acts like a block circulant matrix `A` with 785 shape `[B1,...,Bb, N, N]` for some `b >= 0`. The first `b` indices index a 786 batch member. For every batch index `(i1,...,ib)`, `A[i1,...,ib, : :]` is 787 an `N x N` matrix. This matrix `A` is not materialized, but for 788 purposes of broadcasting this shape will be relevant. 789 790 #### Description in terms of block circulant matrices 791 792 If `A` is block circulant, with block sizes `N0, N1` (`N0 * N1 = N`): 793 `A` has a block circulant structure, composed of `N0 x N0` blocks, with each 794 block an `N1 x N1` circulant matrix. 795 796 For example, with `W`, `X`, `Y`, `Z` each circulant, 797 798 ``` 799 A = |W Z Y X| 800 |X W Z Y| 801 |Y X W Z| 802 |Z Y X W| 803 ``` 804 805 Note that `A` itself will not in general be circulant. 806 807 #### Description in terms of the frequency spectrum 808 809 There is an equivalent description in terms of the [batch] spectrum `H` and 810 Fourier transforms. Here we consider `A.shape = [N, N]` and ignore batch 811 dimensions. 812 813 If `H.shape = [N0, N1]`, (`N0 * N1 = N`): 814 Loosely speaking, matrix multiplication is equal to the action of a 815 Fourier multiplier: `A u = IDFT2[ H DFT2[u] ]`. 816 Precisely speaking, given `[N, R]` matrix `u`, let `DFT2[u]` be the 817 `[N0, N1, R]` `Tensor` defined by re-shaping `u` to `[N0, N1, R]` and taking 818 a two dimensional DFT across the first two dimensions. Let `IDFT2` be the 819 inverse of `DFT2`. Matrix multiplication may be expressed columnwise: 820 821 ```(A u)_r = IDFT2[ H * (DFT2[u])_r ]``` 822 823 #### Operator properties deduced from the spectrum. 824 825 * This operator is positive definite if and only if `Real{H} > 0`. 826 827 A general property of Fourier transforms is the correspondence between 828 Hermitian functions and real valued transforms. 829 830 Suppose `H.shape = [B1,...,Bb, N0, N1]`, we say that `H` is a Hermitian 831 spectrum if, with `%` indicating modulus division, 832 833 ``` 834 H[..., n0 % N0, n1 % N1] = ComplexConjugate[ H[..., (-n0) % N0, (-n1) % N1 ]. 835 ``` 836 837 * This operator corresponds to a real matrix if and only if `H` is Hermitian. 838 * This operator is self-adjoint if and only if `H` is real. 839 840 See e.g. "Discrete-Time Signal Processing", Oppenheim and Schafer. 841 842 ### Example of a self-adjoint positive definite operator 843 844 ```python 845 # spectrum is real ==> operator is self-adjoint 846 # spectrum is positive ==> operator is positive definite 847 spectrum = [[1., 2., 3.], 848 [4., 5., 6.], 849 [7., 8., 9.]] 850 851 operator = LinearOperatorCirculant2D(spectrum) 852 853 # IFFT[spectrum] 854 operator.convolution_kernel() 855 ==> [[5.0+0.0j, -0.5-.3j, -0.5+.3j], 856 [-1.5-.9j, 0, 0], 857 [-1.5+.9j, 0, 0]] 858 859 operator.to_dense() 860 ==> Complex self adjoint 9 x 9 matrix. 861 ``` 862 863 #### Example of defining in terms of a real convolution kernel, 864 865 ```python 866 # convolution_kernel is real ==> spectrum is Hermitian. 867 convolution_kernel = [[1., 2., 1.], [5., -1., 1.]] 868 spectrum = tf.signal.fft2d(tf.cast(convolution_kernel, tf.complex64)) 869 870 # spectrum is shape [2, 3] ==> operator is shape [6, 6] 871 # spectrum is Hermitian ==> operator is real. 872 operator = LinearOperatorCirculant2D(spectrum, input_output_dtype=tf.float32) 873 ``` 874 875 #### Performance 876 877 Suppose `operator` is a `LinearOperatorCirculant` of shape `[N, N]`, 878 and `x.shape = [N, R]`. Then 879 880 * `operator.matmul(x)` is `O(R*N*Log[N])` 881 * `operator.solve(x)` is `O(R*N*Log[N])` 882 * `operator.determinant()` involves a size `N` `reduce_prod`. 883 884 If instead `operator` and `x` have shape `[B1,...,Bb, N, N]` and 885 `[B1,...,Bb, N, R]`, every operation increases in complexity by `B1*...*Bb`. 886 887 #### Matrix property hints 888 889 This `LinearOperator` is initialized with boolean flags of the form `is_X`, 890 for `X = non_singular, self_adjoint, positive_definite, square`. 891 These have the following meaning 892 * If `is_X == True`, callers should expect the operator to have the 893 property `X`. This is a promise that should be fulfilled, but is *not* a 894 runtime assert. For example, finite floating point precision may result 895 in these promises being violated. 896 * If `is_X == False`, callers should expect the operator to not have `X`. 897 * If `is_X == None` (the default), callers should have no expectation either 898 way. 899 """ 900 901 def __init__(self, 902 spectrum, 903 input_output_dtype=dtypes.complex64, 904 is_non_singular=None, 905 is_self_adjoint=None, 906 is_positive_definite=None, 907 is_square=True, 908 name="LinearOperatorCirculant2D"): 909 r"""Initialize an `LinearOperatorCirculant2D`. 910 911 This `LinearOperator` is initialized to have shape `[B1,...,Bb, N, N]` 912 by providing `spectrum`, a `[B1,...,Bb, N0, N1]` `Tensor` with `N0*N1 = N`. 913 914 If `input_output_dtype = DTYPE`: 915 916 * Arguments to methods such as `matmul` or `solve` must be `DTYPE`. 917 * Values returned by all methods, such as `matmul` or `determinant` will be 918 cast to `DTYPE`. 919 920 Note that if the spectrum is not Hermitian, then this operator corresponds 921 to a complex matrix with non-zero imaginary part. In this case, setting 922 `input_output_dtype` to a real type will forcibly cast the output to be 923 real, resulting in incorrect results! 924 925 If on the other hand the spectrum is Hermitian, then this operator 926 corresponds to a real-valued matrix, and setting `input_output_dtype` to 927 a real type is fine. 928 929 Args: 930 spectrum: Shape `[B1,...,Bb, N]` `Tensor`. Allowed dtypes: `float16`, 931 `float32`, `float64`, `complex64`, `complex128`. Type can be different 932 than `input_output_dtype` 933 input_output_dtype: `dtype` for input/output. 934 is_non_singular: Expect that this operator is non-singular. 935 is_self_adjoint: Expect that this operator is equal to its hermitian 936 transpose. If `spectrum` is real, this will always be true. 937 is_positive_definite: Expect that this operator is positive definite, 938 meaning the quadratic form `x^H A x` has positive real part for all 939 nonzero `x`. Note that we do not require the operator to be 940 self-adjoint to be positive-definite. See: 941 https://en.wikipedia.org/wiki/Positive-definite_matrix\ 942 #Extension_for_non_symmetric_matrices 943 is_square: Expect that this operator acts like square [batch] matrices. 944 name: A name to prepend to all ops created by this class. 945 """ 946 parameters = dict( 947 spectrum=spectrum, 948 input_output_dtype=input_output_dtype, 949 is_non_singular=is_non_singular, 950 is_self_adjoint=is_self_adjoint, 951 is_positive_definite=is_positive_definite, 952 is_square=is_square, 953 name=name 954 ) 955 super(LinearOperatorCirculant2D, self).__init__( 956 spectrum, 957 block_depth=2, 958 input_output_dtype=input_output_dtype, 959 is_non_singular=is_non_singular, 960 is_self_adjoint=is_self_adjoint, 961 is_positive_definite=is_positive_definite, 962 is_square=is_square, 963 parameters=parameters, 964 name=name) 965 966 967@tf_export("linalg.LinearOperatorCirculant3D") 968@linear_operator.make_composite_tensor 969class LinearOperatorCirculant3D(_BaseLinearOperatorCirculant): 970 """`LinearOperator` acting like a nested block circulant matrix. 971 972 This operator acts like a block circulant matrix `A` with 973 shape `[B1,...,Bb, N, N]` for some `b >= 0`. The first `b` indices index a 974 batch member. For every batch index `(i1,...,ib)`, `A[i1,...,ib, : :]` is 975 an `N x N` matrix. This matrix `A` is not materialized, but for 976 purposes of broadcasting this shape will be relevant. 977 978 #### Description in terms of block circulant matrices 979 980 If `A` is nested block circulant, with block sizes `N0, N1, N2` 981 (`N0 * N1 * N2 = N`): 982 `A` has a block structure, composed of `N0 x N0` blocks, with each 983 block an `N1 x N1` block circulant matrix. 984 985 For example, with `W`, `X`, `Y`, `Z` each block circulant, 986 987 ``` 988 A = |W Z Y X| 989 |X W Z Y| 990 |Y X W Z| 991 |Z Y X W| 992 ``` 993 994 Note that `A` itself will not in general be circulant. 995 996 #### Description in terms of the frequency spectrum 997 998 There is an equivalent description in terms of the [batch] spectrum `H` and 999 Fourier transforms. Here we consider `A.shape = [N, N]` and ignore batch 1000 dimensions. 1001 1002 If `H.shape = [N0, N1, N2]`, (`N0 * N1 * N2 = N`): 1003 Loosely speaking, matrix multiplication is equal to the action of a 1004 Fourier multiplier: `A u = IDFT3[ H DFT3[u] ]`. 1005 Precisely speaking, given `[N, R]` matrix `u`, let `DFT3[u]` be the 1006 `[N0, N1, N2, R]` `Tensor` defined by re-shaping `u` to `[N0, N1, N2, R]` and 1007 taking a three dimensional DFT across the first three dimensions. Let `IDFT3` 1008 be the inverse of `DFT3`. Matrix multiplication may be expressed columnwise: 1009 1010 ```(A u)_r = IDFT3[ H * (DFT3[u])_r ]``` 1011 1012 #### Operator properties deduced from the spectrum. 1013 1014 * This operator is positive definite if and only if `Real{H} > 0`. 1015 1016 A general property of Fourier transforms is the correspondence between 1017 Hermitian functions and real valued transforms. 1018 1019 Suppose `H.shape = [B1,...,Bb, N0, N1, N2]`, we say that `H` is a Hermitian 1020 spectrum if, with `%` meaning modulus division, 1021 1022 ``` 1023 H[..., n0 % N0, n1 % N1, n2 % N2] 1024 = ComplexConjugate[ H[..., (-n0) % N0, (-n1) % N1, (-n2) % N2] ]. 1025 ``` 1026 1027 * This operator corresponds to a real matrix if and only if `H` is Hermitian. 1028 * This operator is self-adjoint if and only if `H` is real. 1029 1030 See e.g. "Discrete-Time Signal Processing", Oppenheim and Schafer. 1031 1032 ### Examples 1033 1034 See `LinearOperatorCirculant` and `LinearOperatorCirculant2D` for examples. 1035 1036 #### Performance 1037 1038 Suppose `operator` is a `LinearOperatorCirculant` of shape `[N, N]`, 1039 and `x.shape = [N, R]`. Then 1040 1041 * `operator.matmul(x)` is `O(R*N*Log[N])` 1042 * `operator.solve(x)` is `O(R*N*Log[N])` 1043 * `operator.determinant()` involves a size `N` `reduce_prod`. 1044 1045 If instead `operator` and `x` have shape `[B1,...,Bb, N, N]` and 1046 `[B1,...,Bb, N, R]`, every operation increases in complexity by `B1*...*Bb`. 1047 1048 #### Matrix property hints 1049 1050 This `LinearOperator` is initialized with boolean flags of the form `is_X`, 1051 for `X = non_singular, self_adjoint, positive_definite, square`. 1052 These have the following meaning 1053 * If `is_X == True`, callers should expect the operator to have the 1054 property `X`. This is a promise that should be fulfilled, but is *not* a 1055 runtime assert. For example, finite floating point precision may result 1056 in these promises being violated. 1057 * If `is_X == False`, callers should expect the operator to not have `X`. 1058 * If `is_X == None` (the default), callers should have no expectation either 1059 way. 1060 """ 1061 1062 def __init__(self, 1063 spectrum, 1064 input_output_dtype=dtypes.complex64, 1065 is_non_singular=None, 1066 is_self_adjoint=None, 1067 is_positive_definite=None, 1068 is_square=True, 1069 name="LinearOperatorCirculant3D"): 1070 """Initialize an `LinearOperatorCirculant`. 1071 1072 This `LinearOperator` is initialized to have shape `[B1,...,Bb, N, N]` 1073 by providing `spectrum`, a `[B1,...,Bb, N0, N1, N2]` `Tensor` 1074 with `N0*N1*N2 = N`. 1075 1076 If `input_output_dtype = DTYPE`: 1077 1078 * Arguments to methods such as `matmul` or `solve` must be `DTYPE`. 1079 * Values returned by all methods, such as `matmul` or `determinant` will be 1080 cast to `DTYPE`. 1081 1082 Note that if the spectrum is not Hermitian, then this operator corresponds 1083 to a complex matrix with non-zero imaginary part. In this case, setting 1084 `input_output_dtype` to a real type will forcibly cast the output to be 1085 real, resulting in incorrect results! 1086 1087 If on the other hand the spectrum is Hermitian, then this operator 1088 corresponds to a real-valued matrix, and setting `input_output_dtype` to 1089 a real type is fine. 1090 1091 Args: 1092 spectrum: Shape `[B1,...,Bb, N]` `Tensor`. Allowed dtypes: `float16`, 1093 `float32`, `float64`, `complex64`, `complex128`. Type can be different 1094 than `input_output_dtype` 1095 input_output_dtype: `dtype` for input/output. 1096 is_non_singular: Expect that this operator is non-singular. 1097 is_self_adjoint: Expect that this operator is equal to its hermitian 1098 transpose. If `spectrum` is real, this will always be true. 1099 is_positive_definite: Expect that this operator is positive definite, 1100 meaning the real part of all eigenvalues is positive. We do not require 1101 the operator to be self-adjoint to be positive-definite. See: 1102 https://en.wikipedia.org/wiki/Positive-definite_matrix 1103 #Extension_for_non_symmetric_matrices 1104 is_square: Expect that this operator acts like square [batch] matrices. 1105 name: A name to prepend to all ops created by this class. 1106 """ 1107 parameters = dict( 1108 spectrum=spectrum, 1109 input_output_dtype=input_output_dtype, 1110 is_non_singular=is_non_singular, 1111 is_self_adjoint=is_self_adjoint, 1112 is_positive_definite=is_positive_definite, 1113 is_square=is_square, 1114 name=name 1115 ) 1116 super(LinearOperatorCirculant3D, self).__init__( 1117 spectrum, 1118 block_depth=3, 1119 input_output_dtype=input_output_dtype, 1120 is_non_singular=is_non_singular, 1121 is_self_adjoint=is_self_adjoint, 1122 is_positive_definite=is_positive_definite, 1123 is_square=is_square, 1124 parameters=parameters, 1125 name=name) 1126 1127 1128def _to_complex(x): 1129 if x.dtype.is_complex: 1130 return x 1131 dtype = dtypes.complex64 1132 1133 if x.dtype == dtypes.float64: 1134 dtype = dtypes.complex128 1135 return math_ops.cast(x, dtype) 1136