• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
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