From 457bcfcf96b585e25384c5153f48fa1b183a65d5 Mon Sep 17 00:00:00 2001 From: Renato Mello Date: Fri, 4 Oct 2024 10:07:59 +0400 Subject: [PATCH 01/15] numpy --- src/qibo/backends/numpy.py | 21 ++++++++++++++++++++- 1 file changed, 20 insertions(+), 1 deletion(-) diff --git a/src/qibo/backends/numpy.py b/src/qibo/backends/numpy.py index fe1441b7c6..7161086c0c 100644 --- a/src/qibo/backends/numpy.py +++ b/src/qibo/backends/numpy.py @@ -769,12 +769,31 @@ def calculate_matrix_exp(self, a, matrix, eigenvectors=None, eigenvalues=None): ud = self.np.transpose(np.conj(eigenvectors)) return self.np.matmul(eigenvectors, self.np.matmul(expd, ud)) - def calculate_matrix_power(self, matrix, power: Union[float, int]): + def calculate_matrix_power( + self, + matrix, + power: Union[float, int], + precision_singularity: float = 1e-14, + ): if not isinstance(power, (float, int)): raise_error( TypeError, f"``power`` must be either float or int, but it is type {type(power)}.", ) + + if power < 0.0: + # negative powers of singular matrices via SVD + # det with np instead of self.np because of torch + determinant = np.linalg.det(matrix) + if abs(determinant) < precision_singularity: + U, S, Vh = self.np.linalg.svd(matrix) + S_inv = self.np.where( + self.np.abs(S) < precision_singularity, 0.0, S**power + ) + return ( + self.np.linalg.inv(Vh) @ self.np.diag(S_inv) @ self.np.linalg.inv(U) + ) + return fractional_matrix_power(matrix, power) # TODO: remove this method From b7a812933bc6ec16d8a9b4ad420a7acee05c91bf Mon Sep 17 00:00:00 2001 From: Renato Mello Date: Fri, 4 Oct 2024 10:08:03 +0400 Subject: [PATCH 02/15] torch --- src/qibo/backends/pytorch.py | 26 ++++++++++++++++++++++++-- 1 file changed, 24 insertions(+), 2 deletions(-) diff --git a/src/qibo/backends/pytorch.py b/src/qibo/backends/pytorch.py index 22fbbc5640..168a8f3841 100644 --- a/src/qibo/backends/pytorch.py +++ b/src/qibo/backends/pytorch.py @@ -1,5 +1,7 @@ """PyTorch backend.""" +from typing import Union + import numpy as np from qibo import __version__ @@ -199,9 +201,29 @@ def calculate_matrix_exp(self, a, matrix, eigenvectors=None, eigenvalues=None): ud = self.np.conj(eigenvectors).T return self.np.matmul(eigenvectors, self.np.matmul(expd, ud)) - def calculate_matrix_power(self, matrix, power): + def calculate_matrix_power( + self, + matrix, + power: Union[float, int], + precision_singularity: float = 1e-14, + ): + if power < 0.0: + # negative powers of singular matrices via SVD + determinant = self.np.linalg.det(matrix) + if abs(determinant) < precision_singularity: + U, S, Vh = self.np.linalg.svd(matrix) + # cast needed because of different dtypes + S = self.cast(S) + S_inv = self.np.where( + self.np.abs(S) < precision_singularity, 0.0, S**power + ) + + return ( + self.np.linalg.inv(Vh) @ self.np.diag(S_inv) @ self.np.linalg.inv(U) + ) + copied = self.to_numpy(self.np.copy(matrix)) - copied = super().calculate_matrix_power(copied, power) + copied = super().calculate_matrix_power(copied, power, precision_singularity) return self.cast(copied, dtype=copied.dtype) def _test_regressions(self, name): From bd82ddbe8fb556c208ddfae6d5494239cc4282c7 Mon Sep 17 00:00:00 2001 From: Renato Mello Date: Fri, 4 Oct 2024 10:08:09 +0400 Subject: [PATCH 03/15] tensorflow --- src/qibo/backends/tensorflow.py | 26 ++++++++++++++++++++++++++ 1 file changed, 26 insertions(+) diff --git a/src/qibo/backends/tensorflow.py b/src/qibo/backends/tensorflow.py index 2088cd69be..702d993cbe 100644 --- a/src/qibo/backends/tensorflow.py +++ b/src/qibo/backends/tensorflow.py @@ -1,5 +1,6 @@ import collections import os +from typing import Union import numpy as np @@ -192,6 +193,31 @@ def calculate_matrix_exp(self, a, matrix, eigenvectors=None, eigenvalues=None): return self.tf.linalg.expm(-1j * a * matrix) return super().calculate_matrix_exp(a, matrix, eigenvectors, eigenvalues) + def calculate_matrix_power( + self, + matrix, + power: Union[float, int], + precision_singularity: float = 1e-14, + ): + if not isinstance(power, (float, int)): + raise_error( + TypeError, + f"``power`` must be either float or int, but it is type {type(power)}.", + ) + + if power < 0.0: + # negative powers of singular matrices via SVD + determinant = self.tf.linalg.det(matrix) + if abs(determinant) < precision_singularity: + S, U, V = self.tf.linalg.svd(matrix) + S_inv = self.tf.where( + self.tf.abs(S) < precision_singularity, 0.0, S**power + ) + + return V @ self.np.diag(S_inv) @ self.np.linalg.inv(U) + + return super().calculate_matrix_power(matrix, power, precision_singularity) + def calculate_hamiltonian_matrix_product(self, matrix1, matrix2): if self.is_sparse(matrix1) or self.is_sparse(matrix2): raise_error( From bce842416e0b32b51d81d38224bd63bab04b1ccd Mon Sep 17 00:00:00 2001 From: Renato Mello Date: Fri, 4 Oct 2024 10:08:21 +0400 Subject: [PATCH 04/15] quantum_info --- src/qibo/quantum_info/linalg_operations.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/src/qibo/quantum_info/linalg_operations.py b/src/qibo/quantum_info/linalg_operations.py index 6fc48e6e74..23c6020c66 100644 --- a/src/qibo/quantum_info/linalg_operations.py +++ b/src/qibo/quantum_info/linalg_operations.py @@ -198,12 +198,17 @@ def matrix_exponentiation( return backend.calculate_matrix_exp(phase, matrix, eigenvectors, eigenvalues) -def matrix_power(matrix, power: Union[float, int], backend=None): +def matrix_power( + matrix, power: Union[float, int], precision_singularity: float = 1e-14, backend=None +): """Given a ``matrix`` :math:`A` and power :math:`\\alpha`, calculate :math:`A^{\\alpha}`. Args: matrix (ndarray): matrix whose power to calculate. power (float or int): power to raise ``matrix`` to. + precision_singularity (float, optional): If determinant of ``matrix`` is smaller than + ``precision_singularity``, then matrix is considered to be singular. + Used when ``power`` is negative. backend (:class:`qibo.backends.abstract.Backend`, optional): backend to be used in the execution. If ``None``, it uses :class:`qibo.backends.GlobalBackend`. Defaults to ``None``. @@ -213,4 +218,4 @@ def matrix_power(matrix, power: Union[float, int], backend=None): """ backend = _check_backend(backend) - return backend.calculate_matrix_power(matrix, power) + return backend.calculate_matrix_power(matrix, power, precision_singularity) From 99f1f8b7e6a26d3fe1cd40b58610a24464a3abe7 Mon Sep 17 00:00:00 2001 From: Renato Mello Date: Fri, 4 Oct 2024 10:08:29 +0400 Subject: [PATCH 05/15] tests --- tests/test_quantum_info_operations.py | 19 ++++++++++++++++--- 1 file changed, 16 insertions(+), 3 deletions(-) diff --git a/tests/test_quantum_info_operations.py b/tests/test_quantum_info_operations.py index 3d97f1f63a..618ec373b5 100644 --- a/tests/test_quantum_info_operations.py +++ b/tests/test_quantum_info_operations.py @@ -1,5 +1,6 @@ import numpy as np import pytest +from scipy.linalg import sqrtm from qibo import Circuit, gates, matrices from qibo.quantum_info.linalg_operations import ( @@ -113,16 +114,28 @@ def test_partial_trace(backend, density_matrix): backend.assert_allclose(traced, Id) -@pytest.mark.parametrize("power", [2, 2.0, "2"]) -def test_matrix_power(backend, power): +@pytest.mark.parametrize("singular", [False, True]) +@pytest.mark.parametrize("power", [-0.5, 0.5, 2, 2.0, "2"]) +def test_matrix_power(backend, power, singular): nqubits = 2 dims = 2**nqubits - state = random_density_matrix(dims, backend=backend) + state = random_density_matrix(dims, pure=singular, backend=backend) if isinstance(power, str): with pytest.raises(TypeError): test = matrix_power(state, power, backend) + elif power == -0.5 and singular: + # When the singular matrix is a state, this power should be itself + backend.assert_allclose(matrix_power(state, power, backend), state) + elif abs(power) == 0.5 and not singular: + # Should be equal to the (inverse) square root + sqrt = sqrtm(backend.to_numpy(state)).astype(complex) + if power == -0.5: + sqrt = np.linalg.inv(sqrt) + sqrt = backend.cast(sqrt) + + backend.assert_allclose(matrix_power(state, power, backend), sqrt) else: power = matrix_power(state, power, backend) From ecb1642c67ae9fbc000bd5aca034843169e80e56 Mon Sep 17 00:00:00 2001 From: Renato Mello Date: Fri, 4 Oct 2024 10:49:21 +0400 Subject: [PATCH 06/15] fix test --- tests/test_quantum_info_operations.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/tests/test_quantum_info_operations.py b/tests/test_quantum_info_operations.py index 618ec373b5..f61328842e 100644 --- a/tests/test_quantum_info_operations.py +++ b/tests/test_quantum_info_operations.py @@ -124,10 +124,10 @@ def test_matrix_power(backend, power, singular): if isinstance(power, str): with pytest.raises(TypeError): - test = matrix_power(state, power, backend) + test = matrix_power(state, power, backend=backend) elif power == -0.5 and singular: # When the singular matrix is a state, this power should be itself - backend.assert_allclose(matrix_power(state, power, backend), state) + backend.assert_allclose(matrix_power(state, power, backend=backend), state) elif abs(power) == 0.5 and not singular: # Should be equal to the (inverse) square root sqrt = sqrtm(backend.to_numpy(state)).astype(complex) @@ -135,9 +135,9 @@ def test_matrix_power(backend, power, singular): sqrt = np.linalg.inv(sqrt) sqrt = backend.cast(sqrt) - backend.assert_allclose(matrix_power(state, power, backend), sqrt) + backend.assert_allclose(matrix_power(state, power, backend=backend), sqrt) else: - power = matrix_power(state, power, backend) + power = matrix_power(state, power, backend=backend) backend.assert_allclose( float(backend.np.real(backend.np.trace(power))), From 800564ac8da2dc1ee92cd5a4e768d6b11f863323 Mon Sep 17 00:00:00 2001 From: Renato Mello Date: Fri, 4 Oct 2024 10:51:49 +0400 Subject: [PATCH 07/15] fix tests --- src/qibo/quantum_info/entropies.py | 12 ++++++------ tests/test_quantum_info_entropies.py | 10 +++++----- 2 files changed, 11 insertions(+), 11 deletions(-) diff --git a/src/qibo/quantum_info/entropies.py b/src/qibo/quantum_info/entropies.py index 204884c6ef..e0c2596a56 100644 --- a/src/qibo/quantum_info/entropies.py +++ b/src/qibo/quantum_info/entropies.py @@ -772,7 +772,7 @@ def renyi_entropy(state, alpha: Union[float, int], base: float = 2, backend=None / np.log2(base) ) - log = backend.np.log2(backend.np.trace(matrix_power(state, alpha, backend))) + log = backend.np.log2(backend.np.trace(matrix_power(state, alpha, backend=backend))) return (1 / (1 - alpha)) * log / np.log2(base) @@ -871,8 +871,8 @@ def relative_renyi_entropy( return relative_von_neumann_entropy(state, target, base, backend=backend) if alpha == np.inf: - new_state = matrix_power(state, 0.5, backend) - new_target = matrix_power(target, 0.5, backend) + new_state = matrix_power(state, 0.5, backend=backend) + new_target = matrix_power(target, 0.5, backend=backend) log = backend.np.log2( backend.calculate_norm_density_matrix(new_state @ new_target, order=1) @@ -880,8 +880,8 @@ def relative_renyi_entropy( return -2 * log / np.log2(base) - log = matrix_power(state, alpha, backend) - log = log @ matrix_power(target, 1 - alpha, backend) + log = matrix_power(state, alpha, backend=backend) + log = log @ matrix_power(target, 1 - alpha, backend=backend) log = backend.np.log2(backend.np.trace(log)) return (1 / (alpha - 1)) * log / np.log2(base) @@ -939,7 +939,7 @@ def tsallis_entropy(state, alpha: float, base: float = 2, backend=None): return von_neumann_entropy(state, base=base, backend=backend) return (1 / (1 - alpha)) * ( - backend.np.trace(matrix_power(state, alpha, backend)) - 1 + backend.np.trace(matrix_power(state, alpha, backend=backend)) - 1 ) diff --git a/tests/test_quantum_info_entropies.py b/tests/test_quantum_info_entropies.py index 1f19fdb7ed..b955608c39 100644 --- a/tests/test_quantum_info_entropies.py +++ b/tests/test_quantum_info_entropies.py @@ -686,8 +686,8 @@ def test_relative_renyi_entropy(backend, alpha, base, state_flag, target_flag): if target_flag else target ) - new_state = matrix_power(state_outer, 0.5, backend) - new_target = matrix_power(target_outer, 0.5, backend) + new_state = matrix_power(state_outer, 0.5, backend=backend) + new_target = matrix_power(target_outer, 0.5, backend=backend) log = backend.np.log2( backend.calculate_norm_density_matrix( @@ -703,8 +703,8 @@ def test_relative_renyi_entropy(backend, alpha, base, state_flag, target_flag): if len(target.shape) == 1: target = backend.np.outer(target, backend.np.conj(target)) - log = matrix_power(state, alpha, backend) - log = log @ matrix_power(target, 1 - alpha, backend) + log = matrix_power(state, alpha, backend=backend) + log = log @ matrix_power(target, 1 - alpha, backend=backend) log = backend.np.log2(backend.np.trace(log)) log = (1 / (alpha - 1)) * log / np.log2(base) @@ -750,7 +750,7 @@ def test_tsallis_entropy(backend, alpha, base): target = von_neumann_entropy(state, base=base, backend=backend) else: target = (1 / (1 - alpha)) * ( - backend.np.trace(matrix_power(state, alpha, backend)) - 1 + backend.np.trace(matrix_power(state, alpha, backend=backend)) - 1 ) backend.assert_allclose( From dc4d127aacd714457f28c8682be84e121d6e09ca Mon Sep 17 00:00:00 2001 From: Renato Mello Date: Fri, 4 Oct 2024 12:50:30 +0400 Subject: [PATCH 08/15] add arg to abstract method --- src/qibo/backends/abstract.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/qibo/backends/abstract.py b/src/qibo/backends/abstract.py index 468fa89008..9728069a2f 100644 --- a/src/qibo/backends/abstract.py +++ b/src/qibo/backends/abstract.py @@ -358,7 +358,7 @@ def calculate_matrix_exp( @abc.abstractmethod def calculate_matrix_power( - self, matrix, power: Union[float, int] + self, matrix, power: Union[float, int], precision_singularity: float = 1e-14 ): # pragma: no cover """Calculate the (fractional) ``power`` :math:`\\alpha` of ``matrix`` :math:`A`, i.e. :math:`A^{\\alpha}`. From a56b9e253f484a858de8209bc271915c9c9d2792 Mon Sep 17 00:00:00 2001 From: Renato Mello Date: Fri, 4 Oct 2024 15:18:39 +0400 Subject: [PATCH 09/15] cleaner torch implementation --- src/qibo/backends/numpy.py | 4 +++- src/qibo/backends/pytorch.py | 18 ++---------------- 2 files changed, 5 insertions(+), 17 deletions(-) diff --git a/src/qibo/backends/numpy.py b/src/qibo/backends/numpy.py index 7161086c0c..0bf43e6be1 100644 --- a/src/qibo/backends/numpy.py +++ b/src/qibo/backends/numpy.py @@ -784,9 +784,11 @@ def calculate_matrix_power( if power < 0.0: # negative powers of singular matrices via SVD # det with np instead of self.np because of torch - determinant = np.linalg.det(matrix) + determinant = self.np.linalg.det(matrix) if abs(determinant) < precision_singularity: U, S, Vh = self.np.linalg.svd(matrix) + # cast needed because of different dtypes in `torch` + S = self.cast(S) S_inv = self.np.where( self.np.abs(S) < precision_singularity, 0.0, S**power ) diff --git a/src/qibo/backends/pytorch.py b/src/qibo/backends/pytorch.py index 168a8f3841..5392a2583a 100644 --- a/src/qibo/backends/pytorch.py +++ b/src/qibo/backends/pytorch.py @@ -207,22 +207,8 @@ def calculate_matrix_power( power: Union[float, int], precision_singularity: float = 1e-14, ): - if power < 0.0: - # negative powers of singular matrices via SVD - determinant = self.np.linalg.det(matrix) - if abs(determinant) < precision_singularity: - U, S, Vh = self.np.linalg.svd(matrix) - # cast needed because of different dtypes - S = self.cast(S) - S_inv = self.np.where( - self.np.abs(S) < precision_singularity, 0.0, S**power - ) - - return ( - self.np.linalg.inv(Vh) @ self.np.diag(S_inv) @ self.np.linalg.inv(U) - ) - - copied = self.to_numpy(self.np.copy(matrix)) + copied = self.cast(matrix, copy=True) + copied = self.to_numpy(copied) if power >= 0.0 else copied.detach() copied = super().calculate_matrix_power(copied, power, precision_singularity) return self.cast(copied, dtype=copied.dtype) From 791be66f562ba78e52b36ed97b46ac1eb11452eb Mon Sep 17 00:00:00 2001 From: Renato Mello Date: Fri, 4 Oct 2024 15:45:04 +0400 Subject: [PATCH 10/15] remove comment --- src/qibo/backends/numpy.py | 1 - 1 file changed, 1 deletion(-) diff --git a/src/qibo/backends/numpy.py b/src/qibo/backends/numpy.py index 0bf43e6be1..0d1cbd7515 100644 --- a/src/qibo/backends/numpy.py +++ b/src/qibo/backends/numpy.py @@ -783,7 +783,6 @@ def calculate_matrix_power( if power < 0.0: # negative powers of singular matrices via SVD - # det with np instead of self.np because of torch determinant = self.np.linalg.det(matrix) if abs(determinant) < precision_singularity: U, S, Vh = self.np.linalg.svd(matrix) From f913aca22be2b6196f8f9375570b4edd0d995490 Mon Sep 17 00:00:00 2001 From: Renato Mello Date: Tue, 8 Oct 2024 15:02:27 +0400 Subject: [PATCH 11/15] docstring --- src/qibo/quantum_info/linalg_operations.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/qibo/quantum_info/linalg_operations.py b/src/qibo/quantum_info/linalg_operations.py index b380e59d1e..625051a9f3 100644 --- a/src/qibo/quantum_info/linalg_operations.py +++ b/src/qibo/quantum_info/linalg_operations.py @@ -319,7 +319,8 @@ def singular_value_decomposition(matrix, backend=None): :class:`qibo.backends.GlobalBackend`. Defaults to ``None``. Returns: - ndarray, ndarray, ndarray: Singular value decomposition of :math:`A`. + ndarray, ndarray, ndarray: Singular value decomposition of :math:`A`, i.e. + :math:`U`, :math:`S`, and :math:`V^{\\dagger}`, in that order. """ backend = _check_backend(backend) From 879d1a263f4bf54773fe9264c437ddb60e795131 Mon Sep 17 00:00:00 2001 From: Renato Mello Date: Tue, 8 Oct 2024 15:44:23 +0400 Subject: [PATCH 12/15] fix bug after merge --- src/qibo/quantum_info/entanglement.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/qibo/quantum_info/entanglement.py b/src/qibo/quantum_info/entanglement.py index 2eb670dc62..dbd2cb4508 100644 --- a/src/qibo/quantum_info/entanglement.py +++ b/src/qibo/quantum_info/entanglement.py @@ -148,9 +148,9 @@ def negativity(state, bipartition, backend=None): reduced = partial_transpose(state, bipartition, backend) reduced = backend.np.conj(reduced.T) @ reduced - norm = backend.np.trace(matrix_power(reduced, 1 / 2, backend)) + norm = backend.np.trace(matrix_power(reduced, 1 / 2, backend=backend)) - return float(backend.np.real((norm - 1) / 2)) + return backend.np.real((norm - 1) / 2) def entanglement_fidelity( From d1046bab7671ba786349225510e81b7768a31fad Mon Sep 17 00:00:00 2001 From: Renato Mello Date: Wed, 9 Oct 2024 09:09:47 +0400 Subject: [PATCH 13/15] raise_error --- src/qibo/backends/numpy.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/qibo/backends/numpy.py b/src/qibo/backends/numpy.py index 17726cdd49..4c1cdc8a33 100644 --- a/src/qibo/backends/numpy.py +++ b/src/qibo/backends/numpy.py @@ -415,8 +415,12 @@ def execute_circuit(self, circuit, initial_state=None, nshots=1000): if circuit.measurements or circuit.has_collapse: return self.execute_circuit_repeated(circuit, nshots, initial_state) else: - raise RuntimeError( - "Attempting to perform noisy simulation with `density_matrix=False` and no Measurement gate in the Circuit. If you wish to retrieve the statistics of the outcomes please include measurements in the circuit, otherwise set `density_matrix=True` to recover the final state." + raise_error( + RuntimeError, + "Attempting to perform noisy simulation with `density_matrix=False` " + + "and no Measurement gate in the Circuit. If you wish to retrieve the " + + "statistics of the outcomes please include measurements in the circuit, " + + "otherwise set `density_matrix=True` to recover the final state.", ) if circuit.accelerators: # pragma: no cover From 1b8ec29123555e746d3a8098ed7a33c87a18e4e1 Mon Sep 17 00:00:00 2001 From: Renato Mello Date: Wed, 9 Oct 2024 09:10:47 +0400 Subject: [PATCH 14/15] raise_error --- src/qibo/backends/numpy.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/qibo/backends/numpy.py b/src/qibo/backends/numpy.py index 4c1cdc8a33..69bb2f36e1 100644 --- a/src/qibo/backends/numpy.py +++ b/src/qibo/backends/numpy.py @@ -501,8 +501,11 @@ def execute_circuit_repeated(self, circuit, nshots, initial_state=None): and not circuit.measurements and not circuit.density_matrix ): - raise RuntimeError( - "The circuit contains only collapsing measurements (`collapse=True`) but `density_matrix=False`. Please set `density_matrix=True` to retrieve the final state after execution." + raise_error( + RuntimeError, + "The circuit contains only collapsing measurements (`collapse=True`) but " + + "`density_matrix=False`. Please set `density_matrix=True` to retrieve " + + "the final state after execution.", ) results, final_states = [], [] From 2358c3f2b39e13bd50c16402212ac39a92d5903c Mon Sep 17 00:00:00 2001 From: Renato Mello Date: Wed, 9 Oct 2024 09:30:14 +0400 Subject: [PATCH 15/15] Stavros' suggestion --- src/qibo/backends/numpy.py | 22 ++++++++++++++-------- src/qibo/backends/tensorflow.py | 9 +++------ 2 files changed, 17 insertions(+), 14 deletions(-) diff --git a/src/qibo/backends/numpy.py b/src/qibo/backends/numpy.py index 69bb2f36e1..21be02bfce 100644 --- a/src/qibo/backends/numpy.py +++ b/src/qibo/backends/numpy.py @@ -792,14 +792,8 @@ def calculate_matrix_power( # negative powers of singular matrices via SVD determinant = self.np.linalg.det(matrix) if abs(determinant) < precision_singularity: - U, S, Vh = self.np.linalg.svd(matrix) - # cast needed because of different dtypes in `torch` - S = self.cast(S) - S_inv = self.np.where( - self.np.abs(S) < precision_singularity, 0.0, S**power - ) - return ( - self.np.linalg.inv(Vh) @ self.np.diag(S_inv) @ self.np.linalg.inv(U) + return _calculate_negative_power_singular_matrix( + matrix, power, precision_singularity, self.np, self ) return fractional_matrix_power(matrix, power) @@ -847,3 +841,15 @@ def _test_regressions(self, name): {5: 18, 4: 5, 7: 4, 1: 2, 6: 1}, {4: 8, 2: 6, 5: 5, 1: 3, 3: 3, 6: 2, 7: 2, 0: 1}, ] + + +def _calculate_negative_power_singular_matrix( + matrix, power: Union[float, int], precision_singularity: float, engine, backend +): + """Calculate negative power of singular matrix.""" + U, S, Vh = backend.calculate_singular_value_decomposition(matrix) + # cast needed because of different dtypes in `torch` + S = backend.cast(S) + S_inv = engine.where(engine.abs(S) < precision_singularity, 0.0, S**power) + + return engine.linalg.inv(Vh) @ backend.np.diag(S_inv) @ engine.linalg.inv(U) diff --git a/src/qibo/backends/tensorflow.py b/src/qibo/backends/tensorflow.py index 816762ea52..215fbc8062 100644 --- a/src/qibo/backends/tensorflow.py +++ b/src/qibo/backends/tensorflow.py @@ -6,7 +6,7 @@ from qibo import __version__ from qibo.backends.npmatrices import NumpyMatrices -from qibo.backends.numpy import NumpyBackend +from qibo.backends.numpy import NumpyBackend, _calculate_negative_power_singular_matrix from qibo.config import TF_LOG_LEVEL, log, raise_error @@ -209,13 +209,10 @@ def calculate_matrix_power( # negative powers of singular matrices via SVD determinant = self.tf.linalg.det(matrix) if abs(determinant) < precision_singularity: - S, U, V = self.tf.linalg.svd(matrix) - S_inv = self.tf.where( - self.tf.abs(S) < precision_singularity, 0.0, S**power + return _calculate_negative_power_singular_matrix( + matrix, power, precision_singularity, self.tf, self ) - return V @ self.np.diag(S_inv) @ self.np.linalg.inv(U) - return super().calculate_matrix_power(matrix, power, precision_singularity) def calculate_singular_value_decomposition(self, matrix):