An outer-product in Yao?

Specific Domains Quantum
1

Hello,
I am trying to evaluate the following expression:
I3 + ψ*ψ'

First, is there a way in Yao to create an arbitrary sized Identity matrix like I3?

using Yao, YaoPlots
using YaoBlocks: eigenbasis

ψ=ArrayReg(bit"111") |> normalize!
ψ = rand_state(3);
I3 = [
  1 0 0 0 0 0 0 0
  0 1 0 0 0 0 0 0
  0 0 1 0 0 0 0 0
  0 0 0 1 0 0 0 0
  0 0 0 0 1 0 0 0
  0 0 0 0 0 1 0 0
  0 0 0 0 0 0 1 0
  0 0 0 0 0 0 0 1
]

ψ' * ψ

This fails to run, I must be missing something.

print(ψ*adjoint(ψ))

Thanks,

2

We already have similar feature in Yao master branch, which is called projector:
Yao.jl/reflect.jl at master · QuantumBFS/Yao.jl · GitHub. You can try this feature in develop mode already:

(@v1.8) pkg> dev YaoAPI YaoArrayRegister YaoBlocks YaoSym Yao

julia> using Yao

julia> igate(3) + projector(rand_state(3))
nqubits: 3
+
├─ igate(3)
└─ |s⟩⟨s|, nqudits = 3

To stick to the released version, I think we have to use the following implementation should be good enough.

julia> igate(nqubits(ψ)) + matblock(OuterProduct(statevec(ψ), conj.(statevec(ψ))); tag="|ψ⟩⟨ψ|")
nqubits: 12
+
├─ igate(12)
└─ |ψ⟩⟨ψ|

If you want to apply this operator to a register. The block implementation is much faster than computing the matrix explicitly,

julia> ψ = rand_state(12);

julia> @benchmark apply($ψ, igate(nqubits($ψ)) + matblock(OuterProduct(statevec($ψ), conj.(statevec($ψ)))))
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max):  25.104 μs … 998.618 μs  ┊ GC (min … max):  0.00% … 93.86%
 Time  (median):     28.515 μs               ┊ GC (median):     0.00%
 Time  (mean ± σ):   33.041 μs ±  55.170 μs  ┊ GC (mean ± σ):  10.66% ±  6.16%

    ▃▄▆▇████▇▇▆▆▅▄▃▂▁▁              ▁               ▁          ▃
  ▆█████████████████████▇█▆▅▅▅▅▅▆▅▆▇████▇█▆▇▄▆▅▄▆▇▇████▇█▆▆▆▅▄ █
  25.1 μs       Histogram: log(frequency) by time      47.8 μs <

 Memory estimate: 192.95 KiB, allocs estimate: 27.

julia> @benchmark mat(igate(12)) + statevec($ψ) * statevec($ψ')
BenchmarkTools.Trial: 26 samples with 1 evaluation.
 Range (min … max):  158.426 ms … 331.332 ms  ┊ GC (min … max): 17.13% … 52.31%
 Time  (median):     189.718 ms               ┊ GC (median):    14.30%
 Time  (mean ± σ):   192.315 ms ±  40.352 ms  ┊ GC (mean ± σ):  18.93% ± 10.08%

  ▂         ▂█                                                   
  █▁▁▃▁▁▁▁▁▁██▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▃▁▁▁▁▁▁▃ ▁
  158 ms           Histogram: frequency by time          331 ms <

 Memory estimate: 512.00 MiB, allocs estimate: 9.
3

Thank you so much for your prompt reply.
Let me understand whats going on since I have only now started using Julia. Here, in the expression:

ψ=ArrayReg(bit"011") |> normalize!
I3 =igate(3)
ψ1= I3 - 2*projector(ψ)
@show(ψ1)
ψ1 = nqubits: 3
+
├─ igate(3)
└─ [scale: -2] |s⟩⟨s|, nqudits = 3

ψ1 should have resulted in an evaluated numerical value, not an operator right? How do I see that value?

Here:
ψ1=igate(nqubits(ψ)) + matblock(OuterProduct(statevec(ψ), conj.(statevec(ψ))); tag="|ψ⟩⟨ψ|")
The result of the outer product should have been compatible, what does it have to be coerced into a matblock?

I am trying to do something very simple I3 - 2*|ψ⟩⟨ψ| to yield a numerical value.

Thanks,

4

OK I see what I did wrong, I have a state ψ2 = rand_state(3) which I want the operator to act on.

using Yao, YaoPlots
using YaoBlocks: eigenbasis
# state preparation 
cr=(    
    kron(H, H, H))
plot(cr)
ψ0=ArrayReg(bit"000") |> normalize!
ψ0 = (ψ0 |> cr) 
@show (state(ψ0))

#Operator
ψ011=ArrayReg(bit"011") |> normalize!
I3 =igate(3)
ψ_op= I3 - 2*projector(ψ011)

#Apply the operator on the state
print(apply(ψ0,ψ_op))

How do I access the unitary of the operator?
How do I access the resulting value of applying the op to the state?

Thanks,

5

ψ1 should have resulted in an evaluated numerical value, not an operator right? How do I see that value?

It is an operator. If you want to get its matrix representation, you can use

mat(ψ1)

Maybe you will find this notebook helpful: https://giggleliu.github.io/notebooks/notebooks/yaoblocks.html

To obtain the resulting value of applying the op to the state reg, please just use

julia> reg = apply(state, op);

julia> reg.state
6

Perfect thanks.