def testNonSquareMatrix(self):
   with self.assertRaises(ValueError):
     linalg_ops.cholesky(np.array([[1., 2., 3.], [3., 4., 5.]]))
   with self.assertRaises(ValueError):
     linalg_ops.cholesky(
         np.array([[[1., 2., 3.], [3., 4., 5.]], [[1., 2., 3.], [3., 4., 5.]]
                  ]))

  def benchmarkCholeskyOp(self):
    for shape in self.shapes:
      with ops.Graph().as_default(), \
          session.Session() as sess, \
          ops.device("/cpu:0"):
        matrix = variables.Variable(self._GenerateMatrix(shape))
        l = linalg_ops.cholesky(matrix)
        variables.global_variables_initializer().run()
        self.run_op_benchmark(
            sess,
            control_flow_ops.group(
                l,),
            min_iters=25,
            name="cholesky_cpu_{shape}".format(shape=shape))

      if test.is_gpu_available(True):
        with ops.Graph().as_default(), \
            session.Session() as sess, \
            ops.device("/device:GPU:0"):
          matrix = variables.Variable(self._GenerateMatrix(shape))
          l = linalg_ops.cholesky(matrix)
          variables.global_variables_initializer().run()
          self.run_op_benchmark(
              sess,
              control_flow_ops.group(
                  l,),
              min_iters=25,
              name="cholesky_gpu_{shape}".format(shape=shape))

 def testWrongDimensions(self):
   for dtype in self.float_types:
     tensor3 = constant_op.constant([1., 2.], dtype=dtype)
     with self.assertRaises(ValueError):
       linalg_ops.cholesky(tensor3)
     with self.assertRaises(ValueError):
       linalg_ops.cholesky(tensor3)

 def testNonSquareMatrix(self):
   for dtype in self.float_types:
     with self.assertRaises(ValueError):
       linalg_ops.cholesky(np.array([[1., 2., 3.], [3., 4., 5.]], dtype=dtype))
     with self.assertRaises(ValueError):
       linalg_ops.cholesky(
           np.array(
               [[[1., 2., 3.], [3., 4., 5.]], [[1., 2., 3.], [3., 4., 5.]]],
               dtype=dtype))

 def testConcurrentExecutesWithoutError(self):
   with self.test_session(use_gpu=True) as sess:
     matrix1 = random_ops.random_normal([5, 5], seed=42)
     matrix2 = random_ops.random_normal([5, 5], seed=42)
     matrix1 = math_ops.matmul(matrix1, matrix1, adjoint_a=True)
     matrix2 = math_ops.matmul(matrix2, matrix2, adjoint_a=True)
     c1 = linalg_ops.cholesky(matrix1)
     c2 = linalg_ops.cholesky(matrix2)
     c1_val, c2_val = sess.run([c1, c2])
     self.assertAllEqual(c1_val, c2_val)

 def _variance(self):
   x = math_ops.sqrt(self.df) * self.scale_operator_pd.to_dense()
   d = array_ops.expand_dims(array_ops.matrix_diag_part(x), -1)
   v = math_ops.square(x) + math_ops.matmul(d, d, adjoint_b=True)
   if self.cholesky_input_output_matrices:
     return linalg_ops.cholesky(v)
   return v

    def runFiniteDifferences(self, shapes, dtypes=(dtypes_lib.float32, dtypes_lib.float64), scalarTest=False):
        with self.test_session(use_gpu=False):
            for shape in shapes:
                for batch in False, True:
                    for dtype in dtypes:
                        if not scalarTest:
                            x = constant_op.constant(np.random.randn(shape[0], shape[1]), dtype)
                            tensor = math_ops.matmul(x, array_ops.transpose(x)) / shape[0]
                        else:
                            # This is designed to be a faster test for larger matrices.
                            x = constant_op.constant(np.random.randn(), dtype)
                            R = constant_op.constant(np.random.randn(shape[0], shape[1]), dtype)
                            e = math_ops.mul(R, x)
                            tensor = math_ops.matmul(e, array_ops.transpose(e)) / shape[0]

                        # Inner-most matrices in tensor are positive definite.
                        if batch:
                            tensor = array_ops.tile(array_ops.expand_dims(tensor, 0), [4, 1, 1])
                        y = linalg_ops.cholesky(tensor)
                        if scalarTest:
                            y = math_ops.reduce_mean(y)
                        error = gradient_checker.compute_gradient_error(x, x._shape_as_list(), y, y._shape_as_list())
                        tf_logging.info("error = %f", error)
                        if dtype == dtypes_lib.float64:
                            self.assertLess(error, 1e-5)
                        else:
                            self.assertLess(error, 3e-3)

  def _overdetermined(op, grad):
    """Gradients for the overdetermined case of MatrixSolveLs.

    This is the backprop for the solution to the normal equations of the first
    kind:
       X = F(A, B) = (A^T * A + lambda * I)^{-1} * A^T * B
    which solve the least squares problem
       min ||A * X - B||_F^2 + lambda ||X||_F^2.
    """
    a = op.inputs[0]
    b = op.inputs[1]
    l2_regularizer = math_ops.cast(op.inputs[2], a.dtype.base_dtype)
    x = op.outputs[0]
    a_shape = array_ops.shape(a)
    batch_shape = a_shape[:-2]
    n = a_shape[-1]

    identity = linalg_ops.eye(n, batch_shape=batch_shape, dtype=a.dtype)
    gramian = math_ops.matmul(a, a, adjoint_a=True) + l2_regularizer * identity
    chol = linalg_ops.cholesky(gramian)
    # Temporary z = (A^T * A + lambda * I)^{-1} * grad.
    z = linalg_ops.cholesky_solve(chol, grad)
    xzt = math_ops.matmul(x, z, adjoint_b=True)
    zx_sym = xzt + array_ops.matrix_transpose(xzt)
    grad_a = -math_ops.matmul(a, zx_sym) + math_ops.matmul(b, z, adjoint_b=True)
    grad_b = math_ops.matmul(a, z)
    return (grad_a, grad_b, None)

  def _underdetermined(op, grad):
    """Gradients for the underdetermined case of MatrixSolveLs.

    This is the backprop for the solution to the normal equations of the second
    kind:
      X = F(A, B) = A * (A*A^T + lambda*I)^{-1} * B
    that (for lambda=0) solve the least squares problem
      min ||X||_F subject to A*X = B.
    """
    a = op.inputs[0]
    b = op.inputs[1]
    l2_regularizer = math_ops.cast(op.inputs[2], a.dtype.base_dtype)
    a_shape = array_ops.shape(a)
    batch_shape = a_shape[:-2]
    m = a_shape[-2]

    identity = linalg_ops.eye(m, batch_shape=batch_shape, dtype=a.dtype)
    gramian = math_ops.matmul(a, a, adjoint_b=True) + l2_regularizer * identity
    chol = linalg_ops.cholesky(gramian)
    grad_b = linalg_ops.cholesky_solve(chol, math_ops.matmul(a, grad))
    # Temporary tmp = (A * A^T + lambda * I)^{-1} * B.
    tmp = linalg_ops.cholesky_solve(chol, b)
    a1 = math_ops.matmul(tmp, a, adjoint_a=True)
    a1 = -math_ops.matmul(grad_b, a1)
    a2 = grad - math_ops.matmul(a, grad_b, adjoint_a=True)
    a2 = math_ops.matmul(tmp, a2, adjoint_b=True)
    grad_a = a1 + a2
    return (grad_a, grad_b, None)

 def _verifyCholesky(self, x, atol=1e-6):
   # Verify that LL^T == x.
   with self.test_session() as sess:
     placeholder = array_ops.placeholder(
         dtypes.as_dtype(x.dtype), shape=x.shape)
     with self.test_scope():
       chol = linalg_ops.cholesky(placeholder)
     verification = math_ops.matmul(chol, chol, adjoint_b=True)
     self._verifyCholeskyBase(sess, placeholder, x, chol, verification, atol)

 def _log_abs_determinant(self):
   logging.warn(
       "Using (possibly slow) default implementation of determinant."
       "  Requires conversion to a dense matrix and O(N^3) operations.")
   if self._can_use_cholesky():
     diag = array_ops.matrix_diag_part(linalg_ops.cholesky(self.to_dense()))
     return 2 * math_ops.reduce_sum(math_ops.log(diag), axis=[-1])
   _, log_abs_det = linalg.slogdet(self.to_dense())
   return log_abs_det

 def test_cholesky(self):
   with self.test_session(graph=ops.Graph()) as sess:
     sess.graph.seed = random_seed.DEFAULT_GRAPH_SEED
     operator, mat = self.operator_and_matrix(
         shapes_info, dtype, use_placeholder=use_placeholder,
         ensure_self_adjoint_and_pd=True)
     op_chol = operator.cholesky().to_dense()
     mat_chol = linalg_ops.cholesky(mat)
     op_chol_v, mat_chol_v = sess.run([op_chol, mat_chol])
     self.assertAC(mat_chol_v, op_chol_v)

  def benchmarkCholeskyOp(self):
    for size in self.sizes:
      data = self._GenerateData(size)

      with ops.Graph().as_default(), \
          session.Session() as sess, \
          ops.device("/cpu:0"):
        l = linalg_ops.cholesky(data)
        self.run_op_benchmark(
            sess, control_flow_ops.group(l,),
            min_iters=25,
            name="cholesky_cpu_{size}".format(size=size))

      if test.is_gpu_available(True):
        with ops.Graph().as_default(), \
            session.Session() as sess, \
            ops.device("/gpu:0"):
          l = linalg_ops.cholesky(data)
          self.run_op_benchmark(
              sess, l,
              min_iters=25,
              name="cholesky_gpu_{size}".format(size=size))

  def testAgainstSpecialized(self):
    np.random.seed(0)
    data = np.random.randn(33, 33).astype(np.float32)
    data = np.matmul(data, data.T)
    grad_data = np.random.randn(*data.shape).astype(np.float32)

    with ops.Graph().as_default(), self.test_session(use_gpu=False) as s:
      x = constant_op.constant(data, dtypes_lib.float32)
      chol = linalg_ops.cholesky(x)
      composite_grad = gradients_impl.gradients(chol, x, grad_data)[0]
      specialized_grad = SpecializedGrad(chol, grad_data)
      reference, actual = s.run([specialized_grad, composite_grad])
    self.assertAllClose(reference, actual)

  def __init__(self,
               df,
               scale,
               cholesky_input_output_matrices=False,
               validate_args=False,
               allow_nan_stats=True,
               name="WishartFull"):
    """Construct Wishart distributions.

    Args:
      df: `float` or `double` `Tensor`. Degrees of freedom, must be greater than
        or equal to dimension of the scale matrix.
      scale: `float` or `double` `Tensor`. The symmetric positive definite
        scale matrix of the distribution.
      cholesky_input_output_matrices: Python `bool`. Any function which whose
        input or output is a matrix assumes the input is Cholesky and returns a
        Cholesky factored matrix. Example `log_prob` input takes a Cholesky and
        `sample_n` returns a Cholesky when
        `cholesky_input_output_matrices=True`.
      validate_args: Python `bool`, default `False`. When `True` distribution
        parameters are checked for validity despite possibly degrading runtime
        performance. When `False` invalid inputs may silently render incorrect
        outputs.
      allow_nan_stats: Python `bool`, default `True`. When `True`, statistics
        (e.g., mean, mode, variance) use the value "`NaN`" to indicate the
        result is undefined. When `False`, an exception is raised if one or
        more of the statistic's batch members are undefined.
      name: Python `str` name prefixed to Ops created by this class.
    """
    parameters = dict(locals())
    with ops.name_scope(name) as name:
      with ops.name_scope("init", values=[scale]):
        scale = ops.convert_to_tensor(scale)
        if validate_args:
          scale = distribution_util.assert_symmetric(scale)
        chol = linalg_ops.cholesky(scale)
        chol = control_flow_ops.with_dependencies([
            check_ops.assert_positive(array_ops.matrix_diag_part(chol))
        ] if validate_args else [], chol)
    super(WishartFull, self).__init__(
        df=df,
        scale_operator=linalg.LinearOperatorLowerTriangular(
            tril=chol,
            is_non_singular=True,
            is_positive_definite=True,
            is_square=True),
        cholesky_input_output_matrices=cholesky_input_output_matrices,
        validate_args=validate_args,
        allow_nan_stats=allow_nan_stats,
        name=name)
    self._parameters = parameters

 def test_works_with_five_different_random_pos_def_matrices(self):
   for n in range(1, 6):
     for np_type, atol in [(np.float32, 0.05), (np.float64, 1e-5)]:
       with self.session(use_gpu=True):
         # Create 2 x n x n matrix
         array = np.array(
             [_RandomPDMatrix(n, self.rng),
              _RandomPDMatrix(n, self.rng)]).astype(np_type)
         chol = linalg_ops.cholesky(array)
         for k in range(1, 3):
           rhs = self.rng.randn(2, n, k).astype(np_type)
           x = linalg_ops.cholesky_solve(chol, rhs)
           self.assertAllClose(
               rhs, math_ops.matmul(array, x).eval(), atol=atol)

  def __init__(self,
               matrix,
               is_non_singular=None,
               is_self_adjoint=None,
               is_positive_definite=None,
               name="LinearOperatorMatrix"):
    """Initialize a `LinearOperatorMatrix`.

    Args:
      matrix:  Shape `[B1,...,Bb, M, N]` with `b >= 0`, `M, N >= 0`.
        Allowed dtypes: `float32`, `float64`, `complex64`, `complex128`.
      is_non_singular:  Expect that this operator is non-singular.
      is_self_adjoint:  Expect that this operator is equal to its hermitian
        transpose.
      is_positive_definite:  Expect that this operator is positive definite,
        meaning the real part of all eigenvalues is positive.  We do not require
        the operator to be self-adjoint to be positive-definite.  See:
        https://en.wikipedia.org/wiki/Positive-definite_matrix
            #Extension_for_non_symmetric_matrices
      name: A name for this `LinearOperator`.

    Raises:
      TypeError:  If `diag.dtype` is not an allowed type.
    """

    allowed_dtypes = [
        dtypes.float32, dtypes.float64, dtypes.complex64, dtypes.complex128]

    with ops.name_scope(name, values=[matrix]):
      self._matrix = ops.convert_to_tensor(matrix, name="matrix")

      dtype = self._matrix.dtype
      if dtype not in allowed_dtypes:
        raise TypeError(
            "Argument matrix must have dtype in %s.  Found: %s"
            % (allowed_dtypes, dtype))

      # Special treatment for (real) Symmetric Positive Definite.
      self._is_spd = (
          (not dtype.is_complex) and is_self_adjoint and is_positive_definite)
      if self._is_spd:
        self._chol = linalg_ops.cholesky(self._matrix)

      super(LinearOperatorMatrix, self).__init__(
          dtype=self._matrix.dtype,
          graph_parents=[self._matrix],
          is_non_singular=is_non_singular,
          is_self_adjoint=is_self_adjoint,
          is_positive_definite=is_positive_definite,
          name=name)

 def test_cholesky(self):
   self._skip_if_tests_to_skip_contains("cholesky")
   for use_placeholder in self._use_placeholder_options:
     for build_info in self._operator_build_infos:
       for dtype in self._dtypes_to_test:
         with self.test_session(graph=ops.Graph()) as sess:
           sess.graph.seed = random_seed.DEFAULT_GRAPH_SEED
           operator, mat = self._operator_and_matrix(
               build_info, dtype, use_placeholder=use_placeholder,
               ensure_self_adjoint_and_pd=True)
           op_chol = operator.cholesky().to_dense()
           mat_chol = linalg_ops.cholesky(mat)
           op_chol_v, mat_chol_v = sess.run([op_chol, mat_chol])
           self.assertAC(mat_chol_v, op_chol_v)

 def _solve(self, rhs, adjoint=False, adjoint_arg=False):
   """Default implementation of _solve."""
   if self.is_square is False:
     raise NotImplementedError(
         "Solve is not yet implemented for non-square operators.")
   logging.warn(
       "Using (possibly slow) default implementation of solve."
       "  Requires conversion to a dense matrix and O(N^3) operations.")
   rhs = linalg.adjoint(rhs) if adjoint_arg else rhs
   if self._can_use_cholesky():
     return linear_operator_util.cholesky_solve_with_broadcast(
         linalg_ops.cholesky(self.to_dense()), rhs)
   return linear_operator_util.matrix_solve_with_broadcast(
       self.to_dense(), rhs, adjoint=adjoint)

 def _assert_positive_definite(self):
   """Default implementation of _assert_positive_definite."""
   logging.warn(
       "Using (possibly slow) default implementation of "
       "assert_positive_definite."
       "  Requires conversion to a dense matrix and O(N^3) operations.")
   # If the operator is self-adjoint, then checking that
   # Cholesky decomposition succeeds + results in positive diag is necessary
   # and sufficient.
   if self.is_self_adjoint:
     return check_ops.assert_positive(
         array_ops.matrix_diag_part(linalg_ops.cholesky(self.to_dense())),
         message="Matrix was not positive definite.")
   # We have no generic check for positive definite.
   raise NotImplementedError("assert_positive_definite is not implemented.")