References
abbreviations:
import autograd.numpy as anp
from autograd import grad, jacobian
from autograd.extend import primitive, defvjp
def fun(x, y, z):
"""
x : float
y: float
z : bool
"""
return x * y**2 if z else y**2
def d_fun(x, y, z, dwrt):
if z:
if dwrt == 0:
return y**2
elif dwrt == 1:
return 2*x*y
else:
assert False, 'z not differentiable'
else:
if dwrt == 0:
return 0.
elif dwrt == 1:
return 2*y
else:
assert False, 'z not differentiable'
fun_x = grad(fun, 0)
fun_y = grad(fun, 1)
for z in [True,False]:
print(fun_x(.5, .1, z) - d_fun(.5, .1, z, dwrt=0))
print(fun_y(.5, .1, z) - d_fun(.5, .1, z, dwrt=1))
0.0 0.0 0.0 0.0
/home/rrtucci/anaconda3/envs/qubiter/lib/python3.7/site-packages/autograd/tracer.py:14: UserWarning: Output seems independent of input. warnings.warn("Output seems independent of input.")
@primitive
def pfun(x, y, z):
return fun(x, y, z)
defvjp(pfun,
lambda ans, x, y, z: lambda g: g*d_fun(x, y, z, 0), # g and pfun scalars
lambda ans, x, y, z: lambda g: g*d_fun(x, y, z, 1),
argnums=[0, 1])
pfun_x = grad(pfun, 0)
pfun_y = grad(pfun, 1)
for z in [True,False]:
print(pfun_x(.5, .1, z) - d_fun(.5, .1, z, 0))
print(pfun_y(.5, .1, z) - d_fun(.5, .1, z, 1))
0.0 0.0 0.0 0.0
sigx = anp.array([[0, 1],[1,0]])
sigx
array([[0, 1], [1, 0]])
sigy = anp.array([[0, -1j],[1j,0]])
sigy
array([[ 0.+0.j, -0.-1.j], [ 0.+1.j, 0.+0.j]])
sigz = anp.array([[1, 0],[0, -1]])
sigz
array([[ 1, 0], [ 0, -1]])
sig_all = anp.vstack([sigx, sigy, sigz])
sig_all
array([[ 0.+0.j, 1.+0.j], [ 1.+0.j, 0.+0.j], [ 0.+0.j, -0.-1.j], [ 0.+1.j, 0.+0.j], [ 1.+0.j, 0.+0.j], [ 0.+0.j, -1.+0.j]])
sig_all = anp.reshape(sig_all, (3, 2, 2)).transpose(1, 2, 0)
sig_all
array([[[ 0.+0.j, 0.+0.j, 1.+0.j], [ 1.+0.j, -0.-1.j, 0.+0.j]], [[ 1.+0.j, 0.+1.j, 0.+0.j], [ 0.+0.j, 0.+0.j, -1.+0.j]]])
ex = anp.array([1,0,0])
ey = anp.array([0,1,0])
ez = anp.array([0,0,1])
sigx_ = anp.dot(sig_all, ex)
sigy_ = anp.dot(sig_all, ey)
sigz_ = anp.dot(sig_all, ez)
print('sigx_=\n', sigx_)
print('sigy_=\n', sigy_)
print('sigz_=\n', sigz_)
sigx_= [[0.+0.j 1.+0.j] [1.+0.j 0.+0.j]] sigy_= [[0.+0.j 0.-1.j] [0.+1.j 0.+0.j]] sigz_= [[ 1.+0.j 0.+0.j] [ 0.+0.j -1.+0.j]]
$ U = e^{i\sigma_3\theta_3} = C + i\sigma_3 S$
$S = \sin\theta_3, C = \cos \theta_3$
$\frac{dU}{dt} = \dot{\theta}_3(-S + i\sigma_3 C)$
def rot1(t, ax):
assert ax in [1, 2, 3]
return anp.eye(2)*anp.cos(t) + 1j*sig_all[:, :, ax-1]*anp.sin(t)
def rot1_t(t, ax):
assert ax in [1, 2, 3]
return -anp.eye(2)*anp.sin(t) + 1j*sig_all[:, :, ax-1]*anp.cos(t)
def d_auto_rot1(t, ax):
assert ax in [1, 2, 3]
def rot1r(t, ax):
return anp.real(rot1(t, ax))
def rot1i(t, ax):
return anp.imag(rot1(t, ax))
return jacobian(rot1r, 0)(t, ax) + 1j*jacobian(rot1i, 0)(t, ax)
for ax in range(1, 4):
print(d_auto_rot1(.5, 1) - rot1_t(.5, 1))
[[0.+0.j 0.+0.j] [0.+0.j 0.+0.j]] [[0.+0.j 0.+0.j] [0.+0.j 0.+0.j]] [[0.+0.j 0.+0.j] [0.+0.j 0.+0.j]]
$ U = e^{i\sigma_k\theta_k} = C + i\sigma_k \frac{\theta_k}{\theta} S$
$\theta = \sqrt{\theta_k\theta_k}, S = \sin\theta, C = \cos \theta$
$\frac{dU}{dt}=-S \frac{\theta_k}{\theta} \dot{\theta_k}+ i\sigma_k\dot{\theta_r} \left[\frac{\theta_k\theta_r}{\theta^2} C+ \frac{S}{\theta}(-\frac{\theta_k\theta_r}{\theta^2}+\delta_{k, r})\right]$
def u2(*tlist):
assert len(tlist) == 4
t = anp.sqrt(tlist[1]**2 + tlist[2]**2 + tlist[3]**2)
tvec = anp.array([tlist[1]/t, tlist[2]/t, tlist[3]/t])
out = anp.eye(2)*anp.cos(t) + 1j*anp.dot(sig_all, tvec)*anp.sin(t)
return anp.exp(1j*tlist[0])*out
def d_u2(dwrt, *tlist):
assert dwrt in range(4)
assert len(tlist) == 4
if dwrt == 0:
return 1j*u2(*tlist)
dwrt -= 1
t = anp.sqrt(tlist[1]**2 + tlist[2]**2 + tlist[3]**2)
tvec = anp.array([tlist[1]/t, tlist[2]/t, tlist[3]/t])
dotted_vec = tvec*tvec[dwrt]*anp.cos(t) + (anp.sin(t)/t)*(-tvec*tvec[dwrt] + anp.eye(3)[dwrt, :])
out = -anp.sin(t)*tvec[dwrt]*anp.eye(2) + 1j*anp.dot(sig_all, dotted_vec)
return anp.exp(1j*tlist[0])*out
u2(.1, .2, .3, .4)
array([[ 0.81615064+0.46474597j, 0.26526593+0.21804425j], [-0.30329697+0.16099768j, 0.89221274-0.29333789j]])
u2(0., 0., 0., .4) - rot1(.4, ax=3)
array([[0.+0.j, 0.+0.j], [0.+0.j, 0.+0.j]])
def d_auto_u2(dwrt, *tlist):
def u2r(*tlist):
return anp.real(u2(*tlist))
def u2i(*tlist):
return anp.imag(u2(*tlist))
return jacobian(u2r, dwrt)(*tlist) + 1j*jacobian(u2i, dwrt)(*tlist)
for dwrt in range(4):
print(d_auto_u2(dwrt, 0., 0., 0., .5))
[[-0.47942554+0.87758256j 0. +0.j ] [ 0. +0.j 0.47942554+0.87758256j]] [[0.+0.j 0.+0.95885108j] [0.+0.95885108j 0.+0.j ]] [[ 0. +0.j 0.95885108+0.j] [-0.95885108+0.j 0. +0.j]] [[-0.47942554+0.87758256j 0. +0.j ] [ 0. +0.j -0.47942554-0.87758256j]]
/home/rrtucci/anaconda3/envs/qubiter/lib/python3.7/site-packages/numpy/core/_asarray.py:85: ComplexWarning: Casting complex values to real discards the imaginary part return array(a, dtype, copy=False, order=order)
rot1_t(.5, ax=3)
array([[-0.47942554+0.87758256j, 0. +0.j ], [ 0. +0.j , -0.47942554-0.87758256j]])
tlist = [.3, 1.1, .7, .5]
for dwrt in range(4):
print(anp.linalg.norm(d_auto_u2(dwrt, *tlist) - d_u2(dwrt, *tlist)))
0.0 2.5814177514652363e-16 1.4833660915169986e-16 2.786810452962585e-16
@primitive
def pu2r(*tlist):
return anp.real(u2(*tlist))
@primitive
def pu2i(*tlist):
return anp.imag(u2(*tlist))
def pu2(*tlist):
return pu2r(*tlist) + 1j*pu2i(*tlist)
defvjp(pu2r,
lambda ans, *tlist: lambda g: anp.sum(g*anp.real(d_u2(0, *tlist))), # g.shape==pu2r.shape
lambda ans, *tlist: lambda g: anp.sum(g*anp.real(d_u2(1, *tlist))),
lambda ans, *tlist: lambda g: anp.sum(g*anp.real(d_u2(2, *tlist))),
lambda ans, *tlist: lambda g: anp.sum(g*anp.real(d_u2(3, *tlist))),
argnums=range(4))
defvjp(pu2i,
lambda ans, *tlist: lambda g: anp.sum(g*anp.imag(d_u2(0, *tlist))), # g.shape==pu2i.shape
lambda ans, *tlist: lambda g: anp.sum(g*anp.imag(d_u2(1, *tlist))),
lambda ans, *tlist: lambda g: anp.sum(g*anp.imag(d_u2(2, *tlist))),
lambda ans, *tlist: lambda g: anp.sum(g*anp.imag(d_u2(3, *tlist))),
argnums=range(4))
def d_auto_pu2(dwrt, *tlist):
assert dwrt in range(4)
return jacobian(pu2r, dwrt)(*tlist) + 1j*jacobian(pu2i, dwrt)(*tlist)
tlist = [.1, .2, .3, .4]
for dwrt in range(4):
print(anp.linalg.norm(d_auto_pu2(dwrt, *tlist) - d_u2(dwrt, *tlist)))
0.0 0.0 0.0 0.0