Week 11 - 'pencils down'

Sorry for not updating this blog for the last couple of weeks. I forgot about this after the semester began. I'll use this blog post to summarize the project progress in this month. Before that, however, I'd like to note that the current module does indeed perform LLL, LLL with removals and ULLL which was the primary goal of the project.

So, starting off where I left in the previous post, the augmentation of the identity to the basis matrix actually plays a part in bettering the numerical stability of ULLL. Hence, I modified the ULLL function to be similar to the one in flint-1.6.
Also, the default LLL functions (which are essentially wrappers for ULLL) were added. The LLL-reducedness checkers were also changed to implement the certification method introduced by Villard. It will not always work, but if it does the result is guaranteed. It is efficient, and much faster than using the provable precision in the L^2 algorithm or using exact arithmetic for testing.

The Gram matrix version of the provable is_reduced functions was also added.
This was followed by an implementation of a modified LLL reduction algorithm due to Storjohann which has better complexity in terms of the lattice dimension. However, it didn't fare well against the other LLL implementations (specifically classical LLL) during profiling. The mpf versions of the is_reduced functions were removed because they weren't proven and mpfr versions were written.
Adding a function which uses mpf and epsilons to emulate the behaviour of changing the rounding mode is a possibility as mentioned on flint-devel. Finally, the documentation was updated and now should (hopefully) explain about all the strategies and types used in this module (along with the function descriptions of course).

I'd like to thank my mentors, William Hart and Curtis Bright, for their valuable time, support and patience. This project wouldn't have been possible without their help. I'd also like to thank Fredrik Johansson for reviewing some of my code and Burcin Erocal from the lmonade organization for his advice. I'm grateful to the Google OSPO team for running the GSoC program. This has been a great summer and contributing to FLINT has been a tremendous learning experience for me.

Week 10

This week, I removed the redundant store_trans parameter from the context object. Now, the capturing matrix is always updated unless the user passes NULL as the value. Also, in the special case that the user inputs a d x d identity and the number of columns of the input basis is greater than the number of rows (i.e. the embedding dimension is greater than the lattice dimension), I avoid updating the vectors during reduction itself and instead do a matrix multiplication at the end.

Also, as the mpf and wrapper functions of the LLL subroutines are supposed to guarantee that the output is indeed LLL-reduced, we need to check this in the most efficient way. This means that a fp test for reducedness should be used prior to using the exact arithmetic version (which is slower). Thus, I've added is_reduced functions in the module which first test using doubles, then mpfs if the matrix is certified non-reduced in the first test and finally, fmpq.

An initial implementation of the ULLL function was also added. It is different from the one in flint-1.6 because it does not perform any adjoining operations to the matrix. Instead, the option to store the unimodular transformations is utilised here. Also, the original did not use recursion on the truncated data, but the current code does. However, it needs to be tested.

This week, I plan to add test code and documentation for the ULLL function and document the functions for checking LLL-reducedness.

Week 9

This week, I updated the removals code to use a heuristic lower bound, mentioned by Curtis on flint-devel, on the final GS norm while removing the last vector during the execution of the reduction algorithm (after a failed Lovasz test involving kappa = d - 1). This was required because the proof of a theorem in the L^2 paper assumes that the final vector is LLL-reduced while deriving the error bound on the accuracy of the norm. This, of course, doesn't matter in the case of standard LLL but matters here because we need to be sure about the accuracy, lest we may remove something useful.

Besides this, I also unified the code for the cases where the exact Gram matrix is computed or input as mentioned in the todo section of my previous post. This required factoring out the row exponents from the Gram matrix rather than the basis itself, because the latter is not always be available (i.e. fl->rt == GRAM).

Also, I changed the code dealing with the matrix capturing the unimodular transformations to not make any assumptions regarding its column dimension. Earlier, I was assuming U to be a d x d matrix which was to be updated to satisfy the relation B* = UB where B* is the basis obtained by LLL-reducing B, i.e. U was the change-of-basis matrix. However, now U can be any matrix with the same number of rows as B.

I think the main implementation of LLL and LLL with removals is completed now, modulo a few (hopefully minor) changes. So, I plan to at least start working on the ULLL function this week.

Week 8

The eighth week of GSoC involved a lot of debugging to find the reasons for the bad behaviour of the removals function. In theory, the last vector can be removed from the basis if its squared GS length becomes greater than the bound at any point during the execution of the algorithm. However, it is not so straightforward in practice if the GSO is known only approximately. I think this may be a reason why the version in flint-1.6 removed vectors at the end as the algorithm seems to show much better behaviour in this case, i.e. the norms were more accurate.

I added the wrapper function for  LLL with removals optimised for knapsack lattices. Test code for this was also written.  Knapsack LLL differs from the textbook version in the fact that it performs early size reductions on the input basis occasionally which speeds up things in the knapsack case. Speaking of LLL with removals, the code was modified to remove the numerical inaccuracies plaguing it. Earlier, due to a floating point approximation of the Gram Schmidt orthogonalisation being used, the norm was incorrectly flagged as being greater than the bound which led to the removal of some useful vectors. The documentation was also updated and is now up to date with the module. Thus, the prerequisites for ULLL have now been completed.

This week I plan to unify the code for performing LLL on the input matrix in the cases that it is a lattice basis and an exact Gram matrix is to be used for computing the GSO or it is the Gram matrix itself.

Week 7

This week, I added support for LLL with removals on Gram matrix. This finds application in vector rational number reconstruction. Also, the inaccuracy checks (incorrect Lovasz tests and too many Babai loops) were improved to be similar to those used in fpLLL.

Also, the LLL function specialized to knapsack-type lattices, for which the Babai functions were implemented last week was added. It performs early size reductions which tend to make things faster for knapsack problem type lattice bases. It is also a prerequisite for the ULLL with removals function. Another feature implemented was early removal of vectors in LLL with removals. Now, the vectors whose squared GS lengths are greater than the input bound are removed from the basis during the reduction algorithm itself to avoid unnecessary overhead involved in keeping them updated during the computations.

This week, I plan to write the wrapper function for knapsack LLL with removals and fix any loose ends which may be remaining. The plan is to get the existing module to ready before I start working on the actual ULLL algorithm itself.

Week 6

Well, this week's update is late. Sorry for that. I'll try to summarize the preceding week here. I added the heuristic version of LLL with removals, along with it's test code and documentation. In the MPFR version in flint-1.6, the squared GS length is divided by 8 instead of 2 as mentioned in the comments. I don't know if this an overlook or extra precaution. If the latter, I see no reason for this though. As I mentioned in my previous post, however, I avoid division altogether. LLL with removals was completed when its arbitrary precision variant was added and the wrapper function written. Its return value is the new dimension of the basis to be considered for further computation.

This brings us to the third major and perhaps, the most important part of this project: implementing ULLL with removals. This is because of its property of sub-quadratic time complexity in the size of the entries. Also, it is numerically more stable. It isn't mentioned in literature and hence, Bill graciously offered to write a paper on this for reference. Before I start work on the actual ULLL function, however, I need to implement an LLL with removals optimised for knapsack type lattices as it is used in ULLL. This requires a few Babai-like functions as well. These procedures will only reduce the kappa'th vector against the vectors upto cur_kappa (which is an index before which the basis is assumed to be LLL reduced) and not kappa - 1. The Babai functions added differ in their way of computing the dot products like the former versions

Along with the progress reported on the mailing list, one important thing to do now is to update the documentation which is lagging behind. I hope I'll be able to find time to do this task this week.

"The documentation needs documentation."   -- a Bellevue Linux Users Group member, 2005

Week 5

The previous week was quite interesting. I wrote the wrapper function fmpz_lll_wrapper and documented it. Test code for it was also added. Thus, the first milestone of the project was reached. Improvements to the module implemented thus far were made according to mentor comments. In particular, is_reduced functions and fmpz_mat_lll were improved to avoid storing GS vectors as suggested by Curtis on the mailing list. The function fmpz_lll_wrapper should now provide functionality identical to fpLLL and flint-1.6, but an additional feature in this version is that the Gram matrix can also be supplied for reduction instead of the lattice basis. It should be useful because the only other software which I've seen that allows passing a Gram matrix as input is Magma, which isn't open source.

The next step is LLL with removals. There seem to be 2 definitions for LLL_with_removal in literature: Both versions accept a basis (say B) and a bound (say N). Now, the intuitive definition seems to be that B is LLL reduced with bound N if for every vector b in B, norm(b) <= N. However, according to these papers by Mark van Hoeij, et al., it actually deals with the length of the Gram-Schmidt (GS) vectors i.e. the last vector is removed if its GS length is greater than N. Of course, for computational simplicity (and accuracy) I accept the squared norm bounding the GS length. Another point to be observed is that in flint-1.6, the bound is compared with half the squared GS length, probably because a fp approximation of the GS lengths is used. This is done to avoid removing something valuable. Also, the documentation is a bit ambiguous as it mentions that the bound is for the "target vectors". I'm going with van Hoeij's definition because he mentions it in the context of factoring polynomials which applies to a possible use of LLL in FLINT.
I am not sure if LLL with removals requires a version for the Gram matrix as input, as the only mentions of it in literature relate to factoring polynomials which input a lattice basis to the procedure. So, I haven't written a Gram matrix variant for now. Although, I may implement it later, if it's good to have.
My version of LLL_with_removal works even when I directly compare the bound with the squared GS norm because I avoid conversion to doubles and instead compare fmpz_t's. This was ascertained from the test code for fmpz_lll_d_with_removal. Documentation for fmpz_lll_d_with_removal and fmpz_mat_is_reduced_with_removal in the corresponding modules was added.

This week, I plan to write the fmpz_lll_d_heuristic_with_removal function, along with its test code and documentation. If things are smooth, maybe I can even work on the arbitrary precision variant.

Week 4

This last week has been a productive one. I implemented check_babai_heuristic (the Babai version using mpfs) and documented it.  Helper functions for this were also implemented in the fmpz and fmpz_vec modules. The fmpz_lll_mpf2 function was also written. The "2" in the name of the function signifies that it takes the precision to be used for storing the temporary variables and the GSO data (also the approximate Gram matrix if fl->rt == APPROX) as arguments, or in other words mpf_init2 is used for initialising any intermediate mpf_t's used. The wrapper for this is fmpz_lll_mpf which increases the precision until the LLL is successful or God forbid, the precision maxes out. Test code for lll_mpf was added. Also, the test code now uses the fmpz_mat versions of is_reduced and is_reduced_gram which use exact arithmetic to check if a basis is LLL reduced and help cover edge cases.

Also, there's some good news to report. The lll_d and lll_d_heuristic functions now work on all test cases without failure for Z_BASIS as input matrix! With GRAM, they fail sometimes due to inadequate precision of double to store large values. I can confirm that fmpz_lll_d works for the 35 dimensional lattice on the web page for the L^2 paper by Nguyen and Stehle (tested against fplll-4.0.4). I also tested fmpz_lll_mpf with the 55-dimensional lattice which makes NTL's LLL_FP (with delta=0.99) loop forever. It works! :) The output produced is the same as that if the matrix is passed to fpLLL with the "-m proved" option.

This week, I look forward to writing the LLL wrapper fmpz_lll_wrapper and documenting it besides improving the part of the module implemented so far and fixing bugs, if any. I also plan to document those functions which were left undocumented and shift fmpq_mat_lll to the fmpz_mat module and rename it to avoid confusion.

Week 3

Henceforth, starts another exciting week of coding. This week I look forward to implementing arbitrary precision floating point LLL using the mpf_t data type provided by GMP (which is dependency for installing FLINT). According to the discussion on the mailing list, as it isn't a big issue, I've decided to get the code working with mpf for now as some helper functions are already implemented in the mpf_vec and mpf_mat modules (The other option is Arb which will be used eventually).

Last week, I wrote the heuristic version of LLL (which always uses heuristic dot products that attempt to detect cancellation) over doubles and worked on improving the textbook L^2 and Gram matrix variants. I think the textbook L^2 implementation needs more work before it can be proclaimed ready for use. Though it works on most of the test cases, there are still some corner cases to consider. For use in fmpz_lll_check_babai_heuristic, some changes were also made to the mpf_vec, fmpz, fmpz_vec modules to help facilitate conversion between fmpz and mpf and vice versa. Also, _mpf_vec_dot2 now signals possible cancellation through an integer return value.

Week 2

Two weeks have passed since GSoC started and it has been an eventful journey so far. I've implemented the fmpz_lll_d function which performs LLL over doubles on a basis matrix with integer (fmpz_t) entries using either an approximate Gram matrix (having floating point entries) or an exact Gram matrix (with multiprecision integer entries). An LLL context object has been created for supplying LLL parameters and other options to the function based on Bill sir's suggestion on the mailing list.

I wrote test code for fmpz_lll_heuristic_dot and thought the tests pass, I still have to try out the different methods (including using _d_vec_dot_heuristic) to see which one is faster and gives more reliable results. A few helper functions were implemented use in the fmpz_lll_d and its test code. Two such functions are fmpz_mat_get_d_mat for converting an fmpz_mat into a d_mat and fmpz_mat_rq_d for RQ decomposition in the fmpz_mat module. However, these require including d_mat.h in fmpz_mat.h leading to a lot of compiling every time d_mat.h is modified. I don't know if this is an issue or not.

One thing I noticed when testing the LLL function was that when the approximate Gram matrix is used, all tests pass when using randntrulike, randajtai and randsimdioph (with the bitlengths specified in the test code, of course). However, there is a need to shift to arbitrary precision for some cases with randintrel (bits <= 200). However, when the exact Gram matrix of the lattice basis is used, randntrulike fails sometimes as well. Is this to be expected, or is there an error in the implementation?

This week I plan to make changes to the code according to the mentor comments and write fmpz_lll_d_heuristic function, test code and documentation.

GSoC :- Week 1

My proposal "Implementing the LLL algorithm in FLINT" has been accepted for Google Summer of Code 2014. The official coding period began last week (May 19) and so I have started to code. The source code can be found here.

What did you work on this week?

- Implemented the _fmpz_vec_dot function in fmpz_vec module for use in fmpz_lll_heuristic_dot.

- Implemented the fmpz_lll_heuristic_dot function for use in fmpz_lll_check_babai_heuristic_d.

- Implemented the _fmpz_vec_get_d_vec_2exp function in fmpz_vec module for use in fmpz_lll_check_babai.

- Wrote the fmpz_lll_check_babai function.

- Wrote the fmpz_lll_check_babai_heuristic_d function.

- Made changes to d_mat, d_vec modules and fmpq_mat implementation of classical LLL based on mentor comments.

What do you plan on doing next week?

- Change fmpz_lll_heuristic_dot to use the new _d_vec_dot heurstic function and test it.

- Document the Babai functions.

- Write the fmpz_lll_d function.

- Test, profile and document fmpz_lll_d.

- Increase the pace, maybe? :)

Are you stuck on anything?

I am not stuck because of these issues, but would certainly appreciate comments to improve the  code.

- How to check if a double is NaN?: IEEE standard specifies f != f will be true only if f is NaN and so, I use this as the check. However, compilers may optimize this out. This requires checking with different compilers. The condition cannot be changed to use isnan function in math.h as that was introduced in C99.

- Will +-0 cause a problem? -0.0 == 0.0 returns true, so it seems to be okay.

- Hackish way of finding ulps in _d_vec_dot_heurstic. Is there a better way not using functions like nextafter which is since C99?

- y = rint(x); is not ANSI C. I used if(x < 0) y = ceil(x-0.5); else y = floor(x+0.5); instead. Here's the output on my machine with the rounding mode as default.
rint(5.5) = 6.0 rint(-1.5) = -2.0
ceil(5.5-0.5) = 5.0 floor(5.5+0.5) = 6.0 ceil(-1.5-0.5) = -2.0 floor(-1.5+0.5) = -1.0  

Introduction

Hi,
My name is Abhinav Baid. I am a computer science undergraduate at Birla Institute of Technology and Science, Pilani. This is a blog post to go along with my proposal and code sample for the project "Implementing the LLL algorithm in FLINT" as a part of the Google Summer of Code under lmonade. FLINT is a C library for number theory. The aim of this project is to implement a basic LLL in FLINT allowing for parameters to be supplied governing the strength of reduction (delta-eta reduction using floating-point arithmetic), followed by a couple of the more interesting modern versions, including the LLL with removals and ULLL, a version of LLL with better complexity in terms of the size of the entries. The mentors for this project are William Hart and Fredrik Johansson.

A description of the algorithms involved in this project can be found in my proposal. A more comprehensive list of references is available in this thread of the flint-devel mailing list. The proposal also includes an Objective and Deliverables section where I describe how I plan to approach this task. The data types used in this project will be double, mpf_t, and fmpz_t. Unlike flint-1.6, mpf was chosen over mpfr because correct rounding is not required for LLL. The modules which will be added/modified as part of this project are d_mat, d_vec, mpf_mat, mpf_vec. fmpq_mat, fmpz_lll.

• d_mat : Matrices with double precision floating point entries. Functions for operations like memory management, arithmetic, etc. BLAS compatible representation will be used.
• mpf_mat : Matrices with arbitrary precision floating point entries. Functions for operations like memory management, arithmetic, etc.
• d_vec : Vectors with double precision floating point entries. Functions to help modularise the code in the d_mat module.
• mpf_vec : Vectors with arbitrary precision floating point entries. Functions to help modularise the code in the mpf_mat module.
• fmpq_mat : Function to implement the rational algorithm (delta / c reduction).
• fmpz_lll : LLL reduction on the rows of an fmpz_mat_t. Functions related to LLL like the various Babai's, actual LLL (L^2) functions, wrappers, LLL with removals,  ULLL.

I plan to finish the basic helper functions (not related to the algorithm) and the fmpq_mat implementation before the coding period of GSoC. The LLL implementations in FLINT 1.6 do not support passing of parameters and are not compatible with the current FLINT series (FLINT 2.x). Hence, the project aims to build a structured and clean module instead of having all the functions in a single file. The proposal also has a detailed timeline where I provide a schedule with dates and important milestones.

Now, what attracts me to this problem? Well, there are lots of reasons but just for fun, let me enumerate a few:

• It involves C programming.
• It involves linear algebra.
• It has applications in diverse fields.
• It is open source.
• Most importantly, it has excellent mentors - computational number theory experts with lots of experience in implementing efficient algorithms.

If a project had even one of the above features, I would have loved to be a part of it. Now, when I get a (linear (:) combination of all the above plus a chance to flaunt, I'll surely be interested, won't I?