A parallel linear programming solver method based on rlibm

Abstract

The application discloses a parallel linear programming solver method based on RLIBM and belongs to the technical field of computer data processing; the method comprises the following steps: designing a processing link of "convex hull intermediate representation constraint partition-parallel weighted random sampling-maximum consensus solving-parallel verification-iterative optimization", mainly solving the problems in solving large-scale infeasible low-dimensional linear programming (LDLP) encountered in the RLIBM algorithm, and introducing a geometric duality principle to project high-dimensional constraints to a two-dimensional plane and construct a convex hull intermediate representation; in the method, a feasible subset and an infeasible superset are divided by identifying key constraints corresponding to convex hull vertices, the problem scale is effectively reduced, a parallel weighted random sampling strategy is designed, key constraints are efficiently sampled by dynamically adjusting constraint weights and parallel calculation, a maximum consensus solving algorithm is combined to maximize the number of satisfied constraints in the infeasible superset, and a high-quality solution is gradually approached through a predictor-correction iterative mechanism.

Description

Technical Field

[0001] This invention relates to a parallel linear programming solver method, specifically a parallel linear programming solver method based on RLIBM, belonging to the field of computer data processing technology. Background Technology

[0002] With the increasing demand for high-precision mathematical functions in scientific computing and engineering applications, the generation of correctly rounded elementary functions has become an important research area. However, currently available solution methods still have the following problems: 1) For example, the RLIBM project (“RLIBM: A Novel Approach to Generate Correctly RoundedMath Libraries for New Floating-Point Variants”, in Proceedings of the ACM on Programming Languages ​​(PACMPL), OOPSLA 2021, pp. 1–28) aims to generate correctly rounded mathematical functions for 32-bit floating-point inputs. However, in its implementation, it faces the challenge of solving large-scale infeasible low-dimensional linear programming (LDLP) problems. Furthermore, as the authors state, its original method needs to handle billions of rounding interval constraints, resulting in excessively long solution times and huge memory consumption, making it difficult to apply to complex functions or wider floating-point formats. In addition, the method does not effectively utilize the geometric structure between constraints, leading to a large amount of redundant computation.

[0003] 2) More recent research, such as RLIBM-Pro ("RLIBM-Pro: Efficient Polynomial Generation for Correctly Rounded Math Libraries", IEEE Symposium on Computer Arithmetic (ARITH) 2023), attempts to reduce the number of constraints through convex hull pruning, but is still limited by the serial solution strategy, and has limited scalability for large-scale problems.

[0004] 3) In addition, although the random sampling-based algorithm proposed by Clarkson et al. has performed well in feasible LDLP problems, it is difficult to generate high-quality approximate solutions in infeasible problems.

[0005] Therefore, there is an urgent need to design an efficient solution method that is applicable to infeasible LDLP. Summary of the Invention

[0006] The purpose of this invention is to solve the problem of solving large-scale infeasible low-dimensional linear programming (LDLP) in existing methods, and to provide a parallel linear programming solver method based on RLIBM. By designing a processing chain of "convex hull intermediate representation constraint partitioning - parallel weighted random sampling - maximum consensus solution - parallel verification - iterative optimization", this invention solves the problem of solving large-scale infeasible low-dimensional linear programming (LDLP) in the RLIBM algorithm, and thus solves high-quality solutions to infeasible LDLP.

[0007] To achieve the above objectives, the present invention adopts the following technical solution: a parallel linear programming solver method based on RLIBM, comprising the following steps: S1. Constraint Preprocessing: In the RLIBM project, the constraints based on the nonlinear function behavior of the original form are initialized, and the billions of rounding interval constraints generated by RLIBM are transformed into the standard linear programming (LP) form. The rounding interval constraints are dynamically generated by the range of input floating-point values ​​and function characteristics. S2, Convex Hull Projection Segmentation: Data segmentation is achieved by projecting high-dimensional constraints onto a two-dimensional plane using geometric duality. The segmentation method is as follows: S21) Mapping linear constraints to points in a high-dimensional space through geometric duality: for The RLIBM constraints are transformed into points in the dual space: 1) Upper bound constraints: Corresponding hyperplane 2) Lower bound constraint: Corresponding hyperplane ; In the formula, For the variable to be optimized, For constraint variables, and These represent the upper and lower bounds of the constraint relaxation, respectively. It is a single parameter. Let be the unknown quantity to be optimized and solved. T Represents the transpose operation of a vector; S22) The parameterized single parameter has two-dimensional properties Projecting onto a two-dimensional space generates a two-dimensional convex hull. Points on the convex hull represent the "extreme" part of the constraint set. Points falling on the boundary of the convex hull are identified as potential infeasible constraints, forming an infeasible constraint superset V. Points located inside the convex hull are classified as feasible constraint set X. S3, Parallel weighted random sampling: The feasible constraint set X is sampled by parallel weight adjustment. The weight adjustment method is as follows: for constraints that are easily violated in previous solutions during the solution process, the output sampling set S is gradually optimized by increasing the weight of the constraint so that it better represents the key features of the entire constraint set. S4, Maximum Consensus Solution: The process of finding the maximum satisfying subset in the infeasible constraint superset V is as follows: S41) Original constraints: For each constraint j∈S in the sampling set S, the constraint form is: In the formula, The variable to be solved is... a j ∈R k It is the coefficient vector of the constraints. b j ∈R is a constant on the right-hand side of the constraint; For each constraint in the infeasible superset V i ∈V, the constraint form is: In the formula, a i ∈R k It is the coefficient vector of the constraints. is the right-hand constant of the constraint; S42) Introduce slack variables and indicator variables: For each constraint in the infeasible superset V i ∈V, introduce slack variables s i and indicator variables u i Among them, slack variables s i ≥0 indicates a constraint i The degree of violation, if s i =0, then the constraint is satisfied, if s i If the value is greater than 0, the constraint is violated; indicator variable u i ∈{0,1}: when u i =0 indicates a constraint. i When satisfied, u i =1 indicates a constraint. i Violated; The constraints in the infeasible superset V are relaxed based on the introduced variables as follows: Slack variables within this relaxed condition s i Partial violations are permitted, the extent of which is determined by... si The value determines; S43) Optimization objective: The optimization objective is to minimize the number of violated constraints. And the degree of violation is controlled by introducing a penalty. The optimized formula is: In the formula, >0 is a penalty parameter used to balance the number of violations. and degree of violation The trade-off between them; S44) Continuous processing: Indicator variables that are binary variables u i Relaxation as a continuous variable u i ∈[0,1], and gradually approximate the optimal solution that satisfies all constraints in the sampling set S through an iterative method, and finally generate candidate solutions θ*; S5. Parallel verification of the effectiveness of candidate solutions: The evaluation task of a large-scale constraint set A is decomposed into multiple independent subtasks, and efficient aggregation is performed through parallel computation. The global validity of candidate solutions θ* is verified by multiple threads. Specific methods include: S51) Task decomposition: Treat the large-scale constraint set A as a set of points in a high-dimensional constraint space. The large-scale constraint set A contains a feasible constraint set X and an infeasible superset V. Each point represents a constraint condition. Divide the point set into T subsets, where T represents the number of threads. S52) Dynamic Cooperative Parallel Evaluation: Parallel verification is performed across T threads, where each subset is evaluated by a single computation thread. The verification process is equivalent to determining whether a candidate solution θ* falls within the half-space defined by all constraints in that subset. In parallel verification, each thread independently computes the constraint violations in its subset and generates local data, including the number of local constraint violations. and local violation of constraint weights ; S53) Optimize closed-loop feedback drive: Each constraint is checked and it is determined whether it has been violated. If the constraint has not been violated, the local results are then integrated and output as a global metric through a global reduction operation, including the number of violated constraints. and violation of constraint weights The algorithm calculates the final solution θ^* and returns iterative optimization if the constraints are violated. S6. Iterative optimization: In each iteration, a candidate solution θ* is first generated, and its performance on the entire large-scale constraint set A is comprehensively evaluated. If the candidate solution θ* fails to satisfy all constraints, the weights of the optimization constraints are dynamically adjusted based on the verification results. ω i Simultaneously, increase the iteration count by one, and feed back to steps S3-S5 to repeat the above operation until a final solution θ^* satisfying the constraints is found or the maximum iteration count T is reached. max And when the maximum number of iterations T is reached max When the time comes, output the candidate solution θ* of the last iteration as the final solution θ^*.

[0008] In step S1, the initialization steps are as follows: S11) Adaptively discretize the input interval, dynamically divide it into sub-intervals according to the function curvature, ensure that the function changes smoothly in each sub-interval, and avoid constraint explosion caused by overly fine division; S12) Use a linear approximation strategy for each sub-interval to simplify the nonlinear rounding conditions into linear expressions, thereby eliminating higher-order nonlinear terms; S13) Map the linearized constraints to the standard linear programming (LP) form to ensure that all constraints meet the input requirements of the LP solver.

[0009] In step S2, the feasibility of the remaining constraints is verified by iteratively removing the convex hull vertex constraints. The verification process is as follows: 1) Calculate the convex hull vertex set V1; 2) Partition the constraint set A = Х1∪V1, where Х1 = A1-V1; 3) Use the Clarkson method to verify the feasibility of Х1. If it is feasible, continue to optimize V1; otherwise, recursively calculate the new convex hull until Х is feasible.

[0010] In step S3, the parallel weighted random sampling process is as follows: S31) Weight Initialization and Random Number Generation: Weighted Random Sampling (WRS) is used to initialize and generate random numbers based on each constraint. i The importance of dynamically selecting sampling points, each constraint i ∈Х is assigned a weight ω i The larger the weight, the higher the probability that the constraint will be selected. The initial weight is set to... ω i =1, if a constraint is found to be violated during the iteration, its weight is updated and increased; simultaneously, each constraint is processed in parallel. i Generate random numbers u i ~Uniform[0,1]; S32) Parallel computation of weighted values: based on each constraint i weight ω i With random numbers u i Calculate the weighted value The calculation formula is: ; S33) Local top-sszie selection: Based on the constraint set partitioning method in step S2, the feasible constraint set X is divided into T subsets, where T represents the number of threads. After partitioning, a local sorting is performed, and each thread independently processes one subset, calculating all weighted values ​​within its subset. And select the top ssize maximum values ​​in the local area; S34) Global top-ssize merging: First, perform parallel merge sort to merge the local top-ssize of each thread, obtaining the top ssize of the global maximum weighted values. ; S35) Update the sampling set: Record the indices and weights of the selected ssize constraints. The final output sample set S is used for subsequent LP solution.

[0011] In step S3, the specific implementation algorithm for parallel weighted random sampling is as follows: Function: To obtain a representative sample set S from the feasible set; Input: Feasible constraint set Х, sample size ssize, number of threads T; Output: Sample set S; 1. S = empty set / / Initialize the sampling set 2. / / Assign initial weights to constraints 3. for each constraint i in X: 4. w[i]=1 / / Set the initial weight to 1 5. X_subsets = split_into_T_subsets(XT) / / Split X into T subsets for parallel processing 6. / / Parallel sampling between threads 7. parallel for each thread t from 1 to T: 8. local_S_t = empty list / / Local sample set of thread t 9. for each constraint i in X_subsets[t]: / / Process each constraint in the subset 10. u_i = random_uniform(0,1) / / Generate a random number between 0 and 1 11. v_i = u_i ^(1 / w[i]) / / Calculate the weighted value 12. local_S_t = select_top_ssize_by_v(X_subsets[t]) / / Select the first ssize constraints 13. / / Merge the results of all threads 14. S_temp = merge_all_local_S(local_S_t for all t) 15. / / Select the first ssize constraints globally 16. S=select_global_top_ssize(S_temp) 17. return S / / Returns the final sample set.

[0012] In step S4, the specific implementation algorithm for Maximum Consensus (MCS) is as follows: Function: Generate candidate solutions; Input: Sample set S, infeasible constraint superset V; Output: Candidate solution θ*; 1. for each constraint i in V: / / Introduce slack variables for the constraints in V 2. s_i>= 0 / / Non-negative slack variable 3. minimize sum(u_i) + alpha*sum(s_i) / / Optimize the construction problem 4. subject to: 5. a_i^T*θ-b_i<=s_i for all i in V / / For i∈V 6. a_j^T*θ-b_j<=0 for all j in S / / for j∈S 7. s_i>= 0,u_i in [0,1] / / Variable limits 8. / / Initialize binary indicator 9. for each i in V: 10. u_i = 0 / / Initial value of the indicator variable 11. / / Iterative optimization 12. while not converged: / / Continue until convergence 13. solve_LP_for_theta_and_s() / / Fix u_i, solve for s_i and θ 14. for each i in V: / / Update u_i based on the current solution 15. if condition_met: / / Define the condition 16. u_i = 1 17. else: 18. u_i = 0 19. if change < threshold: / / Check if the change in the objective function is small enough 20. break 21.return θ* / / Return the candidate solution.

[0013] In the step S5, the specific parallel verification algorithm is as follows: Function: Quickly verify the feasibility of the candidate solution; Input: Candidate solution θ*, large-scale constraint set A, number of threads T; Output: Number of violated constraints n_v, sum of violated constraint weights ω_v; 1. A_subsets = split_into_T_subsets(A,T) / / Divide A into T subsets 2. n_v = 0 / / Initialize the number of violated constraints 3. ω_v = 0 / / Initialize the sum of violated constraint weights 4. parallel for each thread t from 1 to T: / / Parallel verification among threads 5. local_n_v = 0 / / Local number of violated constraints 6. local_ω_v = 0 / / Local sum of violated constraint weights 7. for each constraint i in A_subsets[t]: / / Check each constraint condition 8. if a_i^T*θ* >b_i: / / If the constraint is violated 9. local_n_v += 1 10. if i in X: / / If the constraint is in the feasible set 11. local_ω_v += ω_i / / Calculate the local weight 12. atomic_add(n_v,local_n_v) / / Add the local to the global 13. atomic_add(ω_v,local_ω_v) / / Add the local to the global 14. / / Adjust the weights 15. if ω_v < (1 / (3k-1)) * ω_s (ω_s is the total weight of the sample set, ω_s = ∑ω_i): 16. for each violated constraint i: 17. ω_i *= 2 / / The weight of a constraint violation is doubled. 18. return ω_v, n_v.

[0014] In step S6, the specific algorithm for iterative optimization is as follows: Function: Iteratively find the optimal solution Input: Large-scale constraint set A, feasible constraint set X, infeasible constraint superset V, maximum number of iterations T_max.

[0015] Output: Final solution θ^*.

[0016] 1.t = 0 / / / Initialize the iteration counter 2. / / Iterate until the maximum number of iterations is reached or the final solution is found. 3. while t <T_max: 4. S_t = PWRS(X,ssize,T) / / Use PWRS to obtain the sample set 5. θ_t^* = MCS(S_t,V) / / Use MCS to obtain candidate solutions 6. n_v,ω_v = Parallel_Validation(θ_t^*,A,T) / / Use a parallel validation algorithm to obtain n_v and ω_v 7. if n_v == 0: / / No constraint is violated 8. return θ_t^* / / Returns the final solution 9. adjust_weights_based_on_ω_v() / / Adjust weights based on validation results 10. t += 1 / / Increase the iteration counter 11. return θ_t^* / / If convergence is not achieved, return the solution from the last iteration.

[0017] The beneficial effects of this invention are: 1) The method of this invention solves the problem of solving large-scale infeasible low-dimensional linear programming (LDLP) encountered in the RLIBM algorithm by designing a processing link of "convex hull intermediate representation constraint segmentation - parallel weighted random sampling - maximum consensus solution - parallel verification - iterative optimization", and introduces the principle of geometric duality to project high-dimensional constraints onto a two-dimensional plane and construct convex hull intermediate representation; finally, it solves the high-quality solution of infeasible LDLP.

[0018] 2) The solution method of this invention identifies the key constraints corresponding to the vertices of the convex hull, and divides the billions of constraints into feasible subsets and infeasible supersets through convex hull intermediate representation and constraint segmentation techniques, effectively reducing the problem size and significantly reducing the solution complexity. On this basis, a parallel weighted random sampling strategy is designed to accelerate the extraction of key constraints and the problem-solving process by dynamically adjusting constraint weights and parallel computing techniques, achieving efficient sampling of key constraints. Combined with the maximum consensus solution algorithm, the number of constraints satisfied in the infeasible superset is maximized under the premise of satisfying the feasible subset constraints, ensuring the generation of high-quality approximate solutions in infeasible problems. The high-quality solution is gradually approximated through the predictor-corrector iterative mechanism. Experiments show that the average speedup is 6.0 times under a 32-core configuration. Attached Figure Description

[0019] Figure 1 This is a flowchart illustrating the parallel solution framework and processing flow of the method of the present invention. Figure 2 This is a comparison chart of the algorithm's execution time in the method of this invention; Figure 3 This is a speedup chart for the algorithm in the method of this invention. Detailed Implementation

[0020] The present invention will be further explained and described below with reference to the accompanying drawings and specific embodiments.

[0021] Example: Figure 1 As shown, this invention provides a parallel linear programming solver method based on RLIBM, characterized by the following steps: S1. Constraint Preprocessing: In the RLIBM project, rounding interval constraints are dynamically generated based on the range of input floating-point values ​​and the function characteristics, reaching billions in number (for example, the e^x function requires handling approximately 500 million constraints within a specific input range). The original form of these constraints is based on nonlinear function behavior, and directly solving them leads to computational inefficiency and difficulty in scaling. Therefore, the constraints, based on nonlinear function behavior in their original form, are initialized, transforming the billions of rounding interval constraints generated by RLIBM into a standard linear programming (LP) form, where the rounding interval constraints are dynamically generated based on the range of input floating-point values ​​and the function characteristics.

[0022] The initialization steps are as follows: S11) Adaptively discretize the input interval, dynamically dividing it into sub-intervals based on the function curvature (such as the exponential growth characteristic of e^x) to ensure that the function changes smoothly within each sub-interval and avoid constraint explosion caused by overly fine division; S12) Use a linear approximation strategy (such as local linearization based on the function derivative) for each sub-interval to simplify the nonlinear rounding conditions into linear expressions, thereby eliminating higher-order nonlinear terms; S13) Map the linearized constraints to the standard linear programming (LP) form (i.e., a linear combination of constraint matrices and vectors) to ensure that all constraints meet the input requirements of the LP solver.

[0023] This preprocessing step significantly improves the structuring of constraints, transforming billions of nonlinear constraints that were previously impossible to process directly into a set of LP constraints that can be efficiently processed in parallel. For example, in solving the exp10 function, the preprocessing compresses 504 million original constraints into an equal number of LP constraints while strictly maintaining rounding accuracy requirements. Experiments show that this step consumes only 5-10% of the total computational cost, yet brings several times the performance improvement to subsequent convex hull partitioning and parallel solving, avoiding solution failures or accuracy losses caused by constraint complexity. Through this transformation, RLIBM can complete large-scale rounding interval verification within a reasonable time, ensuring that the generated elementary functions achieve truly correct rounding.

[0024] S2, Convex Hull Projection Segmentation: Data segmentation is achieved by projecting high-dimensional constraints onto a two-dimensional plane using geometric duality. The segmentation method is as follows: S21) Mapping linear constraints to points in a high-dimensional space through geometric duality: for The RLIBM constraints are transformed into points in the dual space: 1) Upper bound constraints: Corresponding hyperplane 2) Lower bound constraint: Corresponding hyperplane ; In the formula, For the variable to be optimized, For constraint variables, and These represent the upper and lower bounds of the constraint relaxation, respectively. It is a single parameter. Let be the unknown quantity to be optimized and solved. T Represents the transpose operation of a vector; S22) The parameterized single parameter has two-dimensional properties Projected into a two-dimensional space (e.g., selecting) or A two-dimensional convex hull is generated. Points on the convex hull represent the "extreme" part of the constraint set. Points falling on the boundary of the convex hull are identified as potential infeasible constraints, forming an infeasible constraint superset V. Points located inside the convex hull are classified as feasible constraint set X. By using this method of constraint partitioning through convex hull intermediate representation, the size of the infeasible constraint superset V is typically reduced dramatically. For example, it can be reduced from billions of constraints to thousands of constraints. This order-of-magnitude compression makes subsequent analysis and solution feasible.

[0025] The feasibility of the remaining constraints is verified by iteratively removing the convex hull vertex constraints. The verification process is as follows: 1) Calculate the convex hull vertex set V1; 2) Partition the constraint set A = Х1∪V1, where Х1 = A1-V1; 3) Use the Clarkson method to verify the feasibility of Х1. If it is feasible, continue to optimize V1; otherwise, recursively calculate the new convex hull until Х is feasible.

[0026] S3, Parallel weighted random sampling: The feasible constraint set X is sampled by parallel weight adjustment. The weight adjustment method is as follows: for constraints that are easily violated in previous solutions during the solution process, the output sampling set S is gradually optimized by increasing the weight of the constraint so that it better represents the key features of the entire constraint set. The process of parallel weighted random sampling is as follows: S31) Weight Initialization and Random Number Generation: Weighted Random Sampling (WRS) is used to initialize and generate random numbers based on each constraint. i The importance of dynamically selecting sampling points, each constraint i ∈Х is assigned a weight ω i The larger the weight, the higher the probability that the constraint will be selected. The initial weight is set to... ω i =1, if a constraint is found to be violated during the iteration, its weight is updated and increased; simultaneously, each constraint is processed in parallel. i Generate random numbers u i ~Uniform[0,1]; S32) Parallel computation of weighted values: based on each constraint i weight ω i With random numbers u i Calculate the weighted value The calculation formula is: ; S33) Local top-sszie selection: Based on the constraint set partitioning method in step S2, the feasible constraint set X is divided into T subsets, where T represents the number of threads. After partitioning, a local sorting is performed, and each thread independently processes one subset, calculating all weighted values ​​within its subset. And select the top ssize maximum values ​​in the local area; S34) Global top-ssize merging: First, perform parallel merge sort to merge the local top-ssize of each thread, obtaining the top ssize of the global maximum weighted values. ; S35) Update the sampling set: Record the indices and weights of the selected ssize constraints. The final output sample set S is used for subsequent LP solution.

[0027] S4, Maximum Consensus Solution: The process of finding the maximum satisfying subset in the infeasible constraint superset V is as follows: S41) Original constraints: For each constraint j∈S in the sampling set S, the constraint form is: In the formula, The variable to be solved is... a j ∈R k It is the coefficient vector of the constraints. b j ∈R is a constant on the right-hand side of the constraint; For each constraint in the infeasible superset V i ∈V, the constraint form is: In the formula, a i ∈R k It is the coefficient vector of the constraints. is the right-hand constant of the constraint; S42) Introduce slack variables and indicator variables: For each constraint in the infeasible superset V i ∈V, introduce slack variables s i and indicator variables u i Among them, slack variables s i ≥0 indicates a constraint i The degree of violation, if s i =0, then the constraint is satisfied, if s i If the value is greater than 0, the constraint is violated; indicator variable ui ∈{0,1}: when u i =0 indicates a constraint. i When satisfied, u i =1 indicates a constraint. i Violated; The constraints in the infeasible superset V are relaxed based on the introduced variables as follows: Slack variables within this relaxed condition s i Partial violations are permitted, the extent of which is determined by... s i The value determines; S43) Optimization objective: The optimization objective is to minimize the number of violated constraints. And the degree of violation is controlled by introducing a penalty. The optimized formula is: In the formula, >0 is a penalty parameter used to balance the number of violations. and degree of violation The trade-off between them; S44) Continuous processing: Indicator variables that are binary variables u i Relaxation as a continuous variable u i The solution ∈[0,1] is iterated and gradually approximated to the optimal solution that satisfies all constraints in the sample set S, and finally the candidate solution θ* is generated. This continuous approximation can significantly reduce the computational complexity in practical applications.

[0028] S5. Parallel verification of the effectiveness of candidate solutions: The evaluation task of a large-scale constraint set A is decomposed into multiple independent subtasks and efficiently aggregated through parallel computation. The global validity of candidate solutions θ* is verified by multiple threads. This decomposition not only reduces the computational burden of a single thread, but also preserves the overall geometric properties of the constraint set, ensuring the accuracy of the verification results.

[0029] Specific methods include: S51) Task decomposition: Treat the large-scale constraint set A as a set of points in a high-dimensional constraint space. The large-scale constraint set A contains a feasible constraint set X and an infeasible superset V. Each point represents a constraint condition. Divide the point set into T subsets, where T represents the number of threads. S52) Dynamic Cooperative Parallel Evaluation: Parallel verification is performed across T threads, where each subset is evaluated by a single computation thread. The verification process is equivalent to determining whether a candidate solution θ* falls within the half-space defined by all constraints in that subset. In parallel verification, each thread independently computes the constraint violations in its subset and generates local data, including the number of local constraint violations. and local violation of constraint weights ; S53) Optimize closed-loop feedback drive: Each constraint is checked and it is determined whether it has been violated. If the constraint has not been violated, the local results are then integrated and output as a global metric through a global reduction operation, including the number of violated constraints. and violation of constraint weights The algorithm calculates the final solution θ^* and returns iterative optimization if the constraints are violated.

[0030] Distributed collaborative verification mechanism: By deconstructing the data dependencies in the verification process, the originally tightly coupled verification tasks are reconstructed into independent units that can be processed in parallel. This paradigm shift in computing based on the characteristics of multi-core architecture not only achieves an order-of-magnitude improvement in verification efficiency but also theoretically expands the verification boundaries of large-scale optimization problems. Practice shows that this method, through dynamic load balancing strategies and fine-grained thread scheduling, reduces verification time from hours to minutes while ensuring computational accuracy, making it a design in this scheme that balances engineering value and academic innovation.

[0031] S6. Iterative optimization: In each iteration, a candidate solution θ* is first generated, and its performance on the entire large-scale constraint set A is comprehensively evaluated. If the candidate solution θ* fails to satisfy all constraints, the weights of the optimization constraints are dynamically adjusted based on the verification results. ω i Simultaneously, increase the iteration count by one, and feed back to steps S3-S5 to repeat the above operation until a final solution θ^* satisfying the constraints is found or the maximum iteration count T is reached. max And when the maximum number of iterations T is reached max When the time comes, output the candidate solution θ* of the last iteration as the final solution θ^*.

[0032] These mechanisms work together to ensure that the iterative process gradually converges to a high-quality solution, while significantly improving computational efficiency through efficient algorithm design and parallel computing. The theoretical basis of iterative optimization combines the advantages of weighted random sampling and maximum consensus solution. For infeasible problems, although the theoretical guarantee of the global optimal solution is relatively complex, iterative optimization can gradually approach a high-quality approximate solution through the combination of dynamic weight adjustment and MCS.

[0033] Maximum number of iterations Tmax The algorithm's design is based on a trade-off between experimental observations and computational resources. We propose a heuristic rule: the algorithm is considered convergent when the number of constraint violations does not significantly decrease (rate of change < 5%) in k consecutive iterations (k = 20 in this paper). In the experiments, T... max =500, combined with the early stop mechanism, can cover the convergence requirements of 97.5% of test cases.

[0034] Experimental Results and Analysis: 1. Experimental conditions and preparation To evaluate this invention, this paper presents six basic functions for correct rounding (e x ,2 x 10 x ln(x), log2(x), and log 10 (x) creates a new version of the solver that takes a linear programming problem as input and generates a double-precision floating-point solution for all unknown variables.

[0035] The experimental platform for this invention is the Shenwei high-performance server, and the CPU is the domestically produced 28nm Shenwei 1621 processor.

[0036] 2. Analysis of Experimental Results Table 1 compares the capabilities of four solvers: the solver in this paper, the solver in the RLIBM project, the infeasible LDLPs solver, and the MCS solver from the computer vision community, for solving linear programming problems that generate correct rounding implementations.

[0037] Table 1 Table 1 presents the correctness verification results of four solvers under six basic function constraint scenarios. Experimental data shows that the proposed method exhibits the same advantages in infeasible low-dimensional linear programming problems in RLIBM. For test cases with constraints ranging from 7.16 million to 504 million generated by functions such as log and exp, the proposed solver achieved 100% accuracy in infeasible LP scenarios, while the native RLIBM solver failed in all scenarios except for feasible cases with the log2 function. The main reason for its failure is the lack of a mixed-precision constraint handling mechanism. For example, in the 504 million constraint case with the exp10 function, RLIBM still could not meet the precision threshold requirement after 500 iterations. Although the infeasible LDLPs solver can handle infeasible problems, it is limited by its serial architecture and cannot break through the 500 million constraint scale. The MCS solver in the field of computer vision, lacking a floating-point precision compensation strategy, failed to generate valid solutions in all test cases. Its core consensus mechanism misjudged in the non-convex constraint space of the log function.

[0038] Table 2 shows a comparison of solver performance, using the time of infeasible LDLPs solvers as the comparison value. There are 4-core, 8-core, 16-core, and 32-core solvers available.

[0039] Table 2 Experimental results show that the parallel linear programming solver proposed in this invention exhibits good performance advantages on the Shenwei platform. As shown in Table 2, for six basic functions such as exp and log, the solver achieves an average speedup of 6.0 times under a 32-core configuration. Among them, the execution time of the exp10 function under the constraint of 28.61 million is reduced from 1929 seconds to 301 seconds, achieving a speedup of 6.4 times.

[0040] Figure 2 and Figure 3 The comparison curves reveal a strong correlation between parallel efficiency and problem size: the exp series functions, due to their high computational density and strong data locality, exhibit a near-linear speedup of 6.3 times with 16 cores, while the log series functions, limited by branch prediction and memory access latency, show a stable speedup of 4.2-5.1 times with 32 cores. Test data also validates the scalability of the method; as the number of cores increases from 4 to 32, the overall speedup curve conforms to Amdahl's Law prediction, and the communication overhead gradually increases from 18.7% to 26.4%. Further experiments show that, compared to infeasible LDLPs solvers, the proposed method, through constraint partitioning and parallel verification mechanisms, improves the iterative convergence speed of the log2 function by 3.9 times, verifying the effectiveness of the parallel weighted random sampling strategy in conflict constraint identification.

[0041] This invention addresses the problem of solving large-scale infeasible low-dimensional linear programming (LDLP) problems by constructing a processing chain of "convex hull intermediate representation constraint partitioning - parallel weighted random sampling - maximum consensus solution - parallel verification - iterative optimization".

[0042] The above description is only used to illustrate the technical solution of the present invention and is not intended to limit it. Any other modifications or equivalent substitutions made by those skilled in the art to the technical solution of the present invention, as long as they do not depart from the spirit and scope of the technical solution of the present invention, should be covered within the scope of the claims of the present invention.

Claims

1. A parallel linear programming solver method based on RLIBM, characterized in that: Includes the following steps: S1. Constraint Preprocessing: In the RLIBM project, the constraints based on the nonlinear function behavior of the original form are initialized, and the billions of rounding interval constraints generated by RLIBM are transformed into the standard linear programming (LP) form. The rounding interval constraints are dynamically generated by the range of input floating-point values ​​and function characteristics. S2, Convex Hull Projection Segmentation: The data is segmented by projecting high-dimensional constraints onto a two-dimensional plane using geometric duality. The segmentation method is as follows: S21) Mapping linear constraints to points in a high-dimensional space through geometric duality: for The RLIBM constraints are transformed into points in the dual space: 1) Upper bound constraints: Corresponding hyperplane 2) Lower bound constraint: Corresponding hyperplane ; In the formula, For the variable to be optimized, For constraint variables, and These represent the upper and lower bounds of the constraint relaxation, respectively. It is a single parameter. Let be the unknown quantity to be optimized and solved. T Represents the transpose operation of a vector; S22) The parameterized single parameter has two-dimensional properties Projecting onto a two-dimensional space generates a two-dimensional convex hull. Points on the convex hull represent the "extreme" part of the constraint set. Points falling on the boundary of the convex hull are identified as potential infeasible constraints, forming an infeasible constraint superset V. Points located inside the convex hull are classified as feasible constraint set X. S3, Parallel weighted random sampling: The feasible constraint set X is sampled by parallel weight adjustment. The weight adjustment method is as follows: for constraints that are easily violated in previous solutions during the solution process, the output sampling set S is gradually optimized by increasing the weight of the constraint so that it better represents the key features of the entire constraint set. S4, Maximum Consensus Solution: The process of finding the maximum satisfying subset in the infeasible constraint superset V is as follows: S41) Original constraints: For each constraint j∈S in the sampling set S, the constraint form is: In the formula, The variable to be solved is... a j ∈R k It is the coefficient vector of the constraints. b j ∈R is a constant on the right-hand side of the constraint; For each constraint in the infeasible superset V i ∈V, the constraint form is: In the formula, a i ∈R k It is the coefficient vector of the constraints. is the right-hand constant of the constraint; S42) Introduce slack variables and indicator variables: For each constraint in the infeasible superset V i ∈V, introduce slack variables s i and indicator variables u i Among them, slack variables s i ≥0 indicates a constraint i The degree of violation, if s i =0, then the constraint is satisfied, if s i If the value is greater than 0, the constraint is violated; indicator variable u i ∈{0,1}: when u i =0 indicates a constraint. i When satisfied, u i =1 indicates a constraint. i Violated; The constraints in the infeasible superset V are relaxed based on the introduced variables as follows: Slack variables within this relaxed condition s i Partial violations are permitted, the extent of which is determined by... s i The value determines; S43) Optimization objective: The optimization objective is to minimize the number of violated constraints. And the degree of violation is controlled by introducing a penalty. The optimized formula is: In the formula, >0 is a penalty parameter used to balance the number of violations. and degree of violation The trade-off between them; S44) Continuous processing: Indicator variables that are binary variables u i Relaxation as a continuous variable u i ∈[0,1], and gradually approximate the optimal solution that satisfies all constraints in the sampling set S through an iterative method, and finally generate candidate solutions θ*; S5. Parallel verification of the effectiveness of candidate solutions: The evaluation task of a large-scale constraint set A is decomposed into multiple independent subtasks, and efficient aggregation is performed through parallel computation. The global validity of candidate solutions θ* is verified by multiple threads. Specific methods include: S51) Task decomposition: Treat the large-scale constraint set A as a set of points in a high-dimensional constraint space. The large-scale constraint set A contains a feasible constraint set X and an infeasible superset V. Each point represents a constraint condition. Divide the point set into T subsets, where T represents the number of threads. S52) Dynamic Cooperative Parallel Evaluation: Parallel verification is performed across T threads, where each subset is evaluated by a single computation thread. The verification process is equivalent to determining whether a candidate solution θ* falls within the half-space defined by all constraints in that subset. In parallel verification, each thread independently computes the constraint violations in its subset and generates local data, including the number of local constraint violations. and local violation of constraint weights ; S53) Optimize closed-loop feedback drive: Each constraint is checked and it is determined whether it has been violated. If the constraint has not been violated, the local results are then integrated and output as a global metric through a global reduction operation, including the number of violated constraints. and violation of constraint weights The algorithm calculates the final solution θ^* and returns iterative optimization if the constraints are violated. S6. Iterative optimization: In each iteration, a candidate solution θ* is first generated, and its performance on the entire large-scale constraint set A is comprehensively evaluated. If the candidate solution θ* fails to satisfy all constraints, the weights of the optimization constraints are dynamically adjusted based on the verification results. ω i Simultaneously, increase the iteration count by one, and feed back to steps S3-S5 to repeat the above operation until a final solution θ^* satisfying the constraints is found or the maximum iteration count T is reached. max And when the maximum number of iterations T is reached max When the time comes, output the candidate solution θ* of the last iteration as the final solution θ^*.

2. The parallel linear programming solver method based on RLIBM according to claim 1, characterized in that: In step S1, the initialization steps are as follows: S11) Adaptively discretize the input interval, dynamically divide it into sub-intervals according to the function curvature, ensure that the function changes smoothly in each sub-interval, and avoid constraint explosion caused by overly fine division; S12) Use a linear approximation strategy for each sub-interval to simplify the nonlinear rounding conditions into linear expressions, thereby eliminating higher-order nonlinear terms; S13) Map the linearized constraints to the standard linear programming (LP) form to ensure that all constraints meet the input requirements of the LP solver.

3. The parallel linear programming solver method based on RLIBM according to claim 1, characterized in that: In step S2, the feasibility of the remaining constraints is verified by iteratively removing the convex hull vertex constraints. The verification process is as follows: 1) Calculate the convex hull vertex set V1; 2) Partition the constraint set A = Х1∪V1, where Х1 = A1-V1; 3) Use the Clarkson method to verify the feasibility of Х1. If it is feasible, continue to optimize V1; otherwise, recursively calculate the new convex hull until Х is feasible.

4. The parallel linear programming solver method based on RLIBM according to claim 1, characterized in that: In step S3, the parallel weighted random sampling process is as follows: S31) Weight Initialization and Random Number Generation: Weighted Random Sampling (WRS) is used to initialize and generate random numbers based on each constraint. i The importance of dynamically selecting sampling points, each constraint i ∈Х is assigned a weight ω i The larger the weight, the higher the probability that the constraint will be selected. The initial weight is set to... ω i =1, if a constraint is found to be violated during the iteration, update and increase its weight; Simultaneously, each constraint is processed in parallel. i Generate random numbers u i ~Uniform[0,1]; S32) Parallel computation of weighted values: based on each constraint i weight ω i With random numbers u i Calculate the weighted value The calculation formula is: ; S33) Local top-sszie selection: Based on the constraint set partitioning method in step S2, the feasible constraint set X is divided into T subsets, where T represents the number of threads. After partitioning, a local sorting is performed, and each thread independently processes one subset, calculating all weighted values ​​within its subset. And select the top ssize maximum values ​​in the local area; S34) Global top-ssize merging: First, perform parallel merge sort to merge the local top-ssize of each thread, obtaining the top ssize of the global maximum weighted values. ; S35) Update the sampling set: Record the indices and weights of the selected ssize constraints. The final output sample set S is used for subsequent LP solution.

Images