Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension


Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
346 changes: 263 additions & 83 deletions quanestimation/AdaptiveScheme/Adapt.py

Large diffs are not rendered by default.

46 changes: 22 additions & 24 deletions quanestimation/AdaptiveScheme/Adapt_MZI.py
Original file line number Diff line number Diff line change
Expand Up @@ -17,8 +17,8 @@ class Adapt_MZI:
-- Initial state (density matrix).

"""
def __init__(self, x, p, rho0):

def __init__(self, x, p, rho0):
self.x = x
self.p = p
self.rho0 = rho0
Expand All @@ -33,18 +33,16 @@ def online(self, target="sharpness", output="phi"):
Parameters
----------
> **target:** `string`
-- Setting the target function for calculating the tunable phase. Options are:
"sharpness" (default) -- Sharpness.
"MI" -- Mutual information.
-- Setting the target function for calculating the tunable phase. Options are:
"sharpness" (default) -- Sharpness.
"MI" -- Mutual information.

> **output:** `string`
-- The output the class. Options are:
"phi" (default) -- The tunable phase.
"dphi" -- Phase difference.
-- The output the class. Options are:
"phi" (default) -- The tunable phase.
"dphi" -- Phase difference.
"""
phi = QJL.adaptMZI_online(
self.x, self.p, self.rho0, output, target
)
phi = QJL.adaptMZI_online(self.x, self.p, self.rho0, output, target)

def offline(
self,
Expand All @@ -65,13 +63,13 @@ def offline(
Parameters
----------
> **target:** `string`
-- Setting the target function for calculating the tunable phase. Options are:
"sharpness" (default) -- Sharpness.
"MI" -- Mutual information.
-- Setting the target function for calculating the tunable phase. Options are:
"sharpness" (default) -- Sharpness.
"MI" -- Mutual information.

> **method:** `string`
-- The method for the adaptive phase estimation. Options are:
"DE" (default) -- DE algorithm for the adaptive phase estimation.
-- The method for the adaptive phase estimation. Options are:
"DE" (default) -- DE algorithm for the adaptive phase estimation.
"PSO" -- PSO algorithm for the adaptive phase estimation.

If the `method=DE`, the parameters are:
Expand All @@ -83,7 +81,7 @@ def offline(

> **max_episode:** `int`
-- The number of episodes.

> **c:** `float`
-- Mutation constant.

Expand All @@ -95,28 +93,28 @@ def offline(

> **eps:** `float`
-- Machine epsilon.

If the `method=PSO`, the parameters are:

> **deltaphi0:** `list`
-- Initial guesses of phase difference.

> **max_episode:** `int or list`
-- If it is an integer, for example max_episode=1000, it means the
-- If it is an integer, for example max_episode=1000, it means the
program will continuously run 1000 episodes. However, if it is an
array, for example max_episode=[1000,100], the program will run
1000 episodes in total but replace states of all the particles
array, for example max_episode=[1000,100], the program will run
1000 episodes in total but replace states of all the particles
with global best every 100 episodes.

> **c0:** `float`
-- The damping factor that assists convergence, also known as inertia weight.

> **c1:** `float`
-- The exploitation weight that attracts the particle to its best previous
-- The exploitation weight that attracts the particle to its best previous
position, also known as cognitive learning factor.

> **c2:** `float`
-- The exploitation weight that attracts the particle to the best position
-- The exploitation weight that attracts the particle to the best position
in the neighborhood, also known as social learning factor.

> **eps:** `float`
Expand All @@ -127,7 +125,7 @@ def offline(
np.array([int(list(comb_tp[i])[j]) for j in range(self.N)])
for i in range(2**self.N)
]

if method == "DE":
QJL.DE_deltaphiOpt(
self.x,
Expand Down
70 changes: 31 additions & 39 deletions quanestimation/AsymptoticBound/AnalogCramerRao.py
Original file line number Diff line number Diff line change
Expand Up @@ -13,31 +13,31 @@ def HCRB(rho, drho, W, eps=1e-8):
The HCRB is defined as:

$$
\min_{\{X_i\}} \left\{ \mathrm{Tr}(\mathrm{Re}Z) + \mathrm{Tr}(| \mathrm{Im} Z |) \right\},
\min_{\{X_i\}} \left\{ \mathrm{Tr}(\mathrm{Re}Z) + \mathrm{Tr}(| \mathrm{Im} Z |) \right\},
$$

where $Z_{ij} = \mathrm{Tr}(\rho X_i X_j)$ and $V$ is the covariance matrix.

Args:
rho (np.array):
rho (np.array):
Density matrix.
drho (list):
Derivatives of the density matrix with respect to unknown parameters.
drho (list):
Derivatives of the density matrix with respect to unknown parameters.
For example, `drho[0]` is the derivative with respect to the first parameter.
W (np.array):
W (np.array):
Weight matrix for the bound.
eps (float, optional):
eps (float, optional):
Machine epsilon for numerical stability.

Returns:
(float):
Returns:
(float):
The value of the Holevo Cramer-Rao bound.

Raises:
TypeError: If `drho` is not a list.

Notes:
In the single-parameter scenario, the HCRB is equivalent to the QFI.
In the single-parameter scenario, the HCRB is equivalent to the QFI.
For a rank-one weight matrix, the HCRB is equivalent to the inverse of the QFIM.
"""

Expand All @@ -50,7 +50,7 @@ def HCRB(rho, drho, W, eps=1e-8):
"Returning QFI value."
)
return QFIM(rho, drho, eps=eps)

if matrix_rank(W) == 1:
print(
"For rank-one weight matrix, HCRB is equivalent to QFIM. "
Expand All @@ -70,8 +70,7 @@ def HCRB(rho, drho, W, eps=1e-8):
vec_drho = []
for param_idx in range(num_params):
components = [
np.real(np.trace(drho[param_idx] @ basis_mat))
for basis_mat in basis
np.real(np.trace(drho[param_idx] @ basis_mat)) for basis_mat in basis
]
vec_drho.append(np.array(components))

Expand All @@ -91,12 +90,7 @@ def HCRB(rho, drho, W, eps=1e-8):
X = cp.Variable((num, num_params))

# Define constraints
constraints = [
cp.bmat([
[V, X.T @ R.conj().T],
[R @ X, np.identity(num)]
]) >> 0
]
constraints = [cp.bmat([[V, X.T @ R.conj().T], [R @ X, np.identity(num)]]) >> 0]

# Add linear constraints
for i in range(num_params):
Expand All @@ -113,65 +107,63 @@ def HCRB(rho, drho, W, eps=1e-8):

return problem.value


def NHB(rho, drho, W):
r"""
Calculation of the Nagaoka-Hayashi bound (NHB) via the semidefinite program (SDP).

The NHB is defined as:

$$
\min_{X} \left\{ \mathrm{Tr}[W \mathrm{Re}(Z)] + \|\sqrt{W} \mathrm{Im}(Z) \sqrt{W}\|_1 \right\},
\min_{X} \left\{ \mathrm{Tr}[W \mathrm{Re}(Z)] + \|\sqrt{W} \mathrm{Im}(Z) \sqrt{W}\|_1 \right\},
$$

where $Z_{ij} = \mathrm{Tr}(\rho X_i X_j)$ and $V$ is the covariance matrix.

Args:
rho (np.array):
rho (np.array):
Density matrix.
drho (list):
Derivatives of the density matrix with respect to unknown parameters.
drho (list):
Derivatives of the density matrix with respect to unknown parameters.
For example, `drho[0]` is the derivative with respect to the first parameter.
W (np.array):
W (np.array):
Weight matrix for the bound.

Returns:
(float):
Returns:
(float):
The value of the Nagaoka-Hayashi bound.

Raises:
TypeError:
TypeError:
If `drho` is not a list.
"""
if not isinstance(drho, list):
raise TypeError("drho must be a list of derivative matrices")

dim = len(rho)
num_params = len(drho)

# Initialize a temporary matrix for L components
L_temp = [[None for _ in range(num_params)] for _ in range(num_params)]

# Create Hermitian variables for the upper triangle and mirror to lower triangle
for i in range(num_params):
for j in range(i, num_params):
L_temp[i][j] = cp.Variable((dim, dim), hermitian=True)
if i != j:
L_temp[j][i] = L_temp[i][j]

# Construct the block matrix L
L_blocks = [cp.hstack(L_temp[i]) for i in range(num_params)]
L = cp.vstack(L_blocks)

# Create Hermitian variables for X
X = [cp.Variable((dim, dim), hermitian=True) for _ in range(num_params)]

# Construct the block matrix constraint
block_matrix = cp.bmat([
[L, cp.vstack(X)],
[cp.hstack(X), np.identity(dim)]
])
block_matrix = cp.bmat([[L, cp.vstack(X)], [cp.hstack(X), np.identity(dim)]])
constraints = [block_matrix >> 0]

# Add trace constraints
for i in range(num_params):
constraints.append(cp.trace(X[i] @ rho) == 0)
Expand All @@ -180,7 +172,7 @@ def NHB(rho, drho, W):
constraints.append(cp.trace(X[i] @ drho[j]) == 1)
else:
constraints.append(cp.trace(X[i] @ drho[j]) == 0)

# Define and solve the optimization problem
objective = cp.Minimize(cp.real(cp.trace(cp.kron(W, rho) @ L)))
prob = cp.Problem(objective, constraints)
Expand Down
Loading