Change solver API

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,

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)
 

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,33 +96,38 @@ 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
 
     H = Σ Q_ij D_iD_j + Σ l_i D_i
 
 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

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

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?

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

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)