Numerical instability related to squeezing gate with Fock backend when basis=False

[sansiro77]

Examples

for cutoff in [34, 36, 38, 40]:
    r = 2
    cir1 = dq.QumodeCircuit(nmode=2, init_state='vac', cutoff=cutoff, backend='fock', basis=False)
    cir1.s2([0,1], r=r)
    state1 = cir1()

    cir2 = dq.QumodeCircuit(nmode=2, init_state='vac', cutoff=cutoff, backend='fock', basis=False)
    cir2.bs([0,1], [torch.pi/4, 0])
    cir2.s(0, r=r)
    cir2.s(1, r=-r)
    cir2.bs([0,1], [-torch.pi/4, 0])
    state2 = cir2()

    print(f'state1 cutoff={cutoff}', state1.abs().max())
    print(f'state2 cutoff={cutoff}', state2.abs().max())

    print((state1 - state2).abs().max())

state1 cutoff=34 tensor(0.2658)
state2 cutoff=34 tensor(0.2658)
tensor(0.2431)
state1 cutoff=36 tensor(0.2658)
state2 cutoff=36 tensor(0.2658)
tensor(0.2614)
state1 cutoff=38 tensor(0.2658)
state2 cutoff=38 tensor(1.5054)
tensor(1.5149)
state1 cutoff=40 tensor(0.2658)
state2 cutoff=40 tensor(2.2099)
tensor(2.2099)

for cutoff in [40, 60, 100, 140]:
    r = 1
    cir1 = dq.QumodeCircuit(nmode=2, init_state='vac', cutoff=cutoff, backend='fock', basis=False)
    cir1.s(0, r=r)
    cir1.s(0, r=r)
    state1 = cir1()

    cir2 = dq.QumodeCircuit(nmode=2, init_state='vac', cutoff=cutoff, backend='fock', basis=False)
    cir2.s(0, r=2 * r)
    state2 = cir2()

    print(f'state1 cutoff={cutoff}', state1.abs().max())
    print(f'state2 cutoff={cutoff}', state2.abs().max())

    print((state1 - state2).abs().max())

state1 cutoff=40 tensor(0.5156)
state2 cutoff=40 tensor(0.5156)
tensor(0.0001)
state1 cutoff=60 tensor(0.5156)
state2 cutoff=60 tensor(0.5156)
tensor(1.2852e-05)
state1 cutoff=100 tensor(0.5156)
state2 cutoff=100 tensor(0.5156)
tensor(0.0768)
state1 cutoff=140 tensor(97.0402)
state2 cutoff=140 tensor(0.5156)
tensor(97.0294)

for cutoff in [2,4,8,16,32,64]:
    r = 1
    cir1 = dq.QumodeCircuit(nmode=2, init_state='vac', cutoff=cutoff, backend='fock', basis=False)
    cir1.s2([0,1], r=r)
    cir1.s2([0,1], r=r)
    state1 = cir1()

    cir2 = dq.QumodeCircuit(nmode=2, init_state='vac', cutoff=cutoff, backend='fock', basis=False)
    cir2.s2([0,1], r=2 * r)
    state2 = cir2()

    print(f'state1 cutoff={cutoff}', state1.abs().max())
    print(f'state2 cutoff={cutoff}', state2.abs().max())

    print((state1 - state2).abs().max())

state1 cutoff=2 tensor(0.2687)
state2 cutoff=2 tensor(0.2658)
tensor(0.0894)
state1 cutoff=4 tensor(0.2936)
state2 cutoff=4 tensor(0.2658)
tensor(0.0374)
state1 cutoff=8 tensor(0.2680)
state2 cutoff=8 tensor(0.2658)
tensor(0.0117)
state1 cutoff=16 tensor(0.2658)
state2 cutoff=16 tensor(0.2658)
tensor(0.0010)
state1 cutoff=32 tensor(0.2658)
state2 cutoff=32 tensor(0.2658)
tensor(9.9473e-05)
state1 cutoff=64 tensor(91.7321)
state2 cutoff=64 tensor(0.2658)
tensor(91.7057)

Suggestions

• Avoid using a large squeezing rate
• Avoid applying consecutive squeezing gates

[sansiro77]

Possible solution

• Use double instead of float

for cutoff in [34, 36, 38, 40]:
    r = 2
    cir1 = dq.QumodeCircuit(nmode=2, init_state='vac', cutoff=cutoff, backend='fock', basis=False)
    cir1.s2([0,1], r=r)
    cir1.to(torch.double)
    state1 = cir1()

    cir2 = dq.QumodeCircuit(nmode=2, init_state='vac', cutoff=cutoff, backend='fock', basis=False)
    cir2.bs([0,1], [torch.pi/4, 0])
    cir2.s(0, r=r)
    cir2.s(1, r=-r)
    cir2.bs([0,1], [-torch.pi/4, 0])
    cir2.to(torch.double)
    state2 = cir2()

    print(f'state1 cutoff={cutoff}', state1.abs().max())
    print(f'state2 cutoff={cutoff}', state2.abs().max())

    print((state1 - state2).abs().max())

state1 cutoff=34 tensor(0.2658, dtype=torch.float64)
state2 cutoff=34 tensor(0.2658, dtype=torch.float64)
tensor(0.0851, dtype=torch.float64)
state1 cutoff=36 tensor(0.2658, dtype=torch.float64)
state2 cutoff=36 tensor(0.2658, dtype=torch.float64)
tensor(0.0803, dtype=torch.float64)
state1 cutoff=38 tensor(0.2658, dtype=torch.float64)
state2 cutoff=38 tensor(0.2658, dtype=torch.float64)
tensor(0.0758, dtype=torch.float64)
state1 cutoff=40 tensor(0.2658, dtype=torch.float64)
state2 cutoff=40 tensor(0.2658, dtype=torch.float64)
tensor(0.0717, dtype=torch.float64)

• Better, but still unstable

for cutoff in [40, 60, 100, 140]:
    r = 1
    cir1 = dq.QumodeCircuit(nmode=2, init_state='vac', cutoff=cutoff, backend='fock', basis=False)
    cir1.s(0, r=r)
    cir1.s(0, r=r)
    cir1.to(torch.double)
    state1 = cir1()

    cir2 = dq.QumodeCircuit(nmode=2, init_state='vac', cutoff=cutoff, backend='fock', basis=False)
    cir2.s(0, r=2 * r)
    cir2.to(torch.double)
    state2 = cir2()

    print(f'state1 cutoff={cutoff}', state1.abs().max())
    print(f'state2 cutoff={cutoff}', state2.abs().max())

    print((state1 - state2).abs().max())

state1 cutoff=40 tensor(0.5156, dtype=torch.float64)
state2 cutoff=40 tensor(0.5156, dtype=torch.float64)
tensor(0.0001, dtype=torch.float64)
state1 cutoff=60 tensor(0.5156, dtype=torch.float64)
state2 cutoff=60 tensor(0.5156, dtype=torch.float64)
tensor(1.2860e-05, dtype=torch.float64)
state1 cutoff=100 tensor(0.5156, dtype=torch.float64)
state2 cutoff=100 tensor(0.5156, dtype=torch.float64)
tensor(0.0192, dtype=torch.float64)
state1 cutoff=140 tensor(27.8834, dtype=torch.float64)
state2 cutoff=140 tensor(0.5156, dtype=torch.float64)
tensor(27.8941, dtype=torch.float64)

for cutoff in [2,4,8,16,32,64]:
    r = 1
    cir1 = dq.QumodeCircuit(nmode=2, init_state='vac', cutoff=cutoff, backend='fock', basis=False)
    cir1.s2([0,1], r=r)
    cir1.s2([0,1], r=r)
    cir1.to(torch.double)
    state1 = cir1()

    cir2 = dq.QumodeCircuit(nmode=2, init_state='vac', cutoff=cutoff, backend='fock', basis=False)
    cir2.s2([0,1], r=2 * r)
    cir2.to(torch.double)
    state2 = cir2()

    print(f'state1 cutoff={cutoff}', state1.abs().max())
    print(f'state2 cutoff={cutoff}', state2.abs().max())

    print((state1 - state2).abs().max())

state1 cutoff=2 tensor(0.2687, dtype=torch.float64)
state2 cutoff=2 tensor(0.2658, dtype=torch.float64)
tensor(0.0894, dtype=torch.float64)
state1 cutoff=4 tensor(0.2936, dtype=torch.float64)
state2 cutoff=4 tensor(0.2658, dtype=torch.float64)
tensor(0.0374, dtype=torch.float64)
state1 cutoff=8 tensor(0.2680, dtype=torch.float64)
state2 cutoff=8 tensor(0.2658, dtype=torch.float64)
tensor(0.0117, dtype=torch.float64)
state1 cutoff=16 tensor(0.2658, dtype=torch.float64)
state2 cutoff=16 tensor(0.2658, dtype=torch.float64)
tensor(0.0010, dtype=torch.float64)
state1 cutoff=32 tensor(0.2658, dtype=torch.float64)
state2 cutoff=32 tensor(0.2658, dtype=torch.float64)
tensor(0.0002, dtype=torch.float64)
state1 cutoff=64 tensor(77.3696, dtype=torch.float64)
state2 cutoff=64 tensor(0.2658, dtype=torch.float64)
tensor(77.3432, dtype=torch.float64)

[sansiro77]

Reason

• sqrt in get_matrix_state is torch.float