These MATLAB 7.12.0 functions are all the functions required to run the LFS Algorithm as stated in "On Chubanov's Method for Linear Programming". After unzipping the files, if you put them all in the same directory, then you can run Chubanov_Preliminary.m with the proper input (explained below) and determine if the system has a feasible or determine it has no integer solutions. The rest of this README explains each function in detail. These functions are also used in the running of the D&C experiments from the Run_Experiment function. Make sure you have these functions linked to your current MATLAB directory in order to run the D&C as an approximation algorithm. Before proceeding, many of these explanations are written under the assumption you have at least read our paper "On Chubanov's Method for Linear Programming", for which these functions were created. This paper is also available at www.math.ucdavis.edu/~mjunod. ****************************************** LFS.m [sol,infeasible,time_fail] = LFS(A,b,C,d,r) Although this was originally designed to be run as a stand-alone function, it has since been integrated into the function Run_Experiment (see its README for more info). As such, it is easier to run specific examples from Run_Experiment. Given an n-dimensional system Ax = b, Cx <= d, where A has m rows and C has l rows, r must be a radius such that the feasible region of the system is contained in B(0,r), the ball of radius r centered at the origin. Given this input, this function either determines no strictly feasible solution exists, i.e. Ax = b, Cx < d is not feasible, or returns a feasible solution. In the first case, infeasible is set to 1, and either some h and delta are returned and fail is set to 0, or fail is set to 1 and h and delta are set to 0. Either way sol and flag are set to 0. In the second case sol is the feasible solution, and flag is set to 1. INPUT: A = matrix of equality constraints b = right-hand side of equality constraints C = matrix of inequality (<=) constraints d = right-hand side of inequality (<=) constraints r = positive number such that the feasible region is contained in the ball B(0,r) OUTPUT: sol = solution to the strictly feasible system if one exists, otherwise set to 0 infeasible = set to 1 if the system is not strictly feasible, otherwise set to 0 time_fail = 1 if the algorithm reached the time limit, otherwise set to 0 ****************************************** LFS_hoffman.m [sol,infeasible,time_fail] = LFS_hoffman(A,b,C,d,r) This is a version of the LFS algorithm that was used in the Hoffman experiments. It just contains a slight tweak to the original LFS algorithm. Instead of finding the maximum subdeterminant, it simply loads the proper .mat data file, which contains the maximum subdeterminant parameter. ****************************************** LFS_MIPLIB.m [sol,infeasible,time_fail] = LFS_MIPLIB(A,b,C,d,r) This is the version of the LFS algorithm used in our MIPLIB experiments. Rather than trying to find the maximum subdeterminant of these large problems (which could take weeks or years), we loop over progressively larger values of Delta_A. See our paper "On Chubanov's Method for Linear Programming" for more details and the hueristic justification. ****************************************** Max_SubDet.m subdet = Max_SubDet(A) Given an n by m matrix A, with n > m, this function finds the largest determinant of all the m by m submatrices of A. In order to avoid space issues, as even relatively small matrices can have many, many m by m submatrices, this function dos not take advantage of MATLAB vector operations. Instead it loops through all possible submatrices and iteratively stores the largest subdeterminant found so far. ****************************************** DandC_LFS.m [x,h,delta,flag,fail,time_fail] = DandC_LFS(n,m,l,A,b,C,d,z,r,error,c_max) Given an n-dim system Ax = b, Cx <= d, with m equality constraints, l inequality constraints, and a ball of radius r centered at z, for any r this finds either a hyperplane separating the ball and the feasible region, described by h and delta with flag and delta set to 0, or returns an error-approximate solution x with flag set to 1 and fail set to 0. This is called from the function LFS. The parameter c_max is determined by the function find_c_max and is described below. For a full explanation of how this recursive function works, see the paper referenced above. At the beginning of every recursive call of this function, we check if it has exceeded the time limit allotted. If so, time_fail is set to 1 (it is always 0 otherwise), and we immediately quit out of the algorithm with every other returned value set to 0. When called from LFS, the parameters of this function and all the functions that follow are populated automatically. ****************************************** EP_LFS.m [x,h,delta,flag] = EP_LFS(n,m,A,b,C,d,z,r) Given an n-dim system Ax = b, Cx <= d, with m equality constraints, l inequality constraints, and a ball of radius r centered at z, where r is small enough as determined in the D&C, this function either finds a supporting hyperplane of the feasible region that separates the ball from the feasible region, or determines that the ball contains an error-approximate solution, where the error is determined by the user. The D&C determines when r is small enough to enter this function. If a separating hyperplane (h, delta) is found, we return h and delta as well as flag = 0 and the zero vector for x. If a solution x is found, we return x and flag = 1, and the zero vector for h and delta = 0. ****************************************** fail_test.m fail = fail_test(h1,h2) Once the D&C has found two separating hyperplanes, and before it determines the best induced separating hyperplane, we need to check if h1 = -gamma * h2, for gamma > 0. If this happens, then the system is infeasible and we need to quit the recursive D&C. Given hyperplanes h1 and h2 in n-dim space, this function checks for such a gamma. If one is found, we return fail = 1, and otherwise return fail = 0. ****************************************** find_alpha.m [alpha] = find_alpha(h1,delta1,h2,delta2,z,r) This is the final step in each call to the D&C algorithm. Once the two different separating hyperplanes (h1,delta1) and (h2,delta2) are determined, this function calculates the alpha that produces the best induced separating hyperplane (h = alpha*h1 + (1 - alpha)*h2, delta = alpha*delta1 + (1 - alpha)*delta2). ****************************************** slack_version.m [A b C d] = slack_version(A,b,C,d) Given an n-dimensional system of the form Ax = b, Cx <= d, where A has m rows and C has l rows, this function introduces slack variables s to transform it into Ax = b, Cx + s = d, x >= 0, s >= 0. It then returns A as the new equality matrix, obtained by concatenating versions of A and C. The b that is returned is obtained in the same fashion. C and d form the new nonnegative inequalities. ****************************************** equiv_equals.m [A,b,r] = equiv_equals(A,b) Given a linear system Ax = b, this function produces an equivalent full-rank system. It returns the new matrix A, vector b, and r the rank of A. It uses fref.m to determine which rows of A can be kept and which should be discarded in the creation of the full-rank matrix. ****************************************** fref.m [A,v] = fref(A,tol) The algorithm is similar to MATLAB's RREF, but uses vectorized code for increased efficiency, and does not scale the rows to produce leading ones. Instead it simply does the first few steps of the RREF function, creating rank(A) linearly independent rows followed by rows of zeros. If the parameter is not set, it uses MATLAB's default tolerance in the computations. Otherwise it uses the tolerance tol. It is called by equiv_equals.m as [~,v] = fref(A). The purpose of this is not to actually create a row-reduced form of A, but simply to determine which rows need be removed from A in order to create a matrix with full rank. Then last m - rank(A) entries in v precisely tell us these rows that need to be removed. In this manner we create a full-rank matrix with integer entries much faster than the original rref function provided by MATLAB.