generate Haar random gates

ITensor Julia Questions
,
1

Hi,

I try to use itensor to simulate the measurement phase transition. Basically I want both Haar random unitary matrices and random measurements during the time evolution.

So I wrote a short code

using LinearAlgebra
using ITensors
using StatsBase
using Statistics
using StaticArrays
using Random

function rand_haar2(::Val{n}) where n
    M = randn(SMatrix{n,n,ComplexF64}) 
    q = qr(M).Q
    L = cispi.(2 .* @SVector(rand(n)))
    return q*diagm(L)
end

function op(::OpName"haar2",::SiteType"S=1/2",s1::Index,s2::Index)
    h = rand_haar2(Val{4}())
    return itensor(h,s2',s1',s2,s1)
end

"""
psi: wavefunction, n: position where the projector is applyied, p:either up +1 or down -1
"""
function measure_proj(psi::MPS,n,p)
    site_p = (0.5*(op("Id",siteind(psi,n))+2*p*op("Sz",siteind(psi,n))))*psi[n]
    psi[n] = noprime!(site_p)
    psi[n] = psi[n]/(norm(psi[n]))
    return psi
end

function measure_prob(psi::MPS, n)
    orthogonalize!(psi,n)
    sn = siteind(psi,n)
    p = scalar(dag(prime(psi[n],"Site"))*(0.5*(op("Id",sn)+2*op("Sz",sn)))*psi[n])
    return real(p)
end

n = 10
s = siteinds("S=1/2",n)
psi = randomMPS(s)
Gates=[("haar2",3,4)]
apply(ops(Gates,s),psi;cutoff = 1e^-10)

But I got error

ERROR: Older op interface does not support multiple indices with mixed site types. You may want to overload `op(::OpName, ::SiteType..., ::Index...)` or `op!(::ITensor, ::OpName, ::SiteType..., ::Index...) for the operator "haar2" and Index tags ("S=1/2,Site,n=3", "S=1/2,Site,n=4").

I think I should use index rather than integer 3 and 4, but how to fix it?

Also, I use some static arrays since I want to generate random 4x4 matrices frequently. Is this allowed?

Thanks

2

It looks like your first issue is that you didn’t import ITensors.op. The error message is pretty bad, basically it is caused by the fact that it tried to find the op definition for "haar2" but couldn’t. I’ll fix the error message to make that clearer.

To overload op you can write:

import ITensors: op

at the top of your script.

Also, ITensors don’t really handle wrapping StaticArrays right now, so you should change your gate definition to:

function op(::OpName"haar2",::SiteType"S=1/2",s1::Index,s2::Index)
  h = rand_haar2(Val{4}())
   return itensor(Matrix(h),s2',s1',s2,s1)
end

For this case, using a StaticArray wouldn’t help much anyway, you could simplify your code a lot by just using regular Julia Arrays to begin with. You can see the implementation we have in PastaQ.jl here:

It makes use of a function qr_positive from NDTensors.jl which performs a QR decomposition and then rotates Q to make the diagonals of R positive. You should be able to use the same function in your example as well.

1 Like
3

Hi Matthew. Thanks a lot. I will change my haar random function. The reason I try to use StaticArrays is that in my previous work, it turns out that after generating a lot of random matrices (after many time slices), the speed of code will decrease quickly due to allocation of many mutable arrays.

4

Interesting. In general, the tensor contractions and factorizations involved in doing the gate evolution should dominate the cost of the calculation, so you shouldn’t notice the cost of generating the gates. But I could imagine certain limits where generating the gates could be noticeable (for example very small MPS bond dimensions or system sizes).

You should be able to time it easily with:

U = @time ops(Gates,s)
@time apply(ops(U, s),psi;cutoff = 1e^-10)

I would be curious if you have applications where the gate generation is more costly than the gate evolution.

5

Sorry I didn’t make it clear. This phenomenon happened in a related but different problem where the time evolution is replaced by a recursion relation. In that case a large number of sampling is needed to get converged result. I think in the usual MPT problem, the time evolution will be the dominant term. Thanks

1 Like
6

@jisutich in ITensors v0.3.9 there will be a "RandomUnitary" operator defined that makes Haar-random unitaries that you can use for this so you don’t need to define your own. It should be available soon.

See: [ITensors] [ENHANCMENT] Haar random unitary gate by mtfishman · Pull Request #903 · ITensor/ITensors.jl · GitHub

2 Likes