Functions

This part is the functions of the Python-Julia package which written in Python.

Kraus

The parameterization of a state is \begin{align} \rho=\sum_i K_i\rho_0K_i^{\dagger}, \end{align}

where \(\rho\) is the evolved density matrix, \(K_i\) is the Kraus operator.

Parameters

K: list -- Kraus operators.

dK: list -- Derivatives of the Kraus operators with respect to the unknown parameters to be estimated. For example, dK[0] is the derivative vector on the first parameter.

rho0: matrix -- Initial state (density matrix).

Returns

Density matrix and its derivatives on the unknown parameters.

Source code in quanestimation/Parameterization/NonDynamics.py

def Kraus(rho0, K, dK):
    r"""
    The parameterization of a state is
    \begin{align}
    \rho=\sum_i K_i\rho_0K_i^{\dagger},
    \end{align} 

    where $\rho$ is the evolved density matrix, $K_i$ is the Kraus operator.

    Parameters
    ----------
    > **K:** `list`
        -- Kraus operators.

    > **dK:** `list`
        -- Derivatives of the Kraus operators with respect to the unknown parameters to be 
        estimated. For example, dK[0] is the derivative vector on the first 
        parameter.

    > **rho0:** `matrix`
        -- Initial state (density matrix).

    Returns
    ----------
    Density matrix and its derivatives on the unknown parameters.
    """

    k_num = len(K)
    para_num = len(dK[0])
    dK_reshape = [[dK[i][j] for i in range(k_num)] for j in range(para_num)]

    rho = sum([np.dot(Ki, np.dot(rho0, Ki.conj().T)) for Ki in K])
    drho = [sum([(np.dot(dKi, np.dot(rho0, Ki.conj().T))+ np.dot(Ki, np.dot(rho0, dKi.conj().T))) for (Ki, dKi) in zip(K, dKj)]) for dKj in dK_reshape]

    return rho, drho

---

Metrological resources

Calculation of spin squeezing parameter for a density matrix.

Parameters

rho: matrix -- Density matrix.

basis: string -- The basis of the state. Options are:
"Dicke" (default) -- Dicke basis.
"Pauli" -- The original basis of each spin.

output: string -- Types of spin squeezing can be calculated. Options are:
"KU" (default) -- Spin squeezing defined by Kitagawa and Ueda.
"WBIMH" -- Spin squeezing defined by Wineland et al.

Returns

\(\xi\): float -- spin squeezing parameter

Source code in quanestimation/Resource/Resource.py

def SpinSqueezing(rho, basis="Dicke", output="KU"):
    r"""
    Calculation of spin squeezing parameter for a density matrix.

    Parameters
    ----------
    > **rho:** `matrix`
        -- Density matrix.

    > **basis:** `string`
        -- The basis of the state. Options are:  
        "Dicke" (default) -- Dicke basis.  
        "Pauli" -- The original basis of each spin.

    > **output:** `string`
        -- Types of spin squeezing can be calculated. Options are:  
        "KU" (default) -- Spin squeezing defined by Kitagawa and Ueda.  
        "WBIMH" -- Spin squeezing defined by Wineland et al.

    Returns
    ----------
    **$\xi$:** `float`
        -- spin squeezing parameter
    """

    N = len(rho) - 1

    coef = 4.0 / float(N)
    j = N / 2
    if basis == "Pauli":
        sp = np.array([[0.0,1.0],[0.0, 0.0]])
        jp = []
        for i in range(0, N): 
            if i==0 :
                jp_tp = np.kron(sp, np.identity(2**(N-1)))
            elif i==N-1 :
                jp_tp = np.kron(np.identity(2**(N-1)), sp)
            else:
                jp_tp = np.kron(np.identity(2**i), np.kron(sp, np.identity(2**(N-1-i))))
            jp.append(jp_tp)
        Jp = sum(jp)
    else:
        offdiag = [
            np.sqrt(float(j * (j + 1) - m * (m + 1))) for m in np.arange(j, -j - 1, -1)
        ][1:]

        Jp = np.matrix(np.diag(offdiag, 1))
    Jx = 0.5 * (Jp + Jp.H)
    Jy = -0.5 * 1j * (Jp - Jp.H)
    Jz = np.diag(np.arange(j, -j - 1, -1))
    Jx_mean = np.trace(rho * Jx)
    Jy_mean = np.trace(rho * Jy)
    Jz_mean = np.trace(rho * Jz)

    costheta = Jz_mean / np.sqrt(Jx_mean**2 + Jy_mean**2 + Jz_mean**2)
    sintheta = np.sin(np.arccos(costheta))
    cosphi = Jx_mean / np.sqrt(Jx_mean**2 + Jy_mean**2)
    if Jy_mean > 0:
        sinphi = np.sin(np.arccos(cosphi))
    else:
        sinphi = np.sin(2 * np.pi - np.arccos(cosphi))
    Jn1 = -Jx * sinphi + Jy * cosphi
    Jn2 = -Jx * costheta * cosphi - Jy * costheta * sinphi + Jz * sintheta
    A = np.trace(rho * (Jn1 * Jn1 - Jn2 * Jn2))
    B = np.trace(rho * (Jn1 * Jn2 + Jn2 * Jn1))
    C = np.trace(rho * (Jn1 * Jn1 + Jn2 * Jn2))

    V_minus = 0.5 * (C - np.sqrt(A**2 + B**2))
    V_minus = np.real(V_minus)
    Xi = coef * V_minus
    if Xi > 1.0:
        Xi = 1.0

    if output == "KU":
        Xi = Xi
    elif output == "WBIMH":
        Xi = (N / 2) ** 2 * Xi / (Jx_mean**2 + Jy_mean**2 + Jz_mean**2)
    else:
        raise NameError("NameError: output should be choosen in {KU, WBIMH}")

    return Xi

Calculation of the time to reach a given precision limit.

Parameters

f: float -- The given value of the objective function.

tspan: array -- Time length for the evolution.

func: array -- The function for calculating the objective function.

*args: string -- The corresponding input parameter.

**kwargs: string -- Keyword arguments in func.

Returns

time: float -- Time to reach the given target.

Source code in quanestimation/Resource/Resource.py

def TargetTime(f, tspan, func, *args, **kwargs):
    """
    Calculation of the time to reach a given precision limit.

    Parameters
    ----------
    > **f:** `float`
        -- The given value of the objective function.

    > **tspan:** `array`
        -- Time length for the evolution.

    > **func:** `array`
        -- The function for calculating the objective function.

    > ***args:** `string`
        -- The corresponding input parameter.

    > ****kwargs:** `string`
        -- Keyword arguments in `func`.

    Returns
    ----------
    **time:** `float`
        -- Time to reach the given target.
    """

    args = list(zip_broadcast(*args))

    f_last = func(*(args[0]), **kwargs)
    idx = 1
    f_now = func(*(args[1]), **kwargs)

    while (f_now - f) * (f_last - f) > 0 and idx < (len(tspan) - 1):
        f_last = f_now
        idx += 1
        f_now = func(*(args[idx]), **kwargs)

    return tspan[idx]

---

Quantum Cramér-Rao bounds

Calculation of the symmetric logarithmic derivative (SLD) for a density matrix. The SLD operator \(L_a\) is determined by \begin{align} \partial_{a}\rho=\frac{1}{2}(\rho L_{a}+L_{a}\rho) \end{align}

with \(\rho\) the parameterized density matrix. The entries of SLD can be calculated as \begin{align} \langle\lambda_i|L_{a}|\lambda_j\rangle=\frac{2\langle\lambda_i| \partial_{a}\rho |\lambda_j\rangle}{\lambda_i+\lambda_j} \end{align}

for \(\lambda_i~(\lambda_j) \neq 0\). If \(\lambda_i=\lambda_j=0\), the entry of SLD is set to be zero.

Parameters

rho: matrix -- Density matrix.

drho: list -- Derivatives of the density matrix on the unknown parameters to be estimated. For example, drho[0] is the derivative vector on the first parameter.

rep: string -- The basis for the SLDs. Options are:
"original" (default) -- it means the basis is the same with the input density matrix (rho).
"eigen" -- it means the basis is the same with theeigenspace of the density matrix (rho).

eps: float -- Machine epsilon.

Returns

SLD(s): matrix or list --For single parameter estimation (the length of drho is equal to one), the output is a matrix and for multiparameter estimation (the length of drho is more than one), it returns a list.

Source code in quanestimation/AsymptoticBound/CramerRao.py

def SLD(rho, drho, rep="original", eps=1e-8):
    r"""
    Calculation of the symmetric logarithmic derivative (SLD) for a density matrix.
    The SLD operator $L_a$ is determined by
    \begin{align}
    \partial_{a}\rho=\frac{1}{2}(\rho L_{a}+L_{a}\rho)
    \end{align}

    with $\rho$ the parameterized density matrix. The entries of SLD can be calculated
    as 
    \begin{align}
    \langle\lambda_i|L_{a}|\lambda_j\rangle=\frac{2\langle\lambda_i| \partial_{a}\rho |\lambda_j\rangle}{\lambda_i+\lambda_j}
    \end{align}

    for $\lambda_i~(\lambda_j) \neq 0$. If $\lambda_i=\lambda_j=0$, the entry of SLD is set to be zero.

    Parameters
    ----------
    > **rho:** `matrix`
        -- Density matrix.

    > **drho:** `list`
        -- Derivatives of the density matrix on the unknown parameters to be 
        estimated. For example, drho[0] is the derivative vector on the first 
        parameter.

    > **rep:** `string`
        -- The basis for the SLDs. Options are:  
        "original" (default) -- it means the basis is the same with the input density 
        matrix (rho).  
        "eigen" -- it means the basis is the same with theeigenspace of the density
        matrix (rho).

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

    Returns
    ----------
    **SLD(s):** `matrix or list`
        --For single parameter estimation (the length of drho is equal to one), the
        output is a matrix and for multiparameter estimation (the length of drho 
        is more than one), it returns a list.
    """

    if type(drho) != list:
        raise TypeError("Please make sure drho is a list!")

    para_num = len(drho)
    dim = len(rho)
    SLD = [[] for i in range(0, para_num)]

    purity = np.trace(np.dot(rho, rho))

    if np.abs(1 - purity) < eps:
        SLD_org = [[] for i in range(0, para_num)]
        for para_i in range(0, para_num):
            SLD_org[para_i] = 2 * drho[para_i]

            if rep == "original":
                SLD[para_i] = SLD_org[para_i]
            elif rep == "eigen":
                val, vec = np.linalg.eig(rho)
                val = np.real(val)
                SLD[para_i] = np.dot(
                    vec.conj().transpose(), np.dot(SLD_org[para_i], vec)
                )
            else:
                raise ValueError("{!r} is not a valid value for rep, supported values are 'original' and 'eigen'.".format(rep))
        if para_num == 1:
            return SLD[0]
        else:
            return SLD
    else:
        val, vec = np.linalg.eig(rho)
        val = np.real(val)
        for para_i in range(0, para_num):
            SLD_eig = np.array(
                [[0.0 + 0.0 * 1.0j for i in range(0, dim)] for i in range(0, dim)]
            )
            for fi in range(0, dim):
                for fj in range(0, dim):
                    if val[fi] + val[fj] > eps:
                        SLD_eig[fi][fj] = (
                            2
                            * np.dot(
                                vec[:, fi].conj().transpose(),
                                np.dot(drho[para_i], vec[:, fj]),
                            )
                            / (val[fi] + val[fj])
                        )
            SLD_eig[SLD_eig == np.inf] = 0.0

            if rep == "original":
                SLD[para_i] = np.dot(vec, np.dot(SLD_eig, vec.conj().transpose()))
            elif rep == "eigen":
                SLD[para_i] = SLD_eig
            else:
                raise ValueError("{!r} is not a valid value for rep, supported values are 'original' and 'eigen'.".format(rep))

        if para_num == 1:
            return SLD[0]
        else:
            return SLD

Calculation of the right logarithmic derivative (RLD) for a density matrix. The RLD operator defined by \(\partial_{a}\rho=\rho \mathcal{R}_a\) with \(\rho\) the parameterized density matrix. \begin{align} \langle\lambda_i| \mathcal{R}_{a} |\lambda_j\rangle=\frac{1}{\lambda_i}\langle\lambda_i| \partial_a\rho |\lambda_j\rangle \end{align}

for \(\lambda_i\neq 0\) is the \(ij\)th entry of RLD.

Parameters

rho: matrix -- Density matrix.

drho: list -- Derivatives of the density matrix on the unknown parameters to be estimated. For example, drho[0] is the derivative vector on the first parameter.

rep: string -- The basis for the RLD(s). Options are:
"original" (default) -- it means the basis is the same with the input density matrix (rho).
"eigen" -- it means the basis is the same with the eigenspace of the density matrix (rho).

eps: float -- Machine epsilon.

Returns

RLD(s): matrix or list -- For single parameter estimation (the length of drho is equal to one), the output is a matrix and for multiparameter estimation (the length of drho is more than one), it returns a list.

Source code in quanestimation/AsymptoticBound/CramerRao.py

def RLD(rho, drho, rep="original", eps=1e-8):
    r"""
    Calculation of the right logarithmic derivative (RLD) for a density matrix.
    The RLD operator defined by $\partial_{a}\rho=\rho \mathcal{R}_a$
    with $\rho$ the parameterized density matrix. 
    \begin{align}
    \langle\lambda_i| \mathcal{R}_{a} |\lambda_j\rangle=\frac{1}{\lambda_i}\langle\lambda_i| 
    \partial_a\rho |\lambda_j\rangle 
    \end{align}

    for $\lambda_i\neq 0$ is the $ij$th entry of RLD.

    Parameters
    ----------
    > **rho:** `matrix`
        -- Density matrix.

    > **drho:** `list`
        -- Derivatives of the density matrix on the unknown parameters to be 
        estimated. For example, drho[0] is the derivative vector on the first 
        parameter.

    > **rep:** `string`
        -- The basis for the RLD(s). Options are:  
        "original" (default) -- it means the basis is the same with the input density 
        matrix (rho).  
        "eigen" -- it means the basis is the same with the eigenspace of the density 
        matrix (rho).

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

    Returns
    ----------
    **RLD(s):** `matrix or list`
        -- For single parameter estimation (the length of drho is equal to one), the output 
        is a matrix and for multiparameter estimation (the length of drho is more than one), 
        it returns a list.
    """

    if type(drho) != list:
        raise TypeError("Please make sure drho is a list!")

    para_num = len(drho)
    dim = len(rho)
    RLD = [[] for i in range(0, para_num)]

    val, vec = np.linalg.eig(rho)
    val = np.real(val)
    for para_i in range(0, para_num):
        RLD_eig = np.array(
            [[0.0 + 0.0 * 1.0j for i in range(0, dim)] for i in range(0, dim)]
        )
        for fi in range(0, dim):
            for fj in range(0, dim):
                term_tp = np.dot(vec[:, fi].conj().transpose(),np.dot(drho[para_i], vec[:, fj]))
                if np.abs(val[fi]) > eps:
                    RLD_eig[fi][fj] = (term_tp/val[fi])
                else:
                    if np.abs(term_tp) < eps:
                        raise ValueError("The RLD does not exist. It only exist when the support of drho is contained in the support of rho.",
            )
        RLD_eig[RLD_eig == np.inf] = 0.0

        if rep == "original":
            RLD[para_i] = np.dot(vec, np.dot(RLD_eig, vec.conj().transpose()))
        elif rep == "eigen":
            RLD[para_i] = RLD_eig
        else:
            raise ValueError("{!r} is not a valid value for rep, supported values are 'original' and 'eigen'.".format(rep))
    if para_num == 1:
        return RLD[0]
    else:
        return RLD

Calculation of the left logarithmic derivative (LLD) for a density matrix \(\rho\). The LLD operator is defined by \(\partial_{a}\rho=\mathcal{R}_a^{\dagger}\rho\). The entries of LLD can be calculated as \begin{align} \langle\lambda_i| \mathcal{R}_{a}^{\dagger} |\lambda_j\rangle=\frac{1}{\lambda_j}\langle\lambda_i| \partial_a\rho |\lambda_j\rangle \end{align}

for \(\lambda_j\neq 0\).

Parameters

rho: matrix -- Density matrix.

drho: list -- Derivatives of the density matrix on the unknown parameters to be estimated. For example, drho[0] is the derivative vector on the first parameter.

rep: string -- The basis for the LLD(s). Options are:
"original" (default) -- it means the basis is the same with the input density matrix (rho).
"eigen" -- it means the basis is the same with the eigenspace of the density matrix (rho).

eps: float -- Machine epsilon.

Returns

LLD(s): matrix or list -- For single parameter estimation (the length of drho is equal to one), the output is a matrix and for multiparameter estimation (the length of drho is more than one), it returns a list.

Source code in quanestimation/AsymptoticBound/CramerRao.py

def LLD(rho, drho, rep="original", eps=1e-8):
    r"""
    Calculation of the left logarithmic derivative (LLD) for a density matrix $\rho$.
    The LLD operator is defined by $\partial_{a}\rho=\mathcal{R}_a^{\dagger}\rho$. 
    The entries of LLD can be calculated as 
    \begin{align}
    \langle\lambda_i| \mathcal{R}_{a}^{\dagger} |\lambda_j\rangle=\frac{1}{\lambda_j}\langle\lambda_i| 
    \partial_a\rho |\lambda_j\rangle 
    \end{align}

    for $\lambda_j\neq 0$.

    Parameters
    ----------
    > **rho:** `matrix`
        -- Density matrix.

    > **drho:** `list`
        -- Derivatives of the density matrix on the unknown parameters to be 
        estimated. For example, drho[0] is the derivative vector on the first 
        parameter.

    > **rep:** `string`
        -- The basis for the LLD(s). Options are:  
        "original" (default) -- it means the basis is the same with the input density 
        matrix (rho).  
        "eigen" -- it means the basis is the same with the eigenspace of the density 
        matrix (rho).

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

    Returns
    ----------
    **LLD(s):** `matrix or list`
        -- For single parameter estimation (the length of drho is equal to one), the output 
        is a matrix and for multiparameter estimation (the length of drho is more than one), 
        it returns a list.
    """

    if type(drho) != list:
        raise TypeError("Please make sure drho is a list!")

    para_num = len(drho)
    dim = len(rho)
    LLD = [[] for i in range(0, para_num)]

    val, vec = np.linalg.eig(rho)
    val = np.real(val)
    for para_i in range(0, para_num):
        LLD_eig = np.array(
            [[0.0 + 0.0 * 1.0j for i in range(0, dim)] for i in range(0, dim)]
        )
        for fi in range(0, dim):
            for fj in range(0, dim):
                term_tp = np.dot(vec[:, fi].conj().transpose(), np.dot(drho[para_i], vec[:, fj]),)
                if np.abs(val[fj]) > eps:
                    LLD_eig_tp = (term_tp/val[fj])
                    LLD_eig[fj][fi] = LLD_eig_tp.conj()
                else: 
                    if np.abs(term_tp) < eps:
                        raise ValueError("The LLD does not exist. It only exist when the support of drho is contained in the support of rho.",
            )
        LLD_eig[LLD_eig == np.inf] = 0.0

        if rep == "original":
            LLD[para_i] = np.dot(vec, np.dot(LLD_eig, vec.conj().transpose()))
        elif rep == "eigen":
            LLD[para_i] = LLD_eig
        else:
            raise ValueError("{!r} is not a valid value for rep, supported values are 'original' and 'eigen'.".format(rep))

    if para_num == 1:
        return LLD[0]
    else:
        return LLD

Calculation of the quantum Fisher information (QFI) and quantum Fisher information matrix (QFIM) for all types. The entry of QFIM \(\mathcal{F}\) is defined as \begin{align} \mathcal{F}_{ab}=\frac{1}{2}\mathrm{Tr}(\rho{L_a, L_b}) \end{align}

with \(L_a, L_b\) are SLD operators and

and \begin{align} \mathcal{F}_{ab}=\mathrm{Tr}(\rho \mathcal{R}_a \mathcal{R}^{\dagger}_b) \end{align}

with \(\mathcal{R}_a\) the RLD or LLD operator.

Parameters

rho: matrix -- Density matrix.

drho: list Derivatives of the density matrix on the unknown parameters to be estimated. For example, drho[0] is the derivative vector on the first parameter.

LDtype: string -- Types of QFI (QFIM) can be set as the objective function. Options are:
"SLD" (default) -- QFI (QFIM) based on symmetric logarithmic derivative (SLD).
"RLD" -- QFI (QFIM) based on right logarithmic derivative (RLD).
"LLD" -- QFI (QFIM) based on left logarithmic derivative (LLD).

exportLD: bool -- Whether or not to export the values of logarithmic derivatives. If set True then the the values of logarithmic derivatives will be exported.

eps: float -- Machine epsilon.

Returns

QFI or QFIM: float or matrix -- For single parameter estimation (the length of drho is equal to one), the output is QFI and for multiparameter estimation (the length of drho is more than one), it returns QFIM.

Source code in quanestimation/AsymptoticBound/CramerRao.py

def QFIM(rho, drho, LDtype="SLD", exportLD=False, eps=1e-8):
    r"""
    Calculation of the quantum Fisher information (QFI) and quantum Fisher 
    information matrix (QFIM) for all types. The entry of QFIM $\mathcal{F}$
    is defined as
    \begin{align}
    \mathcal{F}_{ab}=\frac{1}{2}\mathrm{Tr}(\rho\{L_a, L_b\})
    \end{align}

    with $L_a, L_b$ are SLD operators and 

    and 
    \begin{align}
    \mathcal{F}_{ab}=\mathrm{Tr}(\rho \mathcal{R}_a \mathcal{R}^{\dagger}_b)
    \end{align}

    with $\mathcal{R}_a$ the RLD or LLD operator.

    Parameters
    ----------
    > **rho:** `matrix`
        -- Density matrix.

    > **drho:** `list`
        Derivatives of the density matrix on the unknown parameters to be 
        estimated. For example, drho[0] is the derivative vector on the first 
        parameter.

    > **LDtype:** `string`
        -- Types of QFI (QFIM) can be set as the objective function. Options are:  
        "SLD" (default) -- QFI (QFIM) based on symmetric logarithmic derivative (SLD).  
        "RLD" -- QFI (QFIM) based on right logarithmic derivative (RLD).  
        "LLD" -- QFI (QFIM) based on left logarithmic derivative (LLD).

    > **exportLD:** `bool`
        -- Whether or not to export the values of logarithmic derivatives. If set True
        then the the values of logarithmic derivatives will be exported.

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

    Returns
    ----------
    **QFI or QFIM:** `float or matrix` 
        -- For single parameter estimation (the length of drho is equal to one), 
        the output is QFI and for multiparameter estimation (the length of drho 
        is more than one), it returns QFIM.
    """

    if type(drho) != list:
        raise TypeError("Please make sure drho is a list")

    para_num = len(drho)

    # single parameter estimation
    if para_num == 1:
        if LDtype == "SLD":
            LD_tp = SLD(rho, drho, eps=eps)
            SLD_ac = np.dot(LD_tp, LD_tp) + np.dot(LD_tp, LD_tp)
            QFIM_res = np.real(0.5 * np.trace(np.dot(rho, SLD_ac)))
        elif LDtype == "RLD":
            LD_tp = RLD(rho, drho, eps=eps)
            QFIM_res = np.real(
                np.trace(np.dot(rho, np.dot(LD_tp, LD_tp.conj().transpose())))
            )
        elif LDtype == "LLD":
            LD_tp = LLD(rho, drho, eps=eps)
            QFIM_res = np.real(
                np.trace(np.dot(rho, np.dot(LD_tp, LD_tp.conj().transpose())))
            )
        else:
            raise ValueError("{!r} is not a valid value for LDtype, supported values are 'SLD', 'RLD' and 'LLD'.".format(LDtype))

    # multiparameter estimation
    else:
        if LDtype == "SLD":
            QFIM_res = np.zeros([para_num, para_num])
            LD_tp = SLD(rho, drho, eps=eps)
            for para_i in range(0, para_num):
                for para_j in range(para_i, para_num):
                    SLD_ac = np.dot(LD_tp[para_i], LD_tp[para_j]) + np.dot(
                        LD_tp[para_j], LD_tp[para_i]
                    )
                    QFIM_res[para_i][para_j] = np.real(
                        0.5 * np.trace(np.dot(rho, SLD_ac))
                    )
                    QFIM_res[para_j][para_i] = QFIM_res[para_i][para_j]
        elif LDtype == "RLD":
            QFIM_res = np.zeros((para_num, para_num), dtype=np.complex128)
            LD_tp = RLD(rho, drho, eps=eps)
            for para_i in range(0, para_num):
                for para_j in range(para_i, para_num):
                    QFIM_res[para_i][para_j] = np.trace(
                            np.dot(
                                rho,
                                np.dot(LD_tp[para_i], LD_tp[para_j].conj().transpose()),
                            )
                        )
                    QFIM_res[para_j][para_i] = QFIM_res[para_i][para_j].conj()
        elif LDtype == "LLD":
            QFIM_res = np.zeros((para_num, para_num), dtype=np.complex128)
            LD_tp = LLD(rho, drho, eps=eps)
            for para_i in range(0, para_num):
                for para_j in range(para_i, para_num):
                    QFIM_res[para_i][para_j] = np.trace(
                            np.dot(
                                rho,
                                np.dot(LD_tp[para_i], LD_tp[para_j].conj().transpose()),
                            )
                        )
                    QFIM_res[para_j][para_i] = QFIM_res[para_i][para_j].conj()
        else:
            raise ValueError("{!r} is not a valid value for LDtype, supported values are 'SLD', 'RLD' and 'LLD'.".format(LDtype))

    if exportLD == False:
        return QFIM_res
    else:
        return QFIM_res, LD_tp

Calculation of the quantum Fisher information (QFI) and quantum Fisher information matrix (QFIM) with Kraus operator(s) for all types.

Parameters

rho0: matrix -- Initial state (density matrix).

K: list -- Kraus operator(s).

dK: list -- Derivatives of the Kraus operator(s) on the unknown parameters to be estimated.

LDtype: string -- Types of QFI (QFIM) can be set as the objective function. Options are:
"SLD" (default) -- QFI (QFIM) based on symmetric logarithmic derivative (SLD).
"RLD" -- QFI (QFIM) based on right logarithmic derivative (RLD).
"LLD" -- QFI (QFIM) based on left logarithmic derivative (LLD).

exportLD: bool -- Whether or not to export the values of logarithmic derivatives. If set True then the the values of logarithmic derivatives will be exported.

eps: float -- Machine epsilon.

Returns

QFI or QFIM: float or matrix -- For single parameter estimation (the length of drho is equal to one), the output is QFI and for multiparameter estimation (the length of drho is more than one), it returns QFIM.

Source code in quanestimation/AsymptoticBound/CramerRao.py

def QFIM_Kraus(rho0, K, dK, LDtype="SLD", exportLD=False, eps=1e-8):
    """
    Calculation of the quantum Fisher information (QFI) and quantum Fisher 
    information matrix (QFIM) with Kraus operator(s) for all types.

    Parameters
    ----------
    > **rho0:** `matrix`
        -- Initial state (density matrix).

    > **K:** `list`
        -- Kraus operator(s).

    > **dK:** `list` 
        -- Derivatives of the Kraus operator(s) on the unknown parameters to be 
        estimated.

    > **LDtype:** `string`
        -- Types of QFI (QFIM) can be set as the objective function. Options are:  
        "SLD" (default) -- QFI (QFIM) based on symmetric logarithmic derivative (SLD).  
        "RLD" -- QFI (QFIM) based on right logarithmic derivative (RLD).  
        "LLD" -- QFI (QFIM) based on left logarithmic derivative (LLD).

    > **exportLD:** `bool`
        -- Whether or not to export the values of logarithmic derivatives. If set True
        then the the values of logarithmic derivatives will be exported.

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

    Returns
    ----------
    **QFI or QFIM:** `float or matrix`
        -- For single parameter estimation (the length of drho is equal to one), 
        the output is QFI and for multiparameter estimation (the length of drho 
        is more than one), it returns QFIM.
    """

    dK = [[dK[i][j] for i in range(len(K))] for j in range(len(dK[0]))]
    rho = sum([np.dot(Ki, np.dot(rho0, Ki.conj().T)) for Ki in K])
    drho = [
                sum(
                    [
                        (
                            np.dot(dKi, np.dot(rho0, Ki.conj().T))
                            + np.dot(Ki, np.dot(rho0, dKi.conj().T))
                        )
                        for (Ki, dKi) in zip(K, dKj)
                    ]
                )
                for dKj in dK
            ]
    return QFIM(rho, drho, LDtype=LDtype, exportLD=exportLD, eps=eps)

Calculation of the classical Fisher information (CFI) and classical Fisher information matrix (CFIM) for a density matrix. The entry of CFIM \(\mathcal{I}\) is defined as \begin{align} \mathcal{I}_{ab}=\sum_y\frac{1}{p(y|\textbf{x})}[\partial_a p(y|\textbf{x})][\partial_b p(y|\textbf{x})], \end{align}

where \(p(y|\textbf{x})=\mathrm{Tr}(\rho\Pi_y)\) with \(\rho\) the parameterized density matrix.

Parameters

rho: matrix -- Density matrix.

drho: list -- Derivatives of the density matrix on the unknown parameters to be estimated. For example, drho[0] is the derivative vector on the first parameter.

M: list of matrices -- A set of positive operator-valued measure (POVM). The default measurement is a set of rank-one symmetric informationally complete POVM (SIC-POVM).

eps: float -- Machine epsilon.

Returns

CFI (CFIM): float or matrix -- For single parameter estimation (the length of drho is equal to one), the output is CFI and for multiparameter estimation (the length of drho is more than one), it returns CFIM.

Note: SIC-POVM is calculated by the Weyl-Heisenberg covariant SIC-POVM fiducial state which can be downloaded from here.

Source code in quanestimation/AsymptoticBound/CramerRao.py

def CFIM(rho, drho, M=[], eps=1e-8):
    r"""
    Calculation of the classical Fisher information (CFI) and classical Fisher 
    information matrix (CFIM) for a density matrix. The entry of CFIM $\mathcal{I}$
    is defined as
    \begin{align}
    \mathcal{I}_{ab}=\sum_y\frac{1}{p(y|\textbf{x})}[\partial_a p(y|\textbf{x})][\partial_b p(y|\textbf{x})],
    \end{align}

    where $p(y|\textbf{x})=\mathrm{Tr}(\rho\Pi_y)$ with $\rho$ the parameterized 
    density matrix.

    Parameters
    ----------
    > **rho:** `matrix`
        -- Density matrix.

    > **drho:** `list`
        -- Derivatives of the density matrix on the unknown parameters to be 
        estimated. For example, drho[0] is the derivative vector on the first 
        parameter.

    > **M:** `list of matrices`
        -- A set of positive operator-valued measure (POVM). The default measurement 
        is a set of rank-one symmetric informationally complete POVM (SIC-POVM).

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

    Returns
    ----------
    **CFI (CFIM):** `float or matrix` 
        -- For single parameter estimation (the length of drho is equal to one), 
        the output is CFI and for multiparameter estimation (the length of drho 
        is more than one), it returns CFIM.

    **Note:** 
        SIC-POVM is calculated by the Weyl-Heisenberg covariant SIC-POVM fiducial state 
        which can be downloaded from [here](http://www.physics.umb.edu/Research/QBism/
        solutions.html).
    """

    if type(drho) != list:
        raise TypeError("Please make sure drho is a list!")

    if M == []:
        M = SIC(len(rho[0]))
    else:
        if type(M) != list:
            raise TypeError("Please make sure M is a list!")

    m_num = len(M)
    para_num = len(drho)
    CFIM_res = np.zeros([para_num, para_num])
    for pi in range(0, m_num):
        Mp = M[pi]
        p = np.real(np.trace(np.dot(rho, Mp)))
        Cadd = np.zeros([para_num, para_num])
        if p > eps:
            for para_i in range(0, para_num):
                drho_i = drho[para_i]
                dp_i = np.real(np.trace(np.dot(drho_i, Mp)))
                for para_j in range(para_i, para_num):
                    drho_j = drho[para_j]
                    dp_j = np.real(np.trace(np.dot(drho_j, Mp)))
                    Cadd[para_i][para_j] = np.real(dp_i * dp_j / p)
                    Cadd[para_j][para_i] = np.real(dp_i * dp_j / p)
        CFIM_res += Cadd

    if para_num == 1:
        return CFIM_res[0][0]
    else:
        return CFIM_res

Calculation of the classical Fisher information (CFI) and classical Fisher information matrix (CFIM) for classical scenarios. The entry of FIM \(I\) is defined as \begin{align} I_{ab}=\sum_{y}\frac{1}{p_y}[\partial_a p_y][\partial_b p_y], \end{align}

where \(\{p_y\}\) is a set of the discrete probability distribution.

Parameters

p: array -- The probability distribution.

dp: list -- Derivatives of the probability distribution on the unknown parameters to be estimated. For example, dp[0] is the derivative vector on the first parameter.

eps: float -- Machine epsilon.

Returns

CFI (CFIM): float or matrix -- For single parameter estimation (the length of drho is equal to one), the output is CFI and for multiparameter estimation (the length of drho is more than one), it returns CFIM.

Source code in quanestimation/AsymptoticBound/CramerRao.py

def FIM(p, dp, eps=1e-8):
    r"""
    Calculation of the classical Fisher information (CFI) and classical Fisher 
    information matrix (CFIM) for classical scenarios. The entry of FIM $I$
    is defined as
    \begin{align}
    I_{ab}=\sum_{y}\frac{1}{p_y}[\partial_a p_y][\partial_b p_y],
    \end{align}

    where $\{p_y\}$ is a set of the discrete probability distribution.

    Parameters
    ----------
    > **p:** `array` 
        -- The probability distribution.

    > **dp:** `list`
        -- Derivatives of the probability distribution on the unknown parameters to 
        be estimated. For example, dp[0] is the derivative vector on the first 
        parameter.

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

    Returns
    ----------
    **CFI (CFIM):** `float or matrix` 
        -- For single parameter estimation (the length of drho is equal to one), 
        the output is CFI and for multiparameter estimation (the length of drho 
        is more than one), it returns CFIM.
    """

    para_num = len(dp[0])
    m_num = len(p)
    FIM_res = np.zeros([para_num, para_num])
    for pi in range(0, m_num):
        p_tp = p[pi]
        Cadd = np.zeros([para_num, para_num])
        if p_tp > eps:
            for para_i in range(0, para_num):
                dp_i = dp[pi][para_i]
                for para_j in range(para_i, para_num):
                    dp_j = dp[pi][para_j]
                    Cadd[para_i][para_j] = np.real(dp_i * dp_j / p_tp)
                    Cadd[para_j][para_i] = np.real(dp_i * dp_j / p_tp)
        FIM_res += Cadd

    if para_num == 1:
        return FIM_res[0][0]
    else:
        return FIM_res

Calculation of the classical Fisher information (CFI) based on the experiment data.

Parameters

y1: array -- Experimental data obtained at the truth value (x).

y2: list -- Experimental data obtained at x+dx.

dx: float -- A known small drift of the parameter.

ftype: string -- The distribution the data follows. Options are:
"norm" (default) -- normal distribution.
"gamma" -- gamma distribution. "rayleigh" -- rayleigh distribution. "poisson" -- poisson distribution.

Returns

CFI: float or matrix

Source code in quanestimation/AsymptoticBound/CramerRao.py

def FI_Expt(y1, y2, dx, ftype="norm"):
    r"""
    Calculation of the classical Fisher information (CFI) based on the experiment data. 

    Parameters
    ----------
    > **y1:** `array` 
        -- Experimental data obtained at the truth value (x).

    > **y2:** `list`
        -- Experimental data obtained at x+dx.

    > **dx:** `float`
        -- A known small drift of the parameter.

    > **ftype:** `string`
        -- The distribution the data follows. Options are:  
        "norm" (default) -- normal distribution.  
        "gamma" -- gamma distribution.
        "rayleigh" -- rayleigh distribution.
        "poisson" -- poisson distribution.

    Returns
    ----------
    **CFI:** `float or matrix` 
    """
    fidelity = 0.0
    if ftype == "norm":
        mu1, std1 = norm.fit(y1)
        mu2, std2 = norm.fit(y2)
        f_func = lambda x: np.sqrt(norm.pdf(x, mu1, std1)*norm.pdf(x, mu2, std2))
        fidelity, err = quad(f_func, -np.inf, np.inf)
    elif ftype == "gamma":
        a1, alpha1, beta1 = gamma.fit(y1)
        a2, alpha2, beta2 = gamma.fit(y2)
        f_func = lambda x: np.sqrt(gamma.pdf(x, a1, alpha1, beta1)*gamma.pdf(x, a2, alpha2, beta2))
        fidelity, err = quad(f_func, 0., np.inf)
    elif ftype == "rayleigh":
        mean1, var1 = rayleigh.fit(y1)
        mean2, var2 = rayleigh.fit(y2)
        f_func = lambda x: np.sqrt(rayleigh.pdf(x, mean1, var1)*rayleigh.pdf(x, mean2, var2))
        fidelity, err = quad(f_func, -np.inf, np.inf)
    elif ftype == "poisson":
        k1 = np.arange(max(y1)+1)
        k2 = np.arange(max(y2)+1)
        p1_pois = poisson.pmf(k1, np.mean(y1))
        p2_pois = poisson.pmf(k2, np.mean(y2))
        p1_pois, p2_pois = p1_pois/sum(p1_pois), p2_pois/sum(p2_pois)
        fidelity = sum([np.sqrt(p1_pois[i]*p2_pois[i]) for i in range(len(p1_pois))])
    else:
        raise ValueError("{!r} is not a valid value for ftype, supported values are 'norm', 'poisson', 'gamma' and 'rayleigh'.".format(ftype))
    Fc = 8*(1-fidelity)/dx**2
    return Fc

Calculation of the SLD based quantum Fisher information (QFI) and quantum
Fisher information matrix (QFIM) in Bloch representation.

Parameters

r: list -- Parameterized Bloch vector.

dr: list -- Derivatives of the Bloch vector on the unknown parameters to be estimated. For example, dr[0] is the derivative vector on the first parameter.

eps: float -- Machine epsilon.

Returns

QFI or QFIM in Bloch representation: float or matrix -- For single parameter estimation (the length of drho is equal to one), the output is QFI and for multiparameter estimation (the length of drho is more than one), it returns QFIM.

Source code in quanestimation/AsymptoticBound/CramerRao.py

def QFIM_Bloch(r, dr, eps=1e-8):
    """
    Calculation of the SLD based quantum Fisher information (QFI) and quantum  
    Fisher information matrix (QFIM) in Bloch representation.

    Parameters
    ----------
    > **r:** `list`
        -- Parameterized Bloch vector.

    > **dr:** `list `
        -- Derivatives of the Bloch vector on the unknown parameters to be 
        estimated. For example, dr[0] is the derivative vector on the first 
        parameter.

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

    Returns
    ----------
    **QFI or QFIM in Bloch representation:** `float or matrix`
        -- For single parameter estimation (the length of drho is equal to one), 
        the output is QFI and for multiparameter estimation (the length of drho 
        is more than one), it returns QFIM.
    """

    if type(dr) != list:
        raise TypeError("Please make sure dr is a list")

    para_num = len(dr)
    QFIM_res = np.zeros([para_num, para_num])

    dim = int(np.sqrt(len(r) + 1))
    Lambda = suN_generator(dim)

    if dim == 2:
        #### single-qubit system ####
        r_norm = np.linalg.norm(r) ** 2
        if np.abs(r_norm - 1.0) < eps:
            for para_i in range(0, para_num):
                for para_j in range(para_i, para_num):
                    QFIM_res[para_i][para_j] = np.real(np.inner(dr[para_i], dr[para_j]))
                    QFIM_res[para_j][para_i] = QFIM_res[para_i][para_j]
        else:
            for para_i in range(0, para_num):
                for para_j in range(para_i, para_num):
                    QFIM_res[para_i][para_j] = np.real(
                        np.inner(dr[para_i], dr[para_j])
                        + np.inner(r, dr[para_i])
                        * np.inner(r, dr[para_j])
                        / (1 - r_norm)
                    )
                    QFIM_res[para_j][para_i] = QFIM_res[para_i][para_j]
    else:
        rho = np.identity(dim, dtype=np.complex128) / dim
        for di in range(dim**2 - 1):
            rho += np.sqrt(dim * (dim - 1) / 2) * r[di] * Lambda[di] / dim

        G = np.zeros((dim**2 - 1, dim**2 - 1), dtype=np.complex128)
        for row_i in range(dim**2 - 1):
            for col_j in range(row_i, dim**2 - 1):
                anti_commu = np.dot(Lambda[row_i], Lambda[col_j]) + np.dot(
                    Lambda[col_j], Lambda[row_i]
                )
                G[row_i][col_j] = 0.5 * np.trace(np.dot(rho, anti_commu))
                G[col_j][row_i] = G[row_i][col_j]

        mat_tp = G * dim / (2 * (dim - 1)) - np.dot(
            np.array(r).reshape(len(r), 1), np.array(r).reshape(1, len(r))
        )
        mat_inv = inv(mat_tp)

        for para_m in range(0, para_num):
            for para_n in range(para_m, para_num):
                QFIM_res[para_m][para_n] = np.real(
                    np.dot(
                        np.array(dr[para_n]).reshape(1, len(r)),
                        np.dot(mat_inv, np.array(dr[para_m]).reshape(len(r), 1)),
                    )[0][0]
                )
                QFIM_res[para_n][para_m] = QFIM_res[para_m][para_n]

    if para_num == 1:
        return QFIM_res[0][0]
    else:
        return QFIM_res

Calculation of the SLD based quantum Fisher information (QFI) and quantum Fisher information matrix (QFIM) with gaussian states.

Parameters

R: array -- First-order moment.

dR: list -- Derivatives of the first-order moment on the unknown parameters to be estimated. For example, dR[0] is the derivative vector on the first parameter.

D: matrix -- Second-order moment.

dD: list -- Derivatives of the second-order moment on the unknown parameters to be estimated. For example, dD[0] is the derivative vector on the first parameter.

eps: float -- Machine epsilon.

Returns

QFI or QFIM with gaussian states: float or matrix -- For single parameter estimation (the length of drho is equal to one), the output is QFI and for multiparameter estimation (the length of drho is more than one), it returns QFIM.

Source code in quanestimation/AsymptoticBound/CramerRao.py

def QFIM_Gauss(R, dR, D, dD):
    """
    Calculation of the SLD based quantum Fisher information (QFI) and quantum 
    Fisher information matrix (QFIM) with gaussian states.

    Parameters
    ----------
    > **R:** `array` 
        -- First-order moment.

    > **dR:** `list`
        -- Derivatives of the first-order moment on the unknown parameters to be 
        estimated. For example, dR[0] is the derivative vector on the first 
        parameter.

    > **D:** `matrix`
        -- Second-order moment.

    > **dD:** `list`
        -- Derivatives of the second-order moment on the unknown parameters to be 
        estimated. For example, dD[0] is the derivative vector on the first 
        parameter.

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

    Returns
    ----------
    **QFI or QFIM with gaussian states:** `float or matrix`
        -- For single parameter estimation (the length of drho is equal to one), 
        the output is QFI and for multiparameter estimation (the length of drho 
        is more than one), it returns QFIM.
    """

    para_num = len(dR)
    m = int(len(R) / 2)
    QFIM_res = np.zeros([para_num, para_num])

    C = np.array(
        [
            [D[i][j] - R[i] * R[j] for j in range(2 * m)]
            for i in range(2 * m)
        ]
    )
    dC = [
        np.array(
            [
                [
                    dD[k][i][j] - dR[k][i] * R[j] - R[i] * dR[k][j]
                    for j in range(2 * m)
                ]
                for i in range(2 * m)
            ]
        )
        for k in range(para_num)
    ]

    C_sqrt = sqrtm(C)
    J = np.kron([[0, 1], [-1, 0]], np.eye(m))
    B = C_sqrt @ J @ C_sqrt
    P = np.eye(2 * m)
    P = np.vstack([P[:][::2], P[:][1::2]])
    T, Q = schur(B)
    vals = eigvals(B)
    c = vals[::2].imag
    Diag = np.diagflat(c**-0.5)
    S = inv(J @ C_sqrt @ Q @ P @ np.kron([[0, 1], [-1, 0]], -Diag)).T @ P.T

    sx = np.array([[0.0, 1.0], [1.0, 0.0]])
    sy = np.array([[0.0, -1.0j], [1.0j, 0.0]])
    sz = np.array([[1.0, 0.0], [0.0, -1.0]])
    a_Gauss = [1j * sy, sz, np.eye(2), sx]

    es = [
        [np.eye(1, m**2, m * i + j).reshape(m, m) for j in range(m)] for i in range(m)
    ]

    As = [[np.kron(s, a_Gauss[i]) / np.sqrt(2) for s in es] for i in range(4)]
    gs = [
        [[[np.trace(inv(S) @ dC @ inv(S.T) @ aa.T) for aa in a] for a in A] for A in As]
        for dC in dC
    ]
    G = [np.zeros((2 * m, 2 * m)).astype(np.longdouble) for _ in range(para_num)]

    for i in range(para_num):
        for j in range(m):
            for k in range(m):
                for l in range(4):
                    G[i] += np.real(
                        gs[i][l][j][k]
                        / (4 * c[j] * c[k] + (-1) ** (l + 1))
                        * inv(S.T)
                        @ As[l][j][k]
                        @ inv(S)
                    )

    QFIM_res += np.real(
        [
            [np.trace(G[i] @ dC[j]) + dR[i] @ inv(C) @ dR[j] for j in range(para_num)]
            for i in range(para_num)
        ]
    )

    if para_num == 1:
        return QFIM_res[0][0]
    else:
        return QFIM_res

---

Holevo Cramér-Rao bound

Calculation of the Holevo Cramer-Rao bound (HCRB) via the semidefinite program (SDP).

Parameters

rho: matrix -- Density matrix.

drho: list -- Derivatives of the density matrix on the unknown parameters to be estimated. For example, drho[0] is the derivative vector on the first parameter.

W: matrix -- Weight matrix.

eps: float -- Machine epsilon.

Returns

HCRB: float -- The value of Holevo Cramer-Rao bound.

Source code in quanestimation/AsymptoticBound/AnalogCramerRao.py

def HCRB(rho, drho, W, eps=1e-8):
    """
    Calculation of the Holevo Cramer-Rao bound (HCRB) via the semidefinite program (SDP).

    Parameters
    ----------
    > **rho:** `matrix`
        -- Density matrix.

    > **drho:** `list`
        -- Derivatives of the density matrix on the unknown parameters to be 
        estimated. For example, drho[0] is the derivative vector on the first 
        parameter.

    > **W:** `matrix`
        -- Weight matrix.

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

    Returns
    ----------
    **HCRB:** `float`
        -- The value of Holevo Cramer-Rao bound.
    """

    if type(drho) != list:
        raise TypeError("Please make sure drho is a list!")

    if len(drho) == 1:
        print(
            "In single parameter scenario, HCRB is equivalent to QFI. This function will return the value of QFI."
        )
        f = QFIM(rho, drho, eps=eps)
        return f
    elif matrix_rank(W) == 1:
        print(
            "For rank-one weight matrix, the HCRB is equivalent to QFIM. This function will return the value of Tr(WF^{-1})."
        )
        F = QFIM(rho, drho, eps=eps)
        return np.trace(np.dot(W, np.linalg.pinv(F)))
    else:
        dim = len(rho)
        num = dim * dim
        para_num = len(drho)

        Lambda = [np.identity(dim)] + suN_generator(dim)
        Lambda = Lambda / np.sqrt(2)

        vec_drho = [[] for i in range(para_num)]
        for pi in range(para_num):
            vec_drho[pi] = np.array(
                [
                    np.real(np.trace(np.dot(drho[pi], Lambda[i])))
                    for i in range(len(Lambda))
                ]
            )

        S = np.zeros((num, num), dtype=np.complex128)
        for a in range(num):
            for b in range(num):
                S[a][b] = np.trace(np.dot(Lambda[a], np.dot(Lambda[b], rho)))

        accu = len(str(int(1 / eps))) - 1
        lu, d, perm = sp.linalg.ldl(S.round(accu))
        R = np.dot(lu, sp.linalg.sqrtm(d)).conj().T
        # ============optimization variables================
        V = cp.Variable((para_num, para_num))
        X = cp.Variable((num, para_num))
        # ================add constraints===================
        constraints = [cp.bmat([[V, X.T @ R.conj().T], [R @ X, np.identity(num)]]) >> 0]
        for i in range(para_num):
            for j in range(para_num):
                if i == j:
                    constraints += [X[:, i].T @ vec_drho[j] == 1]
                else:
                    constraints += [X[:, i].T @ vec_drho[j] == 0]

        prob = cp.Problem(cp.Minimize(cp.trace(W @ V)), constraints)
        prob.solve()

        return prob.value

---

Nagaoka-Hayashi bound

Calculation of the Nagaoka-Hayashi bound (NHB) via the semidefinite program (SDP).

Parameters

rho: matrix -- Density matrix.

drho: list -- Derivatives of the density matrix on the unknown parameters to be estimated. For example, drho[0] is the derivative vector on the first parameter.

W: matrix -- Weight matrix.

Returns

NHB: float -- The value of Nagaoka-Hayashi bound.

Source code in quanestimation/AsymptoticBound/AnalogCramerRao.py

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

    Parameters
    ----------
    > **rho:** `matrix`
        -- Density matrix.

    > **drho:** `list`
        -- Derivatives of the density matrix on the unknown parameters to be 
        estimated. For example, drho[0] is the derivative vector on the first 
        parameter.

    > **W:** `matrix`
        -- Weight matrix.

    Returns
    ----------
    **NHB:** `float`
        -- The value of Nagaoka-Hayashi bound.
    """
    dim = len(rho)
    para_num = len(drho)

    L_tp = [[[] for i in range(para_num)] for j in range(para_num)]
    for para_i in range(para_num):
        for para_j in range(para_i, para_num):
            L_tp[para_i][para_j] = cp.Variable((dim, dim), hermitian=True)
            L_tp[para_j][para_i] = L_tp[para_i][para_j]
    L = cp.vstack([cp.hstack(L_tp[i]) for i in range(para_num)])
    X = [cp.Variable((dim, dim), hermitian=True) for j in range(para_num)]

    constraints = [cp.bmat([[L, cp.vstack(X)], [cp.hstack(X), np.identity(dim)]])  >> 0]

    for i in range(para_num):
        constraints += [cp.trace(X[i] @ rho) == 0]
        for j in range(para_num):
            if i == j:
                constraints += [cp.trace(X[i] @ drho[j]) == 1]
            else:
                constraints += [cp.trace(X[i] @ drho[j]) == 0]
    prob = cp.Problem(cp.Minimize(cp.real(cp.trace(cp.kron(W, rho) @ L))), constraints)
    prob.solve()

    return prob.value

---

Bayesian Cramér-Rao bounds

Calculation of the Bayesian classical Fisher information (BCFI) and the Bayesian classical Fisher information matrix (BCFIM) of the form \begin{align} \mathcal{I}_{\mathrm{Bayes}}=\int p(\textbf{x})\mathcal{I}\mathrm{d}\textbf{x} \end{align}

with \(\mathcal{I}\) the CFIM and \(p(\textbf{x})\) the prior distribution.

Parameters

x: list -- The regimes of the parameters for the integral.

p: multidimensional array -- The prior distribution.

rho: multidimensional list -- Parameterized density matrix.

drho: multidimensional list -- Derivatives of the parameterized density matrix (rho) with respect to the unknown parameters to be estimated.

M: list of matrices -- A set of positive operator-valued measure (POVM). The default measurement is a set of rank-one symmetric informationally complete POVM (SIC-POVM).

eps: float -- Machine epsilon.

Returns

BCFI or BCFIM: float or matrix -- For single parameter estimation (the length of x is equal to one), the output is BCFI and for multiparameter estimation (the length of x is more than one), it returns BCFIM.

Note: SIC-POVM is calculated by the Weyl-Heisenberg covariant SIC-POVM fiducial state which can be downloaded from here.

Source code in quanestimation/BayesianBound/BayesCramerRao.py

def BCFIM(x, p, rho, drho, M=[], eps=1e-8):
    r"""
    Calculation of the Bayesian classical Fisher information (BCFI) and the 
    Bayesian classical Fisher information matrix (BCFIM) of the form
    \begin{align}
    \mathcal{I}_{\mathrm{Bayes}}=\int p(\textbf{x})\mathcal{I}\mathrm{d}\textbf{x}
    \end{align}

    with $\mathcal{I}$ the CFIM and $p(\textbf{x})$ the prior distribution.

    Parameters
    ----------
    > **x:** `list`
        -- The regimes of the parameters for the integral.

    > **p:** `multidimensional array`
        -- The prior distribution.

    > **rho:** `multidimensional list`
        -- Parameterized density matrix.

    > **drho:** `multidimensional list`
        -- Derivatives of the parameterized density matrix (rho) with respect to the unknown
        parameters to be estimated.

    > **M:** `list of matrices`
        -- A set of positive operator-valued measure (POVM). The default measurement 
        is a set of rank-one symmetric informationally complete POVM (SIC-POVM).

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

    Returns
    ----------
    **BCFI or BCFIM:** `float or matrix`
        -- For single parameter estimation (the length of x is equal to one), the output 
        is BCFI and for multiparameter estimation (the length of x is more than one), 
        it returns BCFIM.

    **Note:** 
        SIC-POVM is calculated by the Weyl-Heisenberg covariant SIC-POVM fiducial state 
        which can be downloaded from [here](http://www.physics.umb.edu/Research/QBism/
        solutions.html).
    """

    para_num = len(x)
    if para_num == 1:
        #### single parameter scenario ####
        if M == []:
            M = SIC(len(rho[0]))
        else:
            if type(M) != list:
                raise TypeError("Please make sure M is a list!")

        p_num = len(p)
        if type(drho[0]) == list:
            drho = [drho[i][0] for i in range(p_num)]
        p_num = len(p)
        F_tp = np.zeros(p_num)
        for m in range(p_num):
            F_tp[m] = CFIM(rho[m], [drho[m]], M=M, eps=eps)

        arr = [p[i] * F_tp[i] for i in range(p_num)]
        return simpson(arr, x[0])
    else:
        #### multiparameter scenario ####
        p_shape = np.shape(p)
        p_ext = extract_ele(p, para_num)
        rho_ext = extract_ele(rho, para_num)
        drho_ext = extract_ele(drho, para_num)

        p_list, rho_list, drho_list = [], [], []
        for p_ele, rho_ele, drho_ele in zip(p_ext, rho_ext, drho_ext):
            p_list.append(p_ele)
            rho_list.append(rho_ele)
            drho_list.append(drho_ele)

        dim = len(rho_list[0])
        if M == []:
            M = SIC(dim)
        else:
            if type(M) != list:
                raise TypeError("Please make sure M is a list!")

        F_list = [
            [[0.0 for i in range(len(p_list))] for j in range(para_num)]
            for k in range(para_num)
        ]
        for i in range(len(p_list)):
            F_tp = CFIM(rho_list[i], drho_list[i], M=M, eps=eps)
            for pj in range(para_num):
                for pk in range(para_num):
                    F_list[pj][pk][i] = F_tp[pj][pk]

        BCFIM_res = np.zeros([para_num, para_num])
        for para_i in range(0, para_num):
            for para_j in range(para_i, para_num):
                F_ij = np.array(F_list[para_i][para_j]).reshape(p_shape)
                arr = p * F_ij
                for si in reversed(range(para_num)):
                    arr = simpson(arr, x[si])
                BCFIM_res[para_i][para_j] = arr
                BCFIM_res[para_j][para_i] = arr
        return BCFIM_res

Calculation of the Bayesian quantum Fisher information (BQFI) and the Bayesian quantum Fisher information matrix (BQFIM) of the form \begin{align} \mathcal{F}_{\mathrm{Bayes}}=\int p(\textbf{x})\mathcal{F}\mathrm{d}\textbf{x} \end{align}

with \(\mathcal{F}\) the QFIM of all types and \(p(\textbf{x})\) the prior distribution.

Parameters

x: list -- The regimes of the parameters for the integral.

p: multidimensional array -- The prior distribution.

rho: multidimensional list -- Parameterized density matrix.

drho: multidimensional list -- Derivatives of the parameterized density matrix (rho) with respect to the unknown parameters to be estimated.

LDtype: string -- Types of QFI (QFIM) can be set as the objective function. Options are:
"SLD" (default) -- QFI (QFIM) based on symmetric logarithmic derivative (SLD).
"RLD" -- QFI (QFIM) based on right logarithmic derivative (RLD).
"LLD" -- QFI (QFIM) based on left logarithmic derivative (LLD).

eps: float -- Machine epsilon.

Returns

BQFI or BQFIM: float or matrix -- For single parameter estimation (the length of x is equal to one), the output is BQFI and for multiparameter estimation (the length of x is more than one), it returns BQFIM.

Source code in quanestimation/BayesianBound/BayesCramerRao.py

def BQFIM(x, p, rho, drho, LDtype="SLD", eps=1e-8):
    r"""
    Calculation of the Bayesian quantum Fisher information (BQFI) and the 
    Bayesian quantum Fisher information matrix (BQFIM) of the form
    \begin{align}
    \mathcal{F}_{\mathrm{Bayes}}=\int p(\textbf{x})\mathcal{F}\mathrm{d}\textbf{x}
    \end{align}

    with $\mathcal{F}$ the QFIM of all types and $p(\textbf{x})$ the prior distribution.

    Parameters
    ----------
    > **x:** `list`
        -- The regimes of the parameters for the integral.

    > **p:** `multidimensional array`
        -- The prior distribution.

    > **rho:** `multidimensional list`
        -- Parameterized density matrix.

    > **drho:** `multidimensional list`
        -- Derivatives of the parameterized density matrix (rho) with respect to the unknown
        parameters to be estimated.

    > **LDtype:** `string`
        -- Types of QFI (QFIM) can be set as the objective function. Options are:  
        "SLD" (default) -- QFI (QFIM) based on symmetric logarithmic derivative (SLD).  
        "RLD" -- QFI (QFIM) based on right logarithmic derivative (RLD).  
        "LLD" -- QFI (QFIM) based on left logarithmic derivative (LLD).

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

    Returns
    ----------
    **BQFI or BQFIM:** `float or matrix`
        -- For single parameter estimation (the length of x is equal to one), the output 
        is BQFI and for multiparameter estimation (the length of x is more than one), 
        it returns BQFIM.
    """

    para_num = len(x)
    if para_num == 1:
        #### single parameter scenario ####
        p_num = len(p)
        if type(drho[0]) == list:
            drho = [drho[i][0] for i in range(p_num)]

        F_tp = np.zeros(p_num)
        for m in range(p_num):
            F_tp[m] = QFIM(rho[m], [drho[m]], LDtype=LDtype, eps=eps)
        arr = [p[i] * F_tp[i] for i in range(p_num)]
        return simpson(arr, x[0])
    else:
        #### multiparameter scenario ####
        p_shape = np.shape(p)
        p_ext = extract_ele(p, para_num)
        rho_ext = extract_ele(rho, para_num)
        drho_ext = extract_ele(drho, para_num)

        p_list, rho_list, drho_list = [], [], []
        for p_ele, rho_ele, drho_ele in zip(p_ext, rho_ext, drho_ext):
            p_list.append(p_ele)
            rho_list.append(rho_ele)
            drho_list.append(drho_ele)

        F_list = [
            [[0.0 for i in range(len(p_list))] for j in range(para_num)]
            for k in range(para_num)
        ]
        for i in range(len(p_list)):
            F_tp = QFIM(rho_list[i], drho_list[i], LDtype=LDtype, eps=eps)
            for pj in range(para_num):
                for pk in range(para_num):
                    F_list[pj][pk][i] = F_tp[pj][pk]

        BQFIM_res = np.zeros([para_num, para_num])
        for para_i in range(0, para_num):
            for para_j in range(para_i, para_num):
                F_ij = np.array(F_list[para_i][para_j]).reshape(p_shape)
                arr = p * F_ij
                for si in reversed(range(para_num)):
                    arr = simpson(arr, x[si])
                BQFIM_res[para_i][para_j] = arr
                BQFIM_res[para_j][para_i] = arr
        return BQFIM_res

Calculation of the Bayesian Cramer-Rao bound (BCRB). The covariance matrix with a prior distribution \(p(\textbf{x})\) is defined as \begin{align} \mathrm{cov}(\hat{\textbf{x}},{\Pi_y})=\int p(\textbf{x})\sum_y\mathrm{Tr} (\rho\Pi_y)(\hat{\textbf{x}}-\textbf{x})(\hat{\textbf{x}}-\textbf{x})^{\mathrm{T}} \mathrm{d}\textbf{x} \end{align}

where \(\textbf{x}=(x_0,x_1,\dots)^{\mathrm{T}}\) are the unknown parameters to be estimated and the integral \(\int\mathrm{d}\textbf{x}:=\iiint\mathrm{d}x_0\mathrm{d}x_1\cdots\). \(\{\Pi_y\}\) is a set of positive operator-valued measure (POVM) and \(\rho\) represents the parameterized density matrix.

This function calculates three types BCRB. The first one is \begin{align} \mathrm{cov}(\hat{\textbf{x}},{\Pi_y})\geq \int p(\textbf{x})\left(B\mathcal{I}^{-1}B +\textbf{b}\textbf{b}^{\mathrm{T}}\right)\mathrm{d}\textbf{x}, \end{align}

where \(\textbf{b}\) and \(\textbf{b}'\) are the vectors of biase and its derivatives on parameters. \(B\) is a diagonal matrix with the \(i\)th entry \(B_{ii}=1+[\textbf{b}']_{i}\) and \(\mathcal{I}\) is the CFIM.

The second one is \begin{align} \mathrm{cov}(\hat{\textbf{x}},{\Pi_y})\geq \mathcal{B}\,\mathcal{I}_{\mathrm{Bayes}}^{-1}\, \mathcal{B}+\int p(\textbf{x})\textbf{b}\textbf{b}^{\mathrm{T}}\mathrm{d}\textbf{x}, \end{align}

where \(\mathcal{B}=\int p(\textbf{x})B\mathrm{d}\textbf{x}\) is the average of \(B\) and \(\mathcal{I}_{\mathrm{Bayes}}=\int p(\textbf{x})\mathcal{I}\mathrm{d}\textbf{x}\) is the average CFIM.

The third one is \begin{align} \mathrm{cov}(\hat{\textbf{x}},{\Pi_y})\geq \int p(\textbf{x}) \mathcal{G}\left(\mathcal{I}_p+\mathcal{I}\right)^{-1}\mathcal{G}^{\mathrm{T}}\mathrm{d}\textbf{x} \end{align}

with \([\mathcal{I}_{p}]_{ab}:=[\partial_a \ln p(\textbf{x})][\partial_b \ln p(\textbf{x})]\) and \(\mathcal{G}_{ab}:=[\partial_b\ln p(\textbf{x})][\textbf{b}]_a+B_{aa}\delta_{ab}\).

Parameters

x: list -- The regimes of the parameters for the integral.

p: multidimensional array -- The prior distribution.

rho: multidimensional list -- Parameterized density matrix.

drho: multidimensional list -- Derivatives of the parameterized density matrix (rho) with respect to the unknown parameters to be estimated.

M: list of matrices -- A set of positive operator-valued measure (POVM). The default measurement is a set of rank-one symmetric informationally complete POVM (SIC-POVM).

b: list -- Vector of biases of the form \(\textbf{b}=(b(x_0),b(x_1),\dots)^{\mathrm{T}}\).

db: list -- Derivatives of b with respect to the unknown parameters to be estimated, It should be expressed as \(\textbf{b}'=(\partial_0 b(x_0),\partial_1 b(x_1),\dots)^{\mathrm{T}}\).

btype: int (1, 2, 3) -- Types of the BCRB. Options are:
1 (default) -- It means to calculate the first type of the BCRB.
2 -- It means to calculate the second type of the BCRB. 3 -- It means to calculate the third type of the BCRB.

eps: float -- Machine epsilon.

Returns

BCRB: float or matrix -- For single parameter estimation (the length of x is equal to one), the output is a float and for multiparameter estimation (the length of x is more than one), it returns a matrix.

Note: SIC-POVM is calculated by the Weyl-Heisenberg covariant SIC-POVM fiducial state which can be downloaded from here.

Source code in quanestimation/BayesianBound/BayesCramerRao.py

def BCRB(x, p, dp, rho, drho, M=[], b=[], db=[], btype=1, eps=1e-8):
    r"""
    Calculation of the Bayesian Cramer-Rao bound (BCRB). The covariance matrix 
    with a prior distribution $p(\textbf{x})$ is defined as
    \begin{align}
    \mathrm{cov}(\hat{\textbf{x}},\{\Pi_y\})=\int p(\textbf{x})\sum_y\mathrm{Tr}
    (\rho\Pi_y)(\hat{\textbf{x}}-\textbf{x})(\hat{\textbf{x}}-\textbf{x})^{\mathrm{T}}
    \mathrm{d}\textbf{x}
    \end{align}

    where $\textbf{x}=(x_0,x_1,\dots)^{\mathrm{T}}$ are the unknown parameters to be estimated
    and the integral $\int\mathrm{d}\textbf{x}:=\iiint\mathrm{d}x_0\mathrm{d}x_1\cdots$.
    $\{\Pi_y\}$ is a set of positive operator-valued measure (POVM) and $\rho$ represents 
    the parameterized density matrix.

    This function calculates three types BCRB. The first one is 
    \begin{align}
    \mathrm{cov}(\hat{\textbf{x}},\{\Pi_y\})\geq \int p(\textbf{x})\left(B\mathcal{I}^{-1}B
    +\textbf{b}\textbf{b}^{\mathrm{T}}\right)\mathrm{d}\textbf{x},
    \end{align}

    where $\textbf{b}$ and $\textbf{b}'$ are the vectors of biase and its derivatives on parameters.
    $B$ is a diagonal matrix with the $i$th entry $B_{ii}=1+[\textbf{b}']_{i}$ and $\mathcal{I}$ 
    is the CFIM.

    The second one is
    \begin{align}
    \mathrm{cov}(\hat{\textbf{x}},\{\Pi_y\})\geq \mathcal{B}\,\mathcal{I}_{\mathrm{Bayes}}^{-1}\,
    \mathcal{B}+\int p(\textbf{x})\textbf{b}\textbf{b}^{\mathrm{T}}\mathrm{d}\textbf{x},
    \end{align}

    where $\mathcal{B}=\int p(\textbf{x})B\mathrm{d}\textbf{x}$ is the average of $B$ and 
    $\mathcal{I}_{\mathrm{Bayes}}=\int p(\textbf{x})\mathcal{I}\mathrm{d}\textbf{x}$ is 
    the average CFIM.

    The third one is
    \begin{align}
    \mathrm{cov}(\hat{\textbf{x}},\{\Pi_y\})\geq \int p(\textbf{x})
    \mathcal{G}\left(\mathcal{I}_p+\mathcal{I}\right)^{-1}\mathcal{G}^{\mathrm{T}}\mathrm{d}\textbf{x}
    \end{align}

    with $[\mathcal{I}_{p}]_{ab}:=[\partial_a \ln p(\textbf{x})][\partial_b \ln p(\textbf{x})]$ and
    $\mathcal{G}_{ab}:=[\partial_b\ln p(\textbf{x})][\textbf{b}]_a+B_{aa}\delta_{ab}$.

    Parameters
    ----------
    > **x:** `list`
        -- The regimes of the parameters for the integral.

    > **p:** `multidimensional array`
        -- The prior distribution.

    > **rho:** `multidimensional list`
        -- Parameterized density matrix.

    > **drho:** `multidimensional list`
        -- Derivatives of the parameterized density matrix (rho) with respect to the unknown
        parameters to be estimated.

    > **M:** `list of matrices`
        -- A set of positive operator-valued measure (POVM). The default measurement 
        is a set of rank-one symmetric informationally complete POVM (SIC-POVM).

    > **b:** `list`
        -- Vector of biases of the form $\textbf{b}=(b(x_0),b(x_1),\dots)^{\mathrm{T}}$.

    > **db:** `list`
        -- Derivatives of b with respect to the unknown parameters to be estimated, It should be 
        expressed as $\textbf{b}'=(\partial_0 b(x_0),\partial_1 b(x_1),\dots)^{\mathrm{T}}$.

    > **btype:** `int (1, 2, 3)`
        -- Types of the BCRB. Options are:  
        1 (default) -- It means to calculate the first type of the BCRB.  
        2 -- It means to calculate the second type of the BCRB.
        3 -- It means to calculate the third type of the BCRB.

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

    Returns
    ----------
    **BCRB:** `float or matrix`
        -- For single parameter estimation (the length of x is equal to one), the 
        output is a float and for multiparameter estimation (the length of x is 
        more than one), it returns a matrix.

    **Note:** 
        SIC-POVM is calculated by the Weyl-Heisenberg covariant SIC-POVM fiducial state 
        which can be downloaded from [here](http://www.physics.umb.edu/Research/QBism/
        solutions.html).
    """

    para_num = len(x)
    if para_num == 1:
        #### single parameter scenario ####
        p_num = len(p)
        if not b:
            b = np.zeros(p_num)
            db = np.zeros(p_num)
        elif not db:
            db = np.zeros(p_num)

        if M == []:
            M = SIC(len(rho[0]))
        else:
            if type(M) != list:
                raise TypeError("Please make sure M is a list!")

        if type(drho[0]) == list:
            drho = [drho[i][0] for i in range(p_num)]
        if type(b[0]) == list or type(b[0]) == np.ndarray:
            b = b[0]
        if type(db[0]) == list or type(db[0]) == np.ndarray:
            db = db[0]

        F_tp = np.zeros(p_num)
        for m in range(p_num):
            F_tp[m] = CFIM(rho[m], [drho[m]], M=M, eps=eps)

        if btype == 1:
            arr = [
                p[i] * ((1 + db[i]) ** 2 / F_tp[i] + b[i] ** 2) for i in range(p_num)
            ]
            F = simpson(arr, x[0])
            return F
        elif btype == 2:
            arr = [p[i] * F_tp[i] for i in range(p_num)]
            F1 = simpson(arr, x[0])
            arr2 = [p[j] * (1 + db[j]) for j in range(p_num)]
            B = simpson(arr2, x[0])
            arr3 = [p[k] * b[k] ** 2 for k in range(p_num)]
            bb = simpson(arr3, x[0])
            F = B**2 / F1 + bb
            return F
        elif btype == 3:
            I_tp = [np.real(dp[i] * dp[i] / p[i] ** 2) for i in range(p_num)]
            arr = [p[j]*(dp[j]*b[j]/p[j]+(1 + db[j]))**2 / (I_tp[j] + F_tp[j]) for j in range(p_num)]
            F = simpson(arr, x[0])
            return F
        else:
            raise NameError("NameError: btype should be choosen in {1, 2, 3}.")
    else:
        #### multiparameter scenario ####
        if not b:
            b, db = [], []
            for i in range(para_num):
                b.append(np.zeros(len(x[i])))
                db.append(np.zeros(len(x[i])))
        elif not db:
            db = []
            for i in range(para_num):
                db.append(np.zeros(len(x[i])))

        p_shape = np.shape(p)
        p_ext = extract_ele(p, para_num)
        dp_ext = extract_ele(dp, para_num)
        rho_ext = extract_ele(rho, para_num)
        drho_ext = extract_ele(drho, para_num)
        b_pro = product(*b)
        db_pro = product(*db)

        p_list, rho_list, drho_list = [], [], []
        for p_ele, rho_ele, drho_ele in zip(p_ext, rho_ext, drho_ext):
            p_list.append(p_ele)
            rho_list.append(rho_ele)
            drho_list.append(drho_ele)
        dp_list = [dpi for dpi in dp_ext]

        b_list, db_list = [], []
        for b_ele, db_ele in zip(b_pro, db_pro):
            b_list.append([b_ele[i] for i in range(para_num)])
            db_list.append([db_ele[j] for j in range(para_num)])

        dim = len(rho_list[0])
        if M == []:
            M = SIC(dim)
        else:
            if type(M) != list:
                raise TypeError("Please make sure M is a list!")
        if btype == 1:
            F_list = [
                [[0.0 for i in range(len(p_list))] for j in range(para_num)]
                for k in range(para_num)
            ]
            for i in range(len(p_list)):
                F_tp = CFIM(rho_list[i], drho_list[i], M=M, eps=eps)
                F_inv = np.linalg.pinv(F_tp)
                B = np.diag([(1.0 + db_list[i][j]) for j in range(para_num)])
                term1 = np.dot(B, np.dot(F_inv, B))
                term2 = np.dot(
                    np.array(b_list[i]).reshape(para_num, 1),
                    np.array(b_list[i]).reshape(1, para_num),
                )
                for pj in range(para_num):
                    for pk in range(para_num):
                        F_list[pj][pk][i] = term1[pj][pk] + term2[pj][pk]

            res = np.zeros([para_num, para_num])
            for para_i in range(0, para_num):
                for para_j in range(para_i, para_num):
                    F_ij = np.array(F_list[para_i][para_j]).reshape(p_shape)
                    arr = p * F_ij
                    for si in reversed(range(para_num)):
                        arr = simpson(arr, x[si])
                    res[para_i][para_j] = arr
                    res[para_j][para_i] = arr
            return res
        elif btype == 2:
            F_list = [
                [[0.0 for i in range(len(p_list))] for j in range(para_num)]
                for k in range(para_num)
            ]
            B_list = [
                [[0.0 for i in range(len(p_list))] for j in range(para_num)]
                for k in range(para_num)
            ]
            bb_list = [
                [[0.0 for i in range(len(p_list))] for j in range(para_num)]
                for k in range(para_num)
            ]
            for i in range(len(p_list)):
                F_tp = CFIM(rho_list[i], drho_list[i], M=M, eps=eps)
                B_tp = np.diag([(1.0 + db_list[i][j]) for j in range(para_num)])
                bb_tp = np.dot(
                    np.array(b_list[i]).reshape(para_num, 1),
                    np.array(b_list[i]).reshape(1, para_num),
                )
                for pj in range(para_num):
                    for pk in range(para_num):
                        F_list[pj][pk][i] = F_tp[pj][pk]
                        B_list[pj][pk][i] = B_tp[pj][pk]
                        bb_list[pj][pk][i] = bb_tp[pj][pk]

            F_res = np.zeros([para_num, para_num])
            for para_i in range(0, para_num):
                for para_j in range(para_i, para_num):
                    F_ij = np.array(F_list[para_i][para_j]).reshape(p_shape)
                    arr = p * F_ij
                    for si in reversed(range(para_num)):
                        arr = simpson(arr, x[si])
                    F_res[para_i][para_j] = arr
                    F_res[para_j][para_i] = arr
            B_res = np.zeros([para_num, para_num])
            bb_res = np.zeros([para_num, para_num])
            for para_m in range(para_num):
                for para_n in range(para_num):
                    B_mn = np.array(B_list[para_m][para_n]).reshape(p_shape)
                    bb_mn = np.array(bb_list[para_m][para_n]).reshape(p_shape)
                    arr2 = p * B_mn
                    arr3 = p * bb_mn
                    for sj in reversed(range(para_num)):
                        arr2 = simpson(arr2, x[sj])
                        arr3 = simpson(arr3, x[sj])
                    B_res[para_m][para_n] = arr2
                    bb_res[para_m][para_n] = arr3
            res = np.dot(B_res, np.dot(np.linalg.pinv(F_res), B_res)) + bb_res
            return res
        elif btype == 3:
            F_list = [
                [[0.0 for i in range(len(p_list))] for j in range(para_num)]
                for k in range(para_num)
            ]
            for i in range(len(p_list)):
                F_tp = CFIM(rho_list[i], drho_list[i], M=M, eps=eps)
                I_tp = np.zeros((para_num, para_num))
                G_tp = np.zeros((para_num, para_num))
                for pm in range(para_num):
                    for pn in range(para_num):
                        if pm == pn:
                            G_tp[pm][pn] = dp_list[i][pn]*b_list[i][pm]/p_list[i]+(1.0 + db_list[i][pm])
                        else:
                            G_tp[pm][pn] = dp_list[i][pn]*b_list[i][pm]/p_list[i]
                        I_tp[pm][pn] = dp_list[i][pm] * dp_list[i][pn] / p_list[i] ** 2

                F_tot = np.dot(G_tp, np.dot(np.linalg.pinv(F_tp + I_tp), G_tp.T))
                for pj in range(para_num):
                    for pk in range(para_num):
                        F_list[pj][pk][i] = F_tot[pj][pk]

            res = np.zeros([para_num, para_num])
            for para_i in range(0, para_num):
                for para_j in range(para_i, para_num):
                    F_ij = np.array(F_list[para_i][para_j]).reshape(p_shape)
                    arr = p * F_ij
                    for si in reversed(range(para_num)):
                        arr = simpson(arr, x[si])
                    res[para_i][para_j] = arr
                    res[para_j][para_i] = arr
            return res
        else:
            raise NameError("NameError: btype should be choosen in {1, 2, 3}.")

Calculation of the Bayesian quantum Cramer-Rao bound (BQCRB). The covariance matrix with a prior distribution \(p(\textbf{x})\) is defined as \begin{align} \mathrm{cov}(\hat{\textbf{x}},{\Pi_y})=\int p(\textbf{x})\sum_y\mathrm{Tr} (\rho\Pi_y)(\hat{\textbf{x}}-\textbf{x})(\hat{\textbf{x}}-\textbf{x})^{\mathrm{T}} \mathrm{d}\textbf{x} \end{align}

where \(\textbf{x}=(x_0,x_1,\dots)^{\mathrm{T}}\) are the unknown parameters to be estimated and the integral \(\int\mathrm{d}\textbf{x}:=\iiint\mathrm{d}x_0\mathrm{d}x_1\cdots\). \(\{\Pi_y\}\) is a set of positive operator-valued measure (POVM) and \(\rho\) represent the parameterized density matrix.

This function calculates three types of the BQCRB. The first one is \begin{align} \mathrm{cov}(\hat{\textbf{x}},{\Pi_y})\geq\int p(\textbf{x})\left(B\mathcal{F}^{-1}B +\textbf{b}\textbf{b}^{\mathrm{T}}\right)\mathrm{d}\textbf{x}, \end{align}

where \(\textbf{b}\) and \(\textbf{b}'\) are the vectors of biase and its derivatives on parameters. \(B\) is a diagonal matrix with the \(i\)th entry \(B_{ii}=1+[\textbf{b}']_{i}\) and \(\mathcal{F}\) is the QFIM for all types.

The second one is \begin{align} \mathrm{cov}(\hat{\textbf{x}},{\Pi_y})\geq \mathcal{B}\,\mathcal{F}_{\mathrm{Bayes}}^{-1}\, \mathcal{B}+\int p(\textbf{x})\textbf{b}\textbf{b}^{\mathrm{T}}\mathrm{d}\textbf{x}, \end{align}

where \(\mathcal{B}=\int p(\textbf{x})B\mathrm{d}\textbf{x}\) is the average of \(B\) and \(\mathcal{F}_{\mathrm{Bayes}}=\int p(\textbf{x})\mathcal{F}\mathrm{d}\textbf{x}\) is the average QFIM.

The third one is \begin{align} \mathrm{cov}(\hat{\textbf{x}},{\Pi_y})\geq \int p(\textbf{x}) \mathcal{G}\left(\mathcal{I}_p+\mathcal{F}\right)^{-1}\mathcal{G}^{\mathrm{T}}\mathrm{d}\textbf{x} \end{align}

with \([\mathcal{I}_{p}]_{ab}:=[\partial_a \ln p(\textbf{x})][\partial_b \ln p(\textbf{x})]\) and \(\mathcal{G}_{ab}:=[\partial_b\ln p(\textbf{x})][\textbf{b}]_a+B_{aa}\delta_{ab}\).

Parameters

x: list -- The regimes of the parameters for the integral.

p: multidimensional array -- The prior distribution.

rho: multidimensional list -- Parameterized density matrix.

drho: multidimensional list -- Derivatives of the parameterized density matrix (rho) with respect to the unknown parameters to be estimated.

b: list -- Vector of biases of the form \(\textbf{b}=(b(x_0),b(x_1),\dots)^{\mathrm{T}}\).

db: list -- Derivatives of b with respect to the unknown parameters to be estimated, It should be expressed as \(\textbf{b}'=(\partial_0 b(x_0),\partial_1 b(x_1),\dots)^{\mathrm{T}}\).

btype: int (1, 2, 3) -- Types of the BQCRB. Options are:
1 (default) -- It means to calculate the first type of the BQCRB.
2 -- It means to calculate the second type of the BQCRB. 3 -- It means to calculate the third type of the BCRB.

LDtype: string -- Types of QFI (QFIM) can be set as the objective function. Options are:
"SLD" (default) -- QFI (QFIM) based on symmetric logarithmic derivative (SLD).
"RLD" -- QFI (QFIM) based on right logarithmic derivative (RLD).
"LLD" -- QFI (QFIM) based on left logarithmic derivative (LLD).

eps: float -- Machine epsilon.

Returns

BQCRB: float or matrix -- For single parameter estimation (the length of x is equal to one), the output is a float and for multiparameter estimation (the length of x is more than one), it returns a matrix.

Source code in quanestimation/BayesianBound/BayesCramerRao.py

def BQCRB(x, p, dp, rho, drho, b=[], db=[], btype=1, LDtype="SLD", eps=1e-8):
    r"""
    Calculation of the Bayesian quantum Cramer-Rao bound (BQCRB). The covariance matrix 
    with a prior distribution $p(\textbf{x})$ is defined as
    \begin{align}
    \mathrm{cov}(\hat{\textbf{x}},\{\Pi_y\})=\int p(\textbf{x})\sum_y\mathrm{Tr}
    (\rho\Pi_y)(\hat{\textbf{x}}-\textbf{x})(\hat{\textbf{x}}-\textbf{x})^{\mathrm{T}}
    \mathrm{d}\textbf{x}
    \end{align}

    where $\textbf{x}=(x_0,x_1,\dots)^{\mathrm{T}}$ are the unknown parameters to be estimated
    and the integral $\int\mathrm{d}\textbf{x}:=\iiint\mathrm{d}x_0\mathrm{d}x_1\cdots$.
    $\{\Pi_y\}$ is a set of positive operator-valued measure (POVM) and $\rho$ represent 
    the parameterized density matrix.

    This function calculates three types of the BQCRB. The first one is
    \begin{align}
    \mathrm{cov}(\hat{\textbf{x}},\{\Pi_y\})\geq\int p(\textbf{x})\left(B\mathcal{F}^{-1}B
    +\textbf{b}\textbf{b}^{\mathrm{T}}\right)\mathrm{d}\textbf{x},
    \end{align}

    where $\textbf{b}$ and $\textbf{b}'$ are the vectors of biase and its derivatives on parameters.
    $B$ is a diagonal matrix with the $i$th entry $B_{ii}=1+[\textbf{b}']_{i}$ and $\mathcal{F}$
    is the QFIM for all types.

    The second one is
    \begin{align}
    \mathrm{cov}(\hat{\textbf{x}},\{\Pi_y\})\geq \mathcal{B}\,\mathcal{F}_{\mathrm{Bayes}}^{-1}\,
    \mathcal{B}+\int p(\textbf{x})\textbf{b}\textbf{b}^{\mathrm{T}}\mathrm{d}\textbf{x},
    \end{align}

    where $\mathcal{B}=\int p(\textbf{x})B\mathrm{d}\textbf{x}$ is the average of $B$ and 
    $\mathcal{F}_{\mathrm{Bayes}}=\int p(\textbf{x})\mathcal{F}\mathrm{d}\textbf{x}$ is 
    the average QFIM.

    The third one is
    \begin{align}
    \mathrm{cov}(\hat{\textbf{x}},\{\Pi_y\})\geq \int p(\textbf{x})
    \mathcal{G}\left(\mathcal{I}_p+\mathcal{F}\right)^{-1}\mathcal{G}^{\mathrm{T}}\mathrm{d}\textbf{x}
    \end{align}

    with $[\mathcal{I}_{p}]_{ab}:=[\partial_a \ln p(\textbf{x})][\partial_b \ln p(\textbf{x})]$ and
    $\mathcal{G}_{ab}:=[\partial_b\ln p(\textbf{x})][\textbf{b}]_a+B_{aa}\delta_{ab}$.

    Parameters
    ----------
    > **x:** `list`
        -- The regimes of the parameters for the integral.

    > **p:** `multidimensional array`
        -- The prior distribution.

    > **rho:** `multidimensional list`
        -- Parameterized density matrix.

    > **drho:** `multidimensional list`
        -- Derivatives of the parameterized density matrix (rho) with respect to the unknown
        parameters to be estimated.

    > **b:** `list`
        -- Vector of biases of the form $\textbf{b}=(b(x_0),b(x_1),\dots)^{\mathrm{T}}$.

    > **db:** `list`
        -- Derivatives of b with respect to the unknown parameters to be estimated, It should be 
        expressed as $\textbf{b}'=(\partial_0 b(x_0),\partial_1 b(x_1),\dots)^{\mathrm{T}}$.

    > **btype:** `int (1, 2, 3)`
        -- Types of the BQCRB. Options are:  
        1 (default) -- It means to calculate the first type of the BQCRB.  
        2 -- It means to calculate the second type of the BQCRB.
        3 -- It means to calculate the third type of the BCRB.

    > **LDtype:** `string`
        -- Types of QFI (QFIM) can be set as the objective function. Options are:  
        "SLD" (default) -- QFI (QFIM) based on symmetric logarithmic derivative (SLD).  
        "RLD" -- QFI (QFIM) based on right logarithmic derivative (RLD).  
        "LLD" -- QFI (QFIM) based on left logarithmic derivative (LLD).

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

    Returns
    ----------
    **BQCRB:** `float or matrix`
        -- For single parameter estimation (the length of x is equal to one), the 
        output is a float and for multiparameter estimation (the length of x is 
        more than one), it returns a matrix.
    """

    para_num = len(x)

    if para_num == 1:
        #### single parameter scenario ####
        p_num = len(p)

        if not b:
            b = np.zeros(p_num)
            db = np.zeros(p_num)
        elif not db:
            db = np.zeros(p_num)

        if type(drho[0]) == list:
            drho = [drho[i][0] for i in range(p_num)]
        if type(b[0]) == list or type(b[0]) == np.ndarray:
            b = b[0]
        if type(db[0]) == list or type(db[0]) == np.ndarray:
            db = db[0]

        F_tp = np.zeros(p_num)
        for m in range(p_num):
            F_tp[m] = QFIM(rho[m], [drho[m]], LDtype=LDtype, eps=eps)

        if btype == 1:
            arr = [
                p[i] * ((1 + db[i]) ** 2 / F_tp[i] + b[i] ** 2) for i in range(p_num)
            ]
            F = simpson(arr, x[0])
            return F
        elif btype == 2:
            arr2 = [p[i] * F_tp[i] for i in range(p_num)]
            F2 = simpson(arr2, x[0])
            arr2 = [p[j] * (1 + db[j]) for j in range(p_num)]
            B = simpson(arr2, x[0])
            arr3 = [p[k] * b[k] ** 2 for k in range(p_num)]
            bb = simpson(arr3, x[0])
            F = B**2 / F2 + bb
            return F
        elif btype == 3:
            I_tp = [np.real(dp[i] * dp[i] / p[i] ** 2) for i in range(p_num)]
            arr = [p[j]*(dp[j]*b[j]/p[j]+(1 + db[j]))**2 / (I_tp[j] + F_tp[j]) for j in range(p_num)]
            F = simpson(arr, x[0])
            return F
        else:
            raise NameError("NameError: btype should be choosen in {1, 2, 3}.")
    else:
        #### multiparameter scenario ####
        if not b:
            b, db = [], []
            for i in range(para_num):
                b.append(np.zeros(len(x[i])))
                db.append(np.zeros(len(x[i])))
        elif not db:
            db = []
            for i in range(para_num):
                db.append(np.zeros(len(x[i])))

        p_shape = np.shape(p)
        p_ext = extract_ele(p, para_num)
        dp_ext = extract_ele(dp, para_num)
        rho_ext = extract_ele(rho, para_num)
        drho_ext = extract_ele(drho, para_num)
        b_pro = product(*b)
        db_pro = product(*db)

        p_list, rho_list, drho_list = [], [], []
        for p_ele, rho_ele, drho_ele in zip(p_ext, rho_ext, drho_ext):
            p_list.append(p_ele)
            rho_list.append(rho_ele)
            drho_list.append(drho_ele)
        dp_list = [dpi for dpi in dp_ext]

        b_list, db_list = [], []
        for b_ele, db_ele in zip(b_pro, db_pro):
            b_list.append([b_ele[i] for i in range(para_num)])
            db_list.append([db_ele[j] for j in range(para_num)])

        if btype == 1:
            F_list = [
                [[0.0 for i in range(len(p_list))] for j in range(para_num)]
                for k in range(para_num)
            ]
            for i in range(len(p_list)):
                F_tp = QFIM(rho_list[i], drho_list[i], LDtype=LDtype, eps=eps)
                F_inv = np.linalg.pinv(F_tp)
                B = np.diag([(1.0 + db_list[i][j]) for j in range(para_num)])
                term1 = np.dot(B, np.dot(F_inv, B))
                term2 = np.dot(
                    np.array(b_list[i]).reshape(para_num, 1),
                    np.array(b_list[i]).reshape(1, para_num),
                )
                for pj in range(para_num):
                    for pk in range(para_num):
                        F_list[pj][pk][i] = term1[pj][pk] + term2[pj][pk]

            res = np.zeros([para_num, para_num])
            for para_i in range(0, para_num):
                for para_j in range(para_i, para_num):
                    F_ij = np.array(F_list[para_i][para_j]).reshape(p_shape)
                    arr = p * F_ij
                    for si in reversed(range(para_num)):
                        arr = simpson(arr, x[si])
                    res[para_i][para_j] = arr
                    res[para_j][para_i] = arr
            return res
        elif btype == 2:
            F_list = [
                [[0.0 for i in range(len(p_list))] for j in range(para_num)]
                for k in range(para_num)
            ]
            B_list = [
                [[0.0 for i in range(len(p_list))] for j in range(para_num)]
                for k in range(para_num)
            ]
            bb_list = [
                [[0.0 for i in range(len(p_list))] for j in range(para_num)]
                for k in range(para_num)
            ]
            for i in range(len(p_list)):
                F_tp = QFIM(rho_list[i], drho_list[i], LDtype=LDtype, eps=eps)
                B_tp = np.diag([(1.0 + db_list[i][j]) for j in range(para_num)])
                bb_tp = np.dot(
                    np.array(b_list[i]).reshape(para_num, 1),
                    np.array(b_list[i]).reshape(1, para_num),
                )
                for pj in range(para_num):
                    for pk in range(para_num):
                        F_list[pj][pk][i] = F_tp[pj][pk]
                        B_list[pj][pk][i] = B_tp[pj][pk]
                        bb_list[pj][pk][i] = bb_tp[pj][pk]

            F_res = np.zeros([para_num, para_num])
            for para_i in range(0, para_num):
                for para_j in range(para_i, para_num):
                    F_ij = np.array(F_list[para_i][para_j]).reshape(p_shape)
                    arr = p * F_ij
                    for si in reversed(range(para_num)):
                        arr = simpson(arr, x[si])
                    F_res[para_i][para_j] = arr
                    F_res[para_j][para_i] = arr
            B_res = np.zeros([para_num, para_num])
            bb_res = np.zeros([para_num, para_num])
            for para_m in range(para_num):
                for para_n in range(para_num):
                    B_mn = np.array(B_list[para_m][para_n]).reshape(p_shape)
                    bb_mn = np.array(bb_list[para_m][para_n]).reshape(p_shape)
                    arr2 = p * B_mn
                    arr3 = p * bb_mn
                    for sj in reversed(range(para_num)):
                        arr2 = simpson(arr2, x[sj])
                        arr3 = simpson(arr3, x[sj])
                    B_res[para_m][para_n] = arr2
                    bb_res[para_m][para_n] = arr3
            res = np.dot(B_res, np.dot(np.linalg.pinv(F_res), B_res)) + bb_res
            return res
        elif btype == 3:
            F_list = [
                [[0.0 for i in range(len(p_list))] for j in range(para_num)]
                for k in range(para_num)
            ]
            for i in range(len(p_list)):
                F_tp = QFIM(rho_list[i], drho_list[i], LDtype=LDtype, eps=eps)
                I_tp = np.zeros((para_num, para_num))
                G_tp = np.zeros((para_num, para_num))
                for pm in range(para_num):
                    for pn in range(para_num):
                        if pm == pn:
                            G_tp[pm][pn] = dp_list[i][pn]*b_list[i][pm]/p_list[i]+(1.0 + db_list[i][pm])
                        else:
                            G_tp[pm][pn] = dp_list[i][pn]*b_list[i][pm]/p_list[i]
                        I_tp[pm][pn] = dp_list[i][pm] * dp_list[i][pn] / p_list[i] ** 2

                F_tot = np.dot(G_tp, np.dot(np.linalg.pinv(F_tp + I_tp), G_tp.T))
                for pj in range(para_num):
                    for pk in range(para_num):
                        F_list[pj][pk][i] = F_tot[pj][pk]

            res = np.zeros([para_num, para_num])
            for para_i in range(0, para_num):
                for para_j in range(para_i, para_num):
                    F_ij = np.array(F_list[para_i][para_j]).reshape(p_shape)
                    arr = p * F_ij
                    for si in reversed(range(para_num)):
                        arr = simpson(arr, x[si])
                    res[para_i][para_j] = arr
                    res[para_j][para_i] = arr
            return res
        else:
            raise NameError("NameError: btype should be choosen in {1, 2, 3}.")

Calculation of the optimal biased bound based on the first type of the BQCRB in the case of single parameter estimation. The expression of OBB with a prior distribution \(p(x)\) is \begin{align} \mathrm{var}(\hat{x},{\Pi_y})\geq\int p(x)\left(\frac{(1+b')^2}{F} +b^2\right)\mathrm{d}x, \end{align}

where \(b\) and \(b'\) are the vector of biase and its derivative on \(x\). \(F\) is the QFI for all types.

Parameters

x: list -- The regimes of the parameters for the integral.

p: array -- The prior distribution.

dp: list -- Derivatives of the prior distribution with respect to the unknown parameters to to estimated. For example, dp[0] is the derivative vector with respect to the first parameter.

rho: list -- Parameterized density matrix.

drho: list -- Derivatives of the parameterized density matrix (rho) with respect to the unknown parameters to be estimated.

drho: list -- Second order Derivatives of the parameterized density matrix (rho) with respect to the unknown parameters to be estimated.

LDtype: string -- Types of QFI (QFIM) can be set as the objective function. Options are:
"SLD" (default) -- QFI (QFIM) based on symmetric logarithmic derivative (SLD).
"RLD" -- QFI (QFIM) based on right logarithmic derivative (RLD).
"LLD" -- QFI (QFIM) based on left logarithmic derivative (LLD).

eps: float -- Machine epsilon.

Returns

QVTB: float or matrix -- For single parameter estimation (the length of x is equal to one), the output is a float and for multiparameter estimation (the length of x is more than one), it returns a matrix.

Source code in quanestimation/BayesianBound/BayesCramerRao.py

def OBB(x, p, dp, rho, drho, d2rho, LDtype="SLD", eps=1e-8):
    r"""
    Calculation of the optimal biased bound based on the first type of the BQCRB 
    in the case of single parameter estimation. The expression of OBB with a 
    prior distribution $p(x)$ is
    \begin{align}
    \mathrm{var}(\hat{x},\{\Pi_y\})\geq\int p(x)\left(\frac{(1+b')^2}{F}
    +b^2\right)\mathrm{d}x,
    \end{align}

    where $b$ and $b'$ are the vector of biase and its derivative on $x$. 
    $F$ is the QFI for all types.

    Parameters
    ----------
    > **x:** `list`
        -- The regimes of the parameters for the integral.

    > **p:** `array`
        -- The prior distribution.

    > **dp:** `list`
        -- Derivatives of the prior distribution with respect to the unknown parameters to to
        estimated. For example, dp[0] is the derivative vector with respect to the first 
        parameter.

    > **rho:** `list`
        -- Parameterized density matrix.

    > **drho:** `list`
        -- Derivatives of the parameterized density matrix (rho) with respect to the unknown
        parameters to be estimated.

    > **drho:** `list`
        -- Second order Derivatives of the parameterized density matrix (rho) with respect to the 
        unknown parameters to be estimated.

    > **LDtype:** `string`
        -- Types of QFI (QFIM) can be set as the objective function. Options are:  
        "SLD" (default) -- QFI (QFIM) based on symmetric logarithmic derivative (SLD).  
        "RLD" -- QFI (QFIM) based on right logarithmic derivative (RLD).  
        "LLD" -- QFI (QFIM) based on left logarithmic derivative (LLD).

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

    Returns
    ----------
    **QVTB:** `float or matrix`
        -- For single parameter estimation (the length of x is equal to one), the 
        output is a float and for multiparameter estimation (the length of x is 
        more than one), it returns a matrix.
    """

    #### single parameter scenario ####
    p_num = len(p)

    if type(drho[0]) == list:
        drho = [drho[i][0] for i in range(p_num)]
    if type(d2rho[0]) == list:
        d2rho = [d2rho[i][0] for i in range(p_num)]
    if type(dp[0]) == list or type(dp[0]) == np.ndarray:
        dp = [dp[i][0] for i in range(p_num)]
    if type(x[0]) != float or type(x[0]) != int:
        x = x[0]

    F, J = np.zeros(p_num), np.zeros(p_num)
    bias, dbias = np.zeros(p_num), np.zeros(p_num)
    for m in range(p_num):
        f, LD = QFIM(rho[m], [drho[m]], LDtype=LDtype, exportLD=True, eps=eps)
        F[m] = f
        term1 = np.dot(d2rho[m], LD)
        term2 = np.dot(d2rho[m], LD.conj().T)
        term3 = np.dot(np.dot(LD, LD), drho[m])
        dF = np.real(np.trace(term1 + term2 - term3))
        J[m] = dp[m] / p[m] - dF / f

    y_guess = np.zeros((2, x.size))
    fun = lambda m, n: OBB_func(m, n, x, J, F)
    result = solve_bvp(fun, boundary_condition, x, y_guess)
    res = result.sol(x)
    bias, dbias = res[0], res[1]

    value = [p[i] * ((1 + dbias[i]) ** 2 / F[i] + bias[i] ** 2) for i in range(p_num)]
    return simpson(value, x)

Calculation of the Bayesian version of Cramer-Rao bound introduced by Van Trees (VTB). The covariance matrix with a prior distribution \(p(\textbf{x})\) is defined as \begin{align} \mathrm{cov}(\hat{\textbf{x}},{\Pi_y})=\int p(\textbf{x})\sum_y\mathrm{Tr} (\rho\Pi_y)(\hat{\textbf{x}}-\textbf{x})(\hat{\textbf{x}}-\textbf{x})^{\mathrm{T}} \mathrm{d}\textbf{x} \end{align}

where \(\textbf{x}=(x_0,x_1,\dots)^{\mathrm{T}}\) are the unknown parameters to be estimated and the integral \(\int\mathrm{d}\textbf{x}:=\iiint\mathrm{d}x_0\mathrm{d}x_1\cdots\). \(\{\Pi_y\}\) is a set of positive operator-valued measure (POVM) and \(\rho\) represent the parameterized density matrix.

\[\begin{align} \mathrm{cov}(\hat{\textbf{x}},\{\Pi_y\})\geq \left(\mathcal{I}_{\mathrm{prior}} +\mathcal{I}_{\mathrm{Bayes}}\right)^{-1}, \end{align}\]

where \(\mathcal{I}_{\mathrm{prior}}=\int p(\textbf{x})\mathcal{I}_{p}\mathrm{d}\textbf{x}\) is the CFIM for \(p(\textbf{x})\) and \(\mathcal{I}_{\mathrm{Bayes}}=\int p(\textbf{x})\mathcal{I}\mathrm{d}\textbf{x}\) is the average CFIM.

Parameters

x: list -- The regimes of the parameters for the integral.

p: multidimensional array -- The prior distribution.

dp: list -- Derivatives of the prior distribution with respect to the unknown parameters to be estimated. For example, dp[0] is the derivative vector with respect to the first parameter.

rho: multidimensional list -- Parameterized density matrix.

drho: multidimensional list -- Derivatives of the parameterized density matrix (rho) with respect to the unknown parameters to be estimated.

M: list of matrices -- A set of positive operator-valued measure (POVM). The default measurement is a set of rank-one symmetric informationally complete POVM (SIC-POVM).

eps: float -- Machine epsilon.

Returns

VTB: float or matrix -- For single parameter estimation (the length of x is equal to one), the output is a float and for multiparameter estimation (the length of x is more than one), it returns a matrix.

Note: SIC-POVM is calculated by the Weyl-Heisenberg covariant SIC-POVM fiducial state which can be downloaded from here.

Source code in quanestimation/BayesianBound/BayesCramerRao.py

def VTB(x, p, dp, rho, drho, M=[], eps=1e-8):
    r"""
    Calculation of the Bayesian version of Cramer-Rao bound introduced by
    Van Trees (VTB). The covariance matrix with a prior distribution $p(\textbf{x})$ 
    is defined as
    \begin{align}
    \mathrm{cov}(\hat{\textbf{x}},\{\Pi_y\})=\int p(\textbf{x})\sum_y\mathrm{Tr}
    (\rho\Pi_y)(\hat{\textbf{x}}-\textbf{x})(\hat{\textbf{x}}-\textbf{x})^{\mathrm{T}}
    \mathrm{d}\textbf{x}
    \end{align}

    where $\textbf{x}=(x_0,x_1,\dots)^{\mathrm{T}}$ are the unknown parameters to be estimated
    and the integral $\int\mathrm{d}\textbf{x}:=\iiint\mathrm{d}x_0\mathrm{d}x_1\cdots$.
    $\{\Pi_y\}$ is a set of positive operator-valued measure (POVM) and $\rho$ represent 
    the parameterized density matrix.

    \begin{align}
    \mathrm{cov}(\hat{\textbf{x}},\{\Pi_y\})\geq \left(\mathcal{I}_{\mathrm{prior}}
    +\mathcal{I}_{\mathrm{Bayes}}\right)^{-1},
    \end{align}

    where $\mathcal{I}_{\mathrm{prior}}=\int p(\textbf{x})\mathcal{I}_{p}\mathrm{d}\textbf{x}$ 
    is the CFIM for $p(\textbf{x})$ and 
    $\mathcal{I}_{\mathrm{Bayes}}=\int p(\textbf{x})\mathcal{I}\mathrm{d}\textbf{x}$ is the 
    average CFIM.

    Parameters
    ----------
    > **x:** `list`
        -- The regimes of the parameters for the integral.

    > **p:** `multidimensional array`
        -- The prior distribution.

    > **dp:** `list`
        -- Derivatives of the prior distribution with respect to the unknown parameters 
        to be estimated. For example, dp[0] is the derivative vector with respect to the first 
        parameter.

    > **rho:** `multidimensional list`
        -- Parameterized density matrix.

    > **drho:** `multidimensional list`
        -- Derivatives of the parameterized density matrix (rho) with respect to the 
        unknown parameters to be estimated.

    > **M:** `list of matrices`
        -- A set of positive operator-valued measure (POVM). The default measurement 
        is a set of rank-one symmetric informationally complete POVM (SIC-POVM).

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

    Returns
    ----------
    **VTB:** `float or matrix`
        -- For single parameter estimation (the length of x is equal to one), the 
        output is a float and for multiparameter estimation (the length of x is 
        more than one), it returns a matrix.

    **Note:** 
        SIC-POVM is calculated by the Weyl-Heisenberg covariant SIC-POVM fiducial state 
        which can be downloaded from [here](http://www.physics.umb.edu/Research/QBism/
        solutions.html).
    """

    para_num = len(x)
    p_num = len(p)

    if para_num == 1:
        #### single parameter scenario ####
        if M == []:
            M = SIC(len(rho[0]))
        else:
            if type(M) != list:
                raise TypeError("Please make sure M is a list!")

        if type(drho[0]) == list:
            drho = [drho[i][0] for i in range(p_num)]
        if type(dp[0]) == list or type(dp[0]) == np.ndarray:
            dp = [dp[i][0] for i in range(p_num)]

        F_tp = np.zeros(p_num)
        for m in range(p_num):
            F_tp[m] = CFIM(rho[m], [drho[m]], M=M, eps=eps)


        arr1 = [np.real(dp[i] * dp[i] / p[i]) for i in range(p_num)]
        I = simpson(arr1, x[0])
        arr2 = [np.real(F_tp[j] * p[j]) for j in range(p_num)]
        F = simpson(arr2, x[0])
        return 1.0 / (I + F)
    else:
        #### multiparameter scenario ####
        p_shape = np.shape(p)
        p_ext = extract_ele(p, para_num)
        dp_ext = extract_ele(dp, para_num)
        rho_ext = extract_ele(rho, para_num)
        drho_ext = extract_ele(drho, para_num)

        p_list, rho_list, drho_list = [], [], []
        for p_ele, rho_ele, drho_ele in zip(p_ext, rho_ext, drho_ext):
            p_list.append(p_ele)
            rho_list.append(rho_ele)
            drho_list.append(drho_ele)
        dp_list = [dpi for dpi in dp_ext]

        dim = len(rho_list[0])
        if M == []:
            M = SIC(dim)
        else:
            if type(M) != list:
                raise TypeError("Please make sure M is a list!")

        F_list = [
                [[0.0 for i in range(len(p_list))] for j in range(para_num)]
                for k in range(para_num)
            ]
        I_list = [
                [[0.0 for i in range(len(p_list))] for j in range(para_num)]
                for k in range(para_num)
            ]
        for i in range(len(p_list)):
            F_tp = CFIM(rho_list[i], drho_list[i], M=M, eps=eps)
            for pj in range(para_num):
                for pk in range(para_num):
                    F_list[pj][pk][i] = F_tp[pj][pk]
                    I_list[pj][pk][i] = (
                            dp_list[i][pj] * dp_list[i][pk] / p_list[i] ** 2
                        )

        F_res = np.zeros([para_num, para_num])
        I_res = np.zeros([para_num, para_num])
        for para_i in range(0, para_num):
            for para_j in range(para_i, para_num):
                F_ij = np.array(F_list[para_i][para_j]).reshape(p_shape)
                I_ij = np.array(I_list[para_i][para_j]).reshape(p_shape)
                arr1 = p * F_ij
                arr2 = p * I_ij
                for si in reversed(range(para_num)):
                    arr1 = simpson(arr1, x[si])
                    arr2 = simpson(arr2, x[si])
                F_res[para_i][para_j] = arr1
                F_res[para_j][para_i] = arr1
                I_res[para_i][para_j] = arr2
                I_res[para_j][para_i] = arr2
        return np.linalg.pinv(F_res + I_res)

Calculation of the Bayesian version of quantum Cramer-Rao bound introduced by Van Trees (QVTB). The covariance matrix with a prior distribution p(\textbf{x}) is defined as \begin{align} \mathrm{cov}(\hat{\textbf{x}},{\Pi_y})=\int p(\textbf{x})\sum_y\mathrm{Tr} (\rho\Pi_y)(\hat{\textbf{x}}-\textbf{x})(\hat{\textbf{x}}-\textbf{x})^{\mathrm{T}} \mathrm{d}\textbf{x} \end{align}

where \(\textbf{x}=(x_0,x_1,\dots)^{\mathrm{T}}\) are the unknown parameters to be estimated and the integral \(\int\mathrm{d}\textbf{x}:=\iiint\mathrm{d}x_0\mathrm{d}x_1\cdots\). \(\{\Pi_y\}\) is a set of positive operator-valued measure (POVM) and \(\rho\) represent the parameterized density matrix.

\[\begin{align} \mathrm{cov}(\hat{\textbf{x}},\{\Pi_y\})\geq \left(\mathcal{I}_{\mathrm{prior}} +\mathcal{F}_{\mathrm{Bayes}}\right)^{-1}, \end{align}\]

where \(\mathcal{I}_{\mathrm{prior}}=\int p(\textbf{x})\mathcal{I}_{p}\mathrm{d}\textbf{x}\) is the CFIM for \(p(\textbf{x})\) and \(\mathcal{F}_{\mathrm{Bayes}}=\int p(\textbf{x})\mathcal{F} \mathrm{d}\textbf{x}\) is the average QFIM of all types.

Parameters

x: list -- The regimes of the parameters for the integral.

p: multidimensional array -- The prior distribution.

dp: list -- Derivatives of the prior distribution with respect to the unknown parameters to to estimated. For example, dp[0] is the derivative vector with respect to the first parameter.

rho: multidimensional list -- Parameterized density matrix.

drho: multidimensional list -- Derivatives of the parameterized density matrix (rho) with respect to the unknown parameters to be estimated.

LDtype: string -- Types of QFI (QFIM) can be set as the objective function. Options are:
"SLD" (default) -- QFI (QFIM) based on symmetric logarithmic derivative (SLD).
"RLD" -- QFI (QFIM) based on right logarithmic derivative (RLD).
"LLD" -- QFI (QFIM) based on left logarithmic derivative (LLD).

eps: float -- Machine epsilon.

Returns

QVTB: float or matrix -- For single parameter estimation (the length of x is equal to one), the output is a float and for multiparameter estimation (the length of x is more than one), it returns a matrix.

Source code in quanestimation/BayesianBound/BayesCramerRao.py

def QVTB(x, p, dp, rho, drho, LDtype="SLD", eps=1e-8):
    r"""
    Calculation of the Bayesian version of quantum Cramer-Rao bound introduced 
    by Van Trees (QVTB). The covariance matrix with a prior distribution p(\textbf{x}) 
    is defined as
    \begin{align}
    \mathrm{cov}(\hat{\textbf{x}},\{\Pi_y\})=\int p(\textbf{x})\sum_y\mathrm{Tr}
    (\rho\Pi_y)(\hat{\textbf{x}}-\textbf{x})(\hat{\textbf{x}}-\textbf{x})^{\mathrm{T}}
    \mathrm{d}\textbf{x}
    \end{align}

    where $\textbf{x}=(x_0,x_1,\dots)^{\mathrm{T}}$ are the unknown parameters to be estimated
    and the integral $\int\mathrm{d}\textbf{x}:=\iiint\mathrm{d}x_0\mathrm{d}x_1\cdots$.
    $\{\Pi_y\}$ is a set of positive operator-valued measure (POVM) and $\rho$ represent 
    the parameterized density matrix.

    \begin{align}
    \mathrm{cov}(\hat{\textbf{x}},\{\Pi_y\})\geq \left(\mathcal{I}_{\mathrm{prior}}
    +\mathcal{F}_{\mathrm{Bayes}}\right)^{-1},
    \end{align}

    where $\mathcal{I}_{\mathrm{prior}}=\int p(\textbf{x})\mathcal{I}_{p}\mathrm{d}\textbf{x}$ is 
    the CFIM for $p(\textbf{x})$ and $\mathcal{F}_{\mathrm{Bayes}}=\int p(\textbf{x})\mathcal{F}
    \mathrm{d}\textbf{x}$ is the average QFIM of all types.

    Parameters
    ----------
    > **x:** `list`
        -- The regimes of the parameters for the integral.

    > **p:** multidimensional array
        -- The prior distribution.

    > **dp:** `list`
        -- Derivatives of the prior distribution with respect to the unknown parameters to to
        estimated. For example, dp[0] is the derivative vector with respect to the first 
        parameter.

    > **rho:** `multidimensional list`
        -- Parameterized density matrix.

    > **drho:** `multidimensional list`
        -- Derivatives of the parameterized density matrix (rho) with respect to the unknown
        parameters to be estimated.

    > **LDtype:** `string`
        -- Types of QFI (QFIM) can be set as the objective function. Options are:  
        "SLD" (default) -- QFI (QFIM) based on symmetric logarithmic derivative (SLD).  
        "RLD" -- QFI (QFIM) based on right logarithmic derivative (RLD).  
        "LLD" -- QFI (QFIM) based on left logarithmic derivative (LLD).

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

    Returns
    ----------
    **QVTB:** `float or matrix`
        -- For single parameter estimation (the length of x is equal to one), the 
        output is a float and for multiparameter estimation (the length of x is 
        more than one), it returns a matrix.
    """
    para_num = len(x)
    p_num = len(p)

    if para_num == 1:
        if type(drho[0]) == list:
            drho = [drho[i][0] for i in range(p_num)]
        if type(dp[0]) == list or type(dp[0]) == np.ndarray:
            dp = [dp[i][0] for i in range(p_num)]

        F_tp = np.zeros(p_num)
        for m in range(p_num):
            F_tp[m] = QFIM(rho[m], [drho[m]], LDtype=LDtype, eps=eps)

        arr1 = [np.real(dp[i] * dp[i] / p[i]) for i in range(p_num)]
        I = simpson(arr1, x[0])
        arr2 = [np.real(F_tp[j] * p[j]) for j in range(p_num)]
        F = simpson(arr2, x[0])
        return 1.0 / (I + F)
    else:
        #### multiparameter scenario ####
        p_shape = np.shape(p)
        p_ext = extract_ele(p, para_num)
        dp_ext = extract_ele(dp, para_num)
        rho_ext = extract_ele(rho, para_num)
        drho_ext = extract_ele(drho, para_num)

        p_list, rho_list, drho_list = [], [], []
        for p_ele, rho_ele, drho_ele in zip(p_ext, rho_ext, drho_ext):
            p_list.append(p_ele)
            rho_list.append(rho_ele)
            drho_list.append(drho_ele)
        dp_list = [dpi for dpi in dp_ext]

        F_list = [
                [[0.0 for i in range(len(p_list))] for j in range(para_num)]
                for k in range(para_num)
            ]
        I_list = [
                [[0.0 for i in range(len(p_list))] for j in range(para_num)]
                for k in range(para_num)
            ]
        for i in range(len(p_list)):
            F_tp = QFIM(rho_list[i], drho_list[i], LDtype=LDtype, eps=eps)
            for pj in range(para_num):
                for pk in range(para_num):
                    F_list[pj][pk][i] = F_tp[pj][pk]
                    I_list[pj][pk][i] = (
                            dp_list[i][pj] * dp_list[i][pk] / p_list[i] ** 2
                        )

        F_res = np.zeros([para_num, para_num])
        I_res = np.zeros([para_num, para_num])
        for para_i in range(0, para_num):
            for para_j in range(para_i, para_num):
                F_ij = np.array(F_list[para_i][para_j]).reshape(p_shape)
                I_ij = np.array(I_list[para_i][para_j]).reshape(p_shape)
                arr1 = p * F_ij
                arr2 = p * I_ij
                for si in reversed(range(para_num)):
                    arr1 = simpson(arr1, x[si])
                    arr2 = simpson(arr2, x[si])
                F_res[para_i][para_j] = arr1
                F_res[para_j][para_i] = arr1
                I_res[para_i][para_j] = arr2
                I_res[para_j][para_i] = arr2
        return np.linalg.pinv(F_res + I_res)

---

Quantum Ziv-Zakai bound

Calculation of the quantum Ziv-Zakai bound (QZZB). The expression of QZZB with a prior distribution p(x) in a finite regime \([\alpha,\beta]\) is

\[\begin{eqnarray} \mathrm{var}(\hat{x},\{\Pi_y\}) &\geq & \frac{1}{2}\int_0^\infty \mathrm{d}\tau\tau \mathcal{V}\int_{-\infty}^{\infty} \mathrm{d}x\min\!\left\{p(x), p(x+\tau)\right\} \nonumber \\ & & \times\left(1-\frac{1}{2}||\rho(x)-\rho(x+\tau)||\right), \end{eqnarray}\]

where \(||\cdot||\) represents the trace norm and \(\mathcal{V}\) is the "valley-filling" operator satisfying \(\mathcal{V}f(\tau)=\max_{h\geq 0}f(\tau+h)\). \(\rho(x)\) is the parameterized density matrix.

Parameters

x: list -- The regimes of the parameters for the integral.

p: multidimensional array -- The prior distribution.

rho: multidimensional list -- Parameterized density matrix.

eps: float -- Machine epsilon.

Returns

QZZB: float -- Quantum Ziv-Zakai bound (QZZB).

Source code in quanestimation/BayesianBound/ZivZakai.py

def QZZB(x, p, rho, eps=1e-8):
    r"""
    Calculation of the quantum Ziv-Zakai bound (QZZB). The expression of QZZB with a 
    prior distribution p(x) in a finite regime $[\alpha,\beta]$ is

    \begin{eqnarray}
    \mathrm{var}(\hat{x},\{\Pi_y\}) &\geq & \frac{1}{2}\int_0^\infty \mathrm{d}\tau\tau
    \mathcal{V}\int_{-\infty}^{\infty} \mathrm{d}x\min\!\left\{p(x), p(x+\tau)\right\} \nonumber \\
    & & \times\left(1-\frac{1}{2}||\rho(x)-\rho(x+\tau)||\right),
    \end{eqnarray}

    where $||\cdot||$ represents the trace norm and $\mathcal{V}$ is the "valley-filling" 
    operator satisfying $\mathcal{V}f(\tau)=\max_{h\geq 0}f(\tau+h)$. $\rho(x)$ is the 
    parameterized density matrix. 

    Parameters
    ----------
    > **x:** `list`
        -- The regimes of the parameters for the integral.

    > **p:** `multidimensional array`
        -- The prior distribution.

    > **rho:** `multidimensional list`
        -- Parameterized density matrix.

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

    Returns
    ----------
    **QZZB:** `float`
        -- Quantum Ziv-Zakai bound (QZZB).
    """

    if type(x[0]) == list or type(x[0]) == np.ndarray:
        x = x[0]
    p_num = len(p)
    tau = [xi - x[0] for xi in x]
    f_tau = np.zeros(p_num)
    for i in range(p_num):
        arr = [
            np.real(2 * min(p[j], p[j + i]) * helstrom_dm(rho[j], rho[j + i], eps))
            for j in range(p_num - i)
        ]
        f_tp = simpson(arr, x[0 : p_num - i])
        f_tau[i] = f_tp
    arr2 = [tau[m] * max(f_tau[m:]) for m in range(p_num)]
    I = simpson(arr2, tau)
    return 0.5 * I

---

Bayesian estimation

Bayesian estimation. The prior distribution is updated via the posterior
distribution obtained by the Bayes’ rule and the estimated value of parameters are updated via the expectation value of the distribution or maximum a posteriori probability (MAP).

Parameters

x: list -- The regimes of the parameters for the integral.

p: multidimensional array -- The prior distribution.

rho: multidimensional list -- Parameterized density matrix.

y: array -- The experimental results obtained in practice.

M: list of matrices -- A set of positive operator-valued measure (POVM). The default measurement is a set of rank-one symmetric informationally complete POVM (SIC-POVM).

estimator: string -- Estimators for the bayesian estimation. Options are:
"mean" -- The expectation value of the distribution.
"MAP" -- Maximum a posteriori probability.

savefile: bool -- Whether or not to save all the posterior distributions.
If set True then two files "pout.npy" and "xout.npy" will be generated including the posterior distributions and the estimated values in the iterations. If set False the posterior distribution in the final iteration and the estimated values in all iterations will be saved in "pout.npy" and "xout.npy".

Returns

pout and xout: array and float -- The posterior distribution and the estimated values in the final iteration.

Note: SIC-POVM is calculated by the Weyl-Heisenberg covariant SIC-POVM fiducial state which can be downloaded from here.

Source code in quanestimation/BayesianBound/BayesEstimation.py

def Bayes(x, p, rho, y, M=[], estimator="mean", savefile=False):
    """
    Bayesian estimation. The prior distribution is updated via the posterior  
    distribution obtained by the Bayes’ rule and the estimated value of parameters
    are updated via the expectation value of the distribution or maximum a 
    posteriori probability (MAP).

    Parameters
    ----------
    > **x:** `list`
        -- The regimes of the parameters for the integral.

    > **p:** `multidimensional array`
        -- The prior distribution.

    > **rho:** `multidimensional list`
        -- Parameterized density matrix.

    > **y:** `array`
        -- The experimental results obtained in practice.

    > **M:** `list of matrices`
        -- A set of positive operator-valued measure (POVM). The default measurement 
        is a set of rank-one symmetric informationally complete POVM (SIC-POVM).

    > **estimator:** `string`
        -- Estimators for the bayesian estimation. Options are:  
        "mean" -- The expectation value of the distribution.  
        "MAP" -- Maximum a posteriori probability.

    > **savefile:** `bool`
        -- Whether or not to save all the posterior distributions.  
        If set `True` then two files "pout.npy" and "xout.npy" will be generated including
        the posterior distributions and the estimated values in the iterations. If set 
        `False` the posterior distribution in the final iteration and the estimated values
        in all iterations will be saved in "pout.npy" and "xout.npy". 

    Returns
    ----------
    **pout and xout:** `array and float`
        -- The posterior distribution and the estimated values in the final iteration.

    **Note:** 
        SIC-POVM is calculated by the Weyl-Heisenberg covariant SIC-POVM fiducial state 
        which can be downloaded from [here](http://www.physics.umb.edu/Research/QBism/
        solutions.html).
    """

    para_num = len(x)
    max_episode = len(y)
    if para_num == 1:
        #### single parameter scenario ####
        if M == []:
            M = SIC(len(rho[0]))
        else:
            if type(M) != list:
                raise TypeError("Please make sure M is a list!")
        if savefile == False:
            x_out = []
            if estimator == "mean":
                for mi in range(max_episode):
                    res_exp = int(y[mi])
                    pyx = np.zeros(len(x[0]))
                    for xi in range(len(x[0])):
                        p_tp = np.real(np.trace(np.dot(rho[xi], M[res_exp])))
                        pyx[xi] = p_tp
                    arr = [pyx[m] * p[m] for m in range(len(x[0]))]
                    py = simpson(arr, x[0])
                    p_update = pyx * p / py
                    p = p_update
                    mean = simpson([p[m]*x[0][m] for m in range(len(x[0]))], x[0])
                    x_out.append(mean)
            elif estimator == "MAP":
                for mi in range(max_episode):
                    res_exp = int(y[mi])
                    pyx = np.zeros(len(x[0]))
                    for xi in range(len(x[0])):
                        p_tp = np.real(np.trace(np.dot(rho[xi], M[res_exp])))
                        pyx[xi] = p_tp
                    arr = [pyx[m] * p[m] for m in range(len(x[0]))]
                    py = simpson(arr, x[0])
                    p_update = pyx * p / py
                    p = p_update
                    indx = np.where(p == max(p))[0][0]
                    x_out.append(x[0][indx])
            else:
                raise ValueError(
                "{!r} is not a valid value for estimator, supported values are 'mean' and 'MAP'.".format(estimator))
            np.save("pout", p)
            np.save("xout", x_out)
            return p, x_out[-1]
        else:
            p_out, x_out = [], []
            if estimator == "mean":
                for mi in range(max_episode):
                    res_exp = int(y[mi])
                    pyx = np.zeros(len(x[0]))
                    for xi in range(len(x[0])):
                        p_tp = np.real(np.trace(np.dot(rho[xi], M[res_exp])))
                        pyx[xi] = p_tp
                    arr = [pyx[m] * p[m] for m in range(len(x[0]))]
                    py = simpson(arr, x[0])
                    p_update = pyx * p / py
                    p = p_update
                    mean = simpson([p[m]*x[0][m] for m in range(len(x[0]))], x[0])
                    p_out.append(p)
                    x_out.append(mean)
            elif estimator == "MAP":
                for mi in range(max_episode):
                    res_exp = int(y[mi])
                    pyx = np.zeros(len(x[0]))
                    for xi in range(len(x[0])):
                        p_tp = np.real(np.trace(np.dot(rho[xi], M[res_exp])))
                        pyx[xi] = p_tp
                    arr = [pyx[m] * p[m] for m in range(len(x[0]))]
                    py = simpson(arr, x[0])
                    p_update = pyx * p / py
                    p = p_update
                    indx = np.where(p == max(p))[0][0]
                    p_out.append(p)
                    x_out.append(x[0][indx])
            else:
                raise ValueError(
                "{!r} is not a valid value for estimator, supported values are 'mean' and 'MAP'.".format(estimator))
            np.save("pout", p_out)
            np.save("xout", x_out)
            return p, x_out[-1]
    else:
        #### multiparameter scenario ####
        p_shape = np.shape(p)
        p_ext = extract_ele(p, para_num)
        rho_ext = extract_ele(rho, para_num)

        p_list, rho_list = [], []
        for p_ele, rho_ele in zip(p_ext, rho_ext):
            p_list.append(p_ele)
            rho_list.append(rho_ele)

        dim = len(rho_list[0])
        if M == []:
            M = SIC(dim)
        else:
            if type(M) != list:
                raise TypeError("Please make sure M is a list!")

        if savefile == False:
            x_out = []
            if estimator == "mean":
                for mi in range(max_episode):
                    res_exp = int(y[mi])
                    pyx_list = np.zeros(len(p_list))
                    for xi in range(len(p_list)):
                        p_tp = np.real(np.trace(np.dot(rho_list[xi], M[res_exp])))
                        pyx_list[xi] = p_tp
                    pyx = pyx_list.reshape(p_shape)
                    arr = p * pyx
                    for si in reversed(range(para_num)):
                        arr = simpson(arr, x[si])
                    py = arr
                    p_update = p * pyx / py
                    p = p_update

                    mean = integ(x, p)
                    x_out.append(mean)
            elif estimator == "MAP":
                for mi in range(max_episode):
                    res_exp = int(y[mi])
                    pyx_list = np.zeros(len(p_list))
                    for xi in range(len(p_list)):
                        p_tp = np.real(np.trace(np.dot(rho_list[xi], M[res_exp])))
                        pyx_list[xi] = p_tp
                    pyx = pyx_list.reshape(p_shape)
                    arr = p * pyx
                    for si in reversed(range(para_num)):
                        arr = simpson(arr, x[si])
                    py = arr
                    p_update = p * pyx / py
                    p = p_update

                    indx = np.where(np.array(p) == np.max(np.array(p)))
                    x_out.append([x[i][indx[i][0]] for i in range(para_num)])
            else:
                raise ValueError(
                "{!r} is not a valid value for estimator, supported values are 'mean' and 'MAP'.".format(estimator))
            np.save("pout", p)
            np.save("xout", x_out)
            return p, x_out[-1]
        else:
            p_out, x_out = [], []
            if estimator == "mean":
                for mi in range(max_episode):
                    res_exp = int(y[mi])
                    pyx_list = np.zeros(len(p_list))
                    for xi in range(len(p_list)):
                        p_tp = np.real(np.trace(np.dot(rho_list[xi], M[res_exp])))
                        pyx_list[xi] = p_tp
                    pyx = pyx_list.reshape(p_shape)
                    arr = p * pyx
                    for si in reversed(range(para_num)):
                        arr = simpson(arr, x[si])
                    py = arr
                    p_update = p * pyx / py
                    p = p_update

                    mean = integ(x, p)
                    p_out.append(p)
                    x_out.append(mean)
            elif estimator == "MAP":
                for mi in range(max_episode):
                    res_exp = int(y[mi])
                    pyx_list = np.zeros(len(p_list))
                    for xi in range(len(p_list)):
                        p_tp = np.real(np.trace(np.dot(rho_list[xi], M[res_exp])))
                        pyx_list[xi] = p_tp
                    pyx = pyx_list.reshape(p_shape)
                    arr = p * pyx
                    for si in reversed(range(para_num)):
                        arr = simpson(arr, x[si])
                    py = arr
                    p_update = p * pyx / py
                    p = p_update

                    indx = np.where(np.array(p) == np.max(np.array(p)))
                    p_out.append(p)
                    x_out.append([x[i][indx[i][0]] for i in range(para_num)])
            else:
                raise ValueError(
                "{!r} is not a valid value for estimator, supported values are 'mean' and 'MAP'.".format(estimator))
            np.save("pout", p_out)
            np.save("xout", x_out)
            return p, x_out[-1]

Bayesian estimation. The estimated value of parameters obtained via the maximum likelihood estimation (MLE).

Parameters

x: list -- The regimes of the parameters for the integral.

rho: multidimensional list -- Parameterized density matrix.

y: array -- The experimental results obtained in practice.

M: list of matrices -- A set of positive operator-valued measure (POVM). The default measurement is a set of rank-one symmetric informationally complete POVM (SIC-POVM).

savefile: bool -- Whether or not to save all the likelihood functions.
If set True then two files "Lout.npy" and "xout.npy" will be generated including the likelihood functions and the estimated values in the iterations. If set False the likelihood function in the final iteration and the estimated values in all iterations will be saved in "Lout.npy" and "xout.npy".

Returns

Lout and xout: array and float -- The likelihood function and the estimated values in the final iteration.

Note: SIC-POVM is calculated by the Weyl-Heisenberg covariant SIC-POVM fiducial state which can be downloaded from here.

Source code in quanestimation/BayesianBound/BayesEstimation.py

def MLE(x, rho, y, M=[], savefile=False):
    """
    Bayesian estimation. The estimated value of parameters obtained via the 
    maximum likelihood estimation (MLE).

    Parameters
    ----------
    > **x:** `list`
        -- The regimes of the parameters for the integral.

    > **rho:** `multidimensional list`
        -- Parameterized density matrix.

    > **y:** `array`
        -- The experimental results obtained in practice.

    > **M:** `list of matrices`
        -- A set of positive operator-valued measure (POVM). The default measurement 
        is a set of rank-one symmetric informationally complete POVM (SIC-POVM).

    > **savefile:** `bool`
        -- Whether or not to save all the likelihood functions.  
        If set `True` then two files "Lout.npy" and "xout.npy" will be generated including
        the likelihood functions and the estimated values in the iterations. If set 
        `False` the likelihood function in the final iteration and the estimated values
        in all iterations will be saved in "Lout.npy" and "xout.npy". 

    Returns
    ----------
    **Lout and xout:** `array and float`
        -- The likelihood function and the estimated values in the final iteration.

    **Note:** 
        SIC-POVM is calculated by the Weyl-Heisenberg covariant SIC-POVM fiducial state 
        which can be downloaded from [here](http://www.physics.umb.edu/Research/QBism/
        solutions.html).
    """

    para_num = len(x)
    max_episode = len(y)
    if para_num == 1:
        #### single parameter scenario ####
        if M == []:
            M = SIC(len(rho[0]))
        else:
            if type(M) != list:
                raise TypeError("Please make sure M is a list!")

        if savefile == False:
            x_out = []
            L_out = np.ones(len(x[0]))
            for mi in range(max_episode):
                res_exp = int(y[mi])
                for xi in range(len(x[0])):
                    p_tp = np.real(np.trace(np.dot(rho[xi], M[res_exp])))
                    L_out[xi] = L_out[xi] * p_tp
                indx = np.where(L_out == max(L_out))[0][0]
                x_out.append(x[0][indx])
            np.save("Lout", L_out)
            np.save("xout", x_out)

            return L_out, x_out[-1]
        else:
            L_out, x_out = [], []
            L_tp = np.ones(len(x[0]))
            for mi in range(max_episode):
                res_exp = int(y[mi])
                for xi in range(len(x[0])):
                    p_tp = np.real(np.trace(np.dot(rho[xi], M[res_exp])))
                    L_tp[xi] = L_tp[xi] * p_tp
                indx = np.where(L_tp == max(L_tp))[0][0]
                L_out.append(L_tp)
                x_out.append(x[0][indx])

            np.save("Lout", L_out)
            np.save("xout", x_out)
            return L_tp, x_out[-1]
    else:
        #### multiparameter scenario ####
        p_shape = []
        for i in range(para_num):
            p_shape.append(len(x[i]))
        rho_ext = extract_ele(rho, para_num)

        rho_list = []
        for rho_ele in rho_ext:
            rho_list.append(rho_ele)

        dim = len(rho_list[0])
        if M == []:
            M = SIC(dim)
        else:
            if type(M) != list:
                raise TypeError("Please make sure M is a list!")

        if savefile == False:
            x_out = []
            L_list = np.ones(len(rho_list))
            for mi in range(max_episode):
                res_exp = int(y[mi])
                for xi in range(len(rho_list)):
                    p_tp = np.real(np.trace(np.dot(rho_list[xi], M[res_exp])))
                    L_list[xi] = L_list[xi] * p_tp
                L_out = L_list.reshape(p_shape)
                indx = np.where(L_out == np.max(L_out))
                x_out.append([x[i][indx[i][0]] for i in range(para_num)])
            np.save("Lout", L_out)
            np.save("xout", x_out)

            return L_out, x_out[-1]
        else:
            L_out, x_out = [], []
            L_list = np.ones(len(rho_list))
            for mi in range(max_episode):
                res_exp = int(y[mi])
                for xi in range(len(rho_list)):
                    p_tp = np.real(np.trace(np.dot(rho_list[xi], M[res_exp])))
                    L_list[xi] = L_list[xi] * p_tp
                L_tp = L_list.reshape(p_shape)
                indx = np.where(L_tp == np.max(L_tp))
                L_out.append(L_tp)
                x_out.append([x[i][indx[i][0]] for i in range(para_num)])

            np.save("Lout", L_out)
            np.save("xout", x_out)
            return L_tp, x_out[-1]

Calculation of the average Bayesian cost with a quadratic cost function.

Parameters

x: list -- The regimes of the parameters for the integral.

p: multidimensional array -- The prior distribution.

xest: list -- The estimators.

rho: multidimensional list -- Parameterized density matrix.

M: array -- A set of POVM.

W: array -- Weight matrix.

eps: float -- Machine epsilon.

Returns

The average Bayesian cost: float -- The average Bayesian cost.

Source code in quanestimation/BayesianBound/BayesEstimation.py

def BayesCost(x, p, xest, rho, M, W=[], eps=1e-8):
    """
    Calculation of the average Bayesian cost with a quadratic cost function.

    Parameters
    ----------
    > **x:** `list`
        -- The regimes of the parameters for the integral.

    > **p:** `multidimensional array`
        -- The prior distribution.

    > **xest:** `list`
        -- The estimators.

    > **rho:** `multidimensional list`
        -- Parameterized density matrix.

    > **M:** `array`
        -- A set of POVM.

    > **W:** `array`
        -- Weight matrix.

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

    Returns
    ----------
    **The average Bayesian cost:** `float`
        -- The average Bayesian cost.
    """
    para_num = len(x)
    if para_num == 1:
        # single-parameter scenario
        if M == []:
            M = SIC(len(rho[0]))
        else:
            if type(M) != list:
                raise TypeError("Please make sure M is a list!")
        p_num = len(x[0])
        value = [p[i]*sum([np.trace(np.dot(rho[i], M[mi]))*(x[0][i]-xest[mi][0])**2 for mi in range(len(M))]) for i in range(p_num)]
        C = simpson(value, x[0])
        return np.real(C)
    else:
        # multi-parameter scenario
        p_shape = np.shape(p)
        p_ext = extract_ele(p, para_num)
        rho_ext = extract_ele(rho, para_num)

        p_list, rho_list = [], []
        for p_ele, rho_ele in zip(p_ext, rho_ext):
            p_list.append(p_ele)
            rho_list.append(rho_ele)

        x_pro = product(*x)
        x_list = []
        for x_ele in x_pro:
            x_list.append([x_ele[i] for i in range(para_num)])

        dim = len(rho_list[0])
        p_num = len(p_list)

        if W == []:
            W = np.identity(para_num)

        if M == []:
            M = SIC(dim)
        else:
            if type(M) != list:
                raise TypeError("Please make sure M is a list!")

        value = [0.0 for i in range(p_num)]
        for i in range(p_num):
            x_tp = np.array(x_list[i])
            xCx = 0.0
            for mi in range(len(M)):
                xCx += np.trace(np.dot(rho_list[i], M[mi]))*np.dot((x_tp-xest[mi]).reshape(1, -1), np.dot(W, (x_tp-xest[mi]).reshape(-1, 1)))[0][0]
            value[i] = p_list[i]*xCx
        C = np.array(value).reshape(p_shape)
        for si in reversed(range(para_num)):
            C = simpson(C, x[si])
        return np.real(C)

Calculation of the Bayesian cost bound with a quadratic cost function.

Parameters

x: list -- The regimes of the parameters for the integral.

p: multidimensional array -- The prior distribution.

rho: multidimensional list -- Parameterized density matrix.

W: array -- Weight matrix.

eps: float -- Machine epsilon.

Returns

BCB: float -- The value of the minimum Bayesian cost.

Source code in quanestimation/BayesianBound/BayesEstimation.py

def BCB(x, p, rho, W=[], eps=1e-8):
    """
    Calculation of the Bayesian cost bound with a quadratic cost function.

    Parameters
    ----------
    > **x:** `list`
        -- The regimes of the parameters for the integral.

    > **p:** `multidimensional array`
        -- The prior distribution.

    > **rho:** `multidimensional list`
        -- Parameterized density matrix.

    > **W:** `array`
        -- Weight matrix.

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

    Returns
    ----------
    **BCB:** `float`
        -- The value of the minimum Bayesian cost.
    """
    para_num = len(x)
    if para_num == 1:
        # single-parameter scenario
        dim = len(rho[0])
        p_num = len(x[0])
        value = [p[i]*x[0][i]**2 for i in range(p_num)]
        delta2_x = simpson(value, x[0])
        rho_avg = np.zeros((dim, dim), dtype=np.complex128)
        rho_pri = np.zeros((dim, dim), dtype=np.complex128)
        for di in range(dim):
            for dj in range(dim):
                rho_avg_arr = [p[m]*rho[m][di][dj] for m in range(p_num)]
                rho_pri_arr = [p[n]*x[0][n]*rho[n][di][dj] for n in range(p_num)]
                rho_avg[di][dj] = simpson(rho_avg_arr, x[0])
                rho_pri[di][dj] = simpson(rho_pri_arr, x[0])
        Lambda = Lambda_avg(rho_avg, [rho_pri], eps=eps)
        minBC = delta2_x-np.real(np.trace(np.dot(np.dot(rho_avg, Lambda[0]), Lambda[0])))
        return minBC
    else:
        # multi-parameter scenario
        p_shape = np.shape(p)
        p_ext = extract_ele(p, para_num)
        rho_ext = extract_ele(rho, para_num)

        p_list, rho_list = [], []
        for p_ele, rho_ele in zip(p_ext, rho_ext):
            p_list.append(p_ele)
            rho_list.append(rho_ele)

        dim = len(rho_list[0])
        p_num = len(p_list)

        x_pro = product(*x)
        x_list = []
        for x_ele in x_pro:
            x_list.append([x_ele[i] for i in range(para_num)])

        if W == []:
            W = np.identity(para_num)

        value = [0.0 for i in range(p_num)]
        for i in range(p_num):
            x_tp = np.array(x_list[i])
            xCx = np.dot(x_tp.reshape(1, -1), np.dot(W, x_tp.reshape(-1, 1)))[0][0]
            value[i] = p_list[i]*xCx
        delta2_x = np.array(value).reshape(p_shape)
        for si in reversed(range(para_num)):
            delta2_x = simpson(delta2_x, x[si])
        rho_avg = np.zeros((dim, dim), dtype=np.complex128)
        rho_pri = [np.zeros((dim, dim), dtype=np.complex128) for i in range(para_num)]
        for di in range(dim):
            for dj in range(dim):
                rho_avg_arr = [p_list[m]*rho_list[m][di][dj] for m in range(p_num)]
                rho_avg_tp = np.array(rho_avg_arr).reshape(p_shape)
                for si in reversed(range(para_num)):
                    rho_avg_tp = simpson(rho_avg_tp, x[si])
                rho_avg[di][dj] = rho_avg_tp

                for para_i in range(para_num):
                    rho_pri_arr = [p_list[n]*x_list[n][para_i]*rho_list[n][di][dj] for n in range(p_num)]
                    rho_pri_tp = np.array(rho_pri_arr).reshape(p_shape)
                    for si in reversed(range(para_num)):
                        rho_pri_tp = simpson(rho_pri_tp, x[si])

                    rho_pri[para_i][di][dj] = rho_pri_tp
        Lambda = Lambda_avg(rho_avg, rho_pri, eps=eps)
        Mat = np.zeros((para_num, para_num), dtype=np.complex128)
        for para_m in range(para_num):
            for para_n in range(para_num):
                Mat += W[para_m][para_n]*np.dot(Lambda[para_m], Lambda[para_n])

        minBC = delta2_x-np.real(np.trace(np.dot(rho_avg, Mat)))
        return minBC

---

Common

Generation of the input variables H, dH (or K, dK).

Parameters

x: list -- The regimes of the parameters for the integral.

func: list -- Function defined by the users which returns H or K.

dfunc: list -- Function defined by the users which returns dH or dK.

channel: string -- Seeting the output of this function. Options are:
"dynamics" (default) -- The output of this function is H and dH.
"Kraus" (default) -- The output of this function is K and dHK.

Returns

H, dH (or K, dK).

Source code in quanestimation/Common/Common.py

def BayesInput(x, func, dfunc, channel="dynamics"):
    """
    Generation of the input variables H, dH (or K, dK).

    Parameters
    ----------
    > **x:** `list`
        -- The regimes of the parameters for the integral.

    > **func:** `list`
        -- Function defined by the users which returns H or K.

    > **dfunc:** `list`
        -- Function defined by the users which returns dH or dK.

    > **channel:** `string`
        -- Seeting the output of this function. Options are:  
        "dynamics" (default) --  The output of this function is H and dH.  
        "Kraus" (default) --  The output of this function is K and dHK.

    Returns
    ----------
    H, dH (or K, dK).
    """

    para_num = len(x)
    size = [len(x[i]) for i in range(len(x))]
    x_all = product(*x)
    if channel == "dynamics":
        dim = len(func([0 for i in range(para_num)]))
        H_list, dH_list = [], []
        for xi in x_all:
            H_list.append(func([i for i in xi]))
            dH_list.append(dfunc([i for i in xi]))
        # H_res = np.reshape(H_list, [*size, *[dim, dim]])
        # dH_res = np.reshape(dH_list, [*size, *[para_num, dim, dim]])
        return H_list, dH_list
    elif channel == "Kraus":
        k_num = len(func([0 for i in range(para_num)]))
        dim = len(func([0 for i in range(para_num)])[0])
        K_list, dK_list = [], []
        if para_num == 1:
            for xi in x_all:
                K_list.append(func([i for i in xi]))
            #     dK_list.append(dfunc([i for i in xi]))
            # K_res = np.reshape(K_list, [*size, *[k_num, dim, dim]])
            # dK_res = np.reshape(dK_list, [*size, *[para_num, k_num, dim, dim]])
        else:
            for xi in x_all:
                K_list.append(func([i for i in xi]))
                dK_list.append(dfunc([i for i in xi]))
            # K_res = np.reshape(K_list, [*size, *[k_num, dim, dim]])
            # dK_res = np.reshape(dK_list, [*size, *[k_num, para_num, dim, dim]])
        return K_list, dK_list
    else:
        raise ValueError(
            "{!r} is not a valid value for channel, supported values are 'dynamics' and 'Kraus'.".format(
                channel
            )
        )

Generation of a set of rank-one symmetric informationally complete positive operator-valued measure (SIC-POVM).

Parameters

dim: int -- The dimension of the system.

Returns

A set of SCI-POVM.

Note: SIC-POVM is calculated by the Weyl-Heisenberg covariant SIC-POVM fiducial state which can be downloaded from here.

Source code in quanestimation/Common/Common.py

def SIC(dim):
    """
    Generation of a set of rank-one symmetric informationally complete 
    positive operator-valued measure (SIC-POVM).

    Parameters
    ----------
    > **dim:** `int` 
        -- The dimension of the system.

    Returns
    ----------
    A set of SCI-POVM.

    **Note:** 
        SIC-POVM is calculated by the Weyl-Heisenberg covariant SIC-POVM fiducial state 
        which can be downloaded from [here](http://www.physics.umb.edu/Research/QBism/
        solutions.html).
    """

    if dim <= 151:
        file_path = os.path.join(
            os.path.dirname(os.path.dirname(__file__)),
            "sic_fiducial_vectors/d%d.txt" % (dim),
        )
        data = np.loadtxt(file_path)
        fiducial = data[:, 0] + data[:, 1] * 1.0j
        fiducial = np.array(fiducial).reshape(len(fiducial), 1)
        M = sic_povm(fiducial)
        return M
    else:
        raise ValueError("The dimension of the space should be less or equal to 151.")

Generation of the SU(\(N\)) generators with \(N\) the dimension of the system.

Parameters

n: int -- The dimension of the system.

Returns

SU(\(N\)) generators.

Source code in quanestimation/Common/Common.py

def suN_generator(n):
    r"""
    Generation of the SU($N$) generators with $N$ the dimension of the system.

    Parameters
    ----------
    > **n:** `int` 
        -- The dimension of the system.

    Returns
    ----------
    SU($N$) generators.
    """

    symm, anti_symm, diag = suN_unsorted(n)
    if n == 2:
        return [symm[0], anti_symm[0], diag[0]]
    else:
        Lambda = [0.0 for i in range(len(symm + anti_symm + diag))]

        Lambda[0], Lambda[1], Lambda[2] = symm[0], anti_symm[0], diag[0]

        repeat_times = 2
        m1, n1, k1 = 0, 3, 1
        while True:
            m1 += n1
            j, l = 0, 0
            for i in range(repeat_times):
                Lambda[m1 + j] = symm[k1]
                Lambda[m1 + j + 1] = anti_symm[k1]
                j += 2
                k1 += 1

            repeat_times += 1
            n1 = n1 + 2
            if k1 == len(symm):
                break

        m2, n2, k2 = 2, 5, 1
        while True:
            m2 += n2
            Lambda[m2] = diag[k2]
            n2 = n2 + 2
            k2 = k2 + 1
            if k2 == len(diag):
                break
        return Lambda