diff --git a/README.md b/README.md index 48e2816..dcabd11 100644 --- a/README.md +++ b/README.md @@ -26,7 +26,7 @@ using TenSolver Q = randn(40, 40) -E, psi = TenSolver.solve(Q) +E, psi = TenSolver.minimize(Q) ``` The returned argument `E` is the calculated estimate for the minimum value, @@ -54,7 +54,7 @@ begin @variable(m, x[1:dim], Bin) @objective(m, Min, x'Q*x) - optimize!(m) # <-- Equivalent to running TenSolver.solve(Q) + optimize!(m) # <-- Equivalent to running TenSolver.minimize(Q) end ``` @@ -64,7 +64,7 @@ The solver uses the tensor networks machinery from [ITensors.jl](https://itensor which comes with GPU support for tensor contractions. To run the code in a GPU, all you have to do is passing the appropriate accelerator -to the `solve_qubo` method. +as a keyword to the solver. For example, the code below optimizes the QUBO using `CUDA.jl`. ```julia @@ -73,7 +73,7 @@ import CUDA: cu Q = randn(4, 4) -E, psi = solve(Q; device = CUDA.cu) +E, psi = minimize(Q; device = CUDA.cu) ``` Since ITensor's GPU platform support is always improving, diff --git a/src/TenSolver.jl b/src/TenSolver.jl index b317447..67968cd 100644 --- a/src/TenSolver.jl +++ b/src/TenSolver.jl @@ -9,7 +9,7 @@ include("solution.jl") export sample include("solver.jl") -export solve +export minimize, maximize ## ~:~ Welcome to the QUBOVerse ~:~ ## @@ -30,9 +30,9 @@ function QUBODrivers.sample(sampler::Optimizer{T}) where {T} n, l, Q, a, b = QUBOTools.qubo(sampler, :sparse; sense = :min) # Solve - energy, psi = solve(Q, l) + energy, psi = minimize(Q, l, b) - obj = a * (energy + b) + obj = a * energy x = sample(psi) s = QUBOTools.Sample{T,Int}(x, obj) diff --git a/src/solver.jl b/src/solver.jl index 0c020da..c0dd682 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -42,7 +42,7 @@ function tensorize(sites, Q::AbstractArray{T}, Qs...; cutoff = 1e-8) where T tensorize!(os, sites, x; cutoff) end - return MPO(T, os, sites) + return isempty(os) ? MPO(T,sites) : MPO(T, os, sites) end tensorize!(os, sites, ::Nothing; cutoff = 1e-8) = os @@ -96,13 +96,16 @@ function tensorize!( os:: OpSum end """ - solve(Q [, l, c ; device, cutoff, kwargs...) + minimize(Q::Matrix[, l::Vector[, c::Number ; device, cutoff, kwargs...) -Solve the Quadratic Unconstrained Binary Optimization problem +Solve the Quadratic Unconstrained Binary Optimization (QUBO) problem min b'Qb + l'b + c s.t. b_i in {0, 1} +Return the optimal value `E` and a probability distribution `ψ` over optimal solutions. +You can use [`sample`](@ref) to get an actual bitstring solving the QUBO. + This function uses DMRG with tensor networks to calculate the optimal solution, by finding the ground state (least eigenspace) of the Hamiltonian @@ -110,19 +113,21 @@ by finding the ground state (least eigenspace) of the Hamiltonian where D_i acts locally on the i-th qubit as [0 0; 0 1], i.e, the projection on |1>. +See also [`maximize`](@ref). + The optional keyword `device` controls whether the solver should run on CPU or GPU. For using a GPU, you can import the respective package, e.g. CUDA.jl, and pass its accelerator as argument. ```julia import CUDA -solve(Q; device = CUDA.cu) +minimize(Q; device = CUDA.cu) import Metal -solve(Q; device = Metal.mtl) +minimize(Q; device = Metal.mtl) ``` """ -function solve( Q :: AbstractMatrix{T} +function minimize( Q :: AbstractMatrix{T} , l :: Union{AbstractVector{T}, Nothing} = nothing , c :: T = zero(T) ; cutoff :: Float64 = 1e-8 @@ -153,3 +158,36 @@ function solve( Q :: AbstractMatrix{T} return obj, psi end + +""" +minimize(Q::Matrix, c::Number; kwargs...) + +Solve the Quadratic Unconstrained Binary Optimization problem +without a linear term. + + min b'Qb + c + s.t. b_i in {0, 1} + +See also [`maximize`](@ref). +""" +minimize(Q :: AbstractMatrix{T}, c :: T; kwargs...) where T = minimize(Q, nothing, c; kwargs...) + +""" +maximize(Q::Matrix[, l::Vector[, c::Number; kwargs...) + +Solve the Quadratic Unconstrained Binary Optimization problem +for maximization. + + max b'Qb + l'b + c + s.t. b_i in {0, 1} + +See also [`minimize`](@ref). +""" +function maximize(qs... ; kwargs...) + # Flip the sign of all non-nothing elements + # max p(x) = - min -p(x) + mqs = map(q -> maybe(-, q), qs) + E, psi = minimize(mqs...; kwargs...) + + return -E, psi +end diff --git a/test/cases/vrp.jl b/test/cases/vrp.jl index 2bb4633..88e772f 100644 --- a/test/cases/vrp.jl +++ b/test/cases/vrp.jl @@ -1,6 +1,7 @@ using LinearAlgebra -function vrp_simple() +# See https://secquoia.github.io/QUBOBook/QUBO%20and%20Ising%20Models.html +function vrp_6nodes() A = Float64[ 1 0 0 1 1 1 0 1 1 1 1 0 1 0 1 0 1 1 0 1 1 1 @@ -18,14 +19,13 @@ function vrp_simple() end @testset "Vehicle Routing Problem (VRP)" begin - # See https://secquoia.github.io/QUBOBook/QUBO%20and%20Ising%20Models.html - Q, beta = vrp_simple() + Q, beta = vrp_6nodes() - E, psi = solve(Q, nothing, beta) + E, psi = minimize(Q, beta; iterations = 40, cutoff = 1e-12) x = TenSolver.sample(psi) # Known Solution - @test E ≈ 5.0 atol=1e-4 + @test E ≈ 5.0 @test [0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0] in psi @test [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1] in psi diff --git a/test/qubo.jl b/test/qubo.jl index a930080..bb3c9e1 100644 --- a/test/qubo.jl +++ b/test/qubo.jl @@ -1,11 +1,82 @@ @testset "QUBO Correctness" begin dim = 5 + @testset "Ultra simple sanity checks" begin + @testset "Zero matrix" begin + E, psi = minimize([0.0 0; 0 0]) + + @test E ≈ 0.0 + for x in [[i, j] for i in 0:1, j in 0:1] + @test x in psi + end + end + + @testset "Zero matrix + constant" begin + E, psi = minimize([0.0 0; 0 0], 3.0) + + @test E ≈ 3.0 + for x in [[i, j] for i in 0:1, j in 0:1] + @test x in psi + end + end + + @testset "Zero matrix + zero linear" begin + E, psi = minimize([0.0 0; 0 0], [0.0, 0.0]) + + @test E ≈ 0.0 + for x in [[i, j] for i in 0:1, j in 0:1] + @test x in psi + end + end + + @testset "Zero matrix + linear" begin + E, psi = minimize([0.0 0; 0 0], [1.0, -1.0]) + + @test E ≈ -1.0 + @test TenSolver.sample(psi) == [0, 1] + end + + @testset "Zero matrix + linear + const" begin + E, psi = minimize([0.0 0; 0 0], [1.0, -1.0], 3.0) + + @test E ≈ 2.0 + @test TenSolver.sample(psi) == [0, 1] + end + + @testset "Identity" begin + E, psi = minimize([1.0 0.0; 0.0 1.0]) + + @test E ≈ 0.0 + @test TenSolver.sample(psi) == [0, 0] + end + + @testset "Max: Identity" begin + E, psi = maximize([1.0 0.0; 0.0 1.0]) + + @test E ≈ 2.0 + @test TenSolver.sample(psi) == [1, 1] + end + + @testset "Max: Identity + const" begin + E, psi = maximize([1.0 0.0; 0.0 1.0], 3.0) + + @test E ≈ 5.0 + @test TenSolver.sample(psi) == [1, 1] + end + + @testset "Max: Zero matrix + linear + const" begin + E, psi = maximize([0.0 0.0; 0.0 0.0], [1.0, -1.0], 3.0) + + @test E ≈ 4.0 + @test TenSolver.sample(psi) == [1, 0] + end + end + @testset "Pure quadratic" begin Q = randn(dim, dim) # TenSolver solution - e, psi = TenSolver.solve(Q; cutoff = 1e-16, iterations = 20) + e, psi = TenSolver.minimize(Q; cutoff = 1e-16, iterations = 20) x = TenSolver.sample(psi) # Is the sampled solution part of the ground state? @@ -39,7 +110,7 @@ obj(x) = dot(x, Q, x) + dot(l, x) + c # TenSolver solution - e, psi = TenSolver.solve(Q, l, c; cutoff = 1e-16) + e, psi = TenSolver.minimize(Q, l, c; cutoff = 1e-16) x = TenSolver.sample(psi) # Does the ground energy match solution? diff --git a/test/runtests.jl b/test/runtests.jl index 0ddab6c..bdec08e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -17,7 +17,7 @@ include(filepath("utils.jl")) # Should throw when the matrix is not square Q = randn(dim, dim - 2) - @test_throws BoundsError solve(Q) + @test_throws BoundsError minimize(Q) end # Traditional interface diff --git a/test/utils.jl b/test/utils.jl index c22dc2e..912cd4b 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -6,7 +6,7 @@ """ brute_force(f, n) -A version of `solve` that uses a brute force approach instead of Tensor networks. +A QUBO solver using a brute force approach instead of Tensor networks. Despite being painfully slow, this is useful as a sanity check. """ function brute_force(f::Function, T, n::Int64)