• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
1# Copyright 2017 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.linalg.linalg_impl.matrix_exponential."""
16
17import itertools
18
19import numpy as np
20
21from tensorflow.python.client import session
22from tensorflow.python.framework import constant_op
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 random_ops
28from tensorflow.python.ops import variables
29from tensorflow.python.ops.linalg import linalg_impl
30from tensorflow.python.platform import benchmark
31from tensorflow.python.platform import test
32
33
34def np_expm(x):  # pylint: disable=invalid-name
35  """Slow but accurate Taylor series matrix exponential."""
36  y = np.zeros(x.shape, dtype=x.dtype)
37  xn = np.eye(x.shape[0], dtype=x.dtype)
38  for n in range(40):
39    if n > 0:
40      xn /= float(n)
41    y += xn
42    xn = np.dot(xn, x)
43  return y
44
45
46class ExponentialOpTest(test.TestCase):
47
48  def _verifyExponential(self, x, np_type):
49    inp = x.astype(np_type)
50    with test_util.use_gpu():
51      tf_ans = linalg_impl.matrix_exponential(inp)
52      if x.size == 0:
53        np_ans = np.empty(x.shape, dtype=np_type)
54      else:
55        if x.ndim > 2:
56          np_ans = np.zeros(inp.shape, dtype=np_type)
57          for i in itertools.product(*[range(x) for x in inp.shape[:-2]]):
58            np_ans[i] = np_expm(inp[i])
59        else:
60          np_ans = np_expm(inp)
61      out = self.evaluate(tf_ans)
62      self.assertAllClose(np_ans, out, rtol=1e-3, atol=1e-3)
63
64  def _verifyExponentialReal(self, x):
65    for np_type in [np.float32, np.float64]:
66      self._verifyExponential(x, np_type)
67
68  def _verifyExponentialComplex(self, x):
69    for np_type in [np.complex64, np.complex128]:
70      self._verifyExponential(x, np_type)
71
72  def _makeBatch(self, matrix1, matrix2):
73    matrix_batch = np.concatenate(
74        [np.expand_dims(matrix1, 0),
75         np.expand_dims(matrix2, 0)])
76    matrix_batch = np.tile(matrix_batch, [2, 3, 1, 1])
77    return matrix_batch
78
79  def testNonsymmetricReal(self):
80    # 2x2 matrices
81    matrix1 = np.array([[1., 2.], [3., 4.]])
82    matrix2 = np.array([[1., 3.], [3., 5.]])
83    self._verifyExponentialReal(matrix1)
84    self._verifyExponentialReal(matrix2)
85    # A multidimensional batch of 2x2 matrices
86    self._verifyExponentialReal(self._makeBatch(matrix1, matrix2))
87
88  @test_util.run_deprecated_v1
89  def testNonsymmetricComplex(self):
90    matrix1 = np.array([[1., 2.], [3., 4.]])
91    matrix2 = np.array([[1., 3.], [3., 5.]])
92    matrix1 = matrix1.astype(np.complex64)
93    matrix1 += 1j * matrix1
94    matrix2 = matrix2.astype(np.complex64)
95    matrix2 += 1j * matrix2
96    self._verifyExponentialComplex(matrix1)
97    self._verifyExponentialComplex(matrix2)
98    # Complex batch
99    self._verifyExponentialComplex(self._makeBatch(matrix1, matrix2))
100
101  def testSymmetricPositiveDefiniteReal(self):
102    # 2x2 matrices
103    matrix1 = np.array([[2., 1.], [1., 2.]])
104    matrix2 = np.array([[3., -1.], [-1., 3.]])
105    self._verifyExponentialReal(matrix1)
106    self._verifyExponentialReal(matrix2)
107    # A multidimensional batch of 2x2 matrices
108    self._verifyExponentialReal(self._makeBatch(matrix1, matrix2))
109
110  def testSymmetricPositiveDefiniteComplex(self):
111    matrix1 = np.array([[2., 1.], [1., 2.]])
112    matrix2 = np.array([[3., -1.], [-1., 3.]])
113    matrix1 = matrix1.astype(np.complex64)
114    matrix1 += 1j * matrix1
115    matrix2 = matrix2.astype(np.complex64)
116    matrix2 += 1j * matrix2
117    self._verifyExponentialComplex(matrix1)
118    self._verifyExponentialComplex(matrix2)
119    # Complex batch
120    self._verifyExponentialComplex(self._makeBatch(matrix1, matrix2))
121
122  @test_util.run_deprecated_v1
123  def testNonSquareMatrix(self):
124    # When the exponential of a non-square matrix is attempted we should return
125    # an error
126    with self.assertRaises(ValueError):
127      linalg_impl.matrix_exponential(np.array([[1., 2., 3.], [3., 4., 5.]]))
128
129  @test_util.run_deprecated_v1
130  def testWrongDimensions(self):
131    # The input to the exponential should be at least a 2-dimensional tensor.
132    tensor3 = constant_op.constant([1., 2.])
133    with self.assertRaises(ValueError):
134      linalg_impl.matrix_exponential(tensor3)
135
136  def testInfinite(self):
137    # Check that the op does not loop forever on infinite inputs. (b/158433036)
138    in_tensor = [[np.inf, 1.], [1., 1.]]
139    result = self.evaluate(linalg_impl.matrix_exponential(in_tensor))
140    self.assertTrue(np.all(np.isnan(result)))
141
142  def testEmpty(self):
143    self._verifyExponentialReal(np.empty([0, 2, 2]))
144    self._verifyExponentialReal(np.empty([2, 0, 0]))
145
146  @test_util.run_deprecated_v1
147  def testDynamic(self):
148    with self.session() as sess:
149      inp = array_ops.placeholder(ops.dtypes.float32)
150      expm = linalg_impl.matrix_exponential(inp)
151      matrix = np.array([[1., 2.], [3., 4.]])
152      sess.run(expm, feed_dict={inp: matrix})
153
154  @test_util.run_deprecated_v1
155  def testConcurrentExecutesWithoutError(self):
156    with self.session():
157      matrix1 = random_ops.random_normal([5, 5], seed=42)
158      matrix2 = random_ops.random_normal([5, 5], seed=42)
159      expm1 = linalg_impl.matrix_exponential(matrix1)
160      expm2 = linalg_impl.matrix_exponential(matrix2)
161      expm = self.evaluate([expm1, expm2])
162      self.assertAllEqual(expm[0], expm[1])
163
164
165class MatrixExponentialBenchmark(test.Benchmark):
166
167  shapes = [
168      (4, 4),
169      (10, 10),
170      (16, 16),
171      (101, 101),
172      (256, 256),
173      (1000, 1000),
174      (1024, 1024),
175      (2048, 2048),
176      (513, 4, 4),
177      (513, 16, 16),
178      (513, 256, 256),
179  ]
180
181  def _GenerateMatrix(self, shape):
182    batch_shape = shape[:-2]
183    shape = shape[-2:]
184    assert shape[0] == shape[1]
185    n = shape[0]
186    matrix = np.ones(shape).astype(np.float32) / (2.0 * n) + np.diag(
187        np.ones(n).astype(np.float32))
188    return variables.Variable(np.tile(matrix, batch_shape + (1, 1)))
189
190  def benchmarkMatrixExponentialOp(self):
191    for shape in self.shapes:
192      with ops.Graph().as_default(), \
193          session.Session(config=benchmark.benchmark_config()) as sess, \
194          ops.device("/cpu:0"):
195        matrix = self._GenerateMatrix(shape)
196        expm = linalg_impl.matrix_exponential(matrix)
197        self.evaluate(variables.global_variables_initializer())
198        self.run_op_benchmark(
199            sess,
200            control_flow_ops.group(expm),
201            min_iters=25,
202            name="matrix_exponential_cpu_{shape}".format(shape=shape))
203
204      if test.is_gpu_available(True):
205        with ops.Graph().as_default(), \
206            session.Session(config=benchmark.benchmark_config()) as sess, \
207            ops.device("/gpu:0"):
208          matrix = self._GenerateMatrix(shape)
209          expm = linalg_impl.matrix_exponential(matrix)
210          self.evaluate(variables.global_variables_initializer())
211          self.run_op_benchmark(
212              sess,
213              control_flow_ops.group(expm),
214              min_iters=25,
215              name="matrix_exponential_gpu_{shape}".format(shape=shape))
216
217
218def _TestRandomSmall(dtype, batch_dims, size):
219
220  def Test(self):
221    np.random.seed(42)
222    shape = batch_dims + (size, size)
223    matrix = np.random.uniform(low=-1.0, high=1.0, size=shape).astype(dtype)
224    self._verifyExponentialReal(matrix)
225
226  return Test
227
228
229def _TestL1Norms(dtype, shape, scale):
230
231  def Test(self):
232    np.random.seed(42)
233    matrix = np.random.uniform(
234        low=-1.0, high=1.0, size=np.prod(shape)).reshape(shape).astype(dtype)
235    l1_norm = np.max(np.sum(np.abs(matrix), axis=matrix.ndim - 2))
236    matrix /= l1_norm
237    self._verifyExponentialReal(scale * matrix)
238
239  return Test
240
241
242if __name__ == "__main__":
243  for dtype_ in [np.float32, np.float64, np.complex64, np.complex128]:
244    for batch_ in [(), (2,), (2, 2)]:
245      for size_ in [4, 7]:
246        name = "%s_%d_%d" % (dtype_.__name__, len(batch_), size_)
247        setattr(ExponentialOpTest, "testL1Norms_" + name,
248                _TestRandomSmall(dtype_, batch_, size_))
249
250  for shape_ in [(3, 3), (2, 3, 3)]:
251    for dtype_ in [np.float32, np.complex64]:
252      for scale_ in [0.1, 1.5, 5.0, 20.0]:
253        name = "%s_%d_%d" % (dtype_.__name__, len(shape_), int(scale_ * 10))
254        setattr(ExponentialOpTest, "testL1Norms_" + name,
255                _TestL1Norms(dtype_, shape_, scale_))
256    for dtype_ in [np.float64, np.complex128]:
257      for scale_ in [0.01, 0.2, 0.5, 1.5, 6.0, 25.0]:
258        name = "%s_%d_%d" % (dtype_.__name__, len(shape_), int(scale_ * 100))
259        setattr(ExponentialOpTest, "testL1Norms_" + name,
260                _TestL1Norms(dtype_, shape_, scale_))
261  test.main()
262