Implementing LLL algorithm as been said on Wikipedia, but getting into serious issues

Question

I am not sure my issue is related to programming or related to concept of LLL algorithm and what has been mentioned on Wikipedia.

I decided to implement LLL algorithm as it has been written on Wikipedia (step-by-step / line-by-line) to actually learn the algorithm and make sure it is truly working but I am getting unexpected or invalid results.

So, I used JavaScript (programming language) and node.js (JavaScript engine) to implement it and this is the git repository to get the complete code.

Long story short, value of K gets out of range, for example when we have only 3 vectors (array size is 3, thus maximum value of index would be 2), but k becomes 3 and it is nonsense.

My code is step-by-step (line-by-line) implementation of the algorithm mentioned on Wikipedia and what I did was only implementing it. So I don't what is the issue.

// ** important
//  {b} set of vectors are denoted by this.matrix_before
//  {b*} set of vectors are denoted by this.matrix_after
calculate_LLL() {
    this.matrix_after = new gs(this.matrix_before, false).matrix; // initialize after vectors: perform Gram-Schmidt, but do not normalize
    var flag = false; // invariant
    var k = 1;

    while (k <= this.dimensions && !flag) {
        for (var j = k - 1; j >= 0; j--) {
            if (Math.abs(this.mu(k, j)) > 0.5) {
                var to_subtract = tools.multiply(Math.round(this.mu(k, j)), this.matrix_before[j], this.dimensions);
                this.matrix_before[k] = tools.subtract(this.matrix_before[k], to_subtract, this.dimensions);

                this.matrix_after = new gs(this.matrix_before, false).matrix; // update after vectors: perform Gram-Schmidt, but do not normalize
            }
        }

        if (tools.dot_product(this.matrix_after[k], this.matrix_after[k], this.dimensions) >= (this.delta - Math.pow(this.mu(k, k - 1), 2)) * tools.dot_product(this.matrix_after[k - 1], this.matrix_after[k - 1], this.dimensions)) {
            if (k + 1 >= this.dimensions) { // invariant: there is some issue, something is wrong
                flag = true; // invariant is broken
                console.log("something bad happened ! (1)");
            }

            k++;
            // console.log("if; k, j");
            // console.log(k + ", " + j);
        } else {
            var temp_matrix = this.matrix_before[k];
            this.matrix_before[k] = this.matrix_before[k - 1];
            this.matrix_before[k - 1] = temp_matrix;

            this.matrix_after = new gs(this.matrix_before, false).matrix; // update after vectors: perform Gram-Schmidt, but do not normalize

            if (k === Math.max(k - 1, 1) || k >= this.dimensions || Math.max(k - 1, 1) >= this.dimensions) { // invariant: there is some issue, something is wrong
                flag = true; // invariant is broken
                console.log("something bad happened ! (2)");
            }
            k = Math.max(k - 1, 1);

            // console.log("else; k, j");
            // console.log(k + ", " + j);
        }

        console.log(this.matrix_before);
        console.log("\n");

    } // I added this flag variable to prevent getting exceptions and terminate the loop gracefully

    console.log("final: ");
    console.log(this.matrix_before);
}

// calculated mu as been  mentioned on Wikipedia
// mu(i, j) = <b_i, b*_j> / <b*_j, b*_j>
mu(i, j) {
    var top = tools.dot_product(this.matrix_before[i], this.matrix_after[j], this.dimensions);
    var bottom = tools.dot_product(this.matrix_after[j], this.matrix_after[j], this.dimensions);

    return top / bottom;
}

Here is the screenshot of the algorithm that is on Wikipedia:

---

Update #1: I added more comments to the code to clarify the question hoping that someone would help.

Just in case you are wondering about the already available implementation of the code, you can type: LatticeReduce[{{0,1},{2,0}}] wolfram alpha to see how this code suppose to behave.

Update #2: I cleaned up the code more and added a validate function to make Gram Schmidt code is working correctly, but still code fails and value of k exceeds number of dimensions (or number of vectors) which doesn't make sense.

Answer 1

The algorithm description in Wikipedia uses rather odd notation -- the vectors are numbered 0..n (rather than, say, 0..n-1 or 1..n), so the total number of vectors is n+1.

The code you've posted here treats this.dimensions as if it corresponds to n in the Wikipedia description. Nothing wrong with that so far.

However, the constructor in the full source file on GitHub sets this.dimensions = matrix[0].length. Two things about this look wrong. The first is that surely matrix[0].length is more like m (the dimension of the space) than n (the number of vectors, minus 1 for unclear reasons). The second is that if it's meant to be n then you need to subtract 1 because the number of vectors is n+1, not n.

So if you want to use this.dimensions to mean n, I think you need to initialize it as matrix.length-1. With the square matrix in your test case, using matrix[0].length-1 would work, but I think the code will then break when you feed in a non-square matrix. The name dimensions is kinda misleading, too; maybe just n to match the Wikipedia description?

Or you could call it something like nVectors, let it equal matrix.length, and change the rest of the code appropriately, which just means an adjustment in the termination condition for the main loop.