Using Seysen's Algorithm to Solve Subset Sum Problems

In both [26] and [33] the LLL algorithm was used almost exclusively to perform the required lattice basis reductions. Given the existence of Seysen's basis reduction algorithm and knowledge as to how it performs relative to LLL, it is natural to wonder how Seysen's algorithm could be applied to solving subset sum problems. We know from some of the comparisons performed between LLL and Seysen's algorithm in Chapter 2 that Seysen's algorithm tends to perform fewer row operations than LLL. If we used Seysen's algorithm then as part of an attack on subset sum problems, we should be able to perform more reduction operations (or other work) without increasing the overall execution time of our method. There are also theoretical and empirical reasons to suspect that Seysen's algorithm alone will not perform well on subset sum problems. Seysen's algorithm was originally suggested for simultaneously reducing a lattice and its dual lattice, not for finding the shortest vector in a lattice. If a lattice and its dual are of comparable size, the Seysen algorithm is not likely to perform row operations that generate short vectors in one lattice if that move increases the lengths of vectors in the dual lattice. If the vectors in the dual are significantly larger than those in the prime lattice, Seysen's algorithm may actually produce vectors larger than those originally input. Also, the algorithm's operation is not dependent on the ordering of the basis vectors, which means that the benefits gained by running LLL multiple times on different randomizations of the subset sum lattice disappear. These predicted benefits and deficiencies of the Seysen algorithm suggest that the best place to use Seysen is at the beginning of our attempt to solve a subset sum problem. The initial subset sum lattice L has primal basis vectors which are much larger than those in the dual; this favors the generation of shorter vectors in L and longer vectors in L^*. Also, as Seysen tends to perform many fewer row moves than LLL to reach the same reduced lattice basis, using Seysen's algorithm initially will reduce the total number of multiprecision row moves as compared to using only LLL. Once Seysen's algorithm stops at a local minimum, we can use LLL and other techniques to further reduce just the primal lattice L and look for short vectors.

  

Figure 3-4: Overview of Algorithm SL

Figure 3-4 shows the basic outline of the SL (Seysen-Lovász) algorithm. We start with the initial lattice basis suggested in [26]. If the weights of the subset sum problem are sufficiently large, then the lengths of the basis vectors in L will be much greater than the corresponding lengths of the basis vectors in the dual lattice L^*. We apply the Seysen algorithm to the L,L^* self-dual pair of lattice bases, using a greedy algorithm (Section 2.1.2) to choose at each step the pair of basis vectors $(\bold{b_i},\bold{b_j})$ to be reduced. The measure function S(A) is the one suggested by Seysen [38]:

$$S(A) = \sum_{i=1}^n \Vert\bold{b_i}\Vert^2 \Vert\bold{b_i^*}\Vert^2.$$

When computing the change in S(A) for each pair of vectors $(\bold{b_i},\bold{b_j}), i \neq j$, $\lambda$ is always set to its maximum value of:

$$\lambda = \left\lceil \left( \frac{a_{i,j}^*}{a_{j,j}^*} - \frac{a_{i,j}}{a_{i,i}} \right) \right\rfloor$$

Let L_S1 be the lattice basis which is produced by this (first) application of the Seysen algorithm. Recall that the goal of the SL algorithm is to find a basis vector which describes the solution to the given subset sum problem. In [26,33] the desired solution vector $\bold{e}$ was always of the form (e₁,...,e_n) with e_i equal to either 0 or a fixed constant $\kappa$. We will call all such vectors $\bold{e}$ Type-I solution vectors, to designate them from another class of solution vectors which will appear below. At the completion of the first Seysen reduction in the SL algorithm, it is possible that lattice basis L_S1 contains a Type-I solution vector. If so, then we are finished, and Algorithm SL halts. We assume that L_S1 does not contain a Type-I solution vector, and continue with the next stage of the algorithm. The algorithms of [26,33] searched for short vectors of the form $(e_1,\ldots,e_n,e_{n+1}=0)$ in the LLL-reduced bases, where $e_i \in \{0,\kappa\}$ for $1 \leq i \leq n$ and $\kappa \in {\Bbb Z}$ is any fixed integer. One of the problems with these methods is that often the short vectors produced by LLL reduction had $e_{n+1} \neq 0$; that is, the sum $\sum e_ia_i$ was not of the form s+yt and did not describe any relation involving the target subset sum. One method of reducing the probability that LLL (or Seysen) will include a vector of this form in the reduced basis is to scale the a_i and s by some constant factor N, which increases the length of any vector having a nonzero e_n+1 by about $\sqrt{N}$ if N is sufficiently large. However, this approach has the drawback that it increases the size of numbers which are already quite large and require multi-precision arithmetic. We suggest another method for elimination from consideration all lattice vectors with $e_{n+1} \neq 0$: GCD-reduction. The idea behind GCD-Reduction is to perform row moves on the lattice basis so that the entire $(n+1)^{\text{th}}$ column contains exactly one nonzero element, the GCD of the weights $a_1,\ldots{},a_n$. If (for example) vector $\bold{b_{n+1}}$ contains the GCD, we can remove $\bold{b_{n+1}}$ from the basis and remove the $(n+1)^{\text{st}}$ column completely from the lattice basis. This reduces what was an (n+1)-dimension lattice to an n-dimension lattice, and guarantees that any lattice vector generated by reducing this basis would have had its $(n+1)^{\text{st}}$ component equal to 0 in the (n+1)-dimension space. Implementing GCD-Reduction is easy. The basic algorithm we use was described by Brun [9]. Basis vectors are sorted in order of decreasing b_i,n+1. Then $\bold{b_2}$ is repeatedly subtracted from $\bold{b_1}$ until b_2,n+1 > b_1,n+1. The vectors are then resorted in order of decreasing b_i,n+1 and the process loops. In reality, only vector $\bold{b_1}$ needs to be inserted into the previously generated order. This can be performed very fast by using an auxiliary pointer array; the actual basis vectors don't even need to be moved around. Eventually, the only nonzero element in the $(n+1)^{\text{st}}$ will be b_1,n+1, at which point both the vector $\bold{b_1}$ and all $b_{i,n+1}, 2 \leq i \leq n+1$ can be removed from the lattice. We apply GCD-Reduction to the lattice basis L_S1 output by the first application of Seysen's algorithm. We do not use GCD-Reduction on the original lattice basis L because the resulting lattice basis matrix would be quite dense (which means the Seysen algorithm will take longer to run) and the vectors in the dual lattice would also be much larger. After GCD-Reduction is applied to basis L_S1 (yielding output basis L_G in n dimensions) we again search for the presence of a Type-I solution vector. Assuming we do not find such a vector, the Seysen algorithm is applied to L_G to reduce the lattice to a local minimum of the measure function S(A). Call this lattice L_S2, the output of the second Seysen application. For small values of n ( $n \lesssim 20$), it is possible that L_S2 will contain a Type-I solution vector. For larger values of n and sufficiently high densities d, however, it is unlikely that the Seysen-GCD-Reduction-Seysen portion of the algorithm will have found the desired vector. Thus we begin the second portion of the algorithm: the LLL phase.

  

Figure 3-5: Algorithm SL: LLL-Phase

Figure 3-5 shows the first level of detail in the LLL-phase of Algorithm SL. The input to this phase of the algorithm is the lattice basis L_S2 which was the result of the second application of Seysen reduction. We now take advantage of the theoretical improvements in the density bound given in [12] and Section 3.2 above. (The introduction of the ``one-half'' vector is delayed until after the Seysen stage of Algorithm SL has completed to increase the sparsity of the lattice bases which are Seysen-reduced.) We extend the ${\Bbb R}^n$ lattice described by L_S2 into a lattice in ${\Bbb R}^{n+1}$ by adding an $(n+1)^{\text{st}}$ component to each of the basis vectors $\bold{b_1},\ldots{},\bold{b_n}$:
$$\begin{align*}\text{old } \bold{b_i} & = (b_{i,1},b_{i,2},\ldots{},b_{i,n}) \\ ... ..._{i,1},b_{i,2},\ldots{},b_{i,n},b_{i,n+1} = \sum_{j=1}^n b_{i,j}) \end{align*}$$
The new $(n+1)^{\text{st}}$ component of each basis vector is simply the sum of the first n components. To complete the extension, we need to add one more vector to the lattice basis. New basis vector $\bold{b_{n+1}}$ is defined as follows:

$$\bold{b_{n+1}} = (\underbrace{\tfrac{1}{2},\tfrac{1}{2},\ldots{},\tfrac{1}{2}}_{n\text{ terms}},\tfrac{1}{2}(n-1))$$

For $\bold{b_{n+1}}$, the $(n+1)^{\text{st}}$ component is the sum of the first n terms minus $\tfrac{1}{2}$. The $-\tfrac{1}{2}$ correction is needed because otherwise the $(n+1)^{\text{st}}$ column of the basis would be dependent. Let $L_{\text{in}}$ designate this extended lattice. With the introduction of vector $\bold{b_{n+1}}$ we have also introduced a second class of solution vectors. If $\bold{e}$ is a Type-I solution vector for a subset sum problem with $\tfrac{1}{2}n$ elements in the subset, then there now exists a Type-II vector in the lattice of the form

   $$\bold{\hat{e}} = (\hat{e}_1,\ldots{},\hat{e}_n,\hat{e}_{n+1})... ...e}_i \in \{-\tfrac{1}{2},\tfrac{1}{2}\}, 1 \leq i \leq n+1.$$

Also, if $\bold{e}$ was a solution vector of the n-dimensional lattice, in the extended lattice $L_{\text{in}}$ vector $\bold{e}$ will be transformed into $\bold{e'} = (e_1,\ldots{},e_n,e_{n+1})$ where e_n+1 equals $\kappa$ times the number of elements in the subset which sums to s. Thus, during the LLL phase of Algorithm SL, we search for both extended Type-I solution vectors and Type-II solution vectors. The structure of the LLL phase of Algorithm SL is based upon the successful attacks of [26,33]. Lagarias and Odlyzko showed that using LLL multiple times on random orderings of the input basis improved the chances of finding a solution vector. Thus, Algorithm SL initially tries to reduce lattice $L_{\text{in}}$, and upon failure randomizes $L_{\text{in}}$ and tries again. This process continues until a total of $\pi_1$ reduction attempts have been made ( $L_{\text{in}}$ and $\pi_1 - 1$ random orderings of $L_{\text{in}}$). In [26] $\pi_1 = 5$, but there is nothing special about that particular value.

  

Figure 3-6: Algorithm SL: LLL-Loop Structure

Recall from Section 3.3 above that the Radziszowski-Kreher approach, while setting $\pi_1 = 1$, repeatedly called LLL in conjunction with Weight-Reduction and Sort (Sort just sorts the basis vectors by length.). They showed that repeated iterations of the LLL-Weight-Reduction-Sort loop shown in Figure 3-2 yielded better results than a single application of LLL. Algorithm SL incorporates this approach, applying a number of iterations of LLL-Weight-Reduction- Sort before giving up and trying a new randomization of $L_{\text{in}}$ (Figure 3-5). Each ordering of $L_{\text{in}}$ passes though up to $\pi_2$ iterations of the loop before Algorithm SL gives up. As stated above, in [33] $\pi_2 = 15$, or $\pi_2 = 9$ if all the vectors in the lattice were of length $< \tfrac{1}{2}n$. Again, there is nothing special about these particular values. Theoretically, one could choose $\pi_1$ and $\pi_2$ to be quite large. Any practical algorithm, though, would probably want to use reasonable constants for both $\pi_1$ and $\pi_2$. Although both [26] and [33] run the LLL algorithm with the parameter $y \approx 0.99$, Lenstra, Lenstra, and Lovász show that y may be any value in the range $\tfrac{1}{4}< y < 1$ [27]. For $y \approx \tfrac{1}{4}$, LLL will only exchange vectors in the lattice (Equation 1.2) for relatively large differences between vectors $\bold{b_i}$ and $\bold{b_{i+1}}$. As $y \rightarrow 1$, the amount of improvement required to trigger the LLL exchange step decreases, and LLL will swap vectors and continue running for minimal improvement. Thus, larger values of y will likely lead to improved results, but LLL will also take longer to run. This suggests that instead of running the entire LLL algorithm with y = 0.99, as was done previously, a version of LLL which started with $y \approx \tfrac{1}{4}$ and adjusted it upwards as the algorithm ran might work just as well but require fewer row operations [18].

  

Figure 3-7: Algorithm SL: LLL(x) Internal Structure

Algorithm SL uses a version of the LLL algorithm which varies the parameter y among a finite number of values. Figure 3-6 shows what happens each time SL calls LLL-Loop. For subset sum problems with n < 90, each call to LLL passes through six stages (i.e. six distinct values for y)^3.1. Upon entry to each stage, the y parameter is updated, and LLL-reduction is performed upon the lattice basis using the new y value (see Figure 3-7). Let $L_{y_0,\text{in}}$ and $L_{y_0,\text{out}}$ represent respectively the lattice bases input to and output from an application of the LLL algorithm with y = y₀. We now compute the boolean value of the following expression:

$$\begin{split} \left(\vphantom{\min}\right. & \left. \min\{\V... ...L_{y_0,\text{in}}\}\right)\right). \end{split} \end{split}$$ (3.7)

(This decision point is represented by the diamond in Figure 3-7.) If the boolean value of Equation 3.7 is TRUE, then the Weight-Reduction procedure is applied, the Sort procedure is invoked, and the output is input again into LLL with y=y₀. Thus, we perform LLL-Weight-Reduction-Sort loops until either the length of the shortest vector in the lattice has increased after LLL-reduction, or if the length of the shortest vector has remained constant and the length of the longest vector has not decreased^3.2. By using the recursive structure in Figure 3-7 and the termination condition represented by Equation 3.7, we attempt to have the LLL algorithm perform as many reduction steps as possible with small values of y₀. We delay increasing the value of y until LLL fails to make any improvement in the length of the largest vector in the lattice basis. In this way we restrict the set of row swaps and moves which LLL will consider, which improves the running time of the algorithm. Initially, for y small (say y < 0.75) LLL will only consider row moves which will ``significantly'' reduce the current lattice basis. Later, for values of $y \approx 1$, LLL will consider any row move which reduces the lattice. For n < 90, six specific y values are used: 0.2578125, 0.625, 0.75, 0.875, 0.9375, and 0.9921875. These values were chosen for two reasons. First, all of the numbers have exact fractional binary representations in double-precision floating point. This meant that error would not be introduced into the LLL algorithm by performing arithmetic operations on y. Second, experiments with small subset sum problems ($n \leq 24$) showed that the number of row moves performed by LLL for each of the four middle values were approximately equal. That is, work was evenly distributed across all intermediate stages represented in Figure 3-6. (The endpoint values 0.2578125 and 0.9921875 were fixed). For $n \geq 90$, three additional values were used: 0.34375, 0.4375, and 0.53125. These values evenly divide the interval [0.2578124,0.625] into four equal pieces. As discussed in Section 3.5 below, our initial implementation of the LLL algorithm did not work correctly for lattice bases with $n \geq 90$ because of rounding error introduced into the computations. We were able to reduce the effects of rounding error to acceptable levels by modifying the LLL algorithm itself and by introducing the 0.34375, 0.4375 and 0.53125 stages into the algorithm. We have detailed the operation of Algorithm SL, a combination of the Seysen and the Lenstra-Lenstra-Lovász basis reduction algorithms which also utilizes the GCD-Reduction, Half-Vector, Weight-Reduction, and Sort heuristics. In the next section we present the results of our experiments with this algorithm and show how the combination of these techniques greatly increases the range of subset sum problems which may be solved.