From 662c6d15daa4dfc5dd7bed7f6bbe2edd92a44eb9 Mon Sep 17 00:00:00 2001 From: mganahl Date: Tue, 12 Jan 2021 14:18:54 +0100 Subject: [PATCH 1/4] add eps to backends --- tensornetwork/backends/abstract_backend.py | 14 +++++++++++ tensornetwork/backends/jax/jax_backend.py | 12 ++++++++++ tensornetwork/backends/numpy/numpy_backend.py | 24 ++++++++++++++----- .../backends/pytorch/pytorch_backend.py | 14 ++++++++++- .../backends/symmetric/symmetric_backend.py | 12 ++++++++++ .../backends/tensorflow/tensorflow_backend.py | 12 ++++++++++ 6 files changed, 81 insertions(+), 7 deletions(-) diff --git a/tensornetwork/backends/abstract_backend.py b/tensornetwork/backends/abstract_backend.py index fc28008d3..b9f399354 100644 --- a/tensornetwork/backends/abstract_backend.py +++ b/tensornetwork/backends/abstract_backend.py @@ -1030,3 +1030,17 @@ def cholesky(self, Tuple[Tensor, Tensor]: raise NotImplementedError( f"Backend {self.name} has not implemented cholesky.") + + def eps(self, dtype: Type[np.number]) -> float: + """ + Return machine epsilon for given `dtype` + + Args: + dtype: A dtype. + + Returns: + float: Machine epsilon. + """ + + raise NotImplementedError( + f"Backend {self.name} has not implemented eps.") diff --git a/tensornetwork/backends/jax/jax_backend.py b/tensornetwork/backends/jax/jax_backend.py index 9f5eb9843..f1d0ed0f3 100644 --- a/tensornetwork/backends/jax/jax_backend.py +++ b/tensornetwork/backends/jax/jax_backend.py @@ -888,3 +888,15 @@ def power(self, a: Tensor, b: Union[Tensor, float]) -> Tensor: b: The tensor that contains the exponent or a single scalar. """ return jnp.power(a, b) + + def eps(self, dtype: Type[np.number]) -> float: + """ + Return machine epsilon for given `dtype` + + Args: + dtype: A dtype. + + Returns: + float: Machine epsilon. + """ + return jnp.finfo(dtype).eps diff --git a/tensornetwork/backends/numpy/numpy_backend.py b/tensornetwork/backends/numpy/numpy_backend.py index 40b79fdb1..2d3753911 100644 --- a/tensornetwork/backends/numpy/numpy_backend.py +++ b/tensornetwork/backends/numpy/numpy_backend.py @@ -762,24 +762,36 @@ def deserialize_tensor(self, s: str) -> Tensor: def power(self, a: Tensor, b: Union[Tensor, float]) -> Tensor: """ - Returns the exponentiation of tensor a raised to b. - If b is a tensor, then the exponentiation is element-wise + Returns the exponentiation of tensor a raised to b. + If b is a tensor, then the exponentiation is element-wise between the two tensors, with a as the base and b as the power. - Note that a and b must be broadcastable to the same shape if + Note that a and b must be broadcastable to the same shape if b is a tensor. If b is a scalar, then the exponentiation is each value in a raised to the power of b. - + Args: a: The tensor containing the bases. b: The tensor containing the powers; or a single scalar as the power. Returns: - The tensor that is each element of a raised to the + The tensor that is each element of a raised to the power of b. Note that the shape of the returned tensor is that produced by the broadcast of a and b. """ return np.power(a, b) - + def item(self, tensor): return tensor.item() + + def eps(self, dtype: Type[np.number]) -> float: + """ + Return machine epsilon for given `dtype` + + Args: + dtype: A dtype. + + Returns: + float: Machine epsilon. + """ + return np.finfo(dtype).eps diff --git a/tensornetwork/backends/pytorch/pytorch_backend.py b/tensornetwork/backends/pytorch/pytorch_backend.py index 8467fac39..e0b1dc743 100644 --- a/tensornetwork/backends/pytorch/pytorch_backend.py +++ b/tensornetwork/backends/pytorch/pytorch_backend.py @@ -472,6 +472,18 @@ def sign(self, tensor: Tensor) -> Tensor: tensor: The input tensor. """ return torchlib.sign(tensor) - + def item(self, tensor): return tensor.item() + + def eps(self, dtype: Type[np.number]) -> float: + """ + Return machine epsilon for given `dtype` + + Args: + dtype: A dtype. + + Returns: + float: Machine epsilon. + """ + return torchlib.finfo(dtype).eps diff --git a/tensornetwork/backends/symmetric/symmetric_backend.py b/tensornetwork/backends/symmetric/symmetric_backend.py index 9e55bc694..8cedc6d2d 100644 --- a/tensornetwork/backends/symmetric/symmetric_backend.py +++ b/tensornetwork/backends/symmetric/symmetric_backend.py @@ -689,3 +689,15 @@ def matmul(self, tensor1: Tensor, tensor2: Tensor): if (tensor1.ndim != 2) or (tensor2.ndim != 2): raise ValueError("inputs to `matmul` have to be matrices") return tensor1 @ tensor2 + + def eps(self, dtype: Type[numpy.number]) -> float: + """ + Return machine epsilon for given `dtype` + + Args: + dtype: A dtype. + + Returns: + float: Machine epsilon. + """ + return numpy.finfo(dtype).eps diff --git a/tensornetwork/backends/tensorflow/tensorflow_backend.py b/tensornetwork/backends/tensorflow/tensorflow_backend.py index a46591c37..0cf960838 100644 --- a/tensornetwork/backends/tensorflow/tensorflow_backend.py +++ b/tensornetwork/backends/tensorflow/tensorflow_backend.py @@ -418,3 +418,15 @@ def power(self, a: Tensor, b: Union[Tensor, float]) -> Tensor: is that produced by the broadcast of a and b. """ return tf.math.pow(a, b) + + def eps(self, dtype: Type[np.number]) -> float: + """ + Return machine epsilon for given `dtype` + + Args: + dtype: A dtype. + + Returns: + float: Machine epsilon. + """ + return tf.experimental.numpy.finfo(dtype).eps From 8d680c2c41594b06cc83cb664c5bc3313669b4d1 Mon Sep 17 00:00:00 2001 From: mganahl Date: Tue, 12 Jan 2021 14:32:33 +0100 Subject: [PATCH 2/4] added backend.eps + test --- tensornetwork/backends/jax/jax_backend_test.py | 6 +++++- tensornetwork/backends/numpy/numpy_backend_test.py | 5 +++++ tensornetwork/backends/pytorch/pytorch_backend_test.py | 5 +++++ tensornetwork/backends/symmetric/symmetric_backend_test.py | 5 +++++ .../backends/tensorflow/tensorflow_backend_test.py | 5 +++++ 5 files changed, 25 insertions(+), 1 deletion(-) diff --git a/tensornetwork/backends/jax/jax_backend_test.py b/tensornetwork/backends/jax/jax_backend_test.py index fb87b3224..a3d476710 100644 --- a/tensornetwork/backends/jax/jax_backend_test.py +++ b/tensornetwork/backends/jax/jax_backend_test.py @@ -1254,4 +1254,8 @@ def test_power(dtype): actual = backend.power(base_tensor, power) expected = jax.numpy.power(base_tensor, power) np.testing.assert_allclose(expected, actual) - + +@pytest.mark.parametrize("dtype", np_dtypes) +def test_eps(dtype): + backend = jax_backend.JaxBackend() + assert backend.eps(dtype) == np.finfo(dtype).eps diff --git a/tensornetwork/backends/numpy/numpy_backend_test.py b/tensornetwork/backends/numpy/numpy_backend_test.py index 61533a9d4..990e349bc 100644 --- a/tensornetwork/backends/numpy/numpy_backend_test.py +++ b/tensornetwork/backends/numpy/numpy_backend_test.py @@ -969,3 +969,8 @@ def test_item(dtype): backend = numpy_backend.NumPyBackend() tensor = backend.randn((1,), dtype=dtype, seed=10) assert tensor.item() == backend.item(tensor) + +@pytest.mark.parametrize("dtype", np_dtypes) +def test_eps(dtype): + backend = numpy_backend.NumPyBackend() + assert backend.eps(dtype) == np.finfo(dtype).eps diff --git a/tensornetwork/backends/pytorch/pytorch_backend_test.py b/tensornetwork/backends/pytorch/pytorch_backend_test.py index f44aed5c6..d5fbade5f 100644 --- a/tensornetwork/backends/pytorch/pytorch_backend_test.py +++ b/tensornetwork/backends/pytorch/pytorch_backend_test.py @@ -670,3 +670,8 @@ def test_item(dtype): backend = pytorch_backend.PyTorchBackend() tensor = backend.randn((1,), dtype=dtype, seed=10) assert backend.item(tensor) == tensor.item() + +@pytest.mark.parametrize("dtype", torch_randn_dtypes) +def test_eps(dtype): + backend = pytorch_backend.PyTorchBackend() + assert backend.eps(dtype) == torch.finfo(dtype).eps diff --git a/tensornetwork/backends/symmetric/symmetric_backend_test.py b/tensornetwork/backends/symmetric/symmetric_backend_test.py index a8b78051e..033164236 100644 --- a/tensornetwork/backends/symmetric/symmetric_backend_test.py +++ b/tensornetwork/backends/symmetric/symmetric_backend_test.py @@ -1715,3 +1715,8 @@ def test_matmul_raises(): B = BlockSparseTensor.random(indices=inds2, dtype=dtype) with pytest.raises(ValueError, match="inputs to"): _ = backend.matmul(A, B) + +@pytest.mark.parametrize("dtype", np_dtypes) +def test_eps(dtype): + backend = symmetric_backend.SymmetricBackend() + assert backend.eps(dtype) == np.finfo(dtype).eps diff --git a/tensornetwork/backends/tensorflow/tensorflow_backend_test.py b/tensornetwork/backends/tensorflow/tensorflow_backend_test.py index c9dd1aa43..134c21bfd 100644 --- a/tensornetwork/backends/tensorflow/tensorflow_backend_test.py +++ b/tensornetwork/backends/tensorflow/tensorflow_backend_test.py @@ -624,3 +624,8 @@ def test_power(dtype): actual = backend.power(base_tensor, power) expected = tf.math.pow(base_tensor, power) np.testing.assert_allclose(expected, actual) + +@pytest.mark.parametrize("dtype", tf_dtypes) +def test_eps(dtype): + backend = tensorflow_backend.TensorFlowBackend() + assert backend.eps(dtype) == tf.experimental.numpy.finfo(dtype).eps From 6e216801b9e9b734f40284ea5e93501683c28979 Mon Sep 17 00:00:00 2001 From: mganahl Date: Wed, 20 Jan 2021 07:34:37 +0100 Subject: [PATCH 3/4] bug fix: for simulations with negative small schmidt values due to finite precision arithmetic, sqrt operations resulted in NaNs. This change fixes this. --- .../matrixproductstates/infinite_mps.py | 38 ++++++++++++------- 1 file changed, 24 insertions(+), 14 deletions(-) diff --git a/tensornetwork/matrixproductstates/infinite_mps.py b/tensornetwork/matrixproductstates/infinite_mps.py index 8cfc33a6b..54408f407 100644 --- a/tensornetwork/matrixproductstates/infinite_mps.py +++ b/tensornetwork/matrixproductstates/infinite_mps.py @@ -125,6 +125,7 @@ def transfer_matrix_eigs(self, """ D = self.bond_dimensions[0] + def mv(vector): result = self.unit_cell_transfer_operator( direction, self.backend.reshape(vector, (D, D))) @@ -137,10 +138,20 @@ def mv(vector): initial_state = self.backend.reshape(initial_state, (self.bond_dimensions[0]**2,)) - #note: for real dtype eta and dens are real. - #but scipy.linalg.eigs returns complex dtypes in any case - #since we know that for an MPS transfer matrix the largest - #eigenvalue and corresponding eigenvector are real + if D == 1: + # special case of non boundary entanglement + Z = self.backend.norm(initial_state) + initial_state = initial_state/Z + result = mv(initial_state) + eigval = self.backend.norm(result) + result = self.backend.reshape( + result, (self.bond_dimensions[0], self.bond_dimensions[0])) + return eigval, result + + # note: for real dtype eta and dens are real. + # but scipy.linalg.eigs returns complex dtypes in any case + # since we know that for an MPS transfer matrix the largest + # eigenvalue and corresponding eigenvector are real # we cast them. eta, dens = self.backend.eigs( A=mv, @@ -197,6 +208,9 @@ def canonicalize(self, Returns: None """ + if pseudo_inverse_cutoff is None: + pseudo_inverse_cutoff = self.backend.eps(self.dtype) + if self.center_position is None: self.center_position = 0 @@ -220,12 +234,11 @@ def canonicalize(self, # eigvals_left and u_left are both `Tensor` objects eigvals_left, u_left = self.backend.eigh(l) eigvals_left /= self.backend.norm(eigvals_left) - if pseudo_inverse_cutoff: - mask = eigvals_left <= pseudo_inverse_cutoff + mask = eigvals_left <= pseudo_inverse_cutoff + eigvals_left = self.backend.index_update(eigvals_left, mask, 0.0) inveigvals_left = 1.0 / eigvals_left - if pseudo_inverse_cutoff: - inveigvals_left = self.backend.index_update(inveigvals_left, mask, 0.0) + inveigvals_left = self.backend.index_update(inveigvals_left, mask, 0.0) sqrtl = ncon( [u_left, self.backend.diagflat(self.backend.sqrt(eigvals_left))], @@ -249,18 +262,15 @@ def canonicalize(self, # eigvals_right and u_right are both `Tensor` objects eigvals_right, u_right = self.backend.eigh(r) eigvals_right /= self.backend.norm(eigvals_right) - if pseudo_inverse_cutoff: - mask = eigvals_right <= pseudo_inverse_cutoff + mask = eigvals_right <= pseudo_inverse_cutoff + eigvals_right = self.backend.index_update(eigvals_right, mask, 0.0) inveigvals_right = 1.0 / eigvals_right - if pseudo_inverse_cutoff: - inveigvals_right = self.backend.index_update(inveigvals_right, mask, 0.0) - + inveigvals_right = self.backend.index_update(inveigvals_right, mask, 0.0) sqrtr = ncon( [u_right, self.backend.diagflat(self.backend.sqrt(eigvals_right))], [[-1, 1], [1, -2]], backend=self.backend.name) - inv_sqrtr = ncon([ self.backend.diagflat(self.backend.sqrt(inveigvals_right)), self.backend.conj(u_right) From 657b57037c81ae6d70760377885adf5bce4a95b0 Mon Sep 17 00:00:00 2001 From: mganahl Date: Wed, 20 Jan 2021 07:54:15 +0100 Subject: [PATCH 4/4] linting --- tensornetwork/matrixproductstates/infinite_mps.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tensornetwork/matrixproductstates/infinite_mps.py b/tensornetwork/matrixproductstates/infinite_mps.py index 54408f407..ff7d4ede5 100644 --- a/tensornetwork/matrixproductstates/infinite_mps.py +++ b/tensornetwork/matrixproductstates/infinite_mps.py @@ -141,11 +141,11 @@ def mv(vector): if D == 1: # special case of non boundary entanglement Z = self.backend.norm(initial_state) - initial_state = initial_state/Z + initial_state = initial_state / Z result = mv(initial_state) eigval = self.backend.norm(result) result = self.backend.reshape( - result, (self.bond_dimensions[0], self.bond_dimensions[0])) + result, (self.bond_dimensions[0], self.bond_dimensions[0])) return eigval, result # note: for real dtype eta and dens are real.