1# Copyright 2016 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` acting like a diagonal matrix.""" 16 17from __future__ import absolute_import 18from __future__ import division 19from __future__ import print_function 20 21from tensorflow.python.framework import ops 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.linalg import linalg_impl as linalg 26from tensorflow.python.ops.linalg import linear_operator 27from tensorflow.python.ops.linalg import linear_operator_util 28from tensorflow.python.util.tf_export import tf_export 29 30__all__ = ["LinearOperatorDiag",] 31 32 33@tf_export("linalg.LinearOperatorDiag") 34@linear_operator.make_composite_tensor 35class LinearOperatorDiag(linear_operator.LinearOperator): 36 """`LinearOperator` acting like a [batch] square diagonal matrix. 37 38 This operator acts like a [batch] diagonal matrix `A` with shape 39 `[B1,...,Bb, N, N]` for some `b >= 0`. The first `b` indices index a 40 batch member. For every batch index `(i1,...,ib)`, `A[i1,...,ib, : :]` is 41 an `N x N` matrix. This matrix `A` is not materialized, but for 42 purposes of broadcasting this shape will be relevant. 43 44 `LinearOperatorDiag` is initialized with a (batch) vector. 45 46 ```python 47 # Create a 2 x 2 diagonal linear operator. 48 diag = [1., -1.] 49 operator = LinearOperatorDiag(diag) 50 51 operator.to_dense() 52 ==> [[1., 0.] 53 [0., -1.]] 54 55 operator.shape 56 ==> [2, 2] 57 58 operator.log_abs_determinant() 59 ==> scalar Tensor 60 61 x = ... Shape [2, 4] Tensor 62 operator.matmul(x) 63 ==> Shape [2, 4] Tensor 64 65 # Create a [2, 3] batch of 4 x 4 linear operators. 66 diag = tf.random.normal(shape=[2, 3, 4]) 67 operator = LinearOperatorDiag(diag) 68 69 # Create a shape [2, 1, 4, 2] vector. Note that this shape is compatible 70 # since the batch dimensions, [2, 1], are broadcast to 71 # operator.batch_shape = [2, 3]. 72 y = tf.random.normal(shape=[2, 1, 4, 2]) 73 x = operator.solve(y) 74 ==> operator.matmul(x) = y 75 ``` 76 77 #### Shape compatibility 78 79 This operator acts on [batch] matrix with compatible shape. 80 `x` is a batch matrix with compatible shape for `matmul` and `solve` if 81 82 ``` 83 operator.shape = [B1,...,Bb] + [N, N], with b >= 0 84 x.shape = [C1,...,Cc] + [N, R], 85 and [C1,...,Cc] broadcasts with [B1,...,Bb] to [D1,...,Dd] 86 ``` 87 88 #### Performance 89 90 Suppose `operator` is a `LinearOperatorDiag` of shape `[N, N]`, 91 and `x.shape = [N, R]`. Then 92 93 * `operator.matmul(x)` involves `N * R` multiplications. 94 * `operator.solve(x)` involves `N` divisions and `N * R` multiplications. 95 * `operator.determinant()` involves a size `N` `reduce_prod`. 96 97 If instead `operator` and `x` have shape `[B1,...,Bb, N, N]` and 98 `[B1,...,Bb, N, R]`, every operation increases in complexity by `B1*...*Bb`. 99 100 #### Matrix property hints 101 102 This `LinearOperator` is initialized with boolean flags of the form `is_X`, 103 for `X = non_singular, self_adjoint, positive_definite, square`. 104 These have the following meaning: 105 106 * If `is_X == True`, callers should expect the operator to have the 107 property `X`. This is a promise that should be fulfilled, but is *not* a 108 runtime assert. For example, finite floating point precision may result 109 in these promises being violated. 110 * If `is_X == False`, callers should expect the operator to not have `X`. 111 * If `is_X == None` (the default), callers should have no expectation either 112 way. 113 """ 114 115 def __init__(self, 116 diag, 117 is_non_singular=None, 118 is_self_adjoint=None, 119 is_positive_definite=None, 120 is_square=None, 121 name="LinearOperatorDiag"): 122 r"""Initialize a `LinearOperatorDiag`. 123 124 Args: 125 diag: Shape `[B1,...,Bb, N]` `Tensor` with `b >= 0` `N >= 0`. 126 The diagonal of the operator. Allowed dtypes: `float16`, `float32`, 127 `float64`, `complex64`, `complex128`. 128 is_non_singular: Expect that this operator is non-singular. 129 is_self_adjoint: Expect that this operator is equal to its hermitian 130 transpose. If `diag.dtype` is real, this is auto-set to `True`. 131 is_positive_definite: Expect that this operator is positive definite, 132 meaning the quadratic form `x^H A x` has positive real part for all 133 nonzero `x`. Note that we do not require the operator to be 134 self-adjoint to be positive-definite. See: 135 https://en.wikipedia.org/wiki/Positive-definite_matrix#Extension_for_non-symmetric_matrices 136 is_square: Expect that this operator acts like square [batch] matrices. 137 name: A name for this `LinearOperator`. 138 139 Raises: 140 TypeError: If `diag.dtype` is not an allowed type. 141 ValueError: If `diag.dtype` is real, and `is_self_adjoint` is not `True`. 142 """ 143 parameters = dict( 144 diag=diag, 145 is_non_singular=is_non_singular, 146 is_self_adjoint=is_self_adjoint, 147 is_positive_definite=is_positive_definite, 148 is_square=is_square, 149 name=name 150 ) 151 152 with ops.name_scope(name, values=[diag]): 153 self._diag = linear_operator_util.convert_nonref_to_tensor( 154 diag, name="diag") 155 self._check_diag(self._diag) 156 157 # Check and auto-set hints. 158 if not self._diag.dtype.is_complex: 159 if is_self_adjoint is False: 160 raise ValueError("A real diagonal operator is always self adjoint.") 161 else: 162 is_self_adjoint = True 163 164 if is_square is False: 165 raise ValueError("Only square diagonal operators currently supported.") 166 is_square = True 167 168 super(LinearOperatorDiag, self).__init__( 169 dtype=self._diag.dtype, 170 is_non_singular=is_non_singular, 171 is_self_adjoint=is_self_adjoint, 172 is_positive_definite=is_positive_definite, 173 is_square=is_square, 174 parameters=parameters, 175 name=name) 176 # TODO(b/143910018) Remove graph_parents in V3. 177 self._set_graph_parents([self._diag]) 178 179 def _check_diag(self, diag): 180 """Static check of diag.""" 181 if diag.shape.ndims is not None and diag.shape.ndims < 1: 182 raise ValueError("Argument diag must have at least 1 dimension. " 183 "Found: %s" % diag) 184 185 def _shape(self): 186 # If d_shape = [5, 3], we return [5, 3, 3]. 187 d_shape = self._diag.shape 188 return d_shape.concatenate(d_shape[-1:]) 189 190 def _shape_tensor(self): 191 d_shape = array_ops.shape(self._diag) 192 k = d_shape[-1] 193 return array_ops.concat((d_shape, [k]), 0) 194 195 @property 196 def diag(self): 197 return self._diag 198 199 def _assert_non_singular(self): 200 return linear_operator_util.assert_no_entries_with_modulus_zero( 201 self._diag, 202 message="Singular operator: Diagonal contained zero values.") 203 204 def _assert_positive_definite(self): 205 if self.dtype.is_complex: 206 message = ( 207 "Diagonal operator had diagonal entries with non-positive real part, " 208 "thus was not positive definite.") 209 else: 210 message = ( 211 "Real diagonal operator had non-positive diagonal entries, " 212 "thus was not positive definite.") 213 214 return check_ops.assert_positive( 215 math_ops.real(self._diag), 216 message=message) 217 218 def _assert_self_adjoint(self): 219 return linear_operator_util.assert_zero_imag_part( 220 self._diag, 221 message=( 222 "This diagonal operator contained non-zero imaginary values. " 223 " Thus it was not self-adjoint.")) 224 225 def _matmul(self, x, adjoint=False, adjoint_arg=False): 226 diag_term = math_ops.conj(self._diag) if adjoint else self._diag 227 x = linalg.adjoint(x) if adjoint_arg else x 228 diag_mat = array_ops.expand_dims(diag_term, -1) 229 return diag_mat * x 230 231 def _matvec(self, x, adjoint=False): 232 diag_term = math_ops.conj(self._diag) if adjoint else self._diag 233 return diag_term * x 234 235 def _determinant(self): 236 return math_ops.reduce_prod(self._diag, axis=[-1]) 237 238 def _log_abs_determinant(self): 239 log_det = math_ops.reduce_sum( 240 math_ops.log(math_ops.abs(self._diag)), axis=[-1]) 241 if self.dtype.is_complex: 242 log_det = math_ops.cast(log_det, dtype=self.dtype) 243 return log_det 244 245 def _solve(self, rhs, adjoint=False, adjoint_arg=False): 246 diag_term = math_ops.conj(self._diag) if adjoint else self._diag 247 rhs = linalg.adjoint(rhs) if adjoint_arg else rhs 248 inv_diag_mat = array_ops.expand_dims(1. / diag_term, -1) 249 return rhs * inv_diag_mat 250 251 def _to_dense(self): 252 return array_ops.matrix_diag(self._diag) 253 254 def _diag_part(self): 255 return self.diag 256 257 def _add_to_tensor(self, x): 258 x_diag = array_ops.matrix_diag_part(x) 259 new_diag = self._diag + x_diag 260 return array_ops.matrix_set_diag(x, new_diag) 261 262 def _eigvals(self): 263 return ops.convert_to_tensor_v2_with_dispatch(self.diag) 264 265 def _cond(self): 266 abs_diag = math_ops.abs(self.diag) 267 return (math_ops.reduce_max(abs_diag, axis=-1) / 268 math_ops.reduce_min(abs_diag, axis=-1)) 269 270 @property 271 def _composite_tensor_fields(self): 272 return ("diag",) 273