Added Cholesky Decomposition module in both decompositions.py and decompositions_test.py scripts of both numpy and tensorflow module.

tensornetwork/backends/numpy/decompositions.py

@@ -15,6 +15,7 @@
 
 from typing import Optional, Any, Tuple
 import numpy
+
 Tensor = Any
 
 
@@ -122,3 +123,32 @@ def rq(
   r = np.reshape(r, list(left_dims) + [center_dim])
   q = np.reshape(q, [center_dim] + list(right_dims))
   return r, q
+
+
+def cholesky(
+    np,  # TODO: Typing
+    tensor: Tensor,
+    pivot_axis: int,
+) -> Tuple[Tensor, Tensor]:
+  """Computes the cholesky decomposition of a tensor.
+
+  See tensornetwork.backends.tensorflow.decompositions for details.
+  """
+  left_dims = tensor.shape[:pivot_axis]
+  right_dims = tensor.shape[pivot_axis:]
+  tensor = np.reshape(tensor, [numpy.prod(left_dims), numpy.prod(right_dims)])
+  n = tensor.shape[0]
+  m = tensor.shape[1]
+  if n != m:
+    print("The input must be a square matrix")
+  elif (np.allclose(tensor, np.matrix.getH(tensor)) == False):
+    print("The input must be a Hermitian Matrix")
+  elif (np.all(np.linalg.eigvals(tensor) > 0) == False):
+    print("The input must be a Positive Definite Matrix")
+  else:
+    L = np.linalg.cholesky(tensor)
+    L_trans = np.matrix.getH(L)
+    center_dim = L.shape[1]
+    L = np.reshape(L, list(left_dims) + [center_dim])
+    L_trans = np.reshape(L_trans, [center_dim] + list(right_dims))
+    return L, L_trans

tensornetwork/backends/numpy/decompositions_test.py

@@ -104,6 +104,16 @@ def test_max_truncation_error_relative(self):
     np.testing.assert_almost_equal(trunc_sv_absolute, [0.1])
     np.testing.assert_almost_equal(trunc_sv_relative, [0.2, 0.1])
 
+  def test_expected_shapes_cholesky(self):
+    val = np.array([[[25, 15, -5]], [[15, 18, 0]], [[-5, 0, 11]]])
+    L, L_trans = decompositions.cholesky(np, val, -1)
+    self.assertEqual(L.shape, (3, 1, 3))
+
+  def test_cholesky(self):
+    random_matrix = np.array([[[25, 15, -5]], [[15, 18, 0]], [[-5, 0, 11]]])
+    L, L_trans = decompositions.cholesky(np, random_matrix, -1)
+    self.assertAllClose(L.dot(L_trans), random_matrix)
+
 
 if __name__ == '__main__':
   tf.test.main()

tensornetwork/backends/tensorflow/decompositions.py

@@ -127,8 +127,8 @@ def svd(
 
 
 def qr(
-    tf: Any,
-    tensor: Tensor,
+    tf: Any, 
+    tensor: Tensor, 
     pivot_axis: int,
     non_negative_diagonal: bool
 ) -> Tuple[Tensor, Tensor]:
@@ -177,8 +177,8 @@ def qr(
 
 
 def rq(
-    tf: Any,
-    tensor: Tensor,
+    tf: Any, 
+    tensor: Tensor, 
     pivot_axis: int,
     non_negative_diagonal: bool
 ) -> Tuple[Tensor, Tensor]:
@@ -226,3 +226,61 @@ def rq(
   r = tf.reshape(r, tf.concat([left_dims, [center_dim]], axis=-1))
   q = tf.reshape(q, tf.concat([[center_dim], right_dims], axis=-1))
   return r, q
+
+
+def cholesky(
+    tf: Any,
+    tensor: Tensor,
+    pivot_axis: int,
+) -> Tuple[Tensor, Tensor]:
+  """Computes the Cholesky decomposition of a tensor.
+
+  The Cholesky decomposition is performed by treating the tensor as a matrix,
+  with an effective left (row) index resulting from combining the axes
+  `tensor.shape[:pivot_axis]` and an effective right (column) index
+  resulting from combining the axes `tensor.shape[pivot_axis:]`.
+
+  For example, if `tensor` had a shape (2, 3, 4, 5) and `pivot_axis` was 2,
+  then `r` would have shape (2, 3, 6), and `q` would
+  have shape (6, 4, 5).
+
+  The output consists of two tensors `L, L_trans` such that:
+  ```python
+      L [i1,...,iN, j] * L_trans[j, k1,...,kM] == tensor[i1,...,iN, k1,...,kM]
+  ```
+  Note that the output ordering matches numpy.linalg.svd rather than tf.svd.
+
+  Args:
+    tf: The tensorflow module.
+    tensor: A tensor to be decomposed.
+    pivot_axis: Where to split the tensor's axes before flattening into a
+      matrix.
+
+  Returns:
+    L: Lower Triangular Matrix.
+    L_trans: Conjugate Transpose of L.
+  """
+  left_dims = tf.shape(tensor)[:pivot_axis]
+  right_dims = tf.shape(tensor)[pivot_axis:]
+  tensor = tf.reshape(tensor,
+                      [tf.reduce_prod(left_dims),
+                       tf.reduce_prod(right_dims)])
+  n = tensor.shape[0]
+  m = tensor.shape[1]
+  tensor = tf.cast(tensor, dtype=tf.complex128, name=None)
+  if n != m:
+    print("The input must be a square matrix")
+  elif (tf.experimental.numpy.allclose(tensor,
+                                       tf.linalg.adjoint(tensor)) == False):
+    print("The input must be a Hermitian Matrix")
+  elif (tf.experimental.numpy.all(
+      tf.math.real(tf.linalg.eigvals(tensor)) > 0) == False):
+    print("The input must be a Positive Definite Matrix")
+  else:
+    L = tf.linalg.cholesky(tensor, name=None)
+    L_trans = tf.transpose(L, perm=None, conjugate=True)
+    center_dim = tf.shape(L)[1]
+    L = tf.reshape(L, tf.concat([left_dims, [center_dim]], axis=-1))
+    L_trans = tf.reshape(L_trans, tf.concat([[center_dim], right_dims],
+                                            axis=-1))
+    return L, L_trans

tensornetwork/backends/tensorflow/decompositions_test.py

@@ -122,6 +122,16 @@ def test_max_truncation_error_relative(self):
     np.testing.assert_almost_equal(trunc_sv_absolute, [0.1])
     np.testing.assert_almost_equal(trunc_sv_relative, [0.2, 0.1])
 
+  def test_expected_shapes_cholesky(self):
+    val = tf.constant([[[25, 15, -5]], [[15, 18, 0]], [[-5, 0, 11]]])
+    L, L_trans = decompositions.cholesky(tf, val, -1)
+    self.assertEqual(L.shape, (3, 1, 3))
+
+  def test_cholesky(self):
+    random_matrix = tf.constant([[[25, 15, -5]], [[15, 18, 0]], [[-5, 0, 11]]])
+    L, L_trans = decompositions.cholesky(tf, random_matrix, -1)
+    self.assertAllClose(tf.tensordot(L, L_trans, ([2], [0])), random_matrix)
+
 
 if __name__ == '__main__':
   tf.test.main()