Enumeration Algorithm for the Shortest Vector Problem
Introduction
This document explains how the algorithm described in section 3.1 of [1] works. The explanation is based on the python implementation found in enumeration.py of that algorithm. The algorithm is the Schnorr-Euchner variant of the Kannan-Fincke-Pohst algorithm.
Lattice
A lattice $L$ is a discrete subgroup of $R^n$ seen as a module over $Z$. Thus a lattice can be described by a basis $b \in (R^n)^d$ if d is the dimension of the lattice. Given such a basis $b$ any vector $v \in L$ can be described by the coordinates $x \in Z^d$ in this basis : $v = \sum_i x_i b_i$
Since such a lattice is a discrete subgroup of $R^n$, there are shortest non-zero vectors in terms of the Euclidean norm. This length is called $\lambda$
SVP
Given a lattice, solving the shortest vector problem (SVP), means to find a shortest vector in this lattice, a vector with length $\lambda$.
The Algorithm
The algorithm inputs are a bound $A$ and the Gram-Schmidt coefficients of a lattice basis $b$. The output is a coordinate vector $x \in Z^d$ of a shortest vector in terms of this basis $b$. Note that the basis $b$ is not used by the algorithm. The algorithm calls a function update. The implementation of this function depends on whether the Gram-Schmidt and the computation in the algorithm are carried out exactly or rounded (i.e. floating point arithmetic).
Gram-Schmidt orthogonalization
The Gram-Schmidt coefficients required by the algorithm are:
$\mu_{i,j} =\frac{
The Gram-Schmidt orthogonalization $b^*$ of the basis $b$ is give by:
$b_i^ = b_i - \sum_{j<i} \mu_{i,j} b_j^$
Find Short Vectors
Let $L$ be a lattice of dimension $d$ having basis $b$. The algorithm enumerates the lattice elements in a clever way to find vectors being shorter or equal in length to the bound $\sqrt{A}$. Hence the equation to be solved is:
$| \sum_i x_i b_i |^2 \leq A$
let $b^*$ be the orthogonalisation of the basis $b$, then there is $y \in Z^d$ such that:
$\sum_i x_i b_i = \sum_i y_i b_i^*$
The coordinates in the orthogonal basis can be computed as:
$y_i = x_i + \sum_{j \geq i+1} \mu_{i,j} x_j$
Let $c_i = -\sum_{j \geq i+1} \mu_{i,j} x_j$, then we have $y_i = x_i - c_i$
Now the equation to be solved can be restated as:
$A \geq | \sum_i x_i b_i |^2 = | \sum_i y_i b_i^* |^2 = \sum_i y_i^2 r_i$
Let $l_i = \sum_{j \geq i} y_i^2 r_i$
The algorithm enumerates all coordinates $x$, by first choosing an $x_d$, then an $x_{d-1}$, ... and finally an $x_1$ and at each choice of an $x_i$ shifts it to $y_i = x_i - c_i$ and with the equations above checks if the $x$ construction can still lead to a solution by checking: $A \geq l_i$
Zig-Zag Path
The successive choices of the $x_i$ for a certain $i$ can be made clever by following the Schnorr-Euchner zig-zag path, which ensures that at each step $l_i$ increases. This means if $x_i$ is a choice and $x_i'$ is the next, then $l_i \leq l_i'$. Thus as soon as an $x_i$ leads to $A < l_i$ the subtree given by the choice of $x_{i+1}$ can not lead to more solutions and a new branch can be inspected by choosing a new $x_{i+1}$. Hence coordinates are arranged in a tree and at each node the algorithm can check:
- if $i = 1$ whether this branch is a solution
- else whether this subtree can still contain more solutions
In case 2. if the subtree is finished and $i = d$ the algorithm terminates.
Return Value
The process described above enumerates possible coordinates solving:
$l_1 = \sum_i y_i^2 r_i \leq A$
If now at each found solution, the bound $A$ is substituted with the length $l_1$ of this solution vector, then the last found solution is a shortest vector. The coordinate $x \in Z^d$ of this vector is then returned.
Implementation: enumeration.py
The input arguments are described above. The dimension of the lattice must be positive otherwise an error is raised.
def enumeration(A, m , r):
d = len(r)
if 0 >= d: raise EnumError("dimension of lattice is to small")
The algorithm is started with the coordinate vector $x = (1,0,0,...0)$, this represents the first basis vector $b_1$. A possible choice for the bound to call the algorithm on is $A \geq |b_1|^2$ (or $>$ in case of rounded Gram-Schmidt coefficients). This way this first $x$ will be selected as first possible solution. The zig-zag path is implemented with two vectors in $Z^d$:dx storing the last step width and ddx the last step direction for each coordinate index. The variable sol stores the coordinate of the last found short vector, which as explained above is the smallest up to "now". sol remains the zero vector as long as no vector shorter then $\sqrt{A}$ has been found.
# 1 Initialization
x = [1] + [0] * (d-1)
dx = [1] + [0] * (d-1)
ddx = [1] + [-1] * (d-1)
sol = [0] * d
The variables c and l are as described more above. l has length $d+1$ such that later in the code the general case can also be applied when $i = d$. ys stores the coordinates in respect to the orthogonal basis $b^*$, but squared (thus the "s"), since the $y_i$ are only needed in this form. The values of this vectors are determined by the initialization value of x.
# 2
c = [0] * d
l = [0] * (d+1)
ys = [0] * d
i and x[i:d] represent the state of the computation. In this state x[i:d] stores the chosen coordinates in this computation branch up to "now". loop counts the amount of performed loops for debug purposes and should normally be commented out. while(True) is the main loop of this program. The assert statement checks that x and sol remain integers and are not messed with the floating point values. ys and l are computed following the formulas described more above.
# 3 Main loop
i = 0
loop = 0; solution = False # DEBUG
while(True):
assert int is type(x[0]) and int is type(sol[0])
# 4 Compute next length
ys[i] = (x[i] - c[i])**2
l[i] = l[i+1] + ys[i] * r[i]
If i == 0 the algorithm is at the and of a branch, hence a whole coordinate vector has been constructed. If this coordinates represent a vector of length at most $\sqrt{A}$ this is a candidate to be stored in sol. The update1 does exactly that and updates also the bound A (TODO: more explanations on update1 -2 ... rounding floating points ...). solution is for debug purposes and should be commented out.
# 5 Found possible solution
if i == 0 and l[i] <= A :
(sol, A) = update1(sol, A, x, l[0])
solution = True # DEBUG
The following lines "visualize" the computation tree. This lines are for debug and demonstration purposes and should normally be commented out.
loop += 1 # DEBUG
s = "" # DEBUG
s += '{:>4}'.format(loop)+":" # DEBUG
for j in xrange(i): s+=" " # DEBUG
for j in xrange(i, d): s+='{:>3}'.format(x[j]) # DEBUG
if solution : s += " <"; solution = False # DEBUG
else: s += " " # DEBUG
s += '{:>3}'.format(i) # DEBUG
print s # DEBUG
If i > 0 but the value of l still admits a short vector in this computation branch, the algorithm descends in the tree: i-=1. Then computes the new shift $c_i = y_i - x_i$.
# 6 Build up coordinates
if i > 0 and l[i] <= A :
i -= 1
# 7
c[i] = 0
for j in xrange(d-1, i+1, -1):
c[i] += x[j] * m[j][i]
From $y_i^2 \leq \frac{A - l_{i+1}}{r_i}$ it follows that x_i must remain in the interval $[-\sqrt{\frac{A - l_{i+1}}{r_i}} + c_i , \sqrt{\frac{A - l_{i+1}}{r_i}} + c_i]$. This is an interval centered at $c_i$. The zig zag path of x_i starts at the center of this interval. This guaranties $y_i = 0$ (up to rounding) which leads to the smallest $l_i$. The first zig-zag step is thus zero, the step direction however depends on the rounding of $c_i$.
# 8
x[i] = int(round(c[i]))
dx[i] = 0
if c[i] <= x[i] :
ddx[i] = 1
else :
ddx[i] = -1
The computation tree is finished as soon as a choice of x[d-1] leads to ys[d-1] r[d-1] > A. At this point the last found sol is returned.
# 9 End of computation
elif i == d-1 and l[i] > A :
return sol
If none of the to cases (# 6 and # 9) are fulfilled it's time to change computation branch, the algorithm goes one level up in the computation tree. Then it computes the next coordinate in the zig-zag path by first inverting the step direction, then compute the next step (in opposite direction but one bigger) and finally compute the next x[i].
# 10 Next computation branch
else :
i += 1
# 11
ddx[i] = -ddx[i]
dx[i] = -dx[i] + ddx[i]
x[i] = x[i] + dx[i]
References
[1] X. Pujol, D. Stehle: "Rigorous and Efficient Short Lattice Vectors Enumeration"