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