# The Hyperdrive Protocol: An Automated Market Maker for Fixed and Variable
Rates

Jonny Rhea
<EMAIL_ADDRESS>Alex Towle
<EMAIL_ADDRESS>Mihai Cosma
<EMAIL_ADDRESS>
(April 2024)

###### Abstract

Hyperdrive is a protocol designed to facilitate the trading of fixed and
variable rate assets. The protocol’s unique pricing model consolidates
liquidity into a single pool which addresses the challenges of fragmented
liquidity across terms, eliminates the need for rollovers, and allows terms to
be issued on demand. Its design meaningfully improves trading efficiency,
liquidity provisioning, and user experience over existing fixed and variable
rate protocol models.

##  1 Introduction111This Whitepaper, the underlying code, and related
documents regarding Hyperdrive refer to certain terms for convenience and ease
of reading - such as “bonds,” “shares,” and related terms. These uses are made
for colloquial reference, and do not refer to their official meanings in the
financial industry or in any law, regulation, or rule; they are not intended
to (and do not) have any particular legal or regulatory significance.

At its core, Hyperdrive is an automated market maker (AMM) protocol for fixed
rates and variable rates; however, it has several unique features that set it
apart from existing fixed and variable rate protocols:

  * •

Terms on Demand: Hyperdrive allows for minting to be a part of the AMM, where
the AMM essentially underwrites a new term for the user whenever they open a
position. The user is not constrained to purchasing, selling, or minting into
preexisting terms that are partially matured.

  * •

Continuous Liquidity: Hyperdrive pools never expire and underwrite a variety
of fixed and variable rate terms with differing maturity dates. LPs can
provide liquidity once without needing to roll their liquidity over to new
terms.

  * •

Single-Sided Liquidity: Hyperdrive’s liquidity provision mechanism calls for a
single asset, improving both the UX and capital efficiency of LPing.

A Hyperdrive AMM allows users to open long and short positions to get exposure
to fixed and variable rates, respectively. These rates are generated from
trading activities and its base token, which is an arbitrary yield-bearing
asset underlying the AMM. There are three position types in Hyperdrive:

  1. 1.

Long: A long position represents some amount of fixed rate exposure to the
market. A user can purchase this with a base token and receive a bond token
that will mature to the face value. Long traders receive fixed interest during
the term and surrender any variable interest earned on their investment.

  2. 2.

Short: When a user opens a short, they are selling bonds to the Hyperdrive
pool - this increases the total supply of bonds available to purchase. Short
traders sell long positions, which means they pay fixed interest during the
term and receive all the variable interest collected on the LP’s investment.

  3. 3.

Provide Liquidity: LPs supply liquidity to a single pool, which underwrites a
variety of fixed terms with differing maturity dates, giving LPs continuous
exposure to the fixed rate market and underlying yield source. LPs, by
definition, always take the opposite side of the trade opened by users.

##  2 Preliminaries

A Hyperdrive AMM calculates a market rate to provide users access to a fixed
rate and variable rate. In typical AMMs, LPs deposit pairs of tokens into
liquidity pools. Each pool consists of two different tokens–for example, ETH
and DAI–and a function is then used to price the two assets in terms of each
other for users to trade. One well known function for this, called the
Constant Product Function [1], is defined as:

$\displaystyle x\cdot y=k$

where the pair $(x,y)$ represents the pool reserves and $k$ is some fixed
quantity that is preserved. The following shows how the function preserves $k$
when a user trades $\Delta x$ for $\Delta y$ with the pool:

$\displaystyle(x+\Delta x)\cdot(y-\Delta y)=k$

Hyperdrive works similarly, except that instead of facilitating trades between
ETH/DAI, it allows a user to buy and sell bonds that are backed by an
underlying yield source.

###  2.1 Yield Source

A yield source allows a user to deposit a base token (x) and accrue positive
yield over time. The rate of return of base tokens deposited into a yield
source is not fixed and in many cases changes from block to block. This is
referred to as the “variable interest rate” of a yield source protocol that
holds deposited assets.

Hyperdrive AMMs utilize Yield Sources to back the value of the bonds that are
bought from and sold to it. When an LP provides liquidity to a Hyperdrive AMM,
the assets are immediately deposited into a protocol that accrues variable
rate interest. To support different yield sources without making large changes
to the protocol, Hyperdrive implementations abstract away specifics of how
yield source protocols accrue variable interest by introducing the concept of
shares.

###  2.2 Shares

The balance of assets held in a Hyperdrive pool is tracked in terms of
“shares,” rather than base token amounts. Shares and base token amounts (i.e.
“base”) are related to each other in terms of a “share price,” denoted as $c$,
that changes over time. Given an amount of shares $z$, the corresponding
amount of base $x$ can be calculated as:

$x=c\cdot z$

This share price abstraction provides a consistent way to measure interest
accrual. Suppose that a Hyperdrive AMM’s share reserves $z$ remain constant
from time $t_{0}$ to time $t_{1}$. If $c_{0}$ is the share price at time
$t_{0}$ and $c_{1}$ is the share price at time $t_{1}$, then the base reserves
at $t_{0}$ are $x_{0}=c_{0}\cdot z$ and the base reserves at $t_{1}$ are
$x_{1}=c_{1}\cdot z$. The amount of interest accrued between $t_{0}$ and
$t_{1}$ is simply the difference in the base reserves
$x_{1}-x_{0}=(c_{1}-c_{0})\cdot z$. Note that for the remainder of this paper
the variable, $c$, will represent the current “share price.”

##  3 Pricing Model

Bonds in Hyperdrive are minted and given a maturity date upon purchase. This
means the protocol must be able to price bonds with a variety of maturity
dates from a single liquidity pool. The intuition behind how Hyperdrive
accomplishes this is to re-imagine how a set of minted bonds matures.

Typical intuition would imagine all minted bonds maturing simultaneously:

Figure 1: Bonds maturing simultaneously

In Hyperdrive, bonds mature sequentially:

Figure 2: Bonds maturing sequentially

This allows the AMM to price new bonds that haven’t matured at the market
value and matured bonds at their face value.

###  3.1 Buying and Selling Bonds

Hyperdrive trading is fundamentally comprised of two actions, buying and
selling bonds. The order in which this is done defines the type of position
the user is purchasing and how the reserves are impacted. Since bonds are
minted on demand, they don’t exist until a user opens a position that either
buys them from the pool or sells them to the pool. Bonds in a Hyperdrive pool
can be thought of as a virtual quantity representing how many bonds can be
minted at the pool’s current rate given the amount of share liquidity in the
pool.

In this section (3.1), we will define some notation before addressing the
impact that buying and selling bonds has on the reserves. Since the price of a
new bond is driven by market demand, they are priced as a function of the
share and bond reserve levels (z,y), respectively, using a trading invariant.
We will refer to this trading invariant function as $I$.

$\displaystyle\Delta out=I_{out}^{in}\left(\Delta in,z,y\right)$
$\displaystyle(1)$

For the purposes of this discussion, the only requirement for how $I$ is
defined is that it preserves path-independence. In other words, the end state
of the reserves discussed will always be the same, regardless of the order in
which trades are executed. As mentioned in the previous section, matured bonds
are priced at their face value. It is important to note that, assuming the
pool is solvent, this pricing is also path-independent. We can define this
pricing function as:

$\displaystyle M\left(\Delta y\right)=\frac{\Delta y}{c}=\Delta z$
$\displaystyle(2)$

Combining (1) and (2) together, we define a function $H$ that can calculate a
position’s impact on $z$ given $\Delta y$ at time remaining $t_{r}$. We use
the same “in” and “out” notation defined above for $I$. $H$ is defined as
follows:

$\displaystyle H_{\Delta z}^{\Delta y}\left(\Delta
y,z,y,t_{r}\right)=I_{\Delta z}^{\Delta y}\left(\Delta y\cdot
t_{r},z,y\right)+M\left(\Delta y\cdot(1-t_{r})\right)$ $\displaystyle(3)$

####  3.1.1 Opening Positions

When opening a position, the pricing model uses the trading function $H$ to
calculate the impact on the reserves at $t_{r}=1$ and we always solve for the
$\Delta z$ that keeps the pool’s depth, $k$, constant.

When a user opens a long, they are purchasing fixed-rate exposure at the
current market rate. We can represent the impact that a purchase of $\Delta y$
bonds with $\Delta z$ shares has on the reserves as:

$\displaystyle\Bigl{(}z+H_{\Delta z}^{\Delta y}\left(\Delta
y,z,y,1\right)\,,\,y-\Delta y\Bigr{)}$

When a user opens a short, they are purchasing exposure to the underlying
yield source’s variable rate. In this case, the user is owed the amount of
variable interest accrued over the course of the bond’s term, meaning funds
must be removed from the reserves to accrue interest without backing
additional open positions. We can represent the sale of $\Delta y$ bonds for
$\Delta z$ shares on the reserves as:

$\displaystyle\Bigl{(}z-H_{\Delta z}^{\Delta y}\left(\Delta
y,z,y,1\right)\,,\,y+\Delta y\Bigr{)}$

####  3.1.2 Closing Positions

When closing a position, the pricing model uses the time remaining, $t_{r}$,
to determine the proportion of bonds that Hyperdrive considers new versus
matured. As a result, when closing a position, we use $I$ to price the
proportion of bonds we consider new and $M$ to price the matured ones. Since
matured bonds can just be redeemed for the face value of the bond in terms of
the base token, we know that at maturity $\Delta x=\Delta y$.

When a user closes a long, they are redeeming the $\Delta y\cdot(1-t_{r})$
matured bonds and selling $\Delta y\cdot t_{r}$ new bonds back to Hyperdrive.
This means that the reserves will be updated as follows:

$\displaystyle\Bigl{(}z-H_{\Delta z}^{\Delta y}\left(\Delta
y,z,y,t_{r}\right)\,,\,y+\Delta y\cdot t_{r}\Bigr{)}$

When a user closes a short, they are redeeming the $\Delta y\cdot(1-t_{r})$
matured bonds and purchasing the $\Delta y\cdot t_{r}$ new bonds via
Hyperdrive. Again, since the trader is buying the new bonds from Hyperdrive,
the reserves will be updated as follows:

$\displaystyle\Bigl{(}z+H_{\Delta z}^{\Delta y}\left(\Delta
y,z,y,t_{r}\right)\,,\,y-\Delta y\cdot t_{r}\Bigr{)}$

Unfortunately, these accounting updates are problematic because the impact of
$M$ on the share reserves changes the pool’s depth; this allows a trader to
open a position on one curve and close it on another by sandwiching a closing
trade. As Peter Zeitz’s Curve Vulnerability Report [3] demonstrated, it is not
safe for an AMM to support several different trading curves simultaneously as
it allows malicious traders to extract large amounts of value from the system.

A modification to the accounting scheme is needed to hold the pool’s depth
constant before and after closing matured positions. In the next section, we
will define such a scheme by introducing the zeta adjustment or $\zeta$.

###  3.2 Zeta Adjustment

Simply put, the zeta adjustment replaces the reserve tuple $(z,y)$ used when
pricing new bonds and replaces it with $(z_{e},y)$ where $z_{e}$ is defined
as:

$\displaystyle z_{e}=z-\zeta$

$\zeta$ is introduced as a new state variable to offset any changes to $z$
caused by $M$. When a long is closed at time $t_{r}$, $z_{e}$ is calculated as
follows:

$\displaystyle z_{e}$ $\displaystyle=\biggl{(}z-H_{\Delta z}^{\Delta
y}\left(\Delta y,z,y,t_{r}\right)\biggr{)}-\biggl{(}\zeta-M\bigl{(}\Delta
y\cdot\left(1-t_{r}\right)\bigr{)}\biggr{)}$ $\displaystyle\implies z_{e}$
$\displaystyle=\left(z-\zeta\right)-I_{\Delta z}^{\Delta y}\left(\Delta y\cdot
t_{r},z,y\right)$

Similarly, when a short is closed at time $t_{r}$, $z_{e}$ is calculated as
follows:

$\displaystyle z_{e}$ $\displaystyle=\biggl{(}z-H_{\Delta z}^{\Delta
y}\left(\Delta y,z,y,t_{r}\right)\biggr{)}-\biggl{(}\zeta+M\bigl{(}\Delta
y\cdot\left(1-t_{r}\right)\bigr{)}\biggr{)}$ $\displaystyle\implies z_{e}$
$\displaystyle=\left(z-\zeta\right)+I_{\Delta z}^{\Delta y}\left(\Delta y\cdot
t_{r},z,y\right)$

As a result, when a position is closed at maturity, $z_{e}$ remains constant
because the update to $z$ by $M$ is canceled out by the change made to
$z_{e}$. Using the tuple $(z_{e},y)$, rather than $(z,y)$, to parameterize the
pricing model for new bonds will hold the spot price and depth constant.

The following plot illustrates the impact of the Zeta Adjustment:

$0$$0.2$$0.4$$0.6$$0.8$$0$$0.2$$0.4$$0.6$$0.8$$z$$y$Zeta Adjustment:
$(z,y)$$\zeta=0$$\zeta<0$$\zeta>0$

When zeta is adjusted for a short close, it translates the curve to the left
(red). When zeta is adjusted for a long close, it translates the curve to the
right (blue). As mentioned previously, for the Zeta Adjustment to cancel out
the impact redemptions of matured bonds have on the pricing model, we must use
$z_{e}$ in place of $z$ when calculating trades on $I$. The following plot
illustrates this:

$0$$0.2$$0.4$$0.6$$0.8$$0$$0.2$$0.4$$0.6$$0.8$$z_{e}$$y$Zeta Adjustment:
$(z_{e},y)$$\zeta=0$$\zeta<0$$\zeta>0$

As seen above, all scenarios look identical from the $z_{e}$ axis. Since $I$
is assumed to be path-independent, $M$ is trivially path-independent and has
no impact on $z_{e}$; this implies that the Hyperdrive pricing model as a
whole is path-independent.

##  4 Position Accounting

Since Hyperdrive deposits the capital that underlies bonds into a yield
source, interest accrues during the lifetimes of longs and shorts. In exchange
for the fixed interest that longs receive as their bonds mature, they give up
the variable interest that accrues over time. In contrast, shorts receive the
variable interest that accrues over time to compensate them for paying fixed
interest to the pool.

Fees are charged when traders open and close positions. LPs receive a fee
proportional to the pool’s spot rate when newly-minted bonds are traded that
is parameterized by $\phi_{n}$, $0\leq\phi_{n}\leq 1$. Additionally, LPs
receive a fee proportional to position size when matured bonds are traded that
is parameterized by $\phi_{m}$, $0\leq\phi_{m}\leq 1$. Governance may collect
a portion of trading fees parameterized by $\phi_{g}$, $0\leq\phi_{g}\leq 1$.
For convenience, we will define the following functions that represent the fee
calculations:

$\displaystyle f_{new}(\Delta y,t_{r})$
$\displaystyle=\phi_{n}\cdot\left(1-p\right)\cdot\Delta y\cdot t_{r}$
$\displaystyle f_{mature}(\Delta y,t_{r})$ $\displaystyle=\phi_{m}\cdot\Delta
y\cdot\left(1-t_{r}\right)$ $\displaystyle f(\Delta y,t_{r})$
$\displaystyle=f_{new}(\Delta y,t_{r})+f_{mature}(\Delta y,t_{r})$

These will be referenced as needed in the following sections.

###  4.1 Longs

Traders can open long positions to receive fixed-rate exposure. From the
perspective of the pricing model, opening a long position is equivalent to
purchasing bonds. This has the consequence of pushing the spot rate down and,
conversely, pushing the spot price up towards the face value of the bond.
Assuming the pool has sufficient liquidity, traders can close their existing
long positions at any time before maturity, and regardless of the pool’s idle
liquidity (see Section 5 for more details on idle liquidity), they will always
be able to close their long positions at maturity.

####  4.1.1 Opening Longs

To open a long position, a trader provides $\Delta x$ base to a Hyperdrive AMM
and receives a long position with a face value of $\Delta y$. The long’s face
value is equivalent to the amount of bonds underlying the long and the amount
of base proceeds that can be redeemed when closing the bonds at maturity.
Since longs are opened with a full term remaining before maturity, the size of
long positions for newly minted bonds is priced using the invariant as
follows:

$\Delta y=H_{\Delta y}^{\Delta z}\left(\tfrac{\Delta
x}{c},z,y,1\right)-\phi_{n}\cdot\left(p^{-1}-1\right)\cdot\Delta x$

Since bonds are purchased when a trader opens a long, the share reserves are
increased and the bond reserves are decreased. The purchased bonds are all
newly minted, so the zeta adjustment isn’t changed when updating the reserves.
With all of this in mind, the pool’s state is updated as follows:

$\displaystyle z_{1}$ $\displaystyle=z_{0}+\frac{\Delta
x-\phi_{g}\cdot\phi_{n}\cdot\left(1-p\right)\cdot\Delta x}{c}$ $\displaystyle
y_{1}$ $\displaystyle=y_{0}-\Delta y$ $\displaystyle\zeta_{1}$
$\displaystyle=\zeta_{0}$

####  4.1.2 Closing Longs

When a trader closes a long position, they sell the bonds that underlie the
position back to the pool and receive the base proceeds that were released
from the sale. If the position is closed at or after maturity, the trader
receives the face value of the bonds, and the pool’s spot rate remains
constant. If the position is closed before maturity, some of the bonds
underlying the long are newly minted, and the pool’s spot rate will increase
as a result of the sale.

The base proceeds $\Delta x$ derived from closing a long position of size
$\Delta y$ depend on the time remaining until maturity of $t_{r}$, $0\leq
t_{r}\leq 1$. We decompose the proceeds $\Delta x$ into the proceeds of
selling newly minted bonds $\Delta x_{0}$ and the proceeds of selling matured
bonds $\Delta x_{1}$. We compute the base proceeds $\Delta x$ as follows:

$\displaystyle\Delta x$ $\displaystyle=c\cdot H_{\Delta z}^{\Delta
y}\left(\Delta y,z,y,t_{r}\right)-f(\Delta y,t_{r})$

Since the trader sells bonds when they close a long, the share reserves are
decreased by the amount of base proceeds the trader receives. The bond
reserves are increased by the amount of newly minted bonds that were sold. To
ensure that the pricing model is only affected by the sale of newly minted
bonds, the zeta adjustment is decreased by the proceeds of the sale of matured
bonds. The state updates resulting from closing the long are given by:

$\displaystyle z_{1}$ $\displaystyle=z_{0}-\frac{\Delta x+\phi_{g}\cdot
f(\Delta y,t_{r})}{c}$ $\displaystyle y_{1}$ $\displaystyle=y_{0}+\Delta y$
$\displaystyle\zeta_{1}$ $\displaystyle=\zeta_{0}-M(\Delta
y\cdot(1-t_{r}))+\left(1-\phi_{g}\right)\cdot f_{mature}(\Delta y,t_{r})$

###  4.2 Shorts

With respect to pricing, short positions are the inverse of long positions.
They receive variable-rate exposure in exchange for fixed-rate exposure. From
the perspective of the pricing model, opening a short is equivalent to selling
bonds to the LPs. This has the consequence of decreasing the spot rate and
increasing the spot price, which is the opposite of what happens when longs
are opened. For the system to remain solvent, shorts must deposit the maximum
amount of money that they can lose when the short is opened, which is equal to
the fixed interest that will be owed to LPs at maturity.

When a trader closes their short, they purchase back the partially (or fully)
matured bond that underlies the short position at its current market price. By
locking up only the fixed interest, they’ll earn the current market discount
plus variable interest accrued by the full face value of the bonds. This gives
short positions multiplied exposure to the variable rate.

####  4.2.1 Opening Shorts

To open a short position, a trader specifies the amount of bonds to short
$\Delta y$. The short is responsible for depositing the difference between the
value of the bonds plus any interest that has accrued since the start of the
checkpoint (see section 4.3 for more details on checkpoints),
$\tfrac{c}{c_{0}}\cdot\Delta y$, and the current price of the bonds $\Delta
x$. The short deposit is calculated as:

$\Delta x=\tfrac{c}{c_{0}}\cdot\Delta y-c\cdot H_{\Delta z}^{\Delta y}(\Delta
y,z,y,1)+f_{new}(\Delta y,1)$

Since bonds are sold when a trader opens a short, the share reserves are
decreased by the proceeds from selling the short, and the bond reserves are
increased by the amount of bonds that were short-sold. The bonds that are
short-sold are newly minted, so the zeta adjustment remains constant. With
this in mind, the pool’s state is updated as follows:

$\displaystyle z_{1}$ $\displaystyle=z_{0}-H_{\Delta z}^{\Delta y}(\Delta
y,z,y,1)+f_{new}(\Delta y,1)\cdot(1-\phi_{g})$ $\displaystyle y_{1}$
$\displaystyle=y_{0}+\Delta y$ $\displaystyle\zeta_{1}$
$\displaystyle=\zeta_{0}$

####  4.2.2 Closing Shorts

When a trader closes a short position, they purchase the bonds that were
short-sold when the short was opened. The proceeds that they receive are the
difference between the face value of the underlying bonds $\Delta y$ and the
cost to purchase the bonds plus the variable interest accrued since the short
was opened. Let $c_{0}$ be the share price when the short was opened, and
$c_{1}$ be equal to the current share price, $c$, if the shorts haven’t
matured; otherwise, $c_{1}$ is the share price at the time the shorts were
closed. We can calculate the shorts trading proceeds as:

$\Delta x=\frac{c_{1}}{c_{0}}\cdot\Delta y-c\cdot H_{\Delta z}^{\Delta
y}(\Delta y,z,y,t_{r})-f(\Delta y,t_{r})$

Since the bonds that were short-sold are purchased when the short is closed,
the share reserves are increased by the cost of purchasing the bond, and the
bond reserves are decreased by the amount of bonds that were purchased. To
ensure the pricing model is only impacted by the sale of newly minted bonds,
the zeta adjustment is increased by the cost of the purchase of matured bonds.
The state updates resulting from closing the short are given by:

$\displaystyle z_{1}$ $\displaystyle=z_{0}+H_{\Delta z}^{\Delta y}\left(\Delta
y,z,y,t_{r}\right)+\frac{f(\Delta y,t_{r})}{c}\cdot\left(1-\phi_{g}\right)$
$\displaystyle y_{1}$ $\displaystyle=y_{0}-\Delta y$ $\displaystyle\zeta_{1}$
$\displaystyle=\zeta_{0}+M(\Delta y\cdot(1-t_{r}))+\frac{f_{mature}(\Delta
y,t_{r})}{c}\cdot(1-\phi_{g})$

###  4.3 Checkpoints

Closing positions at maturity automatically increases the system’s capital
efficiency and ensures that positions aren’t open for longer than the
predetermined position duration. While closing positions individually would be
theoretically possible, the costs would be prohibitive. With this in mind,
positions are grouped into checkpoints. Each Hyperdrive pool is parameterized
with a checkpoint duration $d_{c}$ that determines the size of these
checkpoints. The checkpoint duration evenly divides the position duration.
Given the current block time $t$, the start of the latest checkpoint $t_{c}$
is calculated as:

$t_{c}=t-(t\bmod d_{c})$

When positions are opened, they are backdated to the start of the latest
checkpoint. This ensures that all positions opened during a checkpoint mature
simultaneously. The first trade or LP operation in a checkpoint mints the
checkpoint by storing the current share price as the checkpoint’s share price,
which is used to calculate the proceeds of shorts opened in the checkpoint.
Aside from storing the checkpoint’s share price, minting the checkpoint also
applies state updates resulting from closing matured positions in the
checkpoint to the pool’s reserves. Since the closed positions have already
matured, the share reserves and zeta adjustment are updated by the amount of
longs that have matured $y_{l}$ and amount of shorts that have matured $y_{s}$
as:

$\displaystyle z$ $\displaystyle\mathrel{+}=y_{l}-y_{s}$ $\displaystyle\zeta$
$\displaystyle\mathrel{+}=y_{l}-y_{s}$

Traders could hypothetically wait until the end of the checkpoint to open long
positions, allowing them to accrue fixed interest as soon as the checkpoint
advances. To mitigate this arbitrage from negatively impacting LPs, the
checkpoint duration should always be short enough that the fixed interest
accrued is less than the cost of a trade. LPs should verify that the
checkpoint duration is sufficiently small before interacting with a Hyperdrive
pool.

###  4.4 Solvency

Since longs and shorts are entitled to fixed-rate and variable-rate exposure
over the duration of the term, some of a Hyperdrive pool’s capital must be
reserved to ensure the system is able to honor this exposure at maturity. When
shorts are opened, all of the base required to capitalize the shorts is
reserved outside of the share reserves, which ensures that the capital
underlying shorts can’t be consumed by trading or LP actions. In contrast, the
capital that underlies longs is initially held in the share reserves and can
be consumed by trading and LP actions. This means that the system must enforce
a solvency constraint on longs to prevent trading and LPs from removing the
capital reserved for longs.

At maturity, longs receive the face value of their position, and shorts
receive the interest generated on the face value of their position. Since
longs receive the principal and shorts receive the interest, one base asset is
sufficient to capitalize both a long and a short. This means that longs and
shorts can net out. Since shorts already reserve capital equal to their face
value, matching longs and shorts opened in the same checkpoint can use the
same underlying collateral. The solvency constraint $e_{c}$ required to ensure
that positions in checkpoint $c$ are solvent is calculated using the number of
longs in checkpoint $y_{l}$ and shorts in checkpoint $y_{s}$:

$e_{c}=\max(y_{l}-y_{s},0)$

Since longs and shorts in different checkpoints mature at different times,
these positions cannot be matched across checkpoints without compromising the
fixed term of longs and shorts. So while positions within a single checkpoint
can net out, the long exposure across checkpoints must be added up to
determine the true solvency requirement. Assuming the number of checkpoints
per term is $n\geq 1$ and the current checkpoint is $c$, the global solvency
requirement for positions can be calculated as:

$e=\sum_{k=c-(n-1)}^{c}e_{k}$

Since funds for non-netted longs are contained in the share reserves, the
solvency constraint for a Hyperdrive pool can be formulated as:

$z-\frac{e}{c}\geq 0$

To avoid inflation attacks caused by precision errors in fixed-point
arithmetic, the AMM reserves a small amount of capital to ensure its share
reserves never fall below the minimum share reserves, $z_{min}$. For the same
reason, a Hyperdrive pool also burns an amount of LP shares equal to the
minimum share reserves when the pool is initialized. To respect the minimum
share reserves, we can update our solvency constraint to:

$z-\frac{e}{c}\geq z_{min}$

Since longs and shorts can be matched up to offset the solvency constraints
placed on the system, there are situations in which longs and shorts may not
be able to close before maturity. This decision was made because the benefits
of netting positions greatly improve the capital efficiency of the system and
access to liquidity for users who are opening new positions. In economically
rational markets, the risk of being unable to close a position before maturity
should be low.

###  4.5 Zombie Interest

The checkpoint system applies updates to the share and base reserves that
result from closing matured positions. All the proceeds owed to holders of the
matured long and short positions are set aside in the zombie share reserves,
$z_{zombie}$. The interest collected from these matured positions after
maturity is referred to as “zombie interest.” In addition to tracking any
unclaimed proceeds in the share reserves, Hyperdrive also tracks the proceeds
in base with the zombie base reserves quantity, $x_{zombie}$. Whenever long
proceeds, $\Delta x_{long}$, and short proceeds, $\Delta x_{short}$, are set
aside at maturity, $z_{zombie}$ and $x_{zombie}$ are updated in lock-step as
follows:

$\displaystyle z_{zombie,1}$ $\displaystyle=z_{zombie,0}+\frac{\Delta
x_{long}+\Delta x_{short}}{c}$ $\displaystyle x_{zombie,1}$
$\displaystyle=x_{zombie,0}+\Delta x_{long}+\Delta x_{short}$

Long and short positions don’t accrue interest after maturity. Long positions
are opened with a predetermined amount of fixed interest that can be redeemed
in full at maturity. Short positions have paid all of their fixed interest by
maturity, so allowing them to continue to accrue variable interest would be
dangerous for LPs. Whenever a long or short position is closed after maturity,
the system calculates their proceeds, $\Delta x$, identically to how they
would have been calculated at maturity. Instead of updating the share and bond
reserves when the matured position is closed, the system updates the zombie
share and base reserves as follows:

$\displaystyle z_{zombie,1}=z_{zombie,0}-\frac{\Delta x}{c}$ $\displaystyle
x_{zombie,1}=x_{zombie,0}-\Delta x$

Over time, the proceeds that were set aside in the zombie share and base
reserves accrue interest because they are still held in the underlying yield
source. The interest that has been earned, $\Delta x_{interest}$, can be
measured by taking the difference between the current value of $z_{zombie}$
and $x_{zombie}$ as follows:

$\Delta x_{interest}=c\cdot z_{zombie}-x_{zombie}$

Applying this delta to the zombie share reserves results in the zombie share
and base reserves being equal in value after zombie interest is collected.
These collections occur during automated checkpoints and when positions are
closed after maturity. Most of the interest earned on these positions is
distributed to LPs, and a portion is (potentially) distributed to governance,
specified by the parameter $\phi_{g,zombie}$, $0\leq\phi_{g,zombie}<1$. When
zombie interest is collected, the share reserves are increased by the LPs’
portion of the interest, and the zeta adjustment is updated by the same amount
so that interest collection doesn’t impact the pool’s spot price or depth. The
reserves are updated as follows:

$\displaystyle z_{zombie,1}$ $\displaystyle=z_{zombie,0}-\frac{\Delta
x_{interest}}{c}=\frac{x_{zombie}}{c}$ $\displaystyle z_{1}$
$\displaystyle=z_{0}+(1-\phi_{g,zombie})\cdot\frac{\Delta x_{interest}}{c}$
$\displaystyle\zeta_{1}$
$\displaystyle=\zeta_{0}+(1-\phi_{g,zombie})\cdot\frac{\Delta
x_{interest}}{c}$

##  5 Liquidity Provision

For longs and shorts to be opened on a Hyperdrive pool, liquidity must be
provided to back the positions. Traders that add liquidity to Hyperdrive pools
are called liquidity providers, or simply LPs. Instead of opening long or
short positions directly, LPs take the opposing side of new trades by
providing the liquidity needed to open positions. Over time, these passive
traders will collect a portfolio of longs and shorts that back open positions
that cannot be netted out with matching positions. In return for providing
liquidity, LPs are rewarded with trading fees that offset some of the costs
imposed by adverse selection.

Unlike traditional AMMs, immediate liquidity withdrawal in Hyperdrive AMMs is
not always available to LPs because their liquidity may be backing open
positions. To reduce the amount of monitoring needed from LPs to remove their
liquidity, Hyperdrive pools allow LPs to convert their LP shares to withdrawal
shares. Withdrawal shares are LP shares that can be redeemed when liquidity
becomes available after some open positions are closed. Withdrawal shares are
redeemed for base on a first-come-first-served basis, so it is possible for an
LP that removes liquidity after another to redeem their withdrawal shares
first.

Aside from the need for withdrawal shares, Hyperdrive’s LP system also differs
from those of more traditional AMMs in how it prices something called LP
shares (covered in the next section). Since LPs own a portfolio of long and
short positions, share reserves aren’t sufficient for LP share pricing. To
price LP shares safely and fairly, the LP share price and core LP operations
are defined from the idea that the LP share price should be invariant under
all of the trading and LP operations.

###  5.1 LP Shares

The portion of a pool’s liquidity owned by a given liquidity provider is
tracked by the amount of LP shares they own. Liquidity providers receive LP
shares when they add liquidity to a pool. These shares change in value over
time as interest accrues and positions are closed. Since Hyperdrive LPs are
responsible for backing the long and short positions that cannot be netted,
removing liquidity from Hyperdrive is a more complicated process than removing
liquidity from a conventional AMM. To offer liquidity providers fair pricing
and a simple user experience, Hyperdrive’s LP shares pass through three
distinct phases during the LP life cycle.

When liquidity providers add liquidity to the pool, they receive LP shares.
These LP shares are considered to be active LP shares. The proceeds earned by
active LP shares, from fees and closing positions, are continuously reinvested
into the pool to support more trading volume. As active LPs accrue value over
time, the pool becomes more liquid and supports even larger trades. The total
supply of active LP shares is denoted as $l_{a}$. When a liquidity provider
removes some of their liquidity, they receive a combination of base and
withdrawal shares.

Withdrawal shares are LP shares that have been queued for withdrawal. While
the liquidity underlying withdrawal shares still backs open positions, the
proceeds these shares earn are not reinvested into the pool. When solvency
increases after positions are closed, withdrawal shares are moved to the next
phase of the LP share life cycle and marked as ready for withdrawal. The total
supply of withdrawal shares, including withdrawal shares marked as ready for
withdrawal, is denoted as $l_{w}$. Withdrawal shares that haven’t been marked
as ready for withdrawal are priced identically to active LP shares, but those
that have been marked as ready for withdrawal are priced using special
accounting rules.

When capital becomes available to buy back withdrawal shares, the pool will
automatically set aside the proceeds of purchasing withdrawal shares into the
withdrawal pool and mark some withdrawal shares as ready for withdrawal. The
total supply of withdrawal shares that are ready for withdrawal is denoted as
$l_{r}$, and those that are ready for withdrawal are a subset of the total
amount of withdrawal shares. Withdrawal shares can be redeemed on a first-
come-first-serve basis for the proceeds in the withdrawal pool. The capital in
the withdrawal pool is denoted as $z_{r}$, and a liquidity provider that
redeems $\Delta w$ withdraw shares, $\Delta w\leq l_{r}$, will receive the
following amount of base $\Delta x$:

$\Delta x=c\cdot\Delta w\cdot\frac{z_{r}}{l_{r}}$

For the purposes of pricing active LP shares and withdrawal shares, it is
useful to define the total LP supply $l$ as the amount of LP shares that are
not marked as ready for withdrawal. This quantity can be calculated as:

$l=l_{a}+l_{w}-l_{r}$

###  5.2 LP Share Price

Traditional AMMs price LP shares in proportion to the pool’s reserves;
however, this pricing model doesn’t accurately price Hyperdrive LP shares
because it doesn’t account for the open positions owned by the LPs. To price
LP shares, we calculate the LP Share Price of the LP positions. The core
insight used in Hyperdrive’s LP system is that holding the LP share price
invariant during instantaneous trading or the adding or removing of liquidity
prevents sandwich attacks against existing LPs, and ensures that LPs are
treated fairly when adding and removing liquidity. The LP share price $p_{LP}$
is calculated using the LP present value $PV$ and the LP total supply $l$,
which will be defined below, and calculated as the present value per share as
follows:

$p_{LP}=\frac{PV}{l}$

The LP present value is the value of the portfolio of open long and short
positions owned by the pool’s LPs in the current market conditions.
Intuitively, this quantity is calculated by simulating the impact that closing
all the open positions would have on the share reserves. Since the amount of
fees traders pay depends on the spot price when their trade is executed, fees
are not factored into the LP present value. Despite the exclusion of fees, the
present value provides a robust valuation of the assets currently held by LPs.

To calculate the present value efficiently, we take advantage of the fact that
Hyperdrive’s trading function is path-independent by decomposing the
outstanding positions into the net position of new bonds, $y_{new}^{net}$, and
the net position of matured bonds, $y_{mature}^{net}$. To facilitate this
calculation, Hyperdrive tracks the amount of longs outstanding, $y_{l}$, the
average maturity of long positions $t_{l}$, the amount of shorts outstanding
$y_{s}$, and the average maturity of short positions, $t_{s}$. Using this
information, we can calculate the pool’s net position of new bonds,
$y_{new}^{net}$, and the pool’s net position of matured bonds,
$y_{mature}^{net}$, as follows:

$\displaystyle y_{new}^{net}$ $\displaystyle=y_{l}\cdot t_{l}-y_{s}\cdot
t_{s}$ $\displaystyle y_{mature}^{net}$
$\displaystyle=y_{l}\cdot(1-t_{l})-y_{s}\cdot(1-t_{s})$

We define $n_{new}$ as the impact on the share reserves of closing the net
position of new bonds. At the beginning of the term, long positions are closed
by selling the underlying bonds on the trading invariant, which decreases the
share reserves. In contrast, short positions are closed by purchasing the
underlying bonds, which increases the share reserves. To account for the fact
that some trading invariants can run out of inventory above the maximum sell
amount, $y_{sell}^{max}$, or the maximum buy amount, $y_{buy}^{max}$, the part
of the position that can’t be sold is marked to a price of $0$ and the part
that can’t be bought is marked to a price of $1$. The impact of closing the
net position of new bonds is calculated as follows:

$n_{new}=\left\\{\begin{array}[]{lr}\hbox to0.0pt{\hss-}I_{\Delta z}^{\Delta
y}\left(\min(y_{new}^{net},y_{max}^{sell}),z,y\right),&\text{if
}y_{new}^{net}>0\\\ I_{\Delta z}^{\Delta y}\left(\,\,\hbox
to0.0pt{\hss-}\min(|y_{new}^{net}|,y_{max}^{buy}),z,y\right)+\max(|y_{new}^{net}|-y_{max}^{buy},0),&\text{if
}y_{new}^{net}<0\\\ 0,&\text{if }y_{new}^{net}=0\end{array}\right.$

We define $n_{mature}$ as the impact on the share reserves of closing the net
position of matured bonds. At maturity, longs can be redeemed for their full
face value, which reduces the share reserves. On the other hand, redeeming a
matured short increases the share reserves by the full face value since the
underlying bonds are purchased back from the pool. Since all positions can be
closed at maturity, this calculation does not need to account for liquidity
constraints. The impact of closing the net position of mature bonds is
calculated as follows:

$n_{mature}=\left\\{\begin{array}[]{lr}\hbox
to0.0pt{\hss-}M(y_{mature}^{net}),&\text{if }y_{mature}^{net}>0\\\
M(|y_{mature}^{net}|),&\text{if }y_{mature}^{net}<0\\\ 0,&\text{if
}y_{mature}^{net}=0\end{array}\right.$

Now that we’ve defined the impact of closing the net positions of new and
matured bonds, we are in a position to define the present value. This
calculation simply applies the impact of closing each net position to the
share reserves and subtracts $z_{min}$. This calculation simulates the total
amount of liquidity that LPs could remove if all of the outstanding positions
were instantly closed without considering fees. We can calculate the present
value as follows:

$PV=z+net_{new}+net_{mature}-z_{min}$

###  5.3 Adding Liquidity

Traders can open LP positions or add to existing LP positions by adding
liquidity to the pool. When a trader adds a given amount of base $\Delta x$,
the pool’s share reserves are updated to include the additional base as
follows:

$z_{1}=z_{0}+\frac{\Delta x}{c}.$

To ensure that the effective share reserves increase by the same ratio as the
share reserves, the zeta adjustment is updated as follows:

$\zeta_{1}=\zeta_{0}\cdot\frac{z_{1}}{z_{0}}.$

The pool’s spot price should remain invariant when LPs add liquidity. This can
be accomplished by maintaining the same ratio between the effective share
reserves and the bond reserves as we had before, which means that the bond
reserves are updated as follows:

$y_{1}=y_{0}\cdot\frac{z_{1}}{z_{0}}.$

In return for adding liquidity to the pool, LPs are compensated with new LP
shares. To ensure that LPs are fairly rewarded for adding liquidity, we solve
for the amount of LP shares by preserving the LP share price as follows:

$\frac{PV_{0}}{l}=\frac{PV_{1}}{l+\Delta l}\implies\Delta
l=\frac{(PV_{1}-PV_{0})\cdot l}{PV_{0}}.$

###  5.4 Removing Liquidity

When LPs remove liquidity, all of their LP shares are immediately converted
into withdrawal shares. The system will attempt to redeem as many withdrawal
shares as possible to return base to the LP, and any withdrawal shares that
couldn’t be redeemed are minted to the LP. To determine how many of the LP’s
withdrawal shares can be redeemed, the system first attempts to buy back as
many withdrawal shares as possible using the available idle liquidity. Once
the system has distributed the maximum amount of excess idle liquidity
(simplified to “excess idle”), as many of the withdrawal shares as possible
will be paid out from the withdrawal pool.

###  5.5 Distributing Excess Idle

The first time that excess idle is distributed is when LPs initially remove
their liquidity; however, the system attempts to distribute excess idle
whenever checkpoints are minted or when traders and LPs interact with the
system. Distributing excess idle whenever it becomes available ensures that
withdrawal shares are paid out as quickly as liquidity becomes available,
which ensures that they will always be paid out within a single term after
withdrawing.

The idle liquidity in Hyperdrive is defined as the amount of liquidity that
can be removed by LPs without violating the solvency constraints. This idle
liquidity defines an upper limit on the amount of shares that can be used to
buy back withdrawal shares. We can calculate the amount of idle liquidity,
$z_{idle}$, using the solvency equation as follows:

$z_{idle}=z-\tfrac{e}{c}-z_{min}$

The core idea behind the pricing of withdrawal shares is that the LP share
price must be conserved whenever withdrawal shares are bought back. In order
to define the problem formally, we will define $PV(\Delta z)$ to be the LP
present value of the pool after $\Delta z$ shares are removed from the
reserves, $l_{0}$ to be the starting LP total supply, and $\Delta w$ to be the
amount of withdrawal shares that are paid out by the system. Using this
notation, we can express the problem of determining how much idle to
distribute to withdrawal shares as follows:

$\frac{PV(0)}{l_{0}}=\frac{PV(\Delta z)}{l_{0}-\Delta w}$

Unlike the conservation of LP share price formula used when LPs add liquidity,
this formula has two unknowns, $\Delta z$ and $\Delta w$. The fact that we
can’t just solve algebraically means that we must find the solution by solving
an optimization problem. In order to solve this undetermined system, we
formulate our optimization problem to hold the starting and ending LP share
prices equal as follows:

$\displaystyle\max$ $\displaystyle\Delta w$ s.t.
$\displaystyle\frac{PV(0)}{l}=\frac{PV(\Delta z)}{l-\Delta w}$
$\displaystyle\Delta z\leq z_{idle}$ $\displaystyle\Delta w\leq w$
$\displaystyle y_{s}\cdot t_{s}-y_{l}\cdot t_{l}\leq y_{out}^{max}(\Delta z)$

This problem can be solved efficiently using the following algorithm:

  1. 1.

If the pool’s net curve position isn’t net short, i.e. $y_{l}\cdot t_{l}\geq
y_{s}\cdot t_{s}$, or if the maximum amount of bonds that can be purchased
after removing all of the excess idle liquidity still exceeds the net curve
position, i.e. $y_{out}^{max}(z_{idle})\geq y_{s}\cdot t_{s}-y_{l}\cdot
t_{l}$, then we set the upper bound for idle distribution to the total amount
of idle liquidity, $\Delta z_{max}=z_{idle}$, and proceed directly to step (3)
of the algorithm. Otherwise, we proceed to step (2).

  2. 2.

Solve $y_{out}^{max}(\Delta z_{max})=y_{s}\cdot t_{s}-y_{l}\cdot t_{l}$. This
provides an upper bound on the amount of idle that can distributed without
causing the present value to increase. The importance of this step is that it
prevents pathological cases where the present value can begin to increase when
capital is removed from the share reserves. These cases occur when the net
curve position is short and exceeds the maximum amount of bonds that can
purchased from the pool.

  3. 3.

Solve for the amount of withdrawal shares that would be redeemed if $\Delta
z_{max}$ is removed from the share reserves. This can be calculated as
$\Delta w=\left(1-\tfrac{PV(\Delta z_{max})}{PV(0)}\right)\cdot l$. If $\Delta
w\leq w$, then proceed to step (5). Otherwise, set $\Delta w=w$ and continue
to step (4).

  4. 4.

Solve $\frac{PV(0)}{l_{0}}=\frac{PV(\Delta z)}{l-\Delta w}$ for $\Delta z$
using Newton’s method if $y_{l}\cdot t_{l}\neq y_{s}\cdot t_{s}$ or directly
otherwise.

  5. 5.

The amount of idle that can be distributed is $\Delta z$, and the amount of
withdrawal shares that the pool can buy back is $\Delta w$. The algorithm
terminates.

After the optimization problem is solved, the share reserves are debited by
the amount of idle liquidity that could be distributed, $\Delta z$. The zeta
adjustment and bond reserves are updated in a similar fashion to the way that
they are updated when adding liquidity to ensure that the pool’s spot price
remains constant after distributing excess idle. The updated reserves are
calculated as follows:

$\displaystyle z_{1}=z_{0}-\Delta z$
$\displaystyle\zeta_{1}=\zeta_{0}\cdot\frac{z_{1}}{z_{0}}$ $\displaystyle
y_{1}=y_{0}\cdot\frac{z_{1}}{z_{0}}$

Next, the withdrawal shares that the pool can buy back, $\Delta w$, are marked
as ready for withdrawal:

$\displaystyle l_{w}$ $\displaystyle\mathrel{-}=\Delta w$ $\displaystyle
l_{r}$ $\displaystyle\mathrel{+}=\Delta w$

Finally, idle liquidity that could be distributed, $\Delta z$, is added to the
withdrawal pool:

$\Delta z_{r}\mathrel{+}=\Delta z$

## Acknowledgements

We are grateful for the discussions and contributions from Violet Vienhage,
Charles St. Louis, Dylan Paiton, Sheng Lundquist, Giovanni Effio, Daejun Park,
Tim Roughgarden, Allan Niemberg, and Ben Edgington.

## References

  * [1] Alan Lu. 2017. _Building a Decentralized Exchange in Ethereum_ Retrieved from https://medium.com/@gnosisPM/building-a-decentralized-exchange-in-ethereum-ae9df6860b08
  * [2] Allan Niemerg, Dan Robinson, and Lev Livnev. 2020. _YieldSpace: An Automated Liquidity Provider for Fixed Yield Tokens._ Retrieved from https://yield.is/YieldSpace.pdf
  * [3] Peter Zeitz. 2020. _Curve Vulnerability Report_ Retrieved from https://medium.com/@peter_4205/curve-vulnerability-report-a1d7630140ec

## Disclaimer

This paper is for general information purposes only. It does not constitute
investment advice or a recommendation or solicitation to buy or sell any
investment and should not be used in the evaluation of the merits of making
any investment decision. It should not be relied upon for accounting, legal,
or tax advice or investment recommendations. This paper reflects the current
opinions of the authors and is not made on behalf of DELV or their affiliates
and does not necessarily reflect the opinions of DELV, its affiliates, or
individuals or entities associated with DELV. The opinions reflected herein
are subject to change without being updated.


[2]Krishnan Suresh 1]Computer Science, University of Wisconsin, Madison, 1210
W. Dayton Street, Madison, 53706, WI, USA [2]Mechanical Engineering,
University of Wisconsin, Madison, 1513 University Avenue, Madison, 53706, WI,
USA

# Optimal Box Contraction for Solving Linear Systems via Simulated and Quantum
Annealing

Sanjay Suresh<EMAIL_ADDRESS><EMAIL_ADDRESS>[  *

###### Abstract

Solving linear systems of equations is an important problem in science and
engineering. Many quantum algorithms, such as the Harrow-Hassidim-Lloyd (HHL)
algorithm (for quantum-gate computers) and the box algorithm (for quantum-
annealing machines), have been proposed for solving such systems.

The focus of this paper is on improving the efficiency of the box algorithm.
The basic principle behind this algorithm is to transform the linear system
into a series of quadratic unconstrained binary optimization (QUBO) problems,
which are then solved on annealing machines.

The computational efficiency of the box algorithm is entirely determined by
the number of iterations, which, in turn, depends on the box contraction
ratio, typically set to 0.5. Here, we show through theory that a contraction
ratio of 0.5 is sub-optimal and that we can achieve a speed-up with a
contraction ratio of 0.2. This is confirmed through numerical experiments
where a speed-up between $20\%$ to $60\%$ is observed when the optimal
contraction ratio is used.

###### keywords:

QUBO; Linear system of equations; Quantum annealing; Simulated annealing; Box
algorithm; D-WAVE; Quantum computing

##  1 Introduction

Solving least squares problems and linear systems of equations are of utmost
importance in science and engineering. Many algorithms have been proposed to
solve such problems on classical computers. Quantum computers have recently
been proposed as an alternate since they can potentially accelerate the
computation [1, 2]. In particular, the Harrow-Hassidim-Lloyd (HHL) algorithm
is a landmark strategy for solving linear systems of equations on _quantum-
gate computers_. In theory, it offers an exponential speed-up over classical
algorithms [3], and it has been further improved recently [4, 5, 6, 7].
However, due to the accumulation of errors in current noisy intermediate-scale
quantum (NISQ) computers [8], the HHL algorithm and its variants are limited,
in practice, to tiny systems [9]. Furthermore, extracting the target state can
be expensive/impractical [10].

In parallel, _quantum annealing machines_ , such as the D-Wave systems with
several thousand qubits [11], have also been proposed for solving such
problems since they are less susceptible to noise [12, 13]. The basic
principle is to convert the least squares and linear system into a series of
quadratic unconstrained binary optimization (QUBO) problems. For example,
O’Malley and Vesselinov solved the least-squares problem using a finite-
precision qubit representation [14]. Borle and Lomonaco conducted a
theoretical and numerical analysis of this approach [15, 16].

In this paper, we consider solving linear system of equations via the QUBO
formulation. If the matrix is positive-definite, which is often the case, one
can pose this as a energy minimization problem and thereby convert the linear
system into a series of QUBO problems. However, current quantum annealing
machines are only equipped with about 2000 qubits, with additional
restrictions on connectivity [11]. This implies that: (1) the size of the
linear system is somewhat limited, and (2) the solution can only be computed
with limited precision. The first limitation can be addressed through a hybrid
Gauss-Seidel strategy [17]. In contrast, the second limitation typically
relies on the _iterative box algorithm_ [18], the main focus of this paper.
Alternate iterative methods are proposed in [19].

In the box algorithm, the number of iterations strongly depends on the box
contraction ratio, i.e., the ratio by which the box size reduces under certain
conditions (see Section 2). This ratio is typically set to 0.5 [18]. In this
paper, we show, through theoretical analysis, that a contraction ratio of 0.5
is sub-optimal and that we can achieve a speed-up with a contraction ratio of
(approximately) 0.2. This is confirmed through numerical experiments using
simulated and quantum annealing.

##  2 Background

###  2.1 QUBO Formulation

Consider the following linear system of equations:

$\mathbf{A}\mathbf{x}=\mathbf{b}$ (1)

where $\mathbf{A}$ is a $d\times d$ matrix. If $\mathbf{A}$ is positive-
definite (assumed to be true in the remainder of the paper), then solving Eq.
(1) is equivalent to minimizing the potential energy:

$\min_{\mathbf{x}}\Pi=\frac{1}{2}\mathbf{x}^{T}\mathbf{A}\mathbf{x}-\mathbf{x}^{T}\mathbf{b}$
(2)

We represent each real component $x_{j}$ using qubit variables to solve this
on a quantum annealing machine. A well-known strategy is the two’s complement
radix representation [16, 15]; for example, a scalar variable $x$ can be
represented using $m$ qubits as:

$x=-q_{1}2^{m-1}+\sum_{i=2}^{m}q_{i}2^{i-2}$ (3)

Since the number of qubits is often limited in a quantum machine, it is common
to let $m=2$, leading to:

$x=-2q_{1}+q_{2}$ (4)

Of course, this can only capture the numbers $\\{-2,-1,0,1\\}$. However, one
can easily extend this to a wide range of real numbers through scaling $L$ and
offset $c$ via (see Section 2.2 for further explanation):

$x=c+L(-2q_{1}+q_{2})$ (5)

This is often called the _box representation_ [18]. Furthermore, one can
easily generalize this to arbitrary dimension $d$ via:

$\mathbf{x}=\mathbf{c}+L(-2\mathbf{q}_{1}+\mathbf{q}_{2})$ (6)

where $\mathbf{q}_{1}$ and $\mathbf{q}_{2}$ are qubit vectors of length $d$,
i.e., a total of $2d$ qubits is used to capture $\mathbf{x}$. Thus, a
$d$-dimensional system is associated with $4^{d}$ total states. This is
illustrated schematically in Fig. 1 for $d=2$.

Figure 1: The box representation for $d=2$.

Since $\mathbf{x}$ is linear in $\mathbf{q}_{1}$ and $\mathbf{q}_{2}$,
substituting Eq. (6) into Eq. (2) leads to a quadratic unconstrained _binary_
optimization (QUBO) problem:

$\min_{\mathbf{q}=\\{\mathbf{q}_{1},\mathbf{q}_{2}\\}}\Pi=\frac{1}{2}\mathbf{q}^{T}\mathbf{Q}^{\prime}\mathbf{q}+\mathbf{q}^{T}\mathbf{r}$
(7)

Furthermore, since the qubit variables can only take the values $0$ or $1$,
the linear term can be absorbed into the quadratic term [20], resulting in the
standard form:

$\min_{\mathbf{q}=\\{\mathbf{q}_{1},\mathbf{q}_{2}\\}}\Pi=\frac{1}{2}\mathbf{q}^{T}\mathbf{Q}\mathbf{q}$
(8)

where $\mathbf{Q}$ is symmetric.

###  2.2 Box Algorithm

We now describe the box algorithm (see Algorithm 1) that exploits the QUBO
formulation to solve Eq. (1) to a high degree of precision. During each
iteration of the algorithm, the center $\mathbf{c}$ or scale $L$ is updated as
follows. In a particular iteration, when a QUBO problem in Eq. 8 is solved, if
a lower potential energy state than the current state is reached, $\mathbf{c}$
is updated (referred to as a _translation_ ; see Fig 2a), else $L$ is reduced
by a factor of, typically, $0.5$ (referred to as a _contraction_ ; see Fig
2b). The iteration is then continued until $L$ falls below the precision
desired.

Figure 2: Box algorithm: (a) translation, (b) contraction.

We now describe the box algorithm (see Algorithm 1) in detail since it will be
relevant for the remainder of the paper. Observe in the algorithm that:

  * •

Lines 2-7: Various quantities are initialized: (1) the center $\mathbf{c}$ is
initialized to $\mathbf{0}$, (2) the length $L$ is set to $1$ (see remark
below), (3) the qubit vectors $\mathbf{q}_{1}$ and $\mathbf{q}_{2}$ are
created, (4) the number of translations $N_{t}$, and contractions $N_{t}$, are
initialized to $0$, (5) the contraction ratio $\beta$ is initialized to $0.5$
and (5) the potential energy is initialized to 0.

  * •

Line 9: The unknown vector $\mathbf{x}$ is represented via the qubits and the
current $\mathbf{c}$ and $L$ .

  * •

Line 10: The potential energy $\Pi$ is formulated.

  * •

Line 11: The minimum value of $\Pi$ and the corresponding qubit values are
determined either via simulated annealing or quantum annealing.

  * •

Lines 12-15: If a lower energy state is found, the center is translated and
$N_{t}$, the number of translations, is incremented (see remark below).

  * •

Lines 16-18: Else, $N_{c}$, the number of contractions, is incremented and $L$
is reduced by a factor of $\beta$.

  * •

Line 20: The algorithm terminates if $L$ is less than the desired tolerance,
or if the number of total iterations $(N_{t}+N_{c})$ is greater than an
allowable $N_{\text{allowable}}$

Remark on Line 3: Although $L$ is initialized to $1$, the box algorithm is
robust in that it converges for any reasonable value of $L$ [18]. However,
choosing $L$ incorrectly will lead to slower convergence.

Remark on Line 12: We have observed that the box algorithm is more stable and
unnecessary translations can be avoided if we add a small buffer by checking
if $\Pi^{*}<\hat{\Pi}(1+10^{-8}$).

Algorithm 1 Box Algorithm

1:procedure BoxAlg($\mathbf{A}$, $\mathbf{b}$, $\epsilon$,$N_{allowable}$ )

2: $\mathbf{c}\leftarrow\mathbf{0}$ $\triangleright$ Center of length $d$

3: $L\leftarrow 1$ $\triangleright$ Initialize box size

4: $\mathbf{q}_{1},\mathbf{q}_{2}\leftarrow\text{Qubits}(d)$ $\triangleright$
Create qubit arrays of length d

5: $N_{c}=N_{t}=0$ $\triangleright$ Translation and contraction steps set to 0

6: $\beta=0.5$ $\triangleright$ Contraction ratio

7: $\hat{\Pi}=0$ $\triangleright$ Initial potential energy

8: repeat$\triangleright$ Until convergence

9:
$\mathbf{x}\leftarrow\mathbf{c}+L(-2\mathbf{q}_{1}+\mathbf{q}_{2})$$\triangleright$
Symbolic expression

10:
$\Pi\leftarrow\frac{1}{2}\mathbf{x}^{T}\mathbf{A}\mathbf{x}-\mathbf{x}^{T}\mathbf{b}$
$\triangleright$ Construct QUBO

11:
$\Pi^{*},\mathbf{q}_{1}^{*},\mathbf{q}_{2}^{*}\leftarrow\text{minimize}(\Pi)$
$\triangleright$ Solve QUBO

12: if $\Pi^{*}<\hat{\Pi}$ then $\triangleright$ We have found a lower energy
state

13:
$\mathbf{c}\leftarrow\mathbf{c}+L(-2\mathbf{q}_{1}^{*}+\mathbf{q}_{2}^{*})$
$\triangleright$ Translation of box

14: $N_{t}=N_{t}+1$ $\triangleright$ Update translation counter

15: $\hat{\Pi}=\Pi^{*}$ $\triangleright$ Update the lowest potential energy

16: else

17: $L\leftarrow\beta L$ $\triangleright$ Reduce box size

18: $N_{c}=N_{c}+1$ $\triangleright$ Update contraction counter

19: end if

20: until ($L<\epsilon$) or ($N_{c}+N_{t}>N_{allowable}$) $\triangleright$
Termination

21:end procedure$\triangleright$ Output: solution $\mathbf{c}$

A typical convergence of the box algorithm in 2D is illustrated in Fig 3.
Observe that the number of QUBO problems one must solve is equal to the total
number of iterations $(N=N_{t}+N_{c})$. The objective of this paper is to
reduce $N$ by finding an optimal value for $\beta$. In particular, we show, in
the next section, that the default value of $\beta=0.5$ recommended in the
literature is sub-optimal.

Figure 3: Typical box convergence in 2D.

##  3 Box Algorithm Analysis

###  3.1 Contraction steps

Observe that each time the box contracts, the length of the box reduces by a
factor of $\beta$, until $L<\epsilon$. Consequently, the total number of
contractions, independent of the number of translations and the linear system
being solved, is given by

$N_{c}\sim\log_{\beta}\epsilon$ (9)

###  3.2 Translation steps

Although the number of contraction steps is independent of $\beta$, the number
of translations depends on $\beta$. Our objective is to determine an upper
bound $\hat{N}_{t}$ and average estimate $\overline{N}_{t}$ in terms of
$\beta$ and $\epsilon$.

####  3.2.1 1D Box Algorithm

We will first consider one-dimensional problems that requires only $2$ qubits.
Let $L_{i}$ be the length of the box _before_ the $i$-th contraction, i.e.,
$L_{1}=1$, $L_{2}=\beta$, $L_{3}=\beta^{2}$, etc. Further, at the start of the
algorithm, $c=0$, and $c$ gets updated as follows:

$c\leftarrow c+L_{i}(-2q_{1}+q_{2})$ (10)

Let $n_{i}$ be the number of translation steps _before_ the $i^{th}$
contraction. Our objective is to find an upper bound for $n_{i}$.

We will assume that the solution $x$ lies within the range $(-2,1)$, else
arriving at an upper bound for $N_{t}$ is impossible. Consider the possible
sequence of translations, starting at $c=0$, before the first contraction. If
there is no translation, then $n_{1}=0$, else, $c$ gets updated to
$\\{-2,-1,1\\}$, as per Eq (10). Suppose $c\leftarrow-1$, corresponding to
$q_{1}=1$, $q_{2}=1$, then the solution $x$ must lie in the range
$(-1.5,-0.5)$. Similarly, if $c\leftarrow-2$, or $c\leftarrow 1$, $x$ must lie
in the range $(-2,-1.5)$ or $(0.5,1)$ respectively.

No further translation is possible since it would require $c$ to translate to
an inferior solution. In other words, the box must contract in the next
iteration. In summary, $n_{1}\leq 1$.

After this, the box will contract, and $L_{2}=\beta$, and $x$ must lie in the
range $[c-1/2,c+1/2]$, where $c$ is the updated center. Therefore, the next
sequence of translations can move the center by at most $0.5$ units where each
translation is at least $L_{2}=\beta$, as per Eq. (10). Note that the center
can also translate by $-2\beta$ per Eq. (10), However, since we are seeking an
upper bound for $N_{t}$, we consider the worst-case scenario. Therefore, the
maximum number of translations is given by $n_{2}\leq\frac{0.5}{\beta}$. After
this, the box must contract, resulting in $L_{3}=\beta^{2}$.

After the contraction, the solution $x$ must lie in the range
$[c-\frac{\beta}{2},c+\frac{\beta}{2}]$, where $c$ is the updated center. The
next sequence of translations can move the center by at most $\beta/2$ units
and each translation is at least $L_{3}=\beta^{2}$, as per Eq. (10).
Consequently, the maximum number of translations
$n_{3}\leq(\beta/2)/\beta^{2}=\frac{0.5}{\beta}$. After this, the box
contracts to $L_{4}=\beta^{3}$.

Repeating this logic for all $N_{c}$ contractions, we have $n_{1}\leq 1$ and
$n_{2},n_{3},\dots,n_{N_{c}}\leq\frac{1}{2\beta}$. Therefore, the total number
of translation steps has the following upper bound:

$\hat{N}_{t}=\sum_{i=1}^{N_{c}}n_{i}=1+\frac{1}{2\beta}+\frac{1}{2\beta}+\cdots+\frac{1}{2\beta}=1+\frac{N_{c}-1}{2\beta}$
(11)

If we assume $\beta\leq 0.5$ (as we will see below, this is justifiable):

$\hat{N}_{t}=\frac{N_{c}}{2\beta}$ (12)

Combining this with Eq. (9), we have an upper bound on the total iterations

$\hat{N}=\hat{N}_{t}+N_{c}=\left(1+\frac{1}{2\beta}\right)\log_{\beta}\epsilon$
(13)

To find the optimal value of $\beta$, we can the derivative of $\hat{N}$ in
Eq. 13 with respect to $\beta$ and set this equal to $0$, resulting in:

$\ln\beta+2\beta+1=0$ (14)

Solving this numerically, we obtain $\beta^{*}\sim 0.232$, independent of
$\epsilon$. This is illustrated in Fig. 4. For this optimal value, one can
observe a $32\%$ reduction in the maximum box iterations, compared to the
default value of $\beta=0.5$. Further, observe that the number of translations
increases for $\beta>0.5$, justifying our earlier assumption.

Figure 4: Upper bound $\hat{N}$ on the number of box iterations.

Instead of the upper bound, one can also consider the average number of
translations $\overline{N}_{t}$. There are three different scenarios during
each translation ($-1$, $-2$, $1$). If we assume there is an equal probability
of translating in each of these directions (a very simplistic model), then the
expected translation is given by
$(\lvert-1\lvert+\lvert-2\lvert+\lvert+1\lvert)/3=4/3$ (as opposed to $1$ in
the worse case). Consequently

$\overline{N}_{t}=\frac{3N_{c}}{8\beta}$ (15)

Since the number of contractions remains the same, we have:

$\overline{N}=\overline{N}_{t}+N_{c}\approx\left(1+\frac{3}{8\beta}\right)\log_{\beta}\epsilon$
(16)

Fig. 5 illustrates $\overline{N}$ vs $\beta$. Taking the derivative of
$\overline{N}$ with respect to $\beta$ and setting it equal to $0$, we get
$\beta^{*}\approx 0.21$. For this optimal value, one can observe a $44\%$
reduction in the average box iterations, compared to the default value of
$\beta=0.5$.

Figure 5: Average number of box iterations $\overline{N}$.

####  3.2.2 Multi-Dimensional Problems

In this section, we argue that the previous results also holds true in higher
dimensions. For a $d$-dimensional problem, let $n_{i,k}$ represent the number
of translations in the $k^{th}$ dimension before the $i^{th}$ contraction.
Since the dimensions are independent, $n_{i,k}$ are all independent of each
other.

In order to prove this, consider the following hypothetical scenario: Before
the $i^{th}$ contraction, let the center translate in the first dimension
until it cannot translate anymore in this direction. Then let it translate in
dimension $2$, and so on. Once dimension $d$ is reached, this process is
repeated until the box no longer can translate in any dimension. Let $n_{i,k}$
be the total number of times the box contracted in each dimension. By the
premises of the box-algorithm, the box must now contract.

Following the logic from the previous section, we have $\forall i\geq
2,\forall k,n_{i,k}\leq\frac{1}{2\beta}$. However, the maximum number of
translations is dictated by one or more of the dimensions. Therefore, let
$n_{i}=\max_{k}n_{i,k}$. As a result, the upper bound on translation is given
by

$\hat{N}_{t}=\sum_{i=1}^{N_{c}-1}n_{i}=\frac{N_{c}}{2\beta}$ (17)

Consequently,

$\hat{N}=\left(1+\frac{1}{2\beta}\right)\log_{\beta}\epsilon$ (18)

irrespective of the number of dimensions. This is later confirmed in the next
section through numerical experiments. Similar arguments can be made for the
average case.

##  4 Numerical Experiments

We will now carry out numerical experiments to validate the theoretical
analysis. The experiments rely heavily on simulated annealing (SA) since
quantum annealing (QA) is expensive today. However, a limited number of QA
experiments are also carried out. For SA, we rely on D-Wave’s Neal annealer;
for hybrid QA, we rely on D-Wave’s LeapHybridSampler, and for (pure) QA, we
rely on DWaveSampler, with EmbeddingComposite. For all three methods, 20
samples were used.

The QUBO problems were constructed using the pyQUBO package [20]. The Python
code used in generating the results in this section is available from the
github link provided towards the end of the paper.

###  4.1  Positive Definite Matrices

In the first set of experiments, we generate random $d$-dimensional positive
definite matrices $\mathbf{A}$. Further, since we require $\mathbf{x}$ to lie
within $[-2,1]^{d}$, we first generate $\mathbf{x}$, and then construct the
corresponding right-hand-side $\mathbf{b}$. The corresponding Python code is
given below:

⬇

1B = np.random.rand(d, d)

2A = d*np.eye(d)-(B + B.transpose())/2

3x = np.array([random.uniform(-2, 1) for _ in range(d)])

4b = A.dot(xExact)

Listing 1: Generating d-dimensional positive definite matrices and right-hand
side.

In order to capture the average behavior of the box algorithm, we create ten
instances of $\mathbf{A}$ and $\mathbf{b}$, for $d=2,d=10$ and $d=20$.
Finally, for each instance, we use the box algorithm (see Algorithm 1) to
solve for $\mathbf{x}$ for $\epsilon=10^{-6}$ and $\epsilon=10^{-8}$, for
various values of $\beta$. All experiments in this section are carried out
using SA. The results are summarized in Fig. 6. Observe the following:

  * •

All three graphs exhibit a minima around $\beta=0.2$, independent of the
dimension $d$ and desired accuracy $\epsilon$.

  * •

The number of iterations is (nearly) independent of the dimension of the
problem.

  * •

The number of iterations is closer to the theoretical prediction
$\overline{N}$ in Fig. 5 than to the upper bound prediction $\hat{N}$ in Fig.
4.

Figure 6: Observed $N$ vs $\beta$ averaged over ten $d\times d$ problems: (a)
d = 2, (b) d = 10, and (c) d = 20.

###  4.2 1D Poisson Problem

For the next experiment, we construct a $6\times 6$ matrix $\mathbf{A}$ that
arises from a finite difference formulation of 1D Poisson problem [21]:

⬇

1A = np.array([[6,-6,0,0,0,0],[-6,12,-6,0,0,0],

2[0,-6,12,-6,0,0],[0,0,-6,12,-6,0],

3[0,0,0,-6,12,-6],[0,0,0,0,-6,12]])

4xExact = np.array([-np.pi/9, np.pi/11, -np.pi/20,

5np.pi/8, 0.05*np.pi, -np.pi/5 ])

6b = A.dot(xExact)])

Listing 2: Generating a 6-dimensional finite difference matrix and right hand
side.

The results for SA, hybridQA, and QA are summarized in Figure 7. Observe that:

  * •

The overall behavior of $N$ vs. $\beta$ is consistent with the theory.

  * •

The hybridQA results precisely match that of SA for the three sampled points,
suggesting that D-Wave probably relied entirely on CPU for this scenario.

  * •

QA performed poorly compared to SA or hybridQA. This is consistent with the
observations in [16]. However, even in this case, a $50\%$ improvement in
performance can be observed for $\beta=0.2$, compared to $\beta=0.5$.

Figure 7: Observed $N$ versus $\beta$ for a 1D Poisson problem.

##  5 Conclusions

The box algorithm is a popular method for solving linear systems of equations
via the QUBO formulation. In this paper, a theoretical analysis of the box
algorithm was carried out that suggested that a computational speed-up can be
easily achieved by making a simple modification to the algorithm.
Specifically, the theory suggests that a $43\%$ speed-up can be obtained by
reducing the box contraction ratio from 0.5 to 0.2. This was confirmed through
numerical experiments where a speed-up between $20\%$ to $60\%$ was observed.
Additional experiments involving a larger class of linear systems are needed
to corroborate these results.

While the paper focused on linear systems, the strategy can be extended to
least squares systems, and other direct methods for solving linear systems via
the QUBO formulation [19]. Further, the analysis here was restricted to the
case when only 2 qubits are used to represent each scalar variable. The
extension to the more general case needs to be investigated. Finally, as
observed in [16], while it is possible to exploit quantum annealing to solve
linear systems of equations, it is currently limited to small systems and low
precision.

## Compliance with ethical standards

The authors declare that they have no conflict of interest.

## Replication of Results

The Python code pertinent to this paper is available at https://github.com/UW-
ERSL/QUBOBoxContraction.

## Acknowledgments

We would like to thank the Graduate School of UW-Madison for the Vilas
Associate grant.

## References

  *   * [1] G. Tosti Balducci, B. Chen, M. Möller, M. Gerritsma, R. De Breuker, Review and perspectives in quantum computing for partial differential equations in structural mechanics.  Frontiers in Mechanical Engineering p. 75 (2022)
  * [2] Y. Wang, J.E. Kim, K. Suresh, Opportunities and challenges of quantum computing for engineering optimization.  Journal of Computing and Information Science in Engineering 23(6) (2023)
  * [3] A.W. Harrow, A. Hassidim, S. Lloyd, Quantum algorithm for linear systems of equations.  Physical review letters 103(15), 150502 (2009)
  * [4] A. Ambainis, Variable time amplitude amplification and a faster quantum algorithm for solving systems of linear equations.  arXiv preprint arXiv:1010.4458 (2010)
  * [5] A.M. Childs, R. Kothari, R.D. Somma, Quantum algorithm for systems of linear equations with exponentially improved dependence on precision.  SIAM Journal on Computing 46(6), 1920–1950 (2017)
  * [6] X. Liu, H. Xie, Z. Liu, C. Zhao, Survey on the improvement and application of HHL algorithm.  Journal of Physics: Conference Series 2333(1), 012023 (2022)
  * [7] P.C. Costa, D. An, Y.R. Sanders, Y. Su, R. Babbush, D.W. Berry, Optimal scaling quantum linear-systems solver via discrete adiabatic theorem.  PRX Quantum 3(4), 040303 (2022)
  * [8] J. Preskill, Quantum computing in the NISQ era and beyond.  Quantum 2, 79 (2018)
  * [9] W. Ji, X. Meng, in _2020 19th International Symposium on Distributed Computing and Applications for Business Engineering and Science (DCABES)_ (IEEE, 2020), pp. 308–311
  * [10] B.D. Clader, B.C. Jacobs, C.R. Sprouse, Preconditioned quantum linear system algorithm.  Physical review letters 110(25), 250504 (2013)
  * [11] S.W. Shin, G. Smith, J.A. Smolin, U. Vazirani, How quantum is the d-wave machine?  arXiv preprint arXiv:1401.7087 (2014)
  * [12] P. Hauke, H.G. Katzgraber, W. Lechner, H. Nishimori, W.D. Oliver, Perspectives of quantum annealing: Methods and implementations.  Reports on Progress in Physics 83(5), 054401 (2020)
  * [13] S. Yarkoni, E. Raponi, T. Bäck, S. Schmitt, Quantum annealing for industry applications: Introduction and review.  Reports on Progress in Physics (2022)
  * [14] D. O’Malley, V.V. Vesselinov, B.S. Alexandrov, L.B. Alexandrov, Nonnegative/binary matrix factorization with a d-wave quantum annealer.  PloS one 13(12), e0206653 (2018)
  * [15] A. Borle, S.J. Lomonaco, in _WALCOM: Algorithms and Computation: 13th International Conference, WALCOM 2019, Guwahati, India, February 27–March 2, 2019, Proceedings 13_ (Springer, 2019), pp. 289–301
  * [16] A. Borle, S.J. Lomonaco, How viable is quantum annealing for solving linear algebra problems?  arXiv preprint arXiv:2206.10576 (2022)
  * [17] G.G. Pollachini, J.P. Salazar, C.B. Góes, T.O. Maciel, E.I. Duzzioni, Hybrid classical-quantum approach to solve the heat equation using quantum annealers.  Physical Review A 104(3), 032426 (2021)
  * [18] S. Srivastava, V. Sundararaghavan, Box algorithm for the solution of differential equations on a quantum annealer.  Physical Review A 99(5), 052355 (2019)
  * [19] O.M. Raisuddin, S. De, Feqa: Finite element computations on quantum annealers.  Computer Methods in Applied Mechanics and Engineering 395, 115014 (2022)
  * [20] M. Zaman, K. Tanahashi, S. Tanaka, Pyqubo: Python library for mapping combinatorial optimization problems to qubo form.  IEEE Transactions on Computers 71(4), 838–850 (2021)
  * [21] R. Conley, D. Choi, G. Medwig, E. Mroczko, D. Wan, P. Castillo, K. Yu, in _Quantum Computing, Communication, and Simulation III_ , vol. 12446 (SPIE, 2023), pp. 53–63