Click here to go to the documentation for the latest stable release.'; + div.appendChild(closer); + document.body.appendChild(div); +} + +if (document.readyState === "loading") { + document.addEventListener("DOMContentLoaded", maybeAddWarning); +} else { + maybeAddWarning(); +} diff --git a/previews/PR82/booklet/1-intro/index.html b/previews/PR82/booklet/1-intro/index.html new file mode 100644 index 00000000..9b603295 --- /dev/null +++ b/previews/PR82/booklet/1-intro/index.html @@ -0,0 +1,2 @@ + +
ToQUBO.jl Booklet
This booklet aims to gather the theoretical and practical details behind ToQUBO and provide documentation for project internals. The target audience includes, among others, advanced users and those willing to contribute to the project. The latter are advised to read the following sections, as they give a glimpse of the ideas employed up to now.
Table of Contents
Settings
This document was generated with Documenter.jl version 1.1.2 on Sunday 26 November 2023. Using Julia version 1.9.4.
QUBO
Definition
Quadratic Unconstrained Binary Optimization, as the name suggests, refers to the global minimization or maximization of a quadratic polynomial on binary variables. A common presentation, the quadratic matrix form, is written as
\[\begin{array}{rl} + \min_{\mathbf{x}} & \mathbf{x}' Q\,\mathbf{x} \\ + \textrm{s.t.} & \mathbf{x} \in \mathbb{B}^{n} +\end{array}\]
where $Q \in \mathbb{R}^{n \times n}$ is symmetric and $\mathbb{B} = \lbrace{0, 1}\rbrace$. Note that, since $x^{2} = x$ holds for $x \in \mathbb{B}$, the linear terms of the objective function are stored in the main diagonal of $Q$.
OK, but why QUBO?
Mathematically speaking, there is a notorious equivalence between QUBO and Max-Cut problems, e.g. for every QUBO instance there is an information preserving Max-Cut reformulation and vice versa. This statement is followed by two immediate implications:
- In the general case, solving QUBO globally is NP-Hard.
- It is a simple yet expressive mathematical programming framework.
Implication 1. tells us that such problems are computationally intractable and that heuristics and metaheuristics are to be employed instead of exact methods. No 2. relates to the fact that we are able to represent many other optimization models by means of the QUBO formalism.
The Ising Model, on the other hand, is a mathematical abstraction to describe statistical interactions within mechanical systems with interesting properties for encoding combinatorial problems. Its Hamiltonian leads to an optimization formulation in terms of the spin values of their states, given by
\[\begin{array}{rl} + \min_{\mathbf{s}} & \mathbf{h}'\mathbf{s} + \mathbf{s}' J\,\mathbf{s} \\ + \textrm{s.t.} & \mathbf{s} \in \lbrace{-1,+1}\rbrace^{n} +\end{array}\]
with strictly upper triangular $J \in \mathbb{R}^{n \times n}$ and $\mathbf{h} \in \mathbb{R}$.
The Ising reformulation alternative draws the bridge between QUBO problems and devices designed to sample global or approximate ground states of the Ising Hamiltonian with high probability[Mohseni2022]. Some of the paradigms that stand out in this context are quantum gate-based optimization algorithms (QAOA and VQE), quantum annealers, hardware-accelerated platforms (Coherent Ising Machines and Simulated Bifurcation Machines) and physics-inspired methods (Simulated Annealing, Parallel Tempering). The significant advances in these computing systems contributed to the growing popularity of the model across the literature.
- Mohseni2022Mohseni, N., McMahon, P. L. & Byrnes, T. Ising machines as hardware solvers of combinatorial optimization problems. Nat Rev Phys 4, 363–379 (2022). {arXiv}
Settings
This document was generated with Documenter.jl version 1.1.2 on Sunday 26 November 2023. Using Julia version 1.9.4.
Pseudo-Boolean Optimization
Internally, problems are represented through a Pseudo-Boolean Optimization (PBO) framework. The main goal is to represent a given problem using a Pseudo-Boolean Function (PBF) since there is an immediate correspondence between optimization over quadratic PBFs and the QUBO formalism.
Quadratization
In order to successfully achieve a QUBO formulation, sometimes it is needed to quadratize the resulting PBF, i.e., reduce its degree until reaching the quadratic case.
A quadratization is a mapping $\mathcal{Q}: \mathscr{F} \to \mathscr{F}^{2}$ such that
\[\forall f \in \mathscr{F}, \forall x \in \{0, 1\}^{n}, \min_{y} \mathcal{Q}\left\lbrace{}f\right\rbrace{}(x; y) = f(x)\]
There are many quadratization methods available[Dattani2019], and ToQUBO implements two of them for now. However, using Julia's multiple dispatch paradigm, it's possible to extend the quadratization method coverage with your own algorithms.
Quadratization Methods
Stable Quadratization
The quadratization of a PBF does not guarantee that the resulting function will always be the same, as the order of terms can be different each time. This can be an issue in some situations where a deterministic output is required.
With said that, we have introduced the concept of Stable Quadratization, where the terms of the PBF are sorted, guaranteeing that the resulting PBF will be the same every time. We have defined it as an attribute of the compiler, with the ToQUBO.Attributes.StableQuadratization flag.
- Dattani2019Nikesh S. Dattani, Quadratization in discrete optimization and quantum mechanics, ArXiv, 2019 {doi}
Settings
This document was generated with Documenter.jl version 1.1.2 on Sunday 26 November 2023. Using Julia version 1.9.4.
Encoding Methods
ToQUBO.Encoding.encode — Functionencode(var, e::VariableEncodingMethod, x::Union{VI,Nothing}, S)ToQUBO.Encoding.encode! — Functionencode!(target, source...)ToQUBO.Encoding.encodes — Functionencodes(f::AbstractPBF, S::Tuple{T,T}, tol::T) where {T}Variables
As you may already know, QUBO models are comprised only of binary variables. So when we are reformulating general optimization problems, one important step is to encode variables into binary ones.
ToQUBO currently implements 6 encoding techniques. Each method introduces a different number of variables, quadratic terms and linear terms. Also, they differ in the magnitude of their coefficients $\Delta$.
| Encoding | Binary Variables | # Linear terms | # Quadratic terms | $\Delta$ |
|---|---|---|---|---|
| Binary | $O(\log n)$ | $O(\log n)$ | - | $O(n)$ |
| Unary | $O(n)$ | $O(n)$ | - | $O(1)$ |
| One-Hot | $O(n)$ | $O(n)$ | $O(n^2)$ | $O(n)$ |
| Domain-Wall | $O(n)$ | $O(n)$ | $O(n)$ | $O(n)$ |
| Bounded-Coefficient | $O(n)$ | $O(n)$ | - | $O(1)$ |
| Arithmetic Prog | $O(\sqrt{n})$ | $O(\sqrt{n})$ | - | $O(\sqrt{n})$ |
Mirror Encoding
ToQUBO.Encoding.VariableEncodingMethod — TypeVariableEncodingMethodAbstract type for variable encoding methods.
ToQUBO.Encoding.Mirror — TypeMirror()Simply mirrors a binary variable $x \in \mathbb{B}$ with a twin variable $y \in \mathbb{B}$.
Interval Encoding
ToQUBO.Encoding.IntervalVariableEncodingMethod — TypeIntervalVariableEncodingMethodAbstract type for methods that encode variables using a linear function, e.g.,
\[\xi(\mathbf{y}) = \beta + \sum_{i = 1}^{n} \gamma_{i} y_{i}\]
ToQUBO.Encoding.Unary — TypeUnary{T}()Integer
Let $x \in [a, b] \subset \mathbb{Z}$, $n = b - a$ and $\mathbf{y} \in \mathbb{B}^{n}$.
\[\xi{[a, b]}(\mathbf{y}) = a + \sum_{j = 1}^{b - a} y_{j}\]
Real
Given $n \in \mathbb{N}$ for $x \in [a, b] \subset \mathbb{R}$,
\[\xi{[a, b]}(\mathbf{y}) = a + \frac{b - a}{n} \sum_{j = 1}^{n} y_{j}\]
Representation error
Given $\tau > 0$, for the expected encoding error to be less than or equal to $\tau$, at least
\[n \ge 1 + \frac{b - a}{4 \tau}\]
binary variables become necessary.
ToQUBO.Encoding.Binary — TypeBinary{T}()Integer
Let $x \in [a, b] \subset \mathbb{Z}$, $n = \left\lceil \log_{2}(b - a) + 1 \right\rceil$ and $\mathbf{y} \in \mathbb{B}^{n}$.
\[\xi{[a, b]}(\mathbf{y}) = a + \left(b - a - 2^{n - 1} + 1\right) y_{n} + \sum_{j = 1}^{n - 1} 2^{j - 1} y_{j}\]
Real
Given $n \in \mathbb{N}$ for $x \in [a, b] \subset \mathbb{R}$,
\[\xi{[a, b]}(\mathbf{y}) = a + \frac{b - a}{2^{n} - 1} \sum_{j = 1}^{n} 2^{j - 1} y_{j}\]
Representation error
Given $\tau > 0$, for the expected encoding error to be less than or equal to $\tau$, at least
\[n \ge \log_{2} \left[1 + \frac{b - a}{4 \tau}\right]\]
binary variables become necessary.
ToQUBO.Encoding.Arithmetic — TypeArithmetic{T}()Integer
Let $x \in [a, b] \subset \mathbb{Z}$, $n = \left\lceil{ \frac{1}{2} {\sqrt{1 + 8 (b - a)}} - \frac{1}{2} }\right\rceil$ and $\mathbf{y} \in \mathbb{B}^{n}$.
\[\xi{[a, b]}(\mathbf{y}) = a + \left( {b - a - \frac{n (n - 1)}{2}} \right) y_{n} + \sum_{j = 1}^{n - 1} j y_{j}\]
Real
Given $n \in \mathbb{N}$ for $x \in [a, b] \subset \mathbb{R}$,
\[\xi{[a, b]}(\mathbf{y}) = a + \frac{b - a}{n (n + 1)} \sum_{j = 1}^{n} j y_{j}\]
Representation error
Given $\tau > 0$, for the expected encoding error to be less than or equal to $\tau$, at least
\[n \ge \frac{1}{2} \left[ 1 + \sqrt{3 + \frac{(b - a)}{2 \tau})} \right]\]
Bounded Coefficients
ToQUBO.Encoding.Bounded — TypeBounded{E,T}(μ::T) where {E<:Encoding,T}The bounded-coefficient encoding method[Karimi2019] consists in limiting the magnitude of the coefficients in the encoding expansion to a parameter $\mu$. This can be applied to the Unary, Binary and Arithmetic encoding methods.
Let $\xi[a, b] : \mathbb{B}^{n} \to [a, b]$ be an encoding function over the closed interval $[a, b]$. The bounded-coefficient encoding function given $\mu$ is defined as
\[\xi_{\mu}[a, b] = \xi[0, \delta](y_{1}, \dots, y_{k}) + \sum_{j = k + 1}^{r} \mu y{j}\]
```
Arbitrary Set Encoding
ToQUBO.Encoding.SetVariableEncodingMethod — TypeSetVariableEncodingMethodAbstract type for methods that encode variables over an arbitrary set.
ToQUBO.Encoding.OneHot — TypeOneHot{T}()The one-hot encoding is a linear technique used to represent a variable $x \in \set{\gamma_{j}}_{j \in [n]}$.
The associated encoding function is combined with a constraint assuring that only one and exactly one of the expansion's variables $y_{j}$ is activated at a time.
\[\xi[\set{\gamma_{j}}_{j \in [n]}](\mathbf{y}) = \sum_{j = 1}^{n} \gamma_{j} y_{j} ~\textrm{s.t.}~ \sum_{j = 1}^{n} y_{j} = 1\]
When a variable is encoded following this approach, a penalty term of the form
\[\rho \left[ \sum_{j = 1}^{n} y_{j} - 1 \right]^{2}\]
is added to the objective function.
ToQUBO.Encoding.DomainWall — TypeDomainWall{T}()The Domain Wall[Chancellor2019] encoding method is a sequential approach that requires $n - 1$ bits to represent $n$ distinct values.
\[\xi{[\set{\gamma_{j}}_{j \in [n]}]}(\mathbf{y}) = \sum_{j = 1}^{n} \gamma_{j} (y_{j} - y_{j + 1}) ~\textrm{s.t.}~ \sum_{j = 1}^{n} y_{j} \oplus y_{j + 1} = 1, y_{1} = 1, y_{n + 1} = 0\]
where $\mathbf{y} \in \mathbb{B}^{n + 1}$.
Representation Error
ToQUBO.Encoding.encoding_bits — Functionencoding_bits(e::VariableEncodingMethod, S::Tuple{T,T}, tol::T) where {T}ToQUBO.Encoding.encoding_points — Functionencoding_points(e::SetVariableEncodingMethod, S::Tuple{T,T}, tol::T) where {T}Let $\set{x_{i}}_{i \in [k]}$ be the collection of $k$ evenly spaced samples from the discretization of an interval $[a, b] \subseteq \mathbb{R}$.
The representation error for a given point $x$ with respect to $\set{x_{i}}_{i \in [k]}$ is
\[e_{k}(x) = \min_{i \in [k]} \left|x - x_{i}\right|\]
Assuming that $x$ behaves as a uniformly distributed random variable, the expected absolute encoding error is
\[\begin{align*} +\mathbb{E}\left[{e_{k}(x)}\right] &= \frac{1}{b - a} \int_{a}^{b} e_{k}(x) ~\mathrm{d}x \\ + &= \frac{1}{4} \frac{b - a}{k - 1} +\end{align*}\]
Thus, for encoding methods that rely on the regular division of an interval, it is possible to define the number of samples $k$ necessary to limit the expected error according to an upper bound $\tau$, that is,
\[\mathbb{E}\left[{e_{k}(x)}\right] \le \tau \implies k \ge 1 + \frac{b - a}{4 \tau}\]
This allows the compiler to automatically infer the number of bits to allocate for an encoded variable given the tolerance factor.
Constraints
A QUBO model is unconstrained. So when ToQUBO is reformulating a problem, it needs to encode all constraints into the objective function losing as little information as possible.
As constraints are introduced into the objective function, we need to make sure that they won't be violated. In order to do that, ToQUBO multiplies the encoded constraint by a large penalty $\rho$, so that any violation would result in a sub-optimal solution to the problem.
Sometimes, the encoding process might introduce higher-order terms, demanding ToQUBO to reduce the offending polynomials back to a quadratic form.
- Karimi2019Karimi, S. & Ronagh, P. Practical integer-to-binary mapping for quantum annealers. Quantum Inf Process 18, 94 (2019). {doi}
- Chancellor2019Nicholas Chancellor, Domain wall encoding of discrete variables for quantum annealing and QAOA, Quantum Science Technology 4, 2019.
Settings
This document was generated with Documenter.jl version 1.1.2 on Sunday 26 November 2023. Using Julia version 1.9.4.
Virtual Mapping
During reformulation, ToQUBO holds two distinct models, namely the Source Model and the Target Model. The source model is a generic MOI model restricted to the supported constraints. The target one is on the QUBO form used during the solving process. Both lie within a Virtual Model, which provides the necessary API integration and keeps all variable and constraint mapping tied together.
This is done in a transparent fashion for both agents since the user will mostly interact with the presented model, and the solvers will only access the generated one.
Virtual Model
ToQUBO.Virtual.Model — TypeModel{T}(optimizer::Union{Nothing, Type{<:MOI.AbstractOptimizer}} = nothing) where {T}This Virtual Model links the final QUBO formulation to the original one, allowing variable value retrieving and other features.
ToQUBO.Virtual.source — FunctionsourceToQUBO.Virtual.target — FunctiontargetVirtual Variables
Every virtual model stores a collection of virtual variables, intended to provide a link between those in the source and those to be created in the target model. Each virtual variable stores encoding information for later expansion and evaluation.
ToQUBO.Virtual.Variable — TypeVariable{T}ToQUBO.Virtual.encoding — FunctionencodingToQUBO.Virtual.expansion — FunctionexpansionToQUBO.Virtual.penaltyfn — FunctionpenaltyfnSettings
This document was generated with Documenter.jl version 1.1.2 on Sunday 26 November 2023. Using Julia version 1.9.4.
The Compiler
We hope to write this part of the documentation soon. Please come back later!
Compilation Steps
Setup
ToQUBO.Compiler.reset! — Functionreset!(model::Virtual.Model, arch::AbstractArchitecture)ToQUBO.Compiler.setup! — Functionsetup!(model::Virtual.Model, ::AbstractArchitecture)Parsing
ToQUBO.Compiler.parse! — Functionparse!(
+ model::Virtual.Model{T},
+ g::PBO.PBF{VI,T},
+ f::MOI.AbstractFunction,
+ arch::AbstractArchitecture
+) where {T}Parses the given MOI function f into PBF g.
ToQUBO.Compiler._parse — FunctionReformulation
ToQUBO.Compiler.sense! — Functionsense!(model::Virtual.Model, ::AbstractArchitecture)Copies MOI.ObjectiveSense from model.source_model to model.target_model.
ToQUBO.Compiler.variable! — Functionvariable!(model::Virtual.Model{T}) where {T}ToQUBO.Compiler.variables! — Functionvariables!(model::Virtual.Model{T}) where {T}ToQUBO.Compiler.constraint — FunctionconstraintReturns the pseudo-boolean function associated to a given constraint from the source model.
ToQUBO.Compiler.constraints! — Functionconstraints!(model::Virtual.Model, ::AbstractArchitecture)ToQUBO.Compiler.objective! — Functionobjective!(model::Virtual.Model, ::AbstractArchitecture)ToQUBO.Compiler.penalties! — Functionpenalties!(model::Virtual.Model, arch::AbstractArchitecture)Copying
ToQUBO.Compiler.isqubo — Functionisqubo(model::MOI.ModelLike)Tells if a given model is ready to be interpreted as a QUBO model.
For it to be true, a few conditions must be met:
- All variables must be binary (
MOI.VariableIndex ∈ MOI.ZeroOne) - No other constraints are allowed
- The objective function must be of type
MOI.ScalarQuadraticFunction,MOI.ScalarAffineFunctionorMOI.VariableIndex - The objective sense must be either
MOI.MIN_SENSEorMOI.MAX_SENSE
ToQUBO.Compiler.copy! — Functioncopy!(model::Virtual.Model, arch::AbstractArchitecture)Hamiltonian Assembly
ToQUBO.Compiler.build! — Functionbuild!(model::Virtual.Model, arch::AbstractArchitecture)ToQUBO.Compiler.quadratize! — Functionquadratize!(model::Virtual.Model, arch::AbstractArchitecture)Quadratizes the objective function from a model.
Settings
This document was generated with Documenter.jl version 1.1.2 on Sunday 26 November 2023. Using Julia version 1.9.4.
QUBO Solvers
Solvers, Annealers & Samplers
ToQUBO.jl's main goal is to make use of parameterized stochastic optimization solvers, particularly those relying on non-conventional hardware such as Quantum Annealing and other Ising Machines. A few MOI-compliant interfaces for annealers and samplers are bundled within ToQUBO.jl via the QUBODrivers.jl companion package. Some of them are presented below.
Simulated Annealing
Provided by D-Wave's open-source code libraries, this Simulated Annealing engine implements some of the features and configurations you would find using the Quantum API. Its adoption is recommended for basic usage, tests, and research due to its robustness, simplicity and ease of use. The DWave.jl's DWave.Neal module uses QUBODrivers.jl to deliver an interface to this sampler.
Quantum Annealing
Interfacing with D-Wave's quantum annealer is one of the milestones we expect to achieve with this package. Like other proprietary optimization resources such as Gurobi, FICO® Xpress and IBM® CPLEX®, this requires licensing and extra steps are needed to get an access token. In a first moment, for those willing to get started, the DWave.Neal optimizer might be enough to learn the ropes.
Random Sampling
This sampler is implemented for test purposes and simply assigns 0 or 1 to each variable according to a given probability bias $0 \le p \le 1$, which defaults to $p = 0.5$. After running the using QUBODrivers command, RandomSampler.Optimizer will be available.
Exact Solver (Exhaustive Enumeration)
Also made to be used in tests, the ExactSolver.Optimizer interface runs through all possible state configurations, which implies in an exponential time complexity on the number of variables. Thus, only problems with at most $\approxeq 20$ variables should be provided since visiting $2^{20} \approxeq 10^{6}$ states can already take up to a few seconds.
Mixed-Integer Quadratic Programming
The most accessible alternative to the forementioned methods are Mixed-Integer Quadratic Programming (MIQP) solvers such as Gurobi, CPLEX, SCIP and BARON. These are not intended to be of regular use alongside ToQUBO.jl since providing a QUBO reformulation will usually make things harder for non-specialized solvers. Yet, there are still a few cases where they may be suitable, such as tests, benchmarks, or any other situation where global optimality is a must.
Settings
This document was generated with Documenter.jl version 1.1.2 on Sunday 26 November 2023. Using Julia version 1.9.4.
Appendix
ToQUBO.jl's Assets
Logo
The ideia behind ToQUBO.jl's logo comes from a wordplay in Portuguese and Spanish. The package's main purpose is to assemble QUBO Models, which sounds like cubo[1], the translation for cube.
Colors
The colors were picked according to Julia's Reference for logo graphics[2]. Text color matches the innermost shape and renders fairly well in both light and dark background themes.
Typography
The MADETYPE Sunflower[3] font was chosen. It was converted to a SVG path using the Google Font to Svg Path[4] online tool.
Web Icon 
The icon used to decorate the documentation resembles an assembled version of the cube with its blue face making up the background.
- 1IPA: [ˈkubʊ]
- 2github.com/JuliaLang/julia-logo-graphics
- 3Licensed by the authors for use in this project
- 4danmarshall.github.io/google-font-to-svg-path
Settings
This document was generated with Documenter.jl version 1.1.2 on Sunday 26 November 2023. Using Julia version 1.9.4.
Prime Factorization
A central problem in Number Theory and cryptography is to factor $R \in \mathbb{N}$, which is known to be the product of two distinct prime numbers $p, q \in \mathbb{N}$. Shor's Algorithm, designed to address such task is often regarded as one of the major theoretical landmarks in Quantum Computing, being responsible for driving increasingly greater interest to the area.
A naïve approach to model this problem can be stated as a quadratically-constrained integer program:
\[\begin{array}{rl} +\text{s.t.} & p \times q = R \\ + & p, q \ge 0 \\ + & p, q \in \mathbb{Z} +\end{array}\]
From the definition and the basics of number theory, we are able to retrieve a few assumptions about the problem's variables:
- $p$ and $q$ are bounded to the interval $\left[1, R\right]$
- Moreover, it is fine to assume $1 < p \le \sqrt{R} \le q \le R \div 2$.
using JuMP
+using ToQUBO
+using DWave
+
+function factor(R::Integer; optimizer = DWave.Neal.Optimizer)
+ return factor(identity, R; optimizer)
+end
+
+function factor(config!::Function, R::Integer; optimizer = DWave.Neal.Optimizer)
+ model = Model(() -> ToQUBO.Optimizer(optimizer))
+
+ @variable(model, 1 <= p <= √R, Int)
+ @variable(model, √R <= q <= R ÷ 2, Int)
+
+ @constraint(model, p * q == R)
+
+ config!(model)
+
+ optimize!(model)
+
+ p = round(Int, value(model[:p]))
+ q = round(Int, value(model[:q]))
+
+ return (p, q)
+endfactor (generic function with 2 methods)p, q = factor(15) do model
+ set_optimizer_attribute(model, "num_reads", 1_000)
+ set_optimizer_attribute(model, "num_sweeps", 2_000)
+end
+
+print("$p, $q")3, 5Settings
This document was generated with Documenter.jl version 1.1.2 on Sunday 26 November 2023. Using Julia version 1.9.4.
Knapsack
We start with some instances of the discrete Knapsack Problem whose standard formulation is
\[\begin{array}{r l} + \max & \mathbf{c}\, \mathbf{x} \\ + \text{s.t.} & \mathbf{w}\, \mathbf{x} \le C \\ + ~ & \mathbf{x} \in \mathbb{B}^{n} +\end{array}\]
First, consider the following items
| Item ($i$) | Value ($c_{i}$) | Weight ($w_{i}$) |
|---|---|---|
| 1 | 1 | 0.3 |
| 2 | 2 | 0.5 |
| 3 | 3 | 1.0 |
to be carried in a knapsack with capacity $C = 1.6$.
Writing down the data above as a linear program, we have
\[\begin{array}{r l} + \max & x_{1} + 2 x_{2} + 3 x_{3} \\ + \text{s.t.} & 0.3 x_{1} + 0.5 x_{2} + x_{3} \le 1.6 \\ + ~ & \mathbf{x} \in \mathbb{B}^{3} +\end{array}\]
Simple JuMP Model
Writing this in JuMP we end up with
using JuMP
+using ToQUBO
+using DWave
+
+model = Model(() -> ToQUBO.Optimizer(DWave.Neal.Optimizer))
+
+@variable(model, x[1:3], Bin)
+@objective(model, Max, x[1] + 2 * x[2] + 3 * x[3])
+@constraint(model, 0.3 * x[1] + 0.5 * x[2] + x[3] ≤ 1.6)
+
+optimize!(model)
+
+solution_summary(model)* Solver : Virtual QUBO Model
+
+* Status
+ Result count : 23
+ Termination status : LOCALLY_SOLVED
+ Message from the solver:
+ ""
+
+* Candidate solution (result #1)
+ Primal status : FEASIBLE_POINT
+ Dual status : NO_SOLUTION
+ Objective value : 5.00000e+00
+
+* Work counters
+ Solve time (sec) : 1.99027e-01
+The final decision is to take items $2$ and $3$, i.e., $x_{1} = 0, x_{2} = 1, x_{3} = 1$.
value.(x)3-element Vector{Float64}:
+ 0.0
+ 1.0
+ 1.0Using DataFrames
Now, lets fill a few more knapsacks. First, we generate uniform random costs $\mathbf{c}$ and weights $\mathbf{w}$ then set the knapsack's capacity $C$ to be a fraction of the total available weight i.e. $80\%$.
This example was inspired by D-Wave's knapsack example repository.
using CSV
+using DataFrames
+
+df = CSV.read("knapsack.csv", DataFrame)8×2 DataFrame
| Row | cost | weight |
|---|---|---|
| Int64 | Int64 | |
| 1 | 66 | 79 |
| 2 | 85 | 20 |
| 3 | 12 | 40 |
| 4 | 3 | 88 |
| 5 | 7 | 34 |
| 6 | 98 | 66 |
| 7 | 69 | 29 |
| 8 | 7 | 97 |
using JuMP
+using ToQUBO
+using DWave
+
+model = Model(() -> ToQUBO.Optimizer(DWave.Neal.Optimizer))
+
+n = size(df, 1)
+c = collect(Float64, df[!, :cost])
+w = collect(Float64, df[!, :weight])
+C = round(0.8 * sum(w))
+
+@variable(model, x[1:n], Bin)
+@objective(model, Max, c' * x)
+@constraint(model, w' * x <= C)
+
+optimize!(model)
+
+# Add Results as a new column
+df[:,:select] = map(xi -> ifelse(xi > 0, "✅", "❌"), value.(x))
+
+df8×3 DataFrame
| Row | cost | weight | select |
|---|---|---|---|
| Int64 | Int64 | String | |
| 1 | 66 | 79 | ✅ |
| 2 | 85 | 20 | ✅ |
| 3 | 12 | 40 | ✅ |
| 4 | 3 | 88 | ❌ |
| 5 | 7 | 34 | ❌ |
| 6 | 98 | 66 | ✅ |
| 7 | 69 | 29 | ✅ |
| 8 | 7 | 97 | ✅ |
Settings
This document was generated with Documenter.jl version 1.1.2 on Sunday 26 November 2023. Using Julia version 1.9.4.
Portfolio Optimization
In this example, we will be exploring an optimization model for asset distribution where the expected return is maximized while mitigating the financial risk. The following approach was inspired by a JuMP tutorial, where monthly stock prices for three assets are provided, namely IBM, WMT and SEHI.
The modelling presented below aggregates the risk measurement $\mathbf{x}' \Sigma \mathbf{x}$ as a penalty term to the objective function, thus yielding
\[\begin{array}{rll} + \max_{\mathbf{x}} & \mathbf{\mu}'\mathbf{x} - \lambda\, \mathbf{x}' \Sigma \mathbf{x} \\ + \textrm{s.t.} & 0 \le {x}_{i} \le 1 & \forall i \\ + & \sum_{i} {x}_{i} = 1 +\end{array}\]
where $\mu_{i} = \mathbb{E}[r_{i}]$ is the expected return value for each investment $i$; $\Sigma$ is the covariance matrix and $\lambda$ is the risk-aversion penalty factor.
Stock prices
using DataFrames
+using Statistics
+
+assets = [:IBM, :WMT, :SEHI]
+
+df = DataFrames.DataFrame(
+ [
+ 93.043 51.826 1.063
+ 84.585 52.823 0.938
+ 111.453 56.477 1.000
+ 99.525 49.805 0.938
+ 95.819 50.287 1.438
+ 114.708 51.521 1.700
+ 111.515 51.531 2.540
+ 113.211 48.664 2.390
+ 104.942 55.744 3.120
+ 99.827 47.916 2.980
+ 91.607 49.438 1.900
+ 107.937 51.336 1.750
+ 115.590 55.081 1.800
+ ],
+ assets,
+)13×3 DataFrame
| Row | IBM | WMT | SEHI |
|---|---|---|---|
| Float64 | Float64 | Float64 | |
| 1 | 93.043 | 51.826 | 1.063 |
| 2 | 84.585 | 52.823 | 0.938 |
| 3 | 111.453 | 56.477 | 1.0 |
| 4 | 99.525 | 49.805 | 0.938 |
| 5 | 95.819 | 50.287 | 1.438 |
| 6 | 114.708 | 51.521 | 1.7 |
| 7 | 111.515 | 51.531 | 2.54 |
| 8 | 113.211 | 48.664 | 2.39 |
| 9 | 104.942 | 55.744 | 3.12 |
| 10 | 99.827 | 47.916 | 2.98 |
| 11 | 91.607 | 49.438 | 1.9 |
| 12 | 107.937 | 51.336 | 1.75 |
| 13 | 115.59 | 55.081 | 1.8 |
Solving
using JuMP
+using ToQUBO
+using DWave
+
+function solve(
+ config!::Function,
+ df::DataFrame,
+ λ::Float64 = 10.;
+ optimizer = DWave.Neal.Optimizer
+)
+ # Number of assets
+ n = size(df, 2)
+
+ # Relative monthly return
+ r = diff(Matrix(df); dims = 1) ./ Matrix(df[1:end-1, :])
+
+ # Expected monthly return value for each stock
+ μ = vec(Statistics.mean(r; dims = 1))
+
+ # Covariance matrix
+ Σ = Statistics.cov(r)
+
+ # Build model
+ model = Model(() -> ToQUBO.Optimizer(optimizer))
+
+ @variable(model, 0 <= x[1:n] <= 1)
+ @objective(model, Max, μ'x - λ * x' * Σ * x)
+ @constraint(model, sum(x) == 1)
+
+ config!(model)
+
+ optimize!(model)
+
+ return value.(x)
+end
+
+function solve(df::DataFrame, λ::Float64 = 10.; optimizer = DWave.Neal.Optimizer)
+ return solve(identity, df, λ; optimizer)
+endsolve (generic function with 4 methods)solve(df) do model
+ JuMP.set_silent(model)
+ JuMP.set_optimizer_attribute(model, "num_reads", 2_000)
+end3-element Vector{Float64}:
+ 0.0
+ 1.0
+ 0.0Penalty Analysis
To finish our discussion, we are going to sketch some graphics to help our reasoning on how the penalty factor $\lambda$ affects our investments.
using Plots
+
+Λ = collect(0.:5.:50.)
+X = Dict{Symbol,Vector{Float64}}(tag => [] for tag in assets)
+
+for λ = Λ
+ x = solve(df, λ)
+
+ for (i, tag) in enumerate(assets)
+ push!(X[tag], x[i])
+ end
+end
+
+plt = plot(;
+ title="Portfolio Optimization",
+ xlabel=raw"penalty factor ($\lambda$)",
+ ylabel=raw"investment share ($x$)",
+)
+
+for tag in assets
+ plot!(plt, Λ, X[tag]; label=string(tag))
+end
+
+plt
Settings
This document was generated with Documenter.jl version 1.1.2 on Sunday 26 November 2023. Using Julia version 1.9.4.
ToQUBO.jl Documentation
ToQUBO.jl is a Julia Package intended to automatically translate models written in JuMP, into the QUBO mathematical optimization framework.
Quick Start
Installation
julia> import Pkg
+
+julia> Pkg.add("ToQUBO")Example
using JuMP
+using ToQUBO
+using DWave
+
+model = Model(() -> ToQUBO.Optimizer(DWave.Neal.Optimizer))
+
+@variable(model, x[1:3], Bin)
+@objective(model, Max, 1.0*x[1] + 2.0*x[2] + 3.0*x[3])
+@constraint(model, 0.3*x[1] + 0.5*x[2] + 1.0*x[3] <= 1.6)
+
+optimize!(model)
+
+solution_summary(model)* Solver : Virtual QUBO Model
+
+* Status
+ Result count : 23
+ Termination status : LOCALLY_SOLVED
+ Message from the solver:
+ ""
+
+* Candidate solution (result #1)
+ Primal status : FEASIBLE_POINT
+ Dual status : NO_SOLUTION
+ Objective value : 5.00000e+00
+
+* Work counters
+ Solve time (sec) : 1.98995e-01
+Citing ToQUBO.jl
If you use ToQUBO.jl in your work, we kindly ask you to include the following citation:
@software{toqubo:2023,
+ author = {Pedro Maciel Xavier and Pedro Ripper and Tiago Andrade and Joaquim Dias Garcia and David E. Bernal Neira},
+ title = {{ToQUBO.jl}},
+ month = {feb},
+ year = {2023},
+ publisher = {Zenodo},
+ version = {v0.1.5},
+ doi = {10.5281/zenodo.7644291},
+ url = {https://doi.org/10.5281/zenodo.7644291}
+}Settings
This document was generated with Documenter.jl version 1.1.2 on Sunday 26 November 2023. Using Julia version 1.9.4.
Manual
Quick Start Guide
using JuMP
+using ToQUBO
+using DWave
+
+model = Model(() -> ToQUBO.Optimizer(DWave.Neal.Optimizer))
+
+@variable(model, x[1:3], Bin)
+
+@objective(model, Max, 1.0 * x[1] + 2.0 * x[2] + 3.0 * x[3])
+
+@constraint(model, 0.3 * x[1] + 0.5 * x[2] + 1.0 * x[3] <= 1.6)
+
+optimize!(model)
+
+solution_summary(model)* Solver : Virtual QUBO Model
+
+* Status
+ Result count : 22
+ Termination status : LOCALLY_SOLVED
+ Message from the solver:
+ ""
+
+* Candidate solution (result #1)
+ Primal status : FEASIBLE_POINT
+ Dual status : NO_SOLUTION
+ Objective value : 5.00000e+00
+
+* Work counters
+ Solve time (sec) : 1.99414e-01
+Table of Contents
Settings
This document was generated with Documenter.jl version 1.1.2 on Sunday 26 November 2023. Using Julia version 1.9.4.
Running a Model
We hope to write this part of the documentation soon. Please come back later!
Settings
This document was generated with Documenter.jl version 1.1.2 on Sunday 26 November 2023. Using Julia version 1.9.4.
Gathering Results
We hope to write this part of the documentation soon. Please come back later!
Settings
This document was generated with Documenter.jl version 1.1.2 on Sunday 26 November 2023. Using Julia version 1.9.4.
Compiler Settings
ToQUBO.Attributes.StableCompilation — TypeStableCompilation()When set, this boolean flag enables stable reformulation methods, thus yielding predictable results.
Compiler Messages
ToQUBO.Attributes.Warnings — TypeWarnings()Compiler Optimization
ToQUBO.Attributes.Optimization — TypeOptimization()Working with target architectures
ToQUBO.Attributes.Architecture — TypeArchitecture()Selects which solver architecture to use. Defaults to QUBOTools.GenericArchitecture.
Quadratization
ToQUBO.Attributes.Quadratize — TypeQuadratize()Boolean flag to conditionally perform the quadratization step. Is automatically set by the compiler when high-order functions are generated.
ToQUBO.Attributes.QuadratizationMethod — TypeQuadratizationMethod()Defines which quadratization method to use. Available options are defined in the PBO submodule.
ToQUBO.Attributes.StableQuadratization — TypeStableQuadratization()When set, this boolean flag enables stable quadratization methods, thus yielding predictable results. This is intended to be used during tests or other situations where deterministic output is desired. On the other hand, usage in production is not recommended since it requires increased memory and processing resources.
Variable & Constraint Encoding
ToQUBO.Attributes.VariableEncodingBits — TypeVariableEncodingBits()ToQUBO.Attributes.DefaultVariableEncodingBits — TypeDefaultVariableEncodingBits()ToQUBO.Attributes.VariableEncodingATol — TypeVariableEncodingATol()ToQUBO.Attributes.DefaultVariableEncodingATol — TypeDefaultVariableEncodingATol()Fallback value for VariableEncodingATol.
ToQUBO.Attributes.VariableEncodingMethod — TypeVariableEncodingMethod()Available methods are:
Encoding.Binary(default)Encoding.UnaryEncoding.ArithmeticEncoding.OneHotEncoding.DomainWallEncoding.Bounded
The Encoding.Binary, Encoding.Unary and Encoding.Arithmetic encodings can have their expansion coefficients bounded by wrapping them with the Encoding.Bounded method.
ToQUBO.Attributes.DefaultVariableEncodingMethod — TypeDefaultVariableEncodingMethod()Fallback value for VariableEncodingMethod.
ToQUBO.Attributes.VariableEncodingPenalty — TypeVariableEncodingPenalty()Allows the user to set and retrieve the coefficients used for encoding variables when additional constraints are involved.
ToQUBO.Attributes.ConstraintEncodingPenalty — TypeConstraintEncodingPenalty()Allows the user to retrieve the coefficients used for encoding constraints.
Discretization
ToQUBO.Attributes.Discretize — TypeDiscretize()When set, this boolean flag guarantees that every coefficient in the final formulation is an integer.
Settings
This document was generated with Documenter.jl version 1.1.2 on Sunday 26 November 2023. Using Julia version 1.9.4.