Robust Verification of Numerical Software
Abstract
Numerical software are widely used in safetycritical systems such as aircrafts, satellites, car engines and so on, facilitating dynamics control of such systems in real time, it is therefore absolutely necessary to verify their correctness. It is a long standing challenge to guarantee verified properties of numerical software are indeed satisfied by their real behaviours, because most of these verifications are conducted under ideal mathematical models, but their real executions could be influenced essentially by uncertain inputs accounting for roundoff errors and additive perturbations from realworld phenomena to hardware, which are abstracted away in these ideal mathematical models. In this paper, we attempt to address this issue focusing on nontermination analysis of numerical software, where nontermination is often an unexpected behaviour of computer programs and may be problematic for applications such as realtime systems having hard deadlines on transaction execution time, and propose a method for robust conditional nontermination analysis, which can be used to underapproximate the maximal robust nontermination input set for a given program, from which the program never terminates, regardless of the aforementioned disturbances. Finally, several examples are employed to illustrate our approach.
I Introduction
Software is ubiquitous in missioncritical and safetycritical industrial infrastructures since it is, in principle, the most effective way to manipulate complex systems in real time. However, many computer scientist and engineers have experienced costly bugs in embedded software. The failure of the Ariane 5.01 maiden flight (due to an overflow caused by an unprotected data conversion from a too large 64bit floating point to a 16bit signed integer value), the failure of the Patriot missile during the Gulf war (due to an accumulated rounding error), the loss of Mars orbiter(due to a unit error) are a few examples showing that missioncritical and safetycritical software can be far from being safe [11]. It is therefore absolutely necessary to prove the correctness of software by using formal, mathematical techniques, aiding the development of correct and reliable software systems.
The dominant approach to the verification of programs is so called FloydHoareNaur inductive assertion approach [13, 18, 34], basically consisting of pre and postcondition to specify the condition of initial states, and the property that should be satisfied by terminated states, respectively, and Hoare logic to reason about properties of programs. The hardest parts of this approach are invariant generation and termination analysis. It is wellknown that the termination or nontermination problem is undecidable, and even not semidecidable in general. Thus, more practical approaches to termination analysis are either to present some sufficient conditions for termination, or some sufficient conditions for nontermination, or put these two types of conditions in parallel, or prove the decidability for some specific families of programs, e.g., [16, 53, 14, 7, 39, 23, 25, 3].
On the other hand, it is a long standing challenge to guarantee verified properties of numerical programs, which are indeed satisfied by their real behaviours, because most of these verifications are conducted under some ideal mathematical model, but their real executions could be influenced essentially by uncertain inputs accounting for roundoff errors and additive perturbations from realworld phenomena to hardware, which are abstracted away in the ideal mathematical model. Such general remarks on formal verification techniques is certainly applicable to most of existing termination/nontermination analysis approaches. In [54], the authors presented the following example:
Example 1
Consider a simple loop
where , with , and if taking 10 decimal digits of precision.
According to the termination decidability result on simple loops proved in [48], Q1 terminates based on exact computation. But unfortunately, it becomes unterminated in practice as the core decision procedure given in [48] invokes a procedure to compute Jordan normal form based on numeric computation for a given matrix, and thus the representation error has to be taken into account. In order to address this issue, a symbolic decision procedure was proposed in [54]. However, a more interesting and challenging issue is to investigate a systematic way to take all possible disturbances into account during conducting termination and nontermination analysis in practical numerical implementations.
In this paper we attempt to address this challenge, and propose a framework for robust nontermination analysis for numerical software based on control theory as in [40], which proposes a controltheoretic framework based on Lyapunov invariants to conduct verification of numerical software. Nontermination analysis proves that programs, or parts of a program, do not terminate. Nontermination is often an unexpected behaviour of computer programs and exposes a bug in their code, e.g., if a nonterminating computation occurs, this may be problematic for applications such as realtime systems, where there are hard deadlines on transaction execution time and where minimizing work is important. The computer program of interest in this paper is restricted to a class of computer programs composed of a single loop with a complicated switchcase type loop body which are equivalent to constrained piecewise discretetime dynamical systems subjected to timevarying uncertainties. To address the problem of robust conditional nontermination, we by the aid of the constrained piecewise discretetime dynamical system to characterize the maximal robust nontermination input set by means of a value function, which is a solution to a suitable infinite horizon stateconstrained optimal control problem derived via dynamic programming principle. A new finding in disproving termination of computer programs is that if there does not exist a solution to (10), the robust nontermination set is empty. In addition, when the dynamics of the piecewise discretetime system in each mode is polynomial, and the state and uncertain input constraints are semialgebraic, the optimal control problem is relaxed as a semidefinite programming problem, to which its polynomial solution forms an innerapproximation of the maximal robust nontermination input set if it exists. Such relaxation is sound but incomplete. Finally, several examples are employed to illustrate our approach.
Clearly, the concept of robust nontermination input sets is essentially equivalent to the one of maximal robustly positive invariants in control theory, please refer to, e.g., [2, 47, 40, 45]. How to compute the maximal robustly positive invariant of a given dynamical system is still a longstanding and challenging problem not only in the community of computer science but also in control theory. Most of existing works on this subject focus on linear systems, e.g. [38, 21, 47, 45, 51]. Although some methods have been proposed to synthesize positively invariants for nonlinear systems, e.g., the barrier certificate generation method as in [36, 37] and the region of attraction generation method as in [19, 49, 30, 15] could be able to synthesize maximal robustly positive invariants. This, however, leads to bilinear sumofsquares programs, which are notoriously hard to solve. In order to solve the bilinear sumofsquares programs, a commonly used method is to employ some form of alteration (e.g.,[19, 50, 30]) with a feasible initial solution to the bilinear sumofsquares program. Recently,[43, 44] proposed linear programming based methods to synthesize maximal (robustly) positive polyhedral invariants. Contrasting with aforementioned methods, in this paper we propose a linear semidefinite programming based method to compute semialgebraic invariant and our method does not require an initial feasible solution.
Organization of the paper. The structure of this paper is as follows. In Section II, basic notoins used throughout this paper and the problem of interest are introduced. Then we elucidate our approach for performing conditional nontermination analysis in Section III. After demonstrating our approach on sevral illustrative examples in Sectin IV, we discuss related work in Section V and finally conclude this paper in Section VI.
Ii Preliminaries
In this section we describe the programs which are considered in this paper and we explain how to analyze them through their representation as piecewise discretetime dynamical systems.
The following basic notations will be used throughout the rest of this paper: stands for the set of nonnegative integers and for the set of real numbers; denotes the ring of polynomials in variables given by the argument, denotes the vector space of real multivariate polynomials of degree , . Vectors are denoted by boldface letters.
Iia Computer Programs of Interest
In this paper the computer program of interest, as described in Program 1, is composed of a single loop with a possibly complicated switchcase type loop body, in which variables are assigned using parallel assignments , where is the vector of uncertain inputs, of which values are sampled at random from a compact set, i.e. , such as roundoff errors in performing computations. The form of the analyzed program is described in Program 1.
In Program 1, is a compact set in and , is continuous over . stands for the initial condition on inputs; stands for the loop condition, which is a compact set in ; , , stands for the th branch conditions, where . , , , , , are continuous functions over and over respectively. Moreover, forms a complete partition of , i.e. for , where , and .
As described in Program 1, an update of the variable is executed by the th branch if and only if the current value of satisfies the th branch condition .
IiB Piecewise Discretetime Systems
In this subsection we interprete Program 1 as constrained piecewise discretetime dynamical systems with uncertain inputs. Formally,
Definition 1
A constrained piecewise discretetime dynamical system (PS) is a quadruple with

is the condition on initial states;

is the domain constraint, which is a compact set. A path can evolve complying with the discrete dynamics only if its current state is in ;

with as interpreted in Program 1;

is the set of uncertain inputs;

is the family of the continuous functions .
In order to enhance the understanding of PS, we use the following figure, i.e. Fig. 1, to illustrate it further.
From now on, we associate a PS representation to each program of the form Program 1. Since a program may admit several PS representations, we choose one of them, but the choice does not change the results provided in this paper.
Definition 2
An input policy is an ordered sequence , where , and is defined as the set of input policies, i.e. .
If an input policy makes Program 1 nonterminate from an initial state from , the trajectory from following the discrete dynamics is defined by
(1) 
where , , and
with , , representing the indicator function of the set , i.e.
Consequently, Program 1 is said to be robust nontermination starting from an initial state if for any input policy , holds. Formally,
Definition 3
A program of Program 1 is said to be robust nonterminating w.r.t. an initial state , if
(2) 
Now, we define our problem of deciding a set of initial states rendering Program 1 robust nontermination.
Definition 4 (Robust Nontermination Set)
A set of initial states in is a robust nontermination set for a program of the form Program 1 if is robustly nonterminating w.r.t. for any . We call is robustly nonterminating w.r.t. the maximal robust nontermination set, denoted by .
From Definition 4, we observe that is a subset of such that all runs of Program 1 starting from it can not breach it forever, i.e. if , for . Therefore, the set is equivalent to the maximal robust positively invariant for PS (1) in control theory. For the formal concept of maximal robust positively invariant, please refer to, e.g., [2, 47, 45]. In this paper we formulate the problem of robust conditional nontermination for Problem 1 as a constrained optimal control problem in the control theory framework.
Iii Robust NonTermination Set Generation
In this section we elucidate our approach of addressing the problem of robust conditional nontermination for Program 1, i.e. synthesizing robust nontermination sets as presented in Definition 4. For this sake, we firstly in Subsection IIIA characterize the maximal robust nontermination set by means of the value function, which is a solution to a suitable infinite horizon stateconstrained optimal control problem based on PS (1). Any lower semicontinuous solution to this optimal control problem generates a robust nontermination set. Then, in the case that , , is polynomial over and , and the constraint sets over and , i.e. , , and , are of the basic semialgebraic form, the semidefinite program arising from sumofsquares decompositions facilitates the gain of innerapproximations of via solving the relaxation of the derived optimal control problem in Subsection IIIB.
Iiia Characterization of
In this subsection, we firstly introduce the value function to characterize the maximal robust nontermination set and then formulate it as a solution to a constrained optimal control problem.
For , the value function is defined by:
(3) 
Note that may be neither continuous nor semicontinuous. (A function is lower semicontinuous iff for any , is open, e.g., [4].)
The following theorem shows the relation between the value function and the maximal robust nontermination set , that is, the zero sublevel set of is equal to the maximal robust nontermination set .
Theorem 1
, where is the maximal robust nontermination set as in Definition 4.
From Theorem 1, the maximal robust nontermination set could be constructed by computing , which satisfies the dynamic programming principle as presented in Lemma 1 according to standard techniques from optimal control.
Lemma 1
For and , we have:
(5) 
Proof 2
Let
(6) 
We will prove that for , .
According to the definition of , i.e. (3), for any , there exists an input policy such that
We then introduce two infinite uncertain input policies and such that with for and with . Now, let , then we obtain that
(7) 
Therefore,
(8) 
On the other hand, by the definition of , for any , there exists a such that
Also, by the definition of , i.e. (3), for any , there exists a such that
where . We define such that for and for . Then, it follows
(9) 
Based on Lemma 1 stating that the value function complies with the dynamic programming principle (1), we derive a central equation of this paper, to which is a lower semicontinuous solution. The equation is formulated formally in Theorem 2.
Theorem 2
The value function in (3) is a solution to the equation
(10) 
According to Theorem 2, we conclude that if there does not exist a solution to (10), the robust nontermination set is empty. Moreover, according to Theorem 2, as defined in (3) is a solution to (10). Note that the solution to (10) may be not unique, and we do not go deeper into this matter in this paper. However, any solution to (10) forms an innerapproximation of the maximal robust nontermination set, as stated in Corollary 1.
Corollary 1
For any function satisfying (10), is an innerapproximation of the maximal robust nontermination set , i.e. .
Proof 4
Let be a solution to (10). It is evident that satisfies the constraints:
(11) 
From Corollary 1, it is clear that an approximation of from inside, i.e. a robust nontermination set, is able to be constructed by addressing (10). The solution to (10) could be addressed by gridbased numrical methods such as level set methods [12, 32], which are a popular method for interface capturing. Such gridbased methods are prohibitive for systems of dimension greater than four without relying upon specific system structure. Besides, we observe that a robust nontermination set could be searched by solving (11) rather than (10). In the subsection that follows we relax (11) as a sumofsquares decomposition problem in a semidefinite programming formulation, which falls within the convex programming framework and can be efficiently solved by interior point methods when in Program 1, s are polynomials over and , state and uncertain input constraints, i.e. s and s, are restricted to basic semialgebraic sets.
IiiB Semidefinite Programming Implementation
In practice, it is nontrivial to obtain a solution to (2), and thus nontrivial to gain . In this subsection, thanks to (11) and Corollary 1, we present a semidefinite programming based method to solve (10) approximately and construct a robust invariant as presented in Defintion 4 when Assumption 1 holds.
Assumption 1
, , is polynomial over and , and , , are restricted to basic semislgebraic sets in Program 1.
Firstly, (11) has indicator functions on the expression , which is beyond the capability of the solvers we use. We would like to obtain a constraint by removing indicators according to Lemma 2.
Lemma 2 ([8])
Suppose and , where , , and , , . Also, and are respectively disjoint. Then, if and only if (pointwise)
(13) 
Consequently, according to Lemma 2, the equivalent constraint without indicator functions of (11) is equivalently formulated below:
(14) 
Before encoding (14) in sumofsquares programming formulation, we denote the set of sum of squares polynomials over variables by , i.e.
Besides, we define the set of states being reachable from the set within one step computation, i.e.,
(15) 
which can be obtained by semidefinite programming or linear programming methods as in [24, 31]. Herein, we assume that it was already given. Consequently, when Assumption 1 holds and in (14) is constrained to polynomial type and is restricted in a ball , where and , (14) is relaxed as the following sumofsquares programming problem:
(16) 
where , is the vector of the moments of the Lebesgue measure over indexded in the same basis in which the polynomial with coefficients is expressed, , , , , , , are sumofsquares polynomials of appropriate degree. The constraints that polynomials are sumofsquares can be written explicitly as linear matrix inequalities, and the objective is linear in the coefficients of the polynomial ; therefore problem (16) is reformulated as an semidefinite program, which falls within the convex programming framework and can be solved via interiorpoints method in polynomial time (e.g., [52]). Note that the objective of (16) facilitate the gain of the less conservative robust nontermination set.
The implementation based on the sumofsquares program (16) is sound but incomplete. Its soundness is presented in Theorem 3.
Theorem 3 (Soundness)
Proof 5
Since satisfies the constraint in (16), we obtain that satisfies according to procedure in [5]:
(17)  
(18) 
Assume that there exist an initial state and an input policy such that does not hold for . Due to the fact that (18) holds, we have the conclusion that and thus . Let be the first time making violate the constraint , i.e., and for . Also, since , (19) and (18), where is defined in (15), we derive that and , which contradicts (19). Thus, every possible run of Problem 1 initialized in will live in forever while respecting .
Therefore, the conclusion in Theorem 3 is justified.
Iv Experiments
In this section we evaluate the performance of our method built upon the semidefinite program (16). The first two examples, i.e. Examples 2 and 3, are constructed to illustarte the soundness of our method. The third one, i.e. Example 4, is used to evaluate the scalability of our method in dealing with Problem 1 . The parameters that control the performance of our approach in applying (16) to these three examples are presented in Table I. All computations were performed on an i77500U 2.70GHz CPU with 32GB RAM running Windows 10. For numerical implementation, we formulate the sum of squares problem (16) using the MATLAB package YALMIP^{1}^{1}1It can be downloaded from https://yalmip.github.io/. [29] and use Mosek^{2}^{2}2For academic use, the software Mosek can be obtained free from https://www.mosek.com/. [33] as a semidefinite programming solver.
In the following examples, we adapt the following numerical simulation techniques to evaluate the quality of the computed robust nontermination set by solving (16): Given an initial condition , where , the assignments of in the loop body in Program 1 are executed times in total, where is assigned randomly a value from in each iteration of the loop. We repeat this procedure times for the same initial condition. If all values of belong to in this process, we approximately regard the initial condition in the robust nontermination set. We sample multiple initial states randomly from and obtain an estimation of the maximal robust nontermination set by applying the above procedure to every initial state. Although such simulation technique can not give the exact estimation of the set , it is able to provide an insight of the set if the sample initial states are large enough. In the following examples, we regard the estimation obtained via our simulation techniques as . For Examples 2 and 3, we take initial states, and .
Ex.  

1  14  14  14  14  14  11.30 
1  16  16  16  16  16  28.59 
2  6  12  12  12  6  9.06 
2  8  16  16  16  8  65.22 
2  10  20  20  20  10  123.95 
2  12  24  24  24  12  623.95 
4  4  4  4  4  4  58.56 
4  5  4  4  4  4  60.02 
Example 2
This simple example is mainly constructed to illustrate the difference between Program 1 taking uncertain inputs into account and free of disturbances. In both cases, Program 1 is composed of a single loop without switchcase type in loop body, i.e. and .
In case that , and in Program 1, the innerapproximations of the maximal robust nontermination set are illustrated in Fig. 3(Left) when and . By visualizing the results in Fig. 3, the innerapproximation obtained when does not improve the one when a lot. Although there is a gap between the innerapproximations obtained via our method and the set , it is not big.
In the ideal implementation of Program 1, that is, in the loop body is a fixed nominal value, there will exists some initial conditions such that Program 1 in the real implemenation may violate the constraint set , i.e. Program 1 may terminate. We use as an instance to illustrate such situation. The difference between termiantion sets is visualized in Fig. 3(Right). The robust nontermination set in case of is smaller than the nontermination set when . Note that from Fig. 4, we observe that the innerapproximation obtained by our method when can approximate very well.
Example 3
In this example we consider Program 1 with switchcase type in the loop body, where , , , , and . The innerapproximations computed by solving (16) when and respectively are illustrated in Fig. 5. By comparing these results, we observe that polynomials of higher degree facilitate the construction of less conservative estimation of the set .