Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/tensorflow/python/ops/linalg/linear_operator_circulant.py: 28%

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_conversion 

22from tensorflow.python.framework import tensor_shape 

23from tensorflow.python.ops import array_ops 

24from tensorflow.python.ops import check_ops 

25from tensorflow.python.ops import math_ops 

26from tensorflow.python.ops.distributions import util as distribution_util 

27from tensorflow.python.ops.linalg import linalg_impl as linalg 

28from tensorflow.python.ops.linalg import linear_operator 

29from tensorflow.python.ops.linalg import linear_operator_util 

30from tensorflow.python.ops.signal import fft_ops 

31from tensorflow.python.util.tf_export import tf_export 

32 

33__all__ = [ 

34 "LinearOperatorCirculant", 

35 "LinearOperatorCirculant2D", 

36 "LinearOperatorCirculant3D", 

37] 

38 

39# Different FFT Ops will be used for different block depths. 

40_FFT_OP = {1: fft_ops.fft, 2: fft_ops.fft2d, 3: fft_ops.fft3d} 

41_IFFT_OP = {1: fft_ops.ifft, 2: fft_ops.ifft2d, 3: fft_ops.ifft3d} 

42 

43 

44def exponential_power_convolution_kernel( 

45 grid_shape, 

46 length_scale, 

47 power=None, 

48 divisor=None, 

49 zero_inflation=None, 

50): 

51 """Make an exponentiated convolution kernel. 

52 

53 In signal processing, a [kernel] 

54 (https://en.wikipedia.org/wiki/Kernel_(image_processing)) `h` can be convolved 

55 with a signal `x` to filter its spectral content. 

56 

57 This function makes a `d-dimensional` convolution kernel `h` of shape 

58 `grid_shape = [N0, N1, ...]`. For `n` a multi-index with `n[i] < Ni / 2`, 

59 

60 ```h[n] = exp{sum(|n / (length_scale * grid_shape)|**power) / divisor}.``` 

61 

62 For other `n`, `h` is extended to be circularly symmetric. That is 

63 

64 ```h[n0 % N0, ...] = h[(-n0) % N0, ...]``` 

65 

66 Since `h` is circularly symmetric and real valued, `H = FFTd[h]` is the 

67 spectrum of a symmetric (real) circulant operator `A`. 

68 

69 #### Example uses 

70 

71 ``` 

72 # Matern one-half kernel, d=1. 

73 # Will be positive definite without zero_inflation. 

74 h = exponential_power_convolution_kernel( 

75 grid_shape=[10], length_scale=[0.1], power=1) 

76 A = LinearOperatorCirculant( 

77 tf.signal.fft(tf.cast(h, tf.complex64)), 

78 is_self_adjoint=True, is_positive_definite=True) 

79 

80 # Gaussian RBF kernel, d=3. 

81 # Needs zero_inflation since `length_scale` is long enough to cause aliasing. 

82 h = exponential_power_convolution_kernel( 

83 grid_shape=[10, 10, 10], length_scale=[0.1, 0.2, 0.2], power=2, 

84 zero_inflation=0.15) 

85 A = LinearOperatorCirculant3D( 

86 tf.signal.fft3d(tf.cast(h, tf.complex64)), 

87 is_self_adjoint=True, is_positive_definite=True) 

88 ``` 

89 

90 Args: 

91 grid_shape: Length `d` (`d` in {1, 2, 3}) list-like of Python integers. The 

92 shape of the grid on which the convolution kernel is defined. 

93 length_scale: Length `d` `float` `Tensor`. The scale at which the kernel 

94 decays in each direction, as a fraction of `grid_shape`. 

95 power: Scalar `Tensor` of same `dtype` as `length_scale`, default `2`. 

96 Higher (lower) `power` results in nearby points being more (less) 

97 correlated, and far away points being less (more) correlated. 

98 divisor: Scalar `Tensor` of same `dtype` as `length_scale`. The slope of 

99 decay of `log(kernel)` in terms of fractional grid points, along each 

100 axis, at `length_scale`, is `power/divisor`. By default, `divisor` is set 

101 to `power`. This means, by default, `power=2` results in an exponentiated 

102 quadratic (Gaussian) kernel, and `power=1` is a Matern one-half. 

103 zero_inflation: Scalar `Tensor` of same `dtype` as `length_scale`, in 

104 `[0, 1]`. Let `delta` be the Kronecker delta. That is, 

105 `delta[0, ..., 0] = 1` and all other entries are `0`. Then 

106 `zero_inflation` modifies the return value via 

107 `h --> (1 - zero_inflation) * h + zero_inflation * delta`. This may be 

108 needed to ensure a positive definite kernel, especially if `length_scale` 

109 is large enough for aliasing and `power > 1`. 

110 

111 Returns: 

112 `Tensor` of shape `grid_shape` with same `dtype` as `length_scale`. 

113 """ 

114 nd = len(grid_shape) 

115 

116 length_scale = tensor_conversion.convert_to_tensor_v2_with_dispatch( 

117 length_scale, name="length_scale" 

118 ) 

119 dtype = length_scale.dtype 

120 

121 power = 2. if power is None else power 

122 power = tensor_conversion.convert_to_tensor_v2_with_dispatch( 

123 power, name="power", dtype=dtype 

124 ) 

125 divisor = power if divisor is None else divisor 

126 divisor = tensor_conversion.convert_to_tensor_v2_with_dispatch( 

127 divisor, name="divisor", dtype=dtype 

128 ) 

129 

130 # With K = grid_shape[i], we implicitly assume the grid vertices along the 

131 # ith dimension are at: 

132 # 0 = 0 / (K - 1), 1 / (K - 1), 2 / (K - 1), ..., (K - 1) / (K - 1) = 1. 

133 zero = math_ops.cast(0., dtype) 

134 one = math_ops.cast(1., dtype) 

135 ts = [math_ops.linspace(zero, one, num=n) for n in grid_shape] 

136 

137 log_vals = [] 

138 for i, x in enumerate(array_ops.meshgrid(*ts, indexing="ij")): 

139 # midpoint[i] is the vertex just to the left of 1 / 2. 

140 # ifftshift will shift this vertex to position 0. 

141 midpoint = ts[i][math_ops.cast( 

142 math_ops.floor(one / 2. * grid_shape[i]), dtypes.int32)] 

143 log_vals.append(-(math_ops.abs( 

144 (x - midpoint) / length_scale[i]))**power / divisor) 

145 kernel = math_ops.exp( 

146 fft_ops.ifftshift(sum(log_vals), axes=[-i for i in range(1, nd + 1)])) 

147 

148 if zero_inflation: 

149 # delta.shape = grid_shape, delta[0, 0, 0] = 1., all other entries are 0. 

150 zero_inflation = tensor_conversion.convert_to_tensor_v2_with_dispatch( 

151 zero_inflation, name="zero_inflation", dtype=dtype 

152 ) 

153 delta = array_ops.pad( 

154 array_ops.reshape(one, [1] * nd), [[0, dim - 1] for dim in grid_shape]) 

155 kernel = (1. - zero_inflation) * kernel + zero_inflation * delta 

156 

157 return kernel 

158 

159 

160# TODO(langmore) Add transformations that create common spectrums, e.g. 

161# starting with the convolution kernel 

162# start with half a spectrum, and create a Hermitian one. 

163# common filters. 

164# TODO(langmore) Support rectangular Toeplitz matrices. 

165class _BaseLinearOperatorCirculant(linear_operator.LinearOperator): 

166 """Base class for circulant operators. Not user facing. 

167 

168 `LinearOperator` acting like a [batch] [[nested] block] circulant matrix. 

169 """ 

170 

171 def __init__(self, 

172 spectrum, 

173 block_depth, 

174 input_output_dtype=dtypes.complex64, 

175 is_non_singular=None, 

176 is_self_adjoint=None, 

177 is_positive_definite=None, 

178 is_square=True, 

179 parameters=None, 

180 name="LinearOperatorCirculant"): 

181 r"""Initialize an `_BaseLinearOperatorCirculant`. 

182 

183 Args: 

184 spectrum: Shape `[B1,...,Bb] + N` `Tensor`, where `rank(N) in {1, 2, 3}`. 

185 Allowed dtypes: `float16`, `float32`, `float64`, `complex64`, 

186 `complex128`. Type can be different than `input_output_dtype` 

187 block_depth: Python integer, either 1, 2, or 3. Will be 1 for circulant, 

188 2 for block circulant, and 3 for nested block circulant. 

189 input_output_dtype: `dtype` for input/output. 

190 is_non_singular: Expect that this operator is non-singular. 

191 is_self_adjoint: Expect that this operator is equal to its hermitian 

192 transpose. If `spectrum` is real, this will always be true. 

193 is_positive_definite: Expect that this operator is positive definite, 

194 meaning the quadratic form `x^H A x` has positive real part for all 

195 nonzero `x`. Note that we do not require the operator to be 

196 self-adjoint to be positive-definite. See: 

197 https://en.wikipedia.org/wiki/Positive-definite_matrix\ 

198 #Extension_for_non_symmetric_matrices 

199 is_square: Expect that this operator acts like square [batch] matrices. 

200 parameters: Python `dict` of parameters used to instantiate this 

201 `LinearOperator`. 

202 name: A name to prepend to all ops created by this class. 

203 

204 Raises: 

205 ValueError: If `block_depth` is not an allowed value. 

206 TypeError: If `spectrum` is not an allowed type. 

207 """ 

208 

209 allowed_block_depths = [1, 2, 3] 

210 

211 self._name = name 

212 

213 if block_depth not in allowed_block_depths: 

214 raise ValueError( 

215 f"Argument `block_depth` must be one of {allowed_block_depths}. " 

216 f"Received: {block_depth}.") 

217 self._block_depth = block_depth 

218 

219 with ops.name_scope(name, values=[spectrum]): 

220 self._spectrum = self._check_spectrum_and_return_tensor(spectrum) 

221 

222 # Check and auto-set hints. 

223 if not self.spectrum.dtype.is_complex: 

224 if is_self_adjoint is False: 

225 raise ValueError( 

226 f"A real spectrum always corresponds to a self-adjoint operator. " 

227 f"Expected argument `is_self_adjoint` to be True when " 

228 f"`spectrum.dtype.is_complex` = True. " 

229 f"Received: {is_self_adjoint}.") 

230 is_self_adjoint = True 

231 

232 if is_square is False: 

233 raise ValueError( 

234 f"A [[nested] block] circulant operator is always square. " 

235 f"Expected argument `is_square` to be True. Received: {is_square}.") 

236 is_square = True 

237 

238 super(_BaseLinearOperatorCirculant, self).__init__( 

239 dtype=dtypes.as_dtype(input_output_dtype), 

240 is_non_singular=is_non_singular, 

241 is_self_adjoint=is_self_adjoint, 

242 is_positive_definite=is_positive_definite, 

243 is_square=is_square, 

244 parameters=parameters, 

245 name=name) 

246 

247 def _check_spectrum_and_return_tensor(self, spectrum): 

248 """Static check of spectrum. Then return `Tensor` version.""" 

249 spectrum = linear_operator_util.convert_nonref_to_tensor(spectrum, 

250 name="spectrum") 

251 

252 if spectrum.shape.ndims is not None: 

253 if spectrum.shape.ndims < self.block_depth: 

254 raise ValueError( 

255 f"Argument `spectrum` must have at least {self.block_depth} " 

256 f"dimensions. Received: {spectrum}.") 

257 return spectrum 

258 

259 @property 

260 def block_depth(self): 

261 """Depth of recursively defined circulant blocks defining this `Operator`. 

262 

263 With `A` the dense representation of this `Operator`, 

264 

265 `block_depth = 1` means `A` is symmetric circulant. For example, 

266 

267 ``` 

268 A = |w z y x| 

269 |x w z y| 

270 |y x w z| 

271 |z y x w| 

272 ``` 

273 

274 `block_depth = 2` means `A` is block symmetric circulant with symmetric 

275 circulant blocks. For example, with `W`, `X`, `Y`, `Z` symmetric circulant, 

276 

277 ``` 

278 A = |W Z Y X| 

279 |X W Z Y| 

280 |Y X W Z| 

281 |Z Y X W| 

282 ``` 

283 

284 `block_depth = 3` means `A` is block symmetric circulant with block 

285 symmetric circulant blocks. 

286 

287 Returns: 

288 Python `integer`. 

289 """ 

290 return self._block_depth 

291 

292 def block_shape_tensor(self): 

293 """Shape of the block dimensions of `self.spectrum`.""" 

294 # If spectrum.shape = [s0, s1, s2], and block_depth = 2, 

295 # block_shape = [s1, s2] 

296 return self._block_shape_tensor() 

297 

298 def _block_shape_tensor(self, spectrum_shape=None): 

299 if self.block_shape.is_fully_defined(): 

300 return linear_operator_util.shape_tensor( 

301 self.block_shape.as_list(), name="block_shape") 

302 spectrum_shape = ( 

303 array_ops.shape(self.spectrum) 

304 if spectrum_shape is None else spectrum_shape) 

305 return spectrum_shape[-self.block_depth:] 

306 

307 @property 

308 def block_shape(self): 

309 return self.spectrum.shape[-self.block_depth:] 

310 

311 @property 

312 def spectrum(self): 

313 return self._spectrum 

314 

315 def _vectorize_then_blockify(self, matrix): 

316 """Shape batch matrix to batch vector, then blockify trailing dimensions.""" 

317 # Suppose 

318 # matrix.shape = [m0, m1, m2, m3], 

319 # and matrix is a matrix because the final two dimensions are matrix dims. 

320 # self.block_depth = 2, 

321 # self.block_shape = [b0, b1] (note b0 * b1 = m2). 

322 # We will reshape matrix to 

323 # [m3, m0, m1, b0, b1]. 

324 

325 # Vectorize: Reshape to batch vector. 

326 # [m0, m1, m2, m3] --> [m3, m0, m1, m2] 

327 # This is called "vectorize" because we have taken the final two matrix dims 

328 # and turned this into a size m3 batch of vectors. 

329 vec = distribution_util.rotate_transpose(matrix, shift=1) 

330 

331 # Blockify: Blockfy trailing dimensions. 

332 # [m3, m0, m1, m2] --> [m3, m0, m1, b0, b1] 

333 if (vec.shape.is_fully_defined() and 

334 self.block_shape.is_fully_defined()): 

335 # vec_leading_shape = [m3, m0, m1], 

336 # the parts of vec that will not be blockified. 

337 vec_leading_shape = vec.shape[:-1] 

338 final_shape = vec_leading_shape.concatenate(self.block_shape) 

339 else: 

340 vec_leading_shape = array_ops.shape(vec)[:-1] 

341 final_shape = array_ops.concat( 

342 (vec_leading_shape, self.block_shape_tensor()), 0) 

343 return array_ops.reshape(vec, final_shape) 

344 

345 def _unblockify(self, x): 

346 """Flatten the trailing block dimensions.""" 

347 # Suppose 

348 # x.shape = [v0, v1, v2, v3], 

349 # self.block_depth = 2. 

350 # Then 

351 # leading shape = [v0, v1] 

352 # block shape = [v2, v3]. 

353 # We will reshape x to 

354 # [v0, v1, v2*v3]. 

355 if x.shape.is_fully_defined(): 

356 # x_shape = [v0, v1, v2, v3] 

357 x_shape = x.shape.as_list() 

358 # x_leading_shape = [v0, v1] 

359 x_leading_shape = x_shape[:-self.block_depth] 

360 # x_block_shape = [v2, v3] 

361 x_block_shape = x_shape[-self.block_depth:] 

362 # flat_shape = [v0, v1, v2*v3] 

363 flat_shape = x_leading_shape + [np.prod(x_block_shape)] 

364 else: 

365 x_shape = array_ops.shape(x) 

366 x_leading_shape = x_shape[:-self.block_depth] 

367 x_block_shape = x_shape[-self.block_depth:] 

368 flat_shape = array_ops.concat( 

369 (x_leading_shape, [math_ops.reduce_prod(x_block_shape)]), 0) 

370 return array_ops.reshape(x, flat_shape) 

371 

372 def _unblockify_then_matricize(self, vec): 

373 """Flatten the block dimensions then reshape to a batch matrix.""" 

374 # Suppose 

375 # vec.shape = [v0, v1, v2, v3], 

376 # self.block_depth = 2. 

377 # Then 

378 # leading shape = [v0, v1] 

379 # block shape = [v2, v3]. 

380 # We will reshape vec to 

381 # [v1, v2*v3, v0]. 

382 

383 # Un-blockify: Flatten block dimensions. Reshape 

384 # [v0, v1, v2, v3] --> [v0, v1, v2*v3]. 

385 vec_flat = self._unblockify(vec) 

386 

387 # Matricize: Reshape to batch matrix. 

388 # [v0, v1, v2*v3] --> [v1, v2*v3, v0], 

389 # representing a shape [v1] batch of [v2*v3, v0] matrices. 

390 matrix = distribution_util.rotate_transpose(vec_flat, shift=-1) 

391 return matrix 

392 

393 def _fft(self, x): 

394 """FFT along the last self.block_depth dimensions of x. 

395 

396 Args: 

397 x: `Tensor` with floating or complex `dtype`. 

398 Should be in the form returned by self._vectorize_then_blockify. 

399 

400 Returns: 

401 `Tensor` with `dtype` `complex64`. 

402 """ 

403 x_complex = _to_complex(x) 

404 return _FFT_OP[self.block_depth](x_complex) 

405 

406 def _ifft(self, x): 

407 """IFFT along the last self.block_depth dimensions of x. 

408 

409 Args: 

410 x: `Tensor` with floating or complex dtype. Should be in the form 

411 returned by self._vectorize_then_blockify. 

412 

413 Returns: 

414 `Tensor` with `dtype` `complex64`. 

415 """ 

416 x_complex = _to_complex(x) 

417 return _IFFT_OP[self.block_depth](x_complex) 

418 

419 def convolution_kernel(self, name="convolution_kernel"): 

420 """Convolution kernel corresponding to `self.spectrum`. 

421 

422 The `D` dimensional DFT of this kernel is the frequency domain spectrum of 

423 this operator. 

424 

425 Args: 

426 name: A name to give this `Op`. 

427 

428 Returns: 

429 `Tensor` with `dtype` `self.dtype`. 

430 """ 

431 with self._name_scope(name): # pylint: disable=not-callable 

432 h = self._ifft(_to_complex(self.spectrum)) 

433 return math_ops.cast(h, self.dtype) 

434 

435 def _shape(self): 

436 s_shape = self._spectrum.shape 

437 # Suppose spectrum.shape = [a, b, c, d] 

438 # block_depth = 2 

439 # Then: 

440 # batch_shape = [a, b] 

441 # N = c*d 

442 # and we want to return 

443 # [a, b, c*d, c*d] 

444 batch_shape = s_shape[:-self.block_depth] 

445 # trailing_dims = [c, d] 

446 trailing_dims = s_shape[-self.block_depth:] 

447 if trailing_dims.is_fully_defined(): 

448 n = np.prod(trailing_dims.as_list()) 

449 else: 

450 n = None 

451 n_x_n = tensor_shape.TensorShape([n, n]) 

452 return batch_shape.concatenate(n_x_n) 

453 

454 def _shape_tensor(self, spectrum=None): 

455 spectrum = self.spectrum if spectrum is None else spectrum 

456 # See self.shape for explanation of steps 

457 s_shape = array_ops.shape(spectrum) 

458 batch_shape = s_shape[:-self.block_depth] 

459 trailing_dims = s_shape[-self.block_depth:] 

460 n = math_ops.reduce_prod(trailing_dims) 

461 n_x_n = [n, n] 

462 return array_ops.concat((batch_shape, n_x_n), 0) 

463 

464 def assert_hermitian_spectrum(self, name="assert_hermitian_spectrum"): 

465 """Returns an `Op` that asserts this operator has Hermitian spectrum. 

466 

467 This operator corresponds to a real-valued matrix if and only if its 

468 spectrum is Hermitian. 

469 

470 Args: 

471 name: A name to give this `Op`. 

472 

473 Returns: 

474 An `Op` that asserts this operator has Hermitian spectrum. 

475 """ 

476 eps = np.finfo(self.dtype.real_dtype.as_numpy_dtype).eps 

477 with self._name_scope(name): # pylint: disable=not-callable 

478 # Assume linear accumulation of error. 

479 max_err = eps * self.domain_dimension_tensor() 

480 imag_convolution_kernel = math_ops.imag(self.convolution_kernel()) 

481 return check_ops.assert_less( 

482 math_ops.abs(imag_convolution_kernel), 

483 max_err, 

484 message="Spectrum was not Hermitian") 

485 

486 def _assert_non_singular(self): 

487 return linear_operator_util.assert_no_entries_with_modulus_zero( 

488 self.spectrum, 

489 message="Singular operator: Spectrum contained zero values.") 

490 

491 def _assert_positive_definite(self): 

492 # This operator has the action Ax = F^H D F x, 

493 # where D is the diagonal matrix with self.spectrum on the diag. Therefore, 

494 # <x, Ax> = <Fx, DFx>, 

495 # Since F is bijective, the condition for positive definite is the same as 

496 # for a diagonal matrix, i.e. real part of spectrum is positive. 

497 message = ( 

498 "Not positive definite: Real part of spectrum was not all positive.") 

499 return check_ops.assert_positive( 

500 math_ops.real(self.spectrum), message=message) 

501 

502 def _assert_self_adjoint(self): 

503 # Recall correspondence between symmetry and real transforms. See docstring 

504 return linear_operator_util.assert_zero_imag_part( 

505 self.spectrum, 

506 message=( 

507 "Not self-adjoint: The spectrum contained non-zero imaginary part." 

508 )) 

509 

510 def _broadcast_batch_dims(self, x, spectrum): 

511 """Broadcast batch dims of batch matrix `x` and spectrum.""" 

512 spectrum = tensor_conversion.convert_to_tensor_v2_with_dispatch( 

513 spectrum, name="spectrum" 

514 ) 

515 # spectrum.shape = batch_shape + block_shape 

516 # First make spectrum a batch matrix with 

517 # spectrum.shape = batch_shape + [prod(block_shape), 1] 

518 batch_shape = self._batch_shape_tensor( 

519 shape=self._shape_tensor(spectrum=spectrum)) 

520 spec_mat = array_ops.reshape( 

521 spectrum, array_ops.concat((batch_shape, [-1, 1]), axis=0)) 

522 # Second, broadcast, possibly requiring an addition of array of zeros. 

523 x, spec_mat = linear_operator_util.broadcast_matrix_batch_dims((x, 

524 spec_mat)) 

525 # Third, put the block shape back into spectrum. 

526 x_batch_shape = array_ops.shape(x)[:-2] 

527 spectrum_shape = array_ops.shape(spectrum) 

528 spectrum = array_ops.reshape( 

529 spec_mat, 

530 array_ops.concat( 

531 (x_batch_shape, 

532 self._block_shape_tensor(spectrum_shape=spectrum_shape)), 

533 axis=0)) 

534 

535 return x, spectrum 

536 

537 def _cond(self): 

538 # Regardless of whether the operator is real, it is always diagonalizable by 

539 # the Fourier basis F. I.e. A = F S F^H, with S a diagonal matrix 

540 # containing the spectrum. We then have: 

541 # A A^H = F SS^H F^H = F K F^H, 

542 # where K = diag with squared absolute values of the spectrum. 

543 # So in all cases, 

544 abs_singular_values = math_ops.abs(self._unblockify(self.spectrum)) 

545 return (math_ops.reduce_max(abs_singular_values, axis=-1) / 

546 math_ops.reduce_min(abs_singular_values, axis=-1)) 

547 

548 def _eigvals(self): 

549 return tensor_conversion.convert_to_tensor_v2_with_dispatch( 

550 self._unblockify(self.spectrum) 

551 ) 

552 

553 def _matmul(self, x, adjoint=False, adjoint_arg=False): 

554 x = linalg.adjoint(x) if adjoint_arg else x 

555 # With F the matrix of a DFT, and F^{-1}, F^H the inverse and Hermitian 

556 # transpose, one can show that F^{-1} = F^{H} is the IDFT matrix. Therefore 

557 # matmul(x) = F^{-1} diag(spectrum) F x, 

558 # = F^{H} diag(spectrum) F x, 

559 # so that 

560 # matmul(x, adjoint=True) = F^{H} diag(conj(spectrum)) F x. 

561 spectrum = _to_complex(self.spectrum) 

562 if adjoint: 

563 spectrum = math_ops.conj(spectrum) 

564 

565 x = math_ops.cast(x, spectrum.dtype) 

566 

567 x, spectrum = self._broadcast_batch_dims(x, spectrum) 

568 

569 x_vb = self._vectorize_then_blockify(x) 

570 fft_x_vb = self._fft(x_vb) 

571 block_vector_result = self._ifft(spectrum * fft_x_vb) 

572 y = self._unblockify_then_matricize(block_vector_result) 

573 

574 return math_ops.cast(y, self.dtype) 

575 

576 def _determinant(self): 

577 axis = [-(i + 1) for i in range(self.block_depth)] 

578 det = math_ops.reduce_prod(self.spectrum, axis=axis) 

579 return math_ops.cast(det, self.dtype) 

580 

581 def _log_abs_determinant(self): 

582 axis = [-(i + 1) for i in range(self.block_depth)] 

583 lad = math_ops.reduce_sum( 

584 math_ops.log(math_ops.abs(self.spectrum)), axis=axis) 

585 return math_ops.cast(lad, self.dtype) 

586 

587 def _solve(self, rhs, adjoint=False, adjoint_arg=False): 

588 rhs = linalg.adjoint(rhs) if adjoint_arg else rhs 

589 spectrum = _to_complex(self.spectrum) 

590 if adjoint: 

591 spectrum = math_ops.conj(spectrum) 

592 

593 rhs, spectrum = self._broadcast_batch_dims(rhs, spectrum) 

594 

595 rhs_vb = self._vectorize_then_blockify(rhs) 

596 fft_rhs_vb = self._fft(rhs_vb) 

597 solution_vb = self._ifft(fft_rhs_vb / spectrum) 

598 x = self._unblockify_then_matricize(solution_vb) 

599 return math_ops.cast(x, self.dtype) 

600 

601 def _diag_part(self): 

602 # Get ones in shape of diag, which is [B1,...,Bb, N] 

603 # Also get the size of the diag, "N". 

604 if self.shape.is_fully_defined(): 

605 diag_shape = self.shape[:-1] 

606 diag_size = self.domain_dimension.value 

607 else: 

608 diag_shape = self.shape_tensor()[:-1] 

609 diag_size = self.domain_dimension_tensor() 

610 ones_diag = array_ops.ones(diag_shape, dtype=self.dtype) 

611 

612 # As proved in comments in self._trace, the value on the diag is constant, 

613 # repeated N times. This value is the trace divided by N. 

614 

615 # The handling of self.shape = (0, 0) is tricky, and is the reason we choose 

616 # to compute trace and use that to compute diag_part, rather than computing 

617 # the value on the diagonal ("diag_value") directly. Both result in a 0/0, 

618 # but in different places, and the current method gives the right result in 

619 # the end. 

620 

621 # Here, if self.shape = (0, 0), then self.trace() = 0., and then 

622 # diag_value = 0. / 0. = NaN. 

623 diag_value = self.trace() / math_ops.cast(diag_size, self.dtype) 

624 

625 # If self.shape = (0, 0), then ones_diag = [] (empty tensor), and then 

626 # the following line is NaN * [] = [], as needed. 

627 return diag_value[..., array_ops.newaxis] * ones_diag 

628 

629 def _trace(self): 

630 # The diagonal of the [[nested] block] circulant operator is the mean of 

631 # the spectrum. 

632 # Proof: For the [0,...,0] element, this follows from the IDFT formula. 

633 # Then the result follows since all diagonal elements are the same. 

634 

635 # Therefore, the trace is the sum of the spectrum. 

636 

637 # Get shape of diag along with the axis over which to reduce the spectrum. 

638 # We will reduce the spectrum over all block indices. 

639 if self.spectrum.shape.is_fully_defined(): 

640 spec_rank = self.spectrum.shape.ndims 

641 axis = np.arange(spec_rank - self.block_depth, spec_rank, dtype=np.int32) 

642 else: 

643 spec_rank = array_ops.rank(self.spectrum) 

644 axis = math_ops.range(spec_rank - self.block_depth, spec_rank) 

645 

646 # Real diag part "re_d". 

647 # Suppose spectrum.shape = [B1,...,Bb, N1, N2] 

648 # self.shape = [B1,...,Bb, N, N], with N1 * N2 = N. 

649 # re_d_value.shape = [B1,...,Bb] 

650 re_d_value = math_ops.reduce_sum(math_ops.real(self.spectrum), axis=axis) 

651 

652 if not self.dtype.is_complex: 

653 return math_ops.cast(re_d_value, self.dtype) 

654 

655 # Imaginary part, "im_d". 

656 if self.is_self_adjoint: 

657 im_d_value = array_ops.zeros_like(re_d_value) 

658 else: 

659 im_d_value = math_ops.reduce_sum(math_ops.imag(self.spectrum), axis=axis) 

660 

661 return math_ops.cast(math_ops.complex(re_d_value, im_d_value), self.dtype) 

662 

663 @property 

664 def _composite_tensor_fields(self): 

665 return ("spectrum", "input_output_dtype") 

666 

667 @property 

668 def _experimental_parameter_ndims_to_matrix_ndims(self): 

669 return {"spectrum": self.block_depth} 

670 

671 

672@tf_export("linalg.LinearOperatorCirculant") 

673@linear_operator.make_composite_tensor 

674class LinearOperatorCirculant(_BaseLinearOperatorCirculant): 

675 """`LinearOperator` acting like a circulant matrix. 

676 

677 This operator acts like a circulant matrix `A` with 

678 shape `[B1,...,Bb, N, N]` for some `b >= 0`. The first `b` indices index a 

679 batch member. For every batch index `(i1,...,ib)`, `A[i1,...,ib, : :]` is 

680 an `N x N` matrix. This matrix `A` is not materialized, but for 

681 purposes of broadcasting this shape will be relevant. 

682 

683 #### Description in terms of circulant matrices 

684 

685 Circulant means the entries of `A` are generated by a single vector, the 

686 convolution kernel `h`: `A_{mn} := h_{m-n mod N}`. With `h = [w, x, y, z]`, 

687 

688 ``` 

689 A = |w z y x| 

690 |x w z y| 

691 |y x w z| 

692 |z y x w| 

693 ``` 

694 

695 This means that the result of matrix multiplication `v = Au` has `Lth` column 

696 given circular convolution between `h` with the `Lth` column of `u`. 

697 

698 #### Description in terms of the frequency spectrum 

699 

700 There is an equivalent description in terms of the [batch] spectrum `H` and 

701 Fourier transforms. Here we consider `A.shape = [N, N]` and ignore batch 

702 dimensions. Define the discrete Fourier transform (DFT) and its inverse by 

703 

704 ``` 

705 DFT[ h[n] ] = H[k] := sum_{n = 0}^{N - 1} h_n e^{-i 2pi k n / N} 

706 IDFT[ H[k] ] = h[n] = N^{-1} sum_{k = 0}^{N - 1} H_k e^{i 2pi k n / N} 

707 ``` 

708 

709 From these definitions, we see that 

710 

711 ``` 

712 H[0] = sum_{n = 0}^{N - 1} h_n 

713 H[1] = "the first positive frequency" 

714 H[N - 1] = "the first negative frequency" 

715 ``` 

716 

717 Loosely speaking, with `*` element-wise multiplication, matrix multiplication 

718 is equal to the action of a Fourier multiplier: `A u = IDFT[ H * DFT[u] ]`. 

719 Precisely speaking, given `[N, R]` matrix `u`, let `DFT[u]` be the `[N, R]` 

720 matrix with `rth` column equal to the DFT of the `rth` column of `u`. 

721 Define the `IDFT` similarly. 

722 Matrix multiplication may be expressed columnwise: 

723 

724 ```(A u)_r = IDFT[ H * (DFT[u])_r ]``` 

725 

726 #### Operator properties deduced from the spectrum. 

727 

728 Letting `U` be the `kth` Euclidean basis vector, and `U = IDFT[u]`. 

729 The above formulas show that`A U = H_k * U`. We conclude that the elements 

730 of `H` are the eigenvalues of this operator. Therefore 

731 

732 * This operator is positive definite if and only if `Real{H} > 0`. 

733 

734 A general property of Fourier transforms is the correspondence between 

735 Hermitian functions and real valued transforms. 

736 

737 Suppose `H.shape = [B1,...,Bb, N]`. We say that `H` is a Hermitian spectrum 

738 if, with `%` meaning modulus division, 

739 

740 ```H[..., n % N] = ComplexConjugate[ H[..., (-n) % N] ]``` 

741 

742 * This operator corresponds to a real matrix if and only if `H` is Hermitian. 

743 * This operator is self-adjoint if and only if `H` is real. 

744 

745 See e.g. "Discrete-Time Signal Processing", Oppenheim and Schafer. 

746 

747 #### Example of a self-adjoint positive definite operator 

748 

749 ```python 

750 # spectrum is real ==> operator is self-adjoint 

751 # spectrum is positive ==> operator is positive definite 

752 spectrum = [6., 4, 2] 

753 

754 operator = LinearOperatorCirculant(spectrum) 

755 

756 # IFFT[spectrum] 

757 operator.convolution_kernel() 

758 ==> [4 + 0j, 1 + 0.58j, 1 - 0.58j] 

759 

760 operator.to_dense() 

761 ==> [[4 + 0.0j, 1 - 0.6j, 1 + 0.6j], 

762 [1 + 0.6j, 4 + 0.0j, 1 - 0.6j], 

763 [1 - 0.6j, 1 + 0.6j, 4 + 0.0j]] 

764 ``` 

765 

766 #### Example of defining in terms of a real convolution kernel 

767 

768 ```python 

769 # convolution_kernel is real ==> spectrum is Hermitian. 

770 convolution_kernel = [1., 2., 1.]] 

771 spectrum = tf.signal.fft(tf.cast(convolution_kernel, tf.complex64)) 

772 

773 # spectrum is Hermitian ==> operator is real. 

774 # spectrum is shape [3] ==> operator is shape [3, 3] 

775 # We force the input/output type to be real, which allows this to operate 

776 # like a real matrix. 

777 operator = LinearOperatorCirculant(spectrum, input_output_dtype=tf.float32) 

778 

779 operator.to_dense() 

780 ==> [[ 1, 1, 2], 

781 [ 2, 1, 1], 

782 [ 1, 2, 1]] 

783 ``` 

784 

785 #### Example of Hermitian spectrum 

786 

787 ```python 

788 # spectrum is shape [3] ==> operator is shape [3, 3] 

789 # spectrum is Hermitian ==> operator is real. 

790 spectrum = [1, 1j, -1j] 

791 

792 operator = LinearOperatorCirculant(spectrum) 

793 

794 operator.to_dense() 

795 ==> [[ 0.33 + 0j, 0.91 + 0j, -0.24 + 0j], 

796 [-0.24 + 0j, 0.33 + 0j, 0.91 + 0j], 

797 [ 0.91 + 0j, -0.24 + 0j, 0.33 + 0j] 

798 ``` 

799 

800 #### Example of forcing real `dtype` when spectrum is Hermitian 

801 

802 ```python 

803 # spectrum is shape [4] ==> operator is shape [4, 4] 

804 # spectrum is real ==> operator is self-adjoint 

805 # spectrum is Hermitian ==> operator is real 

806 # spectrum has positive real part ==> operator is positive-definite. 

807 spectrum = [6., 4, 2, 4] 

808 

809 # Force the input dtype to be float32. 

810 # Cast the output to float32. This is fine because the operator will be 

811 # real due to Hermitian spectrum. 

812 operator = LinearOperatorCirculant(spectrum, input_output_dtype=tf.float32) 

813 

814 operator.shape 

815 ==> [4, 4] 

816 

817 operator.to_dense() 

818 ==> [[4, 1, 0, 1], 

819 [1, 4, 1, 0], 

820 [0, 1, 4, 1], 

821 [1, 0, 1, 4]] 

822 

823 # convolution_kernel = tf.signal.ifft(spectrum) 

824 operator.convolution_kernel() 

825 ==> [4, 1, 0, 1] 

826 ``` 

827 

828 #### Performance 

829 

830 Suppose `operator` is a `LinearOperatorCirculant` of shape `[N, N]`, 

831 and `x.shape = [N, R]`. Then 

832 

833 * `operator.matmul(x)` is `O(R*N*Log[N])` 

834 * `operator.solve(x)` is `O(R*N*Log[N])` 

835 * `operator.determinant()` involves a size `N` `reduce_prod`. 

836 

837 If instead `operator` and `x` have shape `[B1,...,Bb, N, N]` and 

838 `[B1,...,Bb, N, R]`, every operation increases in complexity by `B1*...*Bb`. 

839 

840 #### Matrix property hints 

841 

842 This `LinearOperator` is initialized with boolean flags of the form `is_X`, 

843 for `X = non_singular, self_adjoint, positive_definite, square`. 

844 These have the following meaning: 

845 

846 * If `is_X == True`, callers should expect the operator to have the 

847 property `X`. This is a promise that should be fulfilled, but is *not* a 

848 runtime assert. For example, finite floating point precision may result 

849 in these promises being violated. 

850 * If `is_X == False`, callers should expect the operator to not have `X`. 

851 * If `is_X == None` (the default), callers should have no expectation either 

852 way. 

853 

854 References: 

855 Toeplitz and Circulant Matrices - A Review: 

856 [Gray, 2006](https://www.nowpublishers.com/article/Details/CIT-006) 

857 ([pdf](https://ee.stanford.edu/~gray/toeplitz.pdf)) 

858 """ 

859 

860 def __init__(self, 

861 spectrum, 

862 input_output_dtype=dtypes.complex64, 

863 is_non_singular=None, 

864 is_self_adjoint=None, 

865 is_positive_definite=None, 

866 is_square=True, 

867 name="LinearOperatorCirculant"): 

868 r"""Initialize an `LinearOperatorCirculant`. 

869 

870 This `LinearOperator` is initialized to have shape `[B1,...,Bb, N, N]` 

871 by providing `spectrum`, a `[B1,...,Bb, N]` `Tensor`. 

872 

873 If `input_output_dtype = DTYPE`: 

874 

875 * Arguments to methods such as `matmul` or `solve` must be `DTYPE`. 

876 * Values returned by all methods, such as `matmul` or `determinant` will be 

877 cast to `DTYPE`. 

878 

879 Note that if the spectrum is not Hermitian, then this operator corresponds 

880 to a complex matrix with non-zero imaginary part. In this case, setting 

881 `input_output_dtype` to a real type will forcibly cast the output to be 

882 real, resulting in incorrect results! 

883 

884 If on the other hand the spectrum is Hermitian, then this operator 

885 corresponds to a real-valued matrix, and setting `input_output_dtype` to 

886 a real type is fine. 

887 

888 Args: 

889 spectrum: Shape `[B1,...,Bb, N]` `Tensor`. Allowed dtypes: `float16`, 

890 `float32`, `float64`, `complex64`, `complex128`. Type can be different 

891 than `input_output_dtype` 

892 input_output_dtype: `dtype` for input/output. 

893 is_non_singular: Expect that this operator is non-singular. 

894 is_self_adjoint: Expect that this operator is equal to its hermitian 

895 transpose. If `spectrum` is real, this will always be true. 

896 is_positive_definite: Expect that this operator is positive definite, 

897 meaning the quadratic form `x^H A x` has positive real part for all 

898 nonzero `x`. Note that we do not require the operator to be 

899 self-adjoint to be positive-definite. See: 

900 https://en.wikipedia.org/wiki/Positive-definite_matrix\ 

901 #Extension_for_non_symmetric_matrices 

902 is_square: Expect that this operator acts like square [batch] matrices. 

903 name: A name to prepend to all ops created by this class. 

904 """ 

905 parameters = dict( 

906 spectrum=spectrum, 

907 input_output_dtype=input_output_dtype, 

908 is_non_singular=is_non_singular, 

909 is_self_adjoint=is_self_adjoint, 

910 is_positive_definite=is_positive_definite, 

911 is_square=is_square, 

912 name=name 

913 ) 

914 super(LinearOperatorCirculant, self).__init__( 

915 spectrum, 

916 block_depth=1, 

917 input_output_dtype=input_output_dtype, 

918 is_non_singular=is_non_singular, 

919 is_self_adjoint=is_self_adjoint, 

920 is_positive_definite=is_positive_definite, 

921 is_square=is_square, 

922 parameters=parameters, 

923 name=name) 

924 

925 

926@tf_export("linalg.LinearOperatorCirculant2D") 

927@linear_operator.make_composite_tensor 

928class LinearOperatorCirculant2D(_BaseLinearOperatorCirculant): 

929 """`LinearOperator` acting like a block circulant matrix. 

930 

931 This operator acts like a block circulant matrix `A` with 

932 shape `[B1,...,Bb, N, N]` for some `b >= 0`. The first `b` indices index a 

933 batch member. For every batch index `(i1,...,ib)`, `A[i1,...,ib, : :]` is 

934 an `N x N` matrix. This matrix `A` is not materialized, but for 

935 purposes of broadcasting this shape will be relevant. 

936 

937 #### Description in terms of block circulant matrices 

938 

939 If `A` is block circulant, with block sizes `N0, N1` (`N0 * N1 = N`): 

940 `A` has a block circulant structure, composed of `N0 x N0` blocks, with each 

941 block an `N1 x N1` circulant matrix. 

942 

943 For example, with `W`, `X`, `Y`, `Z` each circulant, 

944 

945 ``` 

946 A = |W Z Y X| 

947 |X W Z Y| 

948 |Y X W Z| 

949 |Z Y X W| 

950 ``` 

951 

952 Note that `A` itself will not in general be circulant. 

953 

954 #### Description in terms of the frequency spectrum 

955 

956 There is an equivalent description in terms of the [batch] spectrum `H` and 

957 Fourier transforms. Here we consider `A.shape = [N, N]` and ignore batch 

958 dimensions. 

959 

960 If `H.shape = [N0, N1]`, (`N0 * N1 = N`): 

961 Loosely speaking, matrix multiplication is equal to the action of a 

962 Fourier multiplier: `A u = IDFT2[ H DFT2[u] ]`. 

963 Precisely speaking, given `[N, R]` matrix `u`, let `DFT2[u]` be the 

964 `[N0, N1, R]` `Tensor` defined by re-shaping `u` to `[N0, N1, R]` and taking 

965 a two dimensional DFT across the first two dimensions. Let `IDFT2` be the 

966 inverse of `DFT2`. Matrix multiplication may be expressed columnwise: 

967 

968 ```(A u)_r = IDFT2[ H * (DFT2[u])_r ]``` 

969 

970 #### Operator properties deduced from the spectrum. 

971 

972 * This operator is positive definite if and only if `Real{H} > 0`. 

973 

974 A general property of Fourier transforms is the correspondence between 

975 Hermitian functions and real valued transforms. 

976 

977 Suppose `H.shape = [B1,...,Bb, N0, N1]`, we say that `H` is a Hermitian 

978 spectrum if, with `%` indicating modulus division, 

979 

980 ``` 

981 H[..., n0 % N0, n1 % N1] = ComplexConjugate[ H[..., (-n0) % N0, (-n1) % N1 ]. 

982 ``` 

983 

984 * This operator corresponds to a real matrix if and only if `H` is Hermitian. 

985 * This operator is self-adjoint if and only if `H` is real. 

986 

987 See e.g. "Discrete-Time Signal Processing", Oppenheim and Schafer. 

988 

989 ### Example of a self-adjoint positive definite operator 

990 

991 ```python 

992 # spectrum is real ==> operator is self-adjoint 

993 # spectrum is positive ==> operator is positive definite 

994 spectrum = [[1., 2., 3.], 

995 [4., 5., 6.], 

996 [7., 8., 9.]] 

997 

998 operator = LinearOperatorCirculant2D(spectrum) 

999 

1000 # IFFT[spectrum] 

1001 operator.convolution_kernel() 

1002 ==> [[5.0+0.0j, -0.5-.3j, -0.5+.3j], 

1003 [-1.5-.9j, 0, 0], 

1004 [-1.5+.9j, 0, 0]] 

1005 

1006 operator.to_dense() 

1007 ==> Complex self adjoint 9 x 9 matrix. 

1008 ``` 

1009 

1010 #### Example of defining in terms of a real convolution kernel, 

1011 

1012 ```python 

1013 # convolution_kernel is real ==> spectrum is Hermitian.