Home | History | Annotate | Download | only in linalg 
      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 
     16 from __future__ import absolute_import
     17 from __future__ import division
     18 from __future__ import print_function
     19 
     20 import numpy as np
     21 
     22 from tensorflow.python.framework import dtypes
     23 from tensorflow.python.framework import ops
     24 from tensorflow.python.framework import random_seed
     25 from tensorflow.python.ops import array_ops
     26 from tensorflow.python.ops import math_ops
     27 from tensorflow.python.ops.linalg import linalg as linalg_lib
     28 from tensorflow.python.ops.linalg import linear_operator_test_util
     29 from tensorflow.python.platform import test
     30 
     31 linalg = linalg_lib
     32 random_seed.set_random_seed(23)
     33 
     34 
     35 class SquareLinearOperatorFullMatrixTest(
     36     linear_operator_test_util.SquareLinearOperatorDerivedClassTest):
     37   """Most tests done in the base class LinearOperatorDerivedClassTest."""
     38 
     39   def _operator_and_mat_and_feed_dict(self, shape, dtype, use_placeholder):
     40     shape = list(shape)
     41 
     42     matrix = linear_operator_test_util.random_positive_definite_matrix(shape,
     43                                                                        dtype)
     44 
     45     if use_placeholder:
     46       matrix_ph = array_ops.placeholder(dtype=dtype)
     47       # Evaluate here because (i) you cannot feed a tensor, and (ii)
     48       # values are random and we want the same value used for both mat and
     49       # feed_dict.
     50       matrix = matrix.eval()
     51       operator = linalg.LinearOperatorFullMatrix(matrix_ph, is_square=True)
     52       feed_dict = {matrix_ph: matrix}
     53     else:
     54       # is_square should be auto-detected here.
     55       operator = linalg.LinearOperatorFullMatrix(matrix)
     56       feed_dict = None
     57 
     58     # Convert back to Tensor.  Needed if use_placeholder, since then we have
     59     # already evaluated matrix to a numpy array.
     60     mat = ops.convert_to_tensor(matrix)
     61 
     62     return operator, mat, feed_dict
     63 
     64   def test_is_x_flags(self):
     65     # Matrix with two positive eigenvalues.
     66     matrix = [[1., 0.], [1., 11.]]
     67     operator = linalg.LinearOperatorFullMatrix(
     68         matrix,
     69         is_positive_definite=True,
     70         is_non_singular=True,
     71         is_self_adjoint=False)
     72     self.assertTrue(operator.is_positive_definite)
     73     self.assertTrue(operator.is_non_singular)
     74     self.assertFalse(operator.is_self_adjoint)
     75     # Auto-detected.
     76     self.assertTrue(operator.is_square)
     77 
     78   def test_assert_non_singular_raises_if_cond_too_big_but_finite(self):
     79     with self.test_session():
     80       tril = linear_operator_test_util.random_tril_matrix(
     81           shape=(50, 50), dtype=np.float32)
     82       diag = np.logspace(-2, 2, 50).astype(np.float32)
     83       tril = array_ops.matrix_set_diag(tril, diag)
     84       matrix = math_ops.matmul(tril, tril, transpose_b=True).eval()
     85       operator = linalg.LinearOperatorFullMatrix(matrix)
     86       with self.assertRaisesOpError("Singular matrix"):
     87         # Ensure that we have finite condition number...just HUGE.
     88         cond = np.linalg.cond(matrix)
     89         self.assertTrue(np.isfinite(cond))
     90         self.assertGreater(cond, 1e12)
     91         operator.assert_non_singular().run()
     92 
     93   def test_assert_non_singular_raises_if_cond_infinite(self):
     94     with self.test_session():
     95       matrix = [[1., 1.], [1., 1.]]
     96       # We don't pass the is_self_adjoint hint here, which means we take the
     97       # generic code path.
     98       operator = linalg.LinearOperatorFullMatrix(matrix)
     99       with self.assertRaisesOpError("Singular matrix"):
    100         operator.assert_non_singular().run()
    101 
    102   def test_assert_self_adjoint(self):
    103     matrix = [[0., 1.], [0., 1.]]
    104     operator = linalg.LinearOperatorFullMatrix(matrix)
    105     with self.test_session():
    106       with self.assertRaisesOpError("not equal to its adjoint"):
    107         operator.assert_self_adjoint().run()
    108 
    109   def test_assert_positive_definite(self):
    110     matrix = [[1., 1.], [1., 1.]]
    111     operator = linalg.LinearOperatorFullMatrix(matrix, is_self_adjoint=True)
    112     with self.test_session():
    113       with self.assertRaisesOpError("Cholesky decomposition was not success"):
    114         operator.assert_positive_definite().run()
    115 
    116 
    117 class SquareLinearOperatorFullMatrixSymmetricPositiveDefiniteTest(
    118     linear_operator_test_util.SquareLinearOperatorDerivedClassTest):
    119   """Most tests done in the base class LinearOperatorDerivedClassTest.
    120 
    121   In this test, the operator is constructed with hints that invoke the use of
    122   a Cholesky decomposition for solves/determinant.
    123   """
    124 
    125   def setUp(self):
    126     # Increase from 1e-6 to 1e-5.  This reduction in tolerance happens,
    127     # presumably, because we are taking a different code path in the operator
    128     # and the matrix.  The operator uses a Choleksy, the matrix uses standard
    129     # solve.
    130     self._atol[dtypes.float32] = 1e-5
    131     self._rtol[dtypes.float32] = 1e-5
    132     self._atol[dtypes.float64] = 1e-10
    133     self._rtol[dtypes.float64] = 1e-10
    134 
    135   @property
    136   def _dtypes_to_test(self):
    137     return [dtypes.float32, dtypes.float64]
    138 
    139   def _operator_and_mat_and_feed_dict(self, shape, dtype, use_placeholder):
    140     shape = list(shape)
    141 
    142     matrix = linear_operator_test_util.random_positive_definite_matrix(
    143         shape, dtype, force_well_conditioned=True)
    144 
    145     if use_placeholder:
    146       matrix_ph = array_ops.placeholder(dtype=dtype)
    147       # Evaluate here because (i) you cannot feed a tensor, and (ii)
    148       # values are random and we want the same value used for both mat and
    149       # feed_dict.
    150       matrix = matrix.eval()
    151       # is_square is auto-set because of self_adjoint/pd.
    152       operator = linalg.LinearOperatorFullMatrix(
    153           matrix_ph, is_self_adjoint=True, is_positive_definite=True)
    154       feed_dict = {matrix_ph: matrix}
    155     else:
    156       operator = linalg.LinearOperatorFullMatrix(
    157           matrix, is_self_adjoint=True, is_positive_definite=True)
    158       feed_dict = None
    159 
    160     # Convert back to Tensor.  Needed if use_placeholder, since then we have
    161     # already evaluated matrix to a numpy array.
    162     mat = ops.convert_to_tensor(matrix)
    163 
    164     return operator, mat, feed_dict
    165 
    166   def test_is_x_flags(self):
    167     # Matrix with two positive eigenvalues.
    168     matrix = [[1., 0.], [0., 7.]]
    169     operator = linalg.LinearOperatorFullMatrix(
    170         matrix, is_positive_definite=True, is_self_adjoint=True)
    171 
    172     self.assertTrue(operator.is_positive_definite)
    173     self.assertTrue(operator.is_self_adjoint)
    174 
    175     # Should be auto-set
    176     self.assertTrue(operator.is_non_singular)
    177     self.assertTrue(operator._can_use_cholesky)
    178     self.assertTrue(operator.is_square)
    179 
    180   def test_assert_non_singular(self):
    181     matrix = [[1., 1.], [1., 1.]]
    182     operator = linalg.LinearOperatorFullMatrix(
    183         matrix, is_self_adjoint=True, is_positive_definite=True)
    184     with self.test_session():
    185       # Cholesky decomposition may fail, so the error is not specific to
    186       # non-singular.
    187       with self.assertRaisesOpError(""):
    188         operator.assert_non_singular().run()
    189 
    190   def test_assert_self_adjoint(self):
    191     matrix = [[0., 1.], [0., 1.]]
    192     operator = linalg.LinearOperatorFullMatrix(
    193         matrix, is_self_adjoint=True, is_positive_definite=True)
    194     with self.test_session():
    195       with self.assertRaisesOpError("not equal to its adjoint"):
    196         operator.assert_self_adjoint().run()
    197 
    198   def test_assert_positive_definite(self):
    199     matrix = [[1., 1.], [1., 1.]]
    200     operator = linalg.LinearOperatorFullMatrix(
    201         matrix, is_self_adjoint=True, is_positive_definite=True)
    202     with self.test_session():
    203       # Cholesky decomposition may fail, so the error is not specific to
    204       # non-singular.
    205       with self.assertRaisesOpError(""):
    206         operator.assert_positive_definite().run()
    207 
    208 
    209 class NonSquareLinearOperatorFullMatrixTest(
    210     linear_operator_test_util.NonSquareLinearOperatorDerivedClassTest):
    211   """Most tests done in the base class LinearOperatorDerivedClassTest."""
    212 
    213   def _operator_and_mat_and_feed_dict(self, shape, dtype, use_placeholder):
    214     matrix = linear_operator_test_util.random_normal(shape, dtype=dtype)
    215     if use_placeholder:
    216       matrix_ph = array_ops.placeholder(dtype=dtype)
    217       # Evaluate here because (i) you cannot feed a tensor, and (ii)
    218       # values are random and we want the same value used for both mat and
    219       # feed_dict.
    220       matrix = matrix.eval()
    221       operator = linalg.LinearOperatorFullMatrix(matrix_ph)
    222       feed_dict = {matrix_ph: matrix}
    223     else:
    224       operator = linalg.LinearOperatorFullMatrix(matrix)
    225       feed_dict = None
    226 
    227     # Convert back to Tensor.  Needed if use_placeholder, since then we have
    228     # already evaluated matrix to a numpy array.
    229     mat = ops.convert_to_tensor(matrix)
    230 
    231     return operator, mat, feed_dict
    232 
    233   def test_is_x_flags(self):
    234     matrix = [[3., 2., 1.], [1., 1., 1.]]
    235     operator = linalg.LinearOperatorFullMatrix(
    236         matrix,
    237         is_self_adjoint=False)
    238     self.assertEqual(operator.is_positive_definite, None)
    239     self.assertEqual(operator.is_non_singular, None)
    240     self.assertFalse(operator.is_self_adjoint)
    241     self.assertFalse(operator.is_square)
    242 
    243   def test_matrix_must_have_at_least_two_dims_or_raises(self):
    244     with self.assertRaisesRegexp(ValueError, "at least 2 dimensions"):
    245       linalg.LinearOperatorFullMatrix([1.])
    246 
    247 
    248 if __name__ == "__main__":
    249   test.main()
    250