# An outer-product in Yao?

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.

Thanks,

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.
``````

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,

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,

`ψ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
``````

Perfect thanks.