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