Skip to content

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 
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
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 
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
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 
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
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 
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
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(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] = vec.conj().transpose() @ 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 * (vec[:, fi].conj().transpose() @ drho[para_i] @ vec[:, fj]) / (val[fi] + val[fj])

            SLD_eig[SLD_eig == np.inf] = 0.0

            if rep == "original":
                SLD[para_i] = vec @ 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 
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
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 = vec[:, fi].conj().transpose() @ 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] = vec @ 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 
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
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 = vec[:, fi].conj().transpose() @ drho[para_i] @ vec[:, fj]
                if np.abs(val[fj]) > eps:
                    LLD_eig[fi][fj] = term_tp/val[fj]
                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] = vec @ 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 
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
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 = LD_tp @ LD_tp + LD_tp @ LD_tp
            QFIM_res = np.real(0.5 * np.trace(rho @ SLD_ac))
        elif LDtype == "RLD":
            LD_tp = RLD(rho, drho, eps=eps)
            QFIM_res = np.real(np.trace(rho @ LD_tp @ LD_tp.conj().transpose()))
        elif LDtype == "LLD":
            LD_tp = LLD(rho, drho, eps=eps)
            QFIM_res = np.real(np.trace(rho @ 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 = LD_tp[para_i] @ LD_tp[para_j] + LD_tp[para_j] @ LD_tp[para_i]
                    QFIM_res[para_i][para_j] = np.real(0.5 * np.trace(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(rho @ 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(rho @ 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 
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
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([Ki @ rho0 @ Ki.conj().T for Ki in K])
    drho = [
                sum(
                    [
                        (dKi @ rho0 @ Ki.conj().T+ Ki @ 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 
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
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(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(drho_i @ Mp))
                for para_j in range(para_i, para_num):
                    drho_j = drho[para_j]
                    dp_j = np.real(np.trace(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 
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
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 
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
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: np.array -- 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 
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
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:** `np.array`
        -- 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 = Lambda[row_i] @ Lambda[col_j] + Lambda[col_j] @ Lambda[row_i]
                G[row_i][col_j] = 0.5 * np.trace(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 
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
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 
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
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 
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
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 
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
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 
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
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 
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
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 
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
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 
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
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 
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
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 
 925
 926
 927
 928
 929
 930
 931
 932
 933
 934
 935
 936
 937
 938
 939
 940
 941
 942
 943
 944
 945
 946
 947
 948
 949
 950
 951
 952
 953
 954
 955
 956
 957
 958
 959
 960
 961
 962
 963
 964
 965
 966
 967
 968
 969
 970
 971
 972
 973
 974
 975
 976
 977
 978
 979
 980
 981
 982
 983
 984
 985
 986
 987
 988
 989
 990
 991
 992
 993
 994
 995
 996
 997
 998
 999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
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 
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
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 
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
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 
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
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 
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
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 
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
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 
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
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 
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
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 
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
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