Robust Verification of Numerical Software

# Robust Verification of Numerical Software

Bai Xue and Naijun Zhan and Yangjia Li and Qiuye Wang
State Key Lab. of Computer Science, Institute of Software, CAS, China.
Email:
{xuebai,znj,yangjia,wangqye}@ios.ac.cn
###### Abstract

Numerical software are widely used in safety-critical 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 round-off errors and additive perturbations from real-world 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 real-time systems having hard deadlines on transaction execution time, and propose a method for robust conditional nontermination analysis, which can be used to under-approximate 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.

Numerical Software; Nontermination Analysis; Robust Verification

## I Introduction

Software is ubiquitous in mission-critical and safety-critical 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 64-bit floating point to a 16-bit 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 mission-critical and safety-critical 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 Floyd-Hoare-Naur inductive assertion approach [13, 18, 34], basically consisting of pre- and post-condition 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 well-known that the termination or non-termination problem is undecidable, and even not semi-decidable in general. Thus, more practical approaches to termination analysis are either to present some sufficient conditions for termination, or some sufficient conditions for non-termination, 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 round-off errors and additive perturbations from real-world 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/non-termination analysis approaches. In [54], the authors presented the following example:

###### Example 1

Consider a simple loop

 Q1:while (Bx>0) {x:=Ax},

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 non-termination 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 control-theoretic framework based on Lyapunov invariants to conduct verification of numerical software. Non-termination analysis proves that programs, or parts of a program, do not terminate. Non-termination 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 real-time 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 switch-case type loop body which are equivalent to constrained piecewise discrete-time dynamical systems subjected to time-varying uncertainties. To address the problem of robust conditional nontermination, we by the aid of the constrained piecewise discrete-time 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 state-constrained 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 discrete-time system in each mode is polynomial, and the state and uncertain input constraints are semi-algebraic, the optimal control problem is relaxed as a semi-definite programming problem, to which its polynomial solution forms an inner-approximation 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 long-standing 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 sum-of-squares programs, which are notoriously hard to solve. In order to solve the bilinear sum-of-squares 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 sum-of-squares 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 semi-definite programming based method to compute semi-algebraic 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 non-termination 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 discrete-time 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.

### Ii-a 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 switch-case 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 round-off 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 .

### Ii-B Piecewise Discrete-time Systems

In this subsection we interprete Program 1 as constrained piecewise discrete-time dynamical systems with uncertain inputs. Formally,

###### Definition 1

A constrained piecewise discrete-time dynamical system (PS) is a quadruple with

1. is the condition on initial states;

2. 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 ;

3. with as interpreted in Program 1;

4. is the set of uncertain inputs;

5. 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

 xπx0(l+1)=f(xπx0(l),π(l)), (1)

where , , and

 f(x,d)=1X1⋅f1(x,d)+⋯+1Xk⋅fk(x,d)

with , , representing the indicator function of the set , i.e.

 1Xi:={1,if x∈Xi,0,if x∉Xi.

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 non-terminating w.r.t. an initial state , if

 ∀π∈Π. ∀l∈N.  xπx0(l)∈X0. (2)

Now, we define our problem of deciding a set of initial states rendering Program 1 robust non-termination.

###### 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 non-terminating w.r.t. for any . We call is robustly non-terminating w.r.t. the maximal robust non-termination 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 Non-Termination Set Generation

In this section we elucidate our approach of addressing the problem of robust conditional nontermination for Program 1, i.e. synthesizing robust non-termination sets as presented in Definition 4. For this sake, we firstly in Subsection III-A characterize the maximal robust non-termination set by means of the value function, which is a solution to a suitable infinite horizon state-constrained optimal control problem based on PS (1). Any lower semicontinuous solution to this optimal control problem generates a robust non-termination set. Then, in the case that , , is polynomial over and , and the constraint sets over and , i.e. , , and , are of the basic semi-algebraic form, the semi-definite program arising from sum-of-squares decompositions facilitates the gain of inner-approximations of via solving the relaxation of the derived optimal control problem in Subsection III-B.

### Iii-a Characterization of R0

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:

 V(x0):=supπ∈Πsupl∈Nmaxj∈{1,…,n0}{h0,j(xπx0(l))}. (3)

Note that may be neither continuous nor semi-continuous. (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.

###### Proof 1

Let . According to Definition 4, we have that

 ∀i∈N. ∀π∈Π. ∀j∈{1,…,n0}. h0,j(xπy0(i))≤0. (4)

holds. Therefore, and thus .

On the other side, if , then , implying that (4) holds. Therefore, .

This concludes that .

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:

 V(x0)=supπ∈Πmax{V(xπx0(l)),supi∈[0,l)∩Nmaxj∈{1,…,n0}h0,j(xπx0(i))}. (5)
###### Proof 2

Let

 W(l,x0):=supπ∈Πmax{V(xπx0(l)),supi∈[0,l)∩Nmaxj∈{1,…,n0}h0,j(xπx0(i))}. (6)

We will prove that for , .

According to the definition of , i.e. (3), for any , there exists an input policy such that

 V(x0)≤supi∈Nmaxj∈{1,…,n0}{h0,j(xπ′x0(i))}+ϵ1.

We then introduce two infinite uncertain input policies and such that with for and with . Now, let , then we obtain that

 W(l,x0)≥max{V(y),supi∈[0,l)∩Nmaxj∈{1,…,n0}h0,j(xπ1y(i))}≥max{supi∈[l,+∞)∩Nmaxj∈{1,…,n0}{h0,j(xπ2x0(i−l))},supi∈[0,l)∩Nmaxj∈{1,…,n0}{h0,j(xπ1x0(i))}}=max{supi∈[l,+∞)∩Nmaxj∈{1,…,n0}{h0,j(xπ′x0(i))},supi∈[0,l)∩Nmaxj∈{1,…,n0}{h0,j(xπ′x0(i))}}=supi∈Nmaxj∈{1,…,n0}{h0,j(xπ′x0(i))}≥V(x0)−ϵ1. (7)

Therefore,

 V(x0)≤W(l,x0)+ϵ1. (8)

On the other hand, by the definition of , for any , there exists a such that

 W(l,x0)≤max{V(xπ1x0(l)),supi∈[0,l)∩Nmaxj∈{1,…,n0}{h0,j(xπ1x0(i))}}+ϵ1.

Also, by the definition of , i.e. (3), for any , there exists a such that

 V(y)≤supi∈Nmaxj∈{1,…,n0}{h0,j(xπ2y(i))}+ϵ1,

where . We define such that for and for . Then, it follows

 W(l,x0)≤2ϵ1+max{supi∈N∩[l,∞)maxj∈{1,…,n0}{h0,j(xπ2y(i−l))},supi∈[0,l)∩Nmaxj∈{1,…,n0}{h0,j(xπ1x0(i))}}≤supi∈[0,+∞)∩Nmaxj∈{1,…,n0}{h0,j(xπx0(i))}+2ϵ1≤V(x0)+2ϵ1. (9)

Combining (8) and (9), we finally have , implying that since is arbitrary. This completes the proof.

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

 min{infd∈D(V(x0)−V(f(x0,d))),V(x0)−maxj∈{1,…,n0}h0,j(x0)}=0. (10)
###### Proof 3

It is evindent that (10) is derived from (5) when .

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 inner-approximation of the maximal robust nontermination set, as stated in Corollary 1.

###### Corollary 1

For any function satisfying (10), is an inner-approximation of the maximal robust nontermination set , i.e. .

###### Proof 4

Let be a solution to (10). It is evident that satisfies the constraints:

 {u(x0)−u(f(x0,d))≥0,∀x0∈Rn,∀d∈D,u(x0)−h0,j(x0)≥0,∀x0∈Rn,∀j∈{1,…,n0} (11)

Assume . According to (11), we have that for , and ,

 {u(xπx′0(l+1))≤u(xπx′0(l))≤u(x′0)h0,j(xπx′0(l))≤u(xπx′0(l))≤u(x′0). (12)

Therefore, , implying that . Thus, .

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 grid-based numrical methods such as level set methods [12, 32], which are a popular method for interface capturing. Such grid-based 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 sum-of-squares 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 semi-algebraic sets.

### Iii-B Semi-definite Programming Implementation

In practice, it is non-trivial to obtain a solution to (2), and thus non-trivial to gain . In this subsection, thanks to (11) and Corollary 1, we present a semi-definite 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 semi-slgebraic 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)

 k′⋀i=1l′⋀j=1[Fi∧Gj⇒f′i≤g′j]∧k′⋀i=1[Fi∧(l′⋀j=1¬Gj)⇒f′i≤0]∧l′⋀j=1[(k′⋀i=1¬Fi)∧Gj⇒0≤g′j]. (13)

Consequently, according to Lemma 2, the equivalent constraint without indicator functions of (11) is equivalently formulated below:

 k⋀i=1[∀d∈D. ∀x0∈Xi. u(x0)−u(fi(x0,d))≥0]∧n0⋀j=1[∀x0∈Rn. u(x0)−h0,j(x0)≥0]. (14)

Before encoding (14) in sum-of-squares programming formulation, we denote the set of sum of squares polynomials over variables by , i.e.

 SOS(y):={p∈R[y]∣p=r∑i=1q2i,qi∈R[y],i=1,…,r}.

Besides, we define the set of states being reachable from the set within one step computation, i.e.,

 Ω(X0):={x∣x=f(x0,d),x0∈X0,d∈D}∪X0, (15)

which can be obtained by semi-definite 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 sum-of-squares programming problem:

 minu,sXii,l1,sDi,l2,si,l,s′1,jc′⋅wu(x)−u(fi(x,d))+ni∑l1=1sXii,l1hi,l1(x)+nk+1∑l2=1sDi,l2hk+1,l(d)−si,1h(x)∈SOS(x,d),u(x)−h0,j(x)−s′1,jh(x)∈SOS(x),i=1,…,k,j=1,…,n0, (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 sum-of-squares polynomials of appropriate degree. The constraints that polynomials are sum-of-squares 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 semi-definite program, which falls within the convex programming framework and can be solved via interior-points 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 sum-of-squares program (16) is sound but incomplete. Its soundness is presented in Theorem 3.

###### Theorem 3 (Soundness)

Let be solution to (16), then is an inner-approximation of , i.e., every possible run of Program 1 starting from a state in does not terminate.

###### Proof 5

Since satisfies the constraint in (16), we obtain that satisfies according to procedure in [5]:

 k⋀i=1[∀d∈D. ∀x∈Xi∩B. u(x)−u(fi(x,d))≥0]∧ (17) n0⋀j=1[∀x∈B. u(x)−h0,j(x)≥0]. (18)

Due to (17) and the fact that , we obtain that for ,

 ∃i∈{1,…,k}.  ∀d∈D.  u(x0)−u(fi(x0,d))≥0,

implying that

 u(x0)−u(f(x0,d))≥0,∀d∈D. (19)

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 semi-definite 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 i7-7500U 2.70GHz CPU with 32GB RAM running Windows 10. For numerical implementation, we formulate the sum of squares problem (16) using the MATLAB package YALMIP111It can be downloaded from https://yalmip.github.io/. [29] and use Mosek222For academic use, the software Mosek can be obtained free from https://www.mosek.com/. [33] as a semi-definite 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 .

###### 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 switch-case type in loop body, i.e. and .

In case that , and in Program 1, the inner-approximations of the maximal robust nontermination set are illustrated in Fig. 3(Left) when and . By visualizing the results in Fig. 3, the inner-approximation obtained when does not improve the one when a lot. Although there is a gap between the inner-approximations 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 inner-approximation obtained by our method when can approximate very well.

###### Example 3

In this example we consider Program 1 with switch-case type in the loop body, where , , , , and . The inner-approximations 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 .

###### Example 4

In this example, we consider Program 1 with seven variables and illustrate the scalability of our approach. In Program 1, ,