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"""Tests for tensorflow.ops.math_ops.matrix_inverse.""" 16 17import numpy as np 18 19from tensorflow.python.client import session 20from tensorflow.python.eager import context 21from tensorflow.python.framework import constant_op 22from tensorflow.python.framework import errors_impl 23from tensorflow.python.framework import ops 24from tensorflow.python.framework import test_util 25from tensorflow.python.ops import array_ops 26from tensorflow.python.ops import control_flow_ops 27from tensorflow.python.ops import gradient_checker 28from tensorflow.python.ops import gradient_checker_v2 29from tensorflow.python.ops import gradients_impl 30from tensorflow.python.ops import linalg_ops 31from tensorflow.python.ops import math_ops 32from tensorflow.python.ops import stateless_random_ops 33from tensorflow.python.ops import variables 34from tensorflow.python.platform import benchmark 35from tensorflow.python.platform import test 36 37 38def _AddTest(test_class, op_name, testcase_name, fn): 39 test_name = "_".join(["test", op_name, testcase_name]) 40 if hasattr(test_class, test_name): 41 raise RuntimeError("Test %s defined more than once" % test_name) 42 setattr(test_class, test_name, fn) 43 44 45class SvdOpTest(test.TestCase): 46 47 @test_util.run_in_graph_and_eager_modes(use_gpu=True) 48 def testWrongDimensions(self): 49 # The input to svd should be a tensor of at least rank 2. 50 scalar = constant_op.constant(1.) 51 with self.assertRaisesRegex((ValueError, errors_impl.InvalidArgumentError), 52 "rank.* 2.*0"): 53 linalg_ops.svd(scalar) 54 vector = constant_op.constant([1., 2.]) 55 with self.assertRaisesRegex((ValueError, errors_impl.InvalidArgumentError), 56 "rank.* 2.*1"): 57 linalg_ops.svd(vector) 58 59 @test_util.run_in_graph_and_eager_modes(use_gpu=True) 60 def testThrowDeterminismError(self): 61 shape = [6, 5] 62 seed = [42, 24] 63 matrix1 = stateless_random_ops.stateless_random_normal(shape, seed) 64 with test_util.deterministic_ops(): 65 if test_util.is_gpu_available(cuda_only=True): 66 with self.assertRaisesRegex( 67 errors_impl.UnimplementedError, "Determinism is not yet supported " 68 "for Svd."): 69 self.evaluate(linalg_ops.svd(matrix1)) 70 71 @test_util.run_in_graph_and_eager_modes(use_gpu=True) 72 def DISABLED_testBadInputs(self): 73 # TODO(b/185822300): re-enable after the bug is fixed in CUDA-11.x 74 # The input to svd should be a tensor of at least rank 2. 75 for bad_val in [np.nan, np.inf]: 76 matrix = np.array([[1, bad_val], [0, 1]]) 77 s, u, v = linalg_ops.svd(matrix, compute_uv=True) 78 s, u, v = self.evaluate([s, u, v]) 79 for i in range(2): 80 self.assertTrue(np.isnan(s[i])) 81 for j in range(2): 82 self.assertTrue(np.isnan(u[i, j])) 83 self.assertTrue(np.isnan(v[i, j])) 84 85 @test_util.run_in_graph_and_eager_modes(use_gpu=True) 86 def testExecuteMultipleWithoutError(self): 87 all_ops = [] 88 shape = [6, 5] 89 seed = [42, 24] 90 for compute_uv_ in True, False: 91 for full_matrices_ in True, False: 92 matrix1 = stateless_random_ops.stateless_random_normal(shape, seed) 93 matrix2 = stateless_random_ops.stateless_random_normal(shape, seed) 94 self.assertAllEqual(matrix1, matrix2) 95 if compute_uv_: 96 s1, u1, v1 = linalg_ops.svd( 97 matrix1, compute_uv=compute_uv_, full_matrices=full_matrices_) 98 s2, u2, v2 = linalg_ops.svd( 99 matrix2, compute_uv=compute_uv_, full_matrices=full_matrices_) 100 all_ops += [s1, s2, u1, u2, v1, v2] 101 else: 102 s1 = linalg_ops.svd( 103 matrix1, compute_uv=compute_uv_, full_matrices=full_matrices_) 104 s2 = linalg_ops.svd( 105 matrix2, compute_uv=compute_uv_, full_matrices=full_matrices_) 106 all_ops += [s1, s2] 107 val = self.evaluate(all_ops) 108 for i in range(0, len(val), 2): 109 self.assertAllEqual(val[i], val[i + 1]) 110 111 @test_util.run_in_graph_and_eager_modes(use_gpu=True) 112 def testEmptyBatches(self): 113 matrices = constant_op.constant(1.0, shape=[0, 2, 2]) 114 s, u, v = self.evaluate(linalg_ops.svd(matrices)) 115 self.assertAllEqual(s, np.zeros([0, 2])) 116 self.assertAllEqual(u, np.zeros([0, 2, 2])) 117 self.assertAllEqual(v, np.zeros([0, 2, 2])) 118 119 120def _GetSvdOpTest(dtype_, shape_, use_static_shape_, compute_uv_, 121 full_matrices_): 122 123 def CompareSingularValues(self, x, y, tol): 124 atol = (x[0] + y[0]) * tol if len(x) else tol 125 self.assertAllClose(x, y, atol=atol) 126 127 def CompareSingularVectors(self, x, y, rank, tol): 128 # We only compare the first 'rank' singular vectors since the 129 # remainder form an arbitrary orthonormal basis for the 130 # (row- or column-) null space, whose exact value depends on 131 # implementation details. Notice that since we check that the 132 # matrices of singular vectors are unitary elsewhere, we do 133 # implicitly test that the trailing vectors of x and y span the 134 # same space. 135 x = x[..., 0:rank] 136 y = y[..., 0:rank] 137 # Singular vectors are only unique up to sign (complex phase factor for 138 # complex matrices), so we normalize the sign first. 139 sum_of_ratios = np.sum(np.divide(y, x), -2, keepdims=True) 140 phases = np.divide(sum_of_ratios, np.abs(sum_of_ratios)) 141 x *= phases 142 self.assertAllClose(x, y, atol=2 * tol) 143 144 def CheckApproximation(self, a, u, s, v, full_matrices_, tol): 145 # Tests that a ~= u*diag(s)*transpose(v). 146 batch_shape = a.shape[:-2] 147 m = a.shape[-2] 148 n = a.shape[-1] 149 diag_s = math_ops.cast(array_ops.matrix_diag(s), dtype=dtype_) 150 if full_matrices_: 151 if m > n: 152 zeros = array_ops.zeros(batch_shape + (m - n, n), dtype=dtype_) 153 diag_s = array_ops.concat([diag_s, zeros], a.ndim - 2) 154 elif n > m: 155 zeros = array_ops.zeros(batch_shape + (m, n - m), dtype=dtype_) 156 diag_s = array_ops.concat([diag_s, zeros], a.ndim - 1) 157 a_recon = math_ops.matmul(u, diag_s) 158 a_recon = math_ops.matmul(a_recon, v, adjoint_b=True) 159 self.assertAllClose(a_recon, a, rtol=tol, atol=tol) 160 161 def CheckUnitary(self, x, tol): 162 # Tests that x[...,:,:]^H * x[...,:,:] is close to the identity. 163 xx = math_ops.matmul(x, x, adjoint_a=True) 164 identity = array_ops.matrix_band_part(array_ops.ones_like(xx), 0, 0) 165 self.assertAllClose(identity, xx, atol=tol) 166 167 @test_util.run_in_graph_and_eager_modes(use_gpu=True) 168 def Test(self): 169 if not use_static_shape_ and context.executing_eagerly(): 170 return 171 is_complex = dtype_ in (np.complex64, np.complex128) 172 is_single = dtype_ in (np.float32, np.complex64) 173 tol = 3e-4 if is_single else 1e-12 174 if test.is_gpu_available(): 175 # The gpu version returns results that are much less accurate. 176 tol *= 200 177 np.random.seed(42) 178 x_np = np.random.uniform( 179 low=-1.0, high=1.0, size=np.prod(shape_)).reshape(shape_).astype(dtype_) 180 if is_complex: 181 x_np += 1j * np.random.uniform( 182 low=-1.0, high=1.0, 183 size=np.prod(shape_)).reshape(shape_).astype(dtype_) 184 185 if use_static_shape_: 186 x_tf = constant_op.constant(x_np) 187 else: 188 x_tf = array_ops.placeholder(dtype_) 189 190 if compute_uv_: 191 s_tf, u_tf, v_tf = linalg_ops.svd( 192 x_tf, compute_uv=compute_uv_, full_matrices=full_matrices_) 193 if use_static_shape_: 194 s_tf_val, u_tf_val, v_tf_val = self.evaluate([s_tf, u_tf, v_tf]) 195 else: 196 with self.session() as sess: 197 s_tf_val, u_tf_val, v_tf_val = sess.run( 198 [s_tf, u_tf, v_tf], feed_dict={x_tf: x_np}) 199 else: 200 s_tf = linalg_ops.svd( 201 x_tf, compute_uv=compute_uv_, full_matrices=full_matrices_) 202 if use_static_shape_: 203 s_tf_val = self.evaluate(s_tf) 204 else: 205 with self.session() as sess: 206 s_tf_val = sess.run(s_tf, feed_dict={x_tf: x_np}) 207 208 if compute_uv_: 209 u_np, s_np, v_np = np.linalg.svd( 210 x_np, compute_uv=compute_uv_, full_matrices=full_matrices_) 211 else: 212 s_np = np.linalg.svd( 213 x_np, compute_uv=compute_uv_, full_matrices=full_matrices_) 214 # We explicitly avoid the situation where numpy eliminates a first 215 # dimension that is equal to one. 216 s_np = np.reshape(s_np, s_tf_val.shape) 217 218 CompareSingularValues(self, s_np, s_tf_val, tol) 219 if compute_uv_: 220 CompareSingularVectors(self, u_np, u_tf_val, min(shape_[-2:]), tol) 221 CompareSingularVectors(self, np.conj(np.swapaxes(v_np, -2, -1)), v_tf_val, 222 min(shape_[-2:]), tol) 223 CheckApproximation(self, x_np, u_tf_val, s_tf_val, v_tf_val, 224 full_matrices_, tol) 225 CheckUnitary(self, u_tf_val, tol) 226 CheckUnitary(self, v_tf_val, tol) 227 228 return Test 229 230 231class SvdGradOpTest(test.TestCase): 232 pass # Filled in below 233 234 235def _NormalizingSvd(tf_a, full_matrices_): 236 tf_s, tf_u, tf_v = linalg_ops.svd( 237 tf_a, compute_uv=True, full_matrices=full_matrices_) 238 # Singular vectors are only unique up to an arbitrary phase. We normalize 239 # the vectors such that the first component of u (if m >=n) or v (if n > m) 240 # have phase 0. 241 m = tf_a.shape[-2] 242 n = tf_a.shape[-1] 243 if m >= n: 244 top_rows = tf_u[..., 0:1, :] 245 else: 246 top_rows = tf_v[..., 0:1, :] 247 if tf_u.dtype.is_complex: 248 angle = -math_ops.angle(top_rows) 249 phase = math_ops.complex(math_ops.cos(angle), math_ops.sin(angle)) 250 else: 251 phase = math_ops.sign(top_rows) 252 tf_u *= phase[..., :m] 253 tf_v *= phase[..., :n] 254 return tf_s, tf_u, tf_v 255 256 257def _GetSvdGradOpTest(dtype_, shape_, compute_uv_, full_matrices_): 258 259 @test_util.run_in_graph_and_eager_modes(use_gpu=True) 260 def Test(self): 261 262 def RandomInput(): 263 np.random.seed(42) 264 a = np.random.uniform(low=-1.0, high=1.0, size=shape_).astype(dtype_) 265 if dtype_ in [np.complex64, np.complex128]: 266 a += 1j * np.random.uniform( 267 low=-1.0, high=1.0, size=shape_).astype(dtype_) 268 return a 269 270 # Optimal stepsize for central difference is O(epsilon^{1/3}). 271 # See Equation (21) in: 272 # http://www.karenkopecky.net/Teaching/eco613614/Notes_NumericalDifferentiation.pdf 273 # TODO(rmlarsen): Move step size control to gradient checker. 274 epsilon = np.finfo(dtype_).eps 275 delta = 0.25 * epsilon**(1.0 / 3.0) 276 if dtype_ in [np.float32, np.complex64]: 277 tol = 3e-2 278 else: 279 tol = 1e-6 280 if compute_uv_: 281 funcs = [ 282 lambda a: _NormalizingSvd(a, full_matrices_)[0], 283 lambda a: _NormalizingSvd(a, full_matrices_)[1], 284 lambda a: _NormalizingSvd(a, full_matrices_)[2] 285 ] 286 else: 287 funcs = [lambda a: linalg_ops.svd(a, compute_uv=False)] 288 289 for f in funcs: 290 theoretical, numerical = gradient_checker_v2.compute_gradient( 291 f, [RandomInput()], delta=delta) 292 self.assertAllClose(theoretical, numerical, atol=tol, rtol=tol) 293 294 return Test 295 296 297class SvdGradGradOpTest(test.TestCase): 298 pass # Filled in below 299 300 301def _GetSvdGradGradOpTest(dtype_, shape_, compute_uv_, full_matrices_): 302 303 @test_util.run_v1_only("b/120545219") 304 def Test(self): 305 np.random.seed(42) 306 a = np.random.uniform(low=-1.0, high=1.0, size=shape_).astype(dtype_) 307 if dtype_ in [np.complex64, np.complex128]: 308 a += 1j * np.random.uniform( 309 low=-1.0, high=1.0, size=shape_).astype(dtype_) 310 # Optimal stepsize for central difference is O(epsilon^{1/3}). 311 # See Equation (21) in: 312 # http://www.karenkopecky.net/Teaching/eco613614/Notes_NumericalDifferentiation.pdf 313 # TODO(rmlarsen): Move step size control to gradient checker. 314 epsilon = np.finfo(dtype_).eps 315 delta = 0.1 * epsilon**(1.0 / 3.0) 316 tol = 1e-5 317 with self.session(): 318 tf_a = constant_op.constant(a) 319 if compute_uv_: 320 tf_s, tf_u, tf_v = _NormalizingSvd(tf_a, full_matrices_) 321 outputs = [tf_s, tf_u, tf_v] 322 else: 323 tf_s = linalg_ops.svd(tf_a, compute_uv=False) 324 outputs = [tf_s] 325 outputs_sums = [math_ops.reduce_sum(o) for o in outputs] 326 tf_func_outputs = math_ops.add_n(outputs_sums) 327 grad = gradients_impl.gradients(tf_func_outputs, tf_a)[0] 328 x_init = np.random.uniform( 329 low=-1.0, high=1.0, size=shape_).astype(dtype_) 330 if dtype_ in [np.complex64, np.complex128]: 331 x_init += 1j * np.random.uniform( 332 low=-1.0, high=1.0, size=shape_).astype(dtype_) 333 theoretical, numerical = gradient_checker.compute_gradient( 334 tf_a, 335 tf_a.get_shape().as_list(), 336 grad, 337 grad.get_shape().as_list(), 338 x_init_value=x_init, 339 delta=delta) 340 self.assertAllClose(theoretical, numerical, atol=tol, rtol=tol) 341 return Test 342 343 344class SVDBenchmark(test.Benchmark): 345 346 shapes = [ 347 (4, 4), 348 (8, 8), 349 (16, 16), 350 (101, 101), 351 (256, 256), 352 (1024, 1024), 353 (2048, 2048), 354 (1, 8, 8), 355 (10, 8, 8), 356 (100, 8, 8), 357 (1000, 8, 8), 358 (1, 32, 32), 359 (10, 32, 32), 360 (100, 32, 32), 361 (1000, 32, 32), 362 (1, 256, 256), 363 (10, 256, 256), 364 (100, 256, 256), 365 ] 366 367 def benchmarkSVDOp(self): 368 for shape_ in self.shapes: 369 with ops.Graph().as_default(), \ 370 session.Session(config=benchmark.benchmark_config()) as sess, \ 371 ops.device("/cpu:0"): 372 matrix_value = np.random.uniform( 373 low=-1.0, high=1.0, size=shape_).astype(np.float32) 374 matrix = variables.Variable(matrix_value) 375 u, s, v = linalg_ops.svd(matrix) 376 self.evaluate(variables.global_variables_initializer()) 377 self.run_op_benchmark( 378 sess, 379 control_flow_ops.group(u, s, v), 380 min_iters=25, 381 name="SVD_cpu_{shape}".format(shape=shape_)) 382 383 if test.is_gpu_available(True): 384 with ops.Graph().as_default(), \ 385 session.Session(config=benchmark.benchmark_config()) as sess, \ 386 ops.device("/device:GPU:0"): 387 matrix_value = np.random.uniform( 388 low=-1.0, high=1.0, size=shape_).astype(np.float32) 389 matrix = variables.Variable(matrix_value) 390 u, s, v = linalg_ops.svd(matrix) 391 self.evaluate(variables.global_variables_initializer()) 392 self.run_op_benchmark( 393 sess, 394 control_flow_ops.group(u, s, v), 395 min_iters=25, 396 name="SVD_gpu_{shape}".format(shape=shape_)) 397 398 399if __name__ == "__main__": 400 dtypes_to_test = [np.float32, np.float64, np.complex64, np.complex128] 401 for compute_uv in False, True: 402 for full_matrices in False, True: 403 for dtype in dtypes_to_test: 404 for rows in 0, 1, 2, 5, 10, 32, 100: 405 for cols in 0, 1, 2, 5, 10, 32, 100: 406 for batch_dims in [(), (3,)] + [(3, 2)] * (max(rows, cols) < 10): 407 full_shape = batch_dims + (rows, cols) 408 for use_static_shape in set([True, False]): 409 name = "%s_%s_static_shape_%s__compute_uv_%s_full_%s" % ( 410 dtype.__name__, "_".join(map(str, full_shape)), 411 use_static_shape, compute_uv, full_matrices) 412 _AddTest( 413 SvdOpTest, "Svd", name, 414 _GetSvdOpTest(dtype, full_shape, use_static_shape, 415 compute_uv, full_matrices)) 416 for compute_uv in False, True: 417 for full_matrices in False, True: 418 dtypes = ([np.float32, np.float64] + [np.complex64, np.complex128] * 419 (not compute_uv)) 420 for dtype in dtypes: 421 mat_shapes = [(10, 11), (11, 10), (11, 11), (2, 2, 2, 3)] 422 if not full_matrices or not compute_uv: 423 mat_shapes += [(5, 11), (11, 5)] 424 for mat_shape in mat_shapes: 425 for batch_dims in [(), (3,)]: 426 full_shape = batch_dims + mat_shape 427 name = "%s_%s_compute_uv_%s_full_%s" % (dtype.__name__, "_".join( 428 map(str, full_shape)), compute_uv, full_matrices) 429 _AddTest( 430 SvdGradOpTest, "SvdGrad", name, 431 _GetSvdGradOpTest(dtype, full_shape, compute_uv, full_matrices)) 432 # The results are too inaccurate for float32. 433 if dtype in (np.float64, np.complex128): 434 _AddTest( 435 SvdGradGradOpTest, "SvdGradGrad", name, 436 _GetSvdGradGradOpTest(dtype, full_shape, compute_uv, 437 full_matrices)) 438 test.main() 439