Reducing the contraction order optimization time (quantum circuit as an example)

[samuelsonric]

This is a copy of this issue.

original

optimizer = TreeSA(ntrials=10, niters=20, sc_target=52)
time_elapsed = @elapsed net = yao2einsum(c; initial_state=Dict(zip(1:nqubits(c), zeros(Int,nqubits(c)))), final_state=Dict(zip(1:nqubits(c), zeros(Int,nqubits(c)))), optimizer)
@info "contraction complexity: $(contraction_complexity(net)), time cost: $(time_elapsed)s"

The output is

┌ Info: contraction complexity: Time complexity: 2^67.01176482416963
│ Space complexity: 2^52.0
└ Read-write complexity: 2^58.43881862537142, time cost: 52.74190275s

improved

using AbstractTrees, CliqueTrees, Graphs, OMEinsumContractionOrders, OMEinsum

struct TestAlg{A} <: CliqueTrees.EliminationAlgorithm
    alg::A
end

function eincode2order(code::NestedEinsum{L}) where {L}
    elimination_order = Vector{L}()
    OMEinsum.isleaf(code) && return elimination_order
    
    for node in PostOrderDFS(code)
        (node isa OMEinsum.LeafString) && continue
        for id in setdiff(vcat(getixsv(node.eins)...), getiyv(node.eins))
            push!(elimination_order, id)
        end
    end
    
    return elimination_order
end

function CliqueTrees.permutation(weights::AbstractVector, graph, alg::TestAlg)
    # reduce graph
    kernel, stack, label, width = CliqueTrees.pr4(graph, lowerbound(graph))
    
    # feed reduced graph back to OMEinsumContractionOrders
    code = EinCode(maximal_cliques(kernel), Int[])
    sizes = Dict{Int, Int}(v => round(Int, 2^weights[label[v]]) for v in vertices(kernel))
    opt = optimize_code(code, sizes, alg.alg).eins
    
    # compute ordering
    append!(stack, label[eincode2order(opt)])
    return stack, invperm(stack)
end

optimizer = Treewidth(; alg=TestAlg(TreeSA(ntrials=10, niters=20, sc_target=52)))
time_elapsed = @elapsed net = yao2einsum(c; initial_state=Dict(zip(1:nqubits(c), zeros(Int,nqubits(c)))), final_state=Dict(zip(1:nqubits(c), zeros(Int,nqubits(c)))), optimizer)
@info "contraction complexity: $(contraction_complexity(net)), time cost: $(time_elapsed)s"

The output is

┌ Info: contraction complexity: Time complexity: 2^60.20080142774071
│ Space complexity: 2^52.0
└ Read-write complexity: 2^58.2827569012071, time cost: 27.224232333s

[GiggleLiu]

It’s like a magic, do you have a reference about the algorithm of the reduction algorithm pr4? I do not understand why there can be a 2^7 time complexity reduction.

[samuelsonric]

I am using the rule set PR-4 (islet, twig, series, triangle, buddy, cube, simplical) along with the almost-simplicial rule.

• Bodlaender, Hans L., et al. "Pre-processing for triangulation of probabilistic networks." (2001).
• Bodlaender, Hans L., Arie MCA Koster, and Frank van den Eijkhof. "Preprocessing rules for triangulation of probabilistic networks." Computational Intelligence 21.3 (2005): 286-305.

[samuelsonric]

I should add that the eincode2order function was very helpfully pointed out to me by your student @ArrogantGao!

[GiggleLiu]

Thank you for the references!

[ArrogantGao]

It’s like a magic, do you have a reference about the algorithm of the reduction algorithm pr4? I do not understand why there can be a 2^7 time complexity reduction.

Time complexity is given by tree width, so I suggest to show the treewidth of both cases.

[samuelsonric]

treewidth 55 in the reduced case.

Even without reduction, encode -> graph -> eincode (via maximal cliques) seems to get the low-treewidth solution.

using AbstractTrees, CliqueTrees, Graphs, OMEinsumContractionOrders, OMEinsum

struct TestAlg{A} <: CliqueTrees.EliminationAlgorithm
    alg::A
end

struct TestAlg2{A} <: CliqueTrees.EliminationAlgorithm
    alg::A
end

function eincode2order(code::NestedEinsum{L}) where {L}
    elimination_order = Vector{L}()
    OMEinsum.isleaf(code) && return elimination_order
    
    for node in PostOrderDFS(code)
        (node isa OMEinsum.LeafString) && continue
        for id in setdiff(vcat(getixsv(node.eins)...), getiyv(node.eins))
            push!(elimination_order, id)
        end
    end
    
    return elimination_order
end

function CliqueTrees.permutation(weights::AbstractVector, graph, alg::TestAlg)
    # reduce graph
    kernel, stack, label, width = CliqueTrees.pr4(weights, graph, lowerbound(weights, graph))
    
    println(nv(kernel))
    
    # feed reduced graph to OMEinsumContractionOrders
    code = EinCode(maximal_cliques(kernel), Int[])
    sizes = Dict{Int, Int}(v => round(Int, 2^weights[label[v]]) for v in vertices(kernel))
    opt = optimize_code(code, sizes, alg.alg).eins
    
    # compute ordering
    append!(stack, label[eincode2order(opt)])
    
    println(treewidth(graph; alg=stack))
    
    return stack, invperm(stack)
end;

function CliqueTrees.permutation(weights::AbstractVector, graph, alg::TestAlg2)
    # feed graph back to OMEinsumContractionOrders
    code = EinCode(maximal_cliques(graph), Int[])
    sizes = Dict{Int, Int}(v => round(Int, 2^weights[v]) for v in vertices(graph))
    opt = optimize_code(code, sizes, alg.alg).eins
    
    # compute ordering
    stack = eincode2order(opt)

    println(treewidth(graph; alg=stack))
    
    return stack, invperm(stack)
end;

[samuelsonric]

I mean, both TestAlg and TestAlg2 give the same result, but the first is 2x as fast.

[GiggleLiu]

I am considering including the preprocessors like pr4 in your example into OMEinsumContractionOrders, and write a guide for users to speed up their computing. Does it make sense?

If yes, I'd like to repoen this issue ;D

[samuelsonric]

Hello @GiggleLiu! I just saw this comment. I would be happy to help you with this.

Also, I have been working on a new algorithm IND (incomplete nested dissection). It is similar in spirit to your algorithm KaHyParBipartite, but much older. It can achieve a space complexity of 2^53 and time complexity of 2^67 in 0.07 seconds. Not quite optimal...

julia> import Metis

julia> optimizer = Treewidth(; alg=IND(MF(), 100, METIS(; ufactor=150)));

julia> time_elapsed = @elapsed net = yao2einsum(c; initial_state=Dict(zip(1:nqubits(c), zeros(Int,nqubits(c)))), final_state=Dict(zip(1:nqubits(c), zeros(Int,nqubits(c)))), optimizer);

julia> @info "contraction complexity: $(contraction_complexity(net)), time cost: $(time_elapsed)s"
┌ Info: contraction complexity: Time complexity: 2^67.52024268871929
│ Space complexity: 2^53.0
└ Read-write complexity: 2^59.70052832972711, time cost: 0.070303334s

I pushed IND to GitHub, but I have not released a new version, yet.

[GiggleLiu]

Wow, that is amazing! I also tested the N-Queens example. Its performance is significantly improved.

julia> @info "contraction complexity (method: IND+METIS): $(contraction_complexity(optcode, uniformsize(optcode, 2)))"
┌ Info: contraction complexity (method: IND+METIS): Time complexity: 2^140.35769160666345
│ Space complexity: 2^99.0
└ Read-write complexity: 2^101.72695000549584

The previous best space complexity is 105, as recorded in: #62

To smoothen the testing, I have wrapped up the testing cases in this folder:
https://github.com/TensorBFS/OMEinsumContractionOrders.jl/tree/master/examples

I am wondering if it is possible to improve the outcome by giving it more time? Now it only takes 0.17s. Is it possible to find the optimal space complexity 82 in, like, 10min?

[samuelsonric]

Not with METIS, I think. I need a better vertex separator algorithm.

[GiggleLiu]

Not with METIS, I think. I need a better vertex separator algorithm.

It looks very similar to KaHyPar. BTW, this reminds me about another package using flow cutter for optimizing the contraction order of a tensor network: https://github.com/JuliaQX/QXGraphDecompositions.jl/blob/master/src/treewidth.jl

The performance of flow cutter is benchmarked here: JuliaQX/QXTools.jl#32

[samuelsonric]

import FlowCutterPACE17_jll
optimize_code(code, uniformsize(code, 2), Treewidth(alg=FlowCutter(; time=10)))

;)