Frédéric Mazenc
email: frederic.mazenc@l2s.centralesupelec.fr

Saeed Ahmed
email: ahmed@ee.bilkent.edu.tr

Michael Maliso↵
email: malisoff@lsu.edu.

Finite Time Estimation Through a Continuous-Discrete Observer ⇤

Keywords: observer, output feedback, switched system, stabilization

We study two broad classes of nonlinear time-varying continuous time systems with outputs. For the first class, we build an observer in the case where a state dependent disturbance a↵ects the linear approximation. When the disturbances are the zero functions, our observer provides exact values of the state at all times larger than a suitable finite time, and it provides an approximate estimate when there are nonzero disturbances, so our observers are called finite time observers. We use this construction, which is of interest for its own sake, to design a globally exponentially stabilizing dynamic output feedback for a family of nonlinear systems whose outputs are only available on some finite time intervals. Our simulations illustrate the e cacy of our methods.

Introduction

State estimation of nonlinear systems from output measurements is a basic concern in robotics [START_REF] Cedervall | Nonlinear observers for unicycle robots with range sensors[END_REF], chemical and biochemical processes [START_REF] Dochain | State and parameter estimation in chemical and biochemical processes: a tutorial[END_REF], biomedical systems [START_REF] Manchester | A nonlinear observer for on-line estimation of the cerebrospinal fluid outflow resistance[END_REF], communication systems [START_REF] Perruquetti | Finite time observers: application to secure communication[END_REF], automotive systems [START_REF] Phanomchoeng | Nonlinear observer for bounded Jacobian systems, with applications to automotive slip angle estimation[END_REF], networked control systems [START_REF] Wang | Observer design for networked control systems with FlexRay[END_REF], and many other fields. Due to this strong motivation, various techniques to achieve state observation of nonlinear systems have been discussed in the literature. These techniques have included canonical form observers [START_REF] Bestle | Canonical form observer design for non-linear time-variable systems[END_REF], high-gain observers (as in [START_REF] Hammouri | High gain observer based on a triangular structure[END_REF], [START_REF] Khalil | High-gain observers in nonlinear feedback control[END_REF], and [START_REF] Zemouche | On LMI conditions to design observers for Lipschitz nonlinear system[END_REF]), Lyapunov based observers (as in [START_REF] Freidovich | Lyapunov-based switching control of nonlinear systems using high-gain observers[END_REF] and [START_REF] Kazantzis | Nonlinear observer design using Lyapunov's auxiliary theorem[END_REF]), and extended Kalman and Luenberger observers (as in [START_REF] Wang | High-gain observers with limited gain power for systems with observability canonical form[END_REF] and [START_REF] Zeitz | The extended Luenberger observer for nonlinear systems[END_REF]).

The above mentioned observer design techniques have the common disadvantage that they only ensure asymptotic convergence of the estimation error to zero, whereas finite time convergence of estimation errors to zero is often desirable for control and supervision purposes; see [START_REF] Du | Recursive design of finite-time convergent observers for a class of time-varying nonlinear systems[END_REF], [START_REF] Menhold | Finite time convergent observers for linear time-varying systems[END_REF] and [START_REF] Menhold | Finite time convergent observers for nonlinear systems[END_REF], and see [START_REF] Song | Time-varying feedback for finite-time robust regulation of normal-form nonlinear systems[END_REF] for more motivation for finite time control. One can distinguish between two broad families of finite time converging observers, namely, the families composed of nonsmooth observers without delay and the family of the smooth observers with delays. Non-smooth finite time converging observers have been proposed for instance in [START_REF] Du | Recursive design of finite-time convergent observers for a class of time-varying nonlinear systems[END_REF], [START_REF] Levant | Robust exact di↵erentiation via sliding mode technique[END_REF], [START_REF] Menard | A global high-gain finite time observer[END_REF], [START_REF] Perruquetti | Finite time observers: application to secure communication[END_REF], and [START_REF] Shen | Global finite time observers for Lipschitz nonlinear systems[END_REF]. The main drawback of nonsmooth finite time observers is that their lack of smoothness may generate poor robustness performance, but they have no delay, which is an advantage because the presence of a delay may complicate the implementation of an observer. Another possible drawback is due to the fact that the time of convergence of each trajectory depends on its initial condition. This is not the case for the observers that use artificial delays, for which their instant of convergence of the solutions is independent of the initial conditions and can be rendered as small as desired by the selection of a parameter (called the artificial delay). Finite time converging smooth observer designs have been introduced more recently. They were first presented in [START_REF] Engel | A continuous time observer which converges in finite time[END_REF], which was only applicable to linear time invariant (LTI) systems; see also [START_REF] Ra↵ | An impulsive observer that estimates the exact state of a linear continuous time system in predetermined finite time[END_REF] for finite time observers for LTI systems. An extension of the design presented in [START_REF] Engel | A continuous time observer which converges in finite time[END_REF] was carried out in [START_REF] Menhold | Finite time convergent observers for linear time-varying systems[END_REF] and [START_REF] Menhold | Finite time convergent observers for nonlinear systems[END_REF] for linear time varying (LTV) systems and nonlinear systems in observer canonical (normal) form, respectively. A generalized finite time converging observer design technique for nonlinear systems was proposed in [START_REF] Sauvage | Design of a nonlinear finite time converging observer for a class of nonlinear systems[END_REF] and was applicable to noncanonical form nonlinear systems as well. This approach was developed through a Lyapunov based observer from [START_REF] Kazantzis | Nonlinear observer design using Lyapunov's auxiliary theorem[END_REF].

The preceding finite time observer design approaches were carried out without considering disturbances in the measurements or dynamics. Such disturbances are usually present in practical applications and they a↵ect measurements and can be state dependent. Motivated by this fact, a finite time state estimation algorithm for nonlinear systems with bounded and time-varying disturbances in the dynamics and measurements was recently proposed in [START_REF] Mazenc | Estimation of solutions of observable nonlinear systems with disturbances[END_REF]. The design approach used in [START_REF] Mazenc | Estimation of solutions of observable nonlinear systems with disturbances[END_REF] is similar to that of [START_REF] Sauvage | Design of a nonlinear finite time converging observer for a class of nonlinear systems[END_REF]. Other issues pertaining to finite time converging observer design for nonlinear systems have been discussed in [START_REF] Menard | A global high-gain finite time observer[END_REF], [START_REF] Shen | Uniformly observable and globally Lipschitzian nonlinear systems admit global finite time observers[END_REF], and [START_REF] Shen | Semi-global finite time observers for nonlinear systems[END_REF].

The present paper has two main objectives. Our first aim is to complement [START_REF] Mazenc | Estimation of solutions of observable nonlinear systems with disturbances[END_REF] and [START_REF] Sauvage | Design of a nonlinear finite time converging observer for a class of nonlinear systems[END_REF], by proposing a finite time converging observer design for Lipschitz nonlinear systems of the form

⇢ ẋ(t) = [A + ✏(t, x(t))]x(t) + f (t, y(t), u(t)) y(t) = Cx(t) (1)
where A 2 R n⇥n and C 2 R q⇥n are constant matrices, the state x is valued in R n , the input u is valued in R p , the output y is valued in R q , f is a nonlinear function which is assumed to be locally Lipschitz with respect to y and piecewise continuous in its other arguments, the locally Lipschitz function ✏ : [0, +1)⇥R n ! R n⇥n can represent a disturbance, and the dimensions are arbitrary. Systems of the form (1) arise in many engineering contexts, e.g., in the modeling of vibrations of elastic membranes; see our examples section below.

The key di↵erence between the nonlinear systems in [START_REF] Mazenc | Estimation of solutions of observable nonlinear systems with disturbances[END_REF] and [START_REF] Sauvage | Design of a nonlinear finite time converging observer for a class of nonlinear systems[END_REF] and ( 1) is the presence of the function ✏. This disturbance significantly increases the di culty of constructing a finite time observer, since it makes it impossible to apply [START_REF] Mazenc | Estimation of solutions of observable nonlinear systems with disturbances[END_REF] or [START_REF] Sauvage | Design of a nonlinear finite time converging observer for a class of nonlinear systems[END_REF] to (1) and it seems that there is no direct way to extend them to [START_REF] Bestle | Canonical form observer design for non-linear time-variable systems[END_REF]. The nonlinear term ✏ is worth considering because (i) disturbances of this type often a↵ect systems and (ii) this term will enable us to use our finite time observer to construct dynamic output feedback for a broad family of nonlinear time-varying systems with temporary loss of output measurements. Very few works design finite time observers for Lipschitz nonlinear systems; see, e.g., [START_REF] Shen | Global finite time observers for Lipschitz nonlinear systems[END_REF] for a finite time observer for Lipschitz nonlinear systems under homogeneity conditions that are not required here. An advantage of our observer approach is that the finite convergence time in our method for (1) is independent of the initial state. Our dynamic output feedback design for our class of nonlinear systems with temporary loss of output measurements is motivated by the fact that in many engineering applications, the state is not available for measurement, and the output measurements are only available intermittently, meaning there may be intervals during which there is no output measurement, e.g., due to communication failures in GPS-denied environments. Our strategy has the following steps. We use an assumption that is inspired by [START_REF] Yan | Stability analysis of switched time-delay systems[END_REF] and [START_REF] Yan | Robust stabilization of parameter varying time delay systems by switched controllers[END_REF] and our finite time observer design to construct a dynamic output feedback through a continuous-discrete observer; see, e.g., [START_REF] Mazenc | Design of continuous-discrete observers for time-varying nonlinear systems[END_REF] for continuous-discrete observers when the output values are available at all times instead of being intermittent. Our dynamic output feedback globally exponentially stabilizes the origin of the nonlinear systems with a temporary loss of measurements. To the best of our knowledge, the stabilization problem we described with temporary loss of output measurements had remained unsolved in the literature, even in the case of linear systems, so our proposed tools are of considerable independent interest.

Throughout the sequel, the notation will be simplified whenever no confusion can arise from the context. The dimensions of our Euclidean spaces are arbitrary unless otherwise noted. The Euclidean norm in R a , and the induced norm of matrices, are denoted by | • |, and we assume that the initial times for our solutions are 0. We define X t by X t (s) = X(t + s) for all choices of X, s  0, and t 0 for which the equality is defined. Let I denote an identity matrix of any dimension. A matrix is called Schur stable provided its spectral radius satisfies < 1. A constant matrix is called positive provided all of its entries are positive. Let | • | 1 denote the sup norm of any matrix valued function over its entire domain.

Three Lemmas

The following lemmas will be used to prove our main results. For their proofs, see Appendices B and D. Lemma 1. Let M 2 R n⇥n be an invertible matrix. Let N 2 R n⇥n be a matrix. Let n and m be two constants such that |M 1 |  m and |N |  n. Assume that mn < 1 .

(

) 2
Then the matrix M + N is invertible and

(M + N ) 1 M 1  m2 n 1 mn (3)
is satisfied.

Lemma 2. Let A 2 R n⇥n be a constant matrix. Consider the system

⇣(t) = [A + E(t)] ⇣(t) ( 4
)
where ⇣ is valued in R n and E : [0, +1) ! R n⇥n is a bounded locally Lipschitz function. Let denote the fundamental solution of the system (4). Then for all t 1 2 R and t 2 2 R such that t 1 t 2 , the inequality

(t 1 , t 2 ) e A(t1 t2)  |E| 1 (t 1 t 2 )e (|A|+|E|1)(t1 t2) (5)
is satisfied.

Lemma 3. Let M =
 ↵ ! be a positive Schur stable matrix. Then for each constant N 0, we can find constants c i > 0 for i = 1 and 2 and k 2 N such that for all piecewise continuous functions z i : [0, +1) ! [0, +1) for i = 1, 2 that satisfy

z 1 (t)  ↵ sup s2[t N ,t] z 1 (s) + sup s2[t N ,t] z 2 (s) z 2 (t)  sup s2[t N ,t] z 1 (s) + ! sup s2[t N ,t] z 2 (s) (6)
for all t N and that have finite left limits z i (t ) at each point t 0 for i = 1 and 2, we have

|z(t)|  c 1 e c2t sup s2[0,kN ] |z(s)| for all t kN .
3 Finite Time Observer

Statement of Result

In this section, we complement the papers [START_REF] Mazenc | Estimation of solutions of observable nonlinear systems with disturbances[END_REF] and [START_REF] Sauvage | Design of a nonlinear finite time converging observer for a class of nonlinear systems[END_REF], where a finite time observer is provided, by allowing the more general class of systems (1) for general choices of the state dependent uncertainties ✏. We assume:

Assumption 1. The function f in (1) is locally Lipschitz with respect to y and u and piecewise continuous with respect to t. The pair (A, C) is observable and ( 1) is forward complete. The function ✏ is bounded and locally Lipschitz. Finally, u is piecewise continuous and locally bounded. ⇤

When Assumption 1 is satisfied, we can use [START_REF] Mazenc | Estimation of solutions of observable nonlinear systems with disturbances[END_REF]Lemma 1] to show that the observability of (A, C) implies that we can select a matrix L 2 R n⇥q with the following property: For each constant ⌧ > 0, there is a constant ⌧ 2 (0, ⌧ ) such that with the choice F = A + LC, the matrix

E(⌧ ) = e A⌧ e F ⌧ (7)
is invertible. This follows from an analytic continuity argument by first using [15, Lemma 1] to find a constant ⌧ a > 0 and a matrix L such that ( 7) is invertible with the choices F = A + LC and ⌧ = ⌧ a , and then noting that if there were a ⌧ 2 (0, ⌧ a ) such that det(E(r)) = 0 for all r 2 (0, ⌧ ), then real analyticity of det(E(r)) as a function of r would also give the contradiction det(E(⌧ a )) = 0. Our final assumption is:

Assumption 2.
There exist a positive constant ⌧ and a constant matrix L such that (i) the matrix E(⌧ ) defined in [START_REF] Hammouri | High gain observer based on a triangular structure[END_REF] with the choice F = A + LC is invertible and (ii) the bound

|✏| 1 ⌧ h e (|A|+|✏|1)⌧ + e (|F |+|✏|1)⌧ i |E 1 (⌧ )| < 1 ( 8
)
is satisfied. ⇤

We are ready to state and prove the following result:

Theorem 1. Let the system (1) satisfy Assumptions 1-2. Let

x(t) = E 1 (⌧ ) Z t t ⌧ h⇣ e (t ⌧ `)A e (t ⌧ `)F ⌘ f (`, y(`), u(`)) + e (t ⌧ `)F Ly(`) i d`(9)
for all t ⌧ , where E, F and ⌧ are from Assumption 2. Then, in terms of the constants

c(⌧ ) = |✏| 1 ⇣ e (|A|+|✏|1)⌧ + e (|F |+|✏|1)⌧ ⌘ ⇥ J(⌧ ) + |✏| 1 J(⌧ ) + |E 1 (⌧ )| ⌧ ⇤ ( 10
)
and

J(⌧ ) = ⌧ ⇥ e (|A|+|✏|1)⌧ + e (|F |+|✏|1)⌧ ⇤ |E 1 (⌧ )| 2 1 |✏| 1 ⌧ ⇥ e (|A|+|✏|1)⌧ + e (|F |+|✏|1)⌧ ⇤ |E 1 (⌧ )| , (11)
we have

|x(t) x(t)|  c(⌧ ) Z t t ⌧ |f (`, y(`), u(`))| + |L||y(`)| d`( 12
)
for all t ⌧ along all maximal solutions of the system (1).

Remark 1. In general, ⌧ ⇥ e |A|⌧ + e |F |⌧ ⇤ |E 1 (⌧ )|
does not converge to zero when ⌧ converges to zero, so Assumption 2 is a constraint on |✏| 1 . Since the inequality ( 12) holds for all t ⌧ , the function ( 9) is a finite time observer which gives an approximate value of the solution in finite time, and which agrees with the true state variable for all t ⌧ when ✏ = 0. To simplify, and in contrast with [START_REF] Mazenc | Estimation of solutions of observable nonlinear systems with disturbances[END_REF], we do not assume that there are disturbances in f and y. However, an extension to this case can be proved by combining the proof of Theorem 1 with the key ideas of [START_REF] Mazenc | Estimation of solutions of observable nonlinear systems with disturbances[END_REF].

Proof of Theorem 1

Since the system (1) is forward complete, it follows that for any initial condition, both x(t) and x(t) are defined over [⌧, +1). Fix any maximal solution x(t) of (1). The proof has two parts. In the first part, we consider the case where ✏ does not depend on x. In the second part, we use the result of the first part to show how the case where ✏ depends on x can also be handled. To simplify the notation, throughout the proof, we write f (t, y(t), u(t)) and g(t, y(t), u(t)) as f (t) and g(t) respectively, where g(t, y(t), u(t)) = f (t, y(t), u(t)) Ly(t) and L is from Assumption 2.

Let 1 (t, s) and 2 (t, s) denote the fundamental solutions of the systems

⇠1 (t) = [A + ✏(t)]⇠ 1 (t) and ⇠2 (t) = [F + ✏(t)]⇠ 2 (t) respectively. Let 1 (t, s) = 1 (t, s) 1 and 2 (t, s) = 2 (t, s) 1 . It is well known that @ > 1 @t (t, 0) = [A + ✏(t)] > 1 (t, 0) > (13)
and

@ > 2 @t (t, 0) = [F + ✏(t)] > 2 (t, 0) > (14)
hold for all t 0; see [START_REF] Sontag | Deterministic Finite-Dimensional Systems[END_REF]Appendix C.4]. Let

z i (t) = i (t, 0)x(t) for i = 1, 2. Then ż1 (t) = 1 (t, 0)[A + ✏(t)]x(t) + 1 (t, 0)f (t) + @ 1 @t (t, 0)x(t) = 1 (t, 0)f (t) (15)
for all t 0. By integrating this equality, we obtain

z 1 (t) = z 1 (t ⌧ ) + R t t ⌧ 1 (`, 0)f (`)d`. ( 16
)
Here and in the sequel, all inequality and equalities are for all t ⌧ unless otherwise noted. Then the semigroup property of the flow map 1 gives 1 (t, t ⌧ ) = 1 (t, 0) 1 1 (t ⌧, 0) and 1 (t, 0) 1 1 (`, 0) = 1 (t, `) for all `2 [t ⌧, t], and therefore also

x(t) = 1 (t, 0) 1 1 (t ⌧, 0)x(t ⌧ ) + R t t ⌧ 1 (t, 0) 1 1 (`, 0)f (`)d= 1 (t, t ⌧ )x(t ⌧ ) + R t t ⌧ 1 (t, `)f (`)d`. (17)
Observing that ẋ

(t) = [F + ✏(t)]x(t) + g(t, y(t), u(t)) (18)
and using variation of parameters gives

x(t) = 2 (t, t ⌧ )x(t ⌧ ) + Z t t ⌧ 2 (t, `)g(`)d`. (19)
Left multiplying the second line of ( 17) and ( 19) by 1 (t, t ⌧ ) and 2 (t, t ⌧ ) respectively and computing the di↵erence of the results gives

[ 1 (t, t ⌧ ) 2 (t, t ⌧ )] x(t) = R t t ⌧ 1 (t, t ⌧ ) 1 (t, `)f (`)d` R t t ⌧ 2 (t, t ⌧ ) 2 (t, `)g(`)d= R t t ⌧ 1 (`, t ⌧ )f (`)d` R t t ⌧ 2 (`, t ⌧ )g(`)d`, (20)
where the second equality used the semigroup property of the i 's. Therefore,

[E(⌧ ) + Z(t)] x(t) = R t t ⌧ 1 (`, t ⌧ )f (`)d` R t t ⌧ 2 (`, t ⌧ )g(`)d`(21) holds with Z(t) = 1 (t, t ⌧ ) 2 (t, t ⌧ ) E(⌧
) and E defined in [START_REF] Hammouri | High gain observer based on a triangular structure[END_REF]. From ( 13)-( 14), and by Lemma 2 with the choice E = ✏, we deduce that

1 (`, t ⌧ ) e A(` t+⌧ )  ✏(` t + ⌧ )e (|A|+✏)(` t+⌧ ) 2 (`, t ⌧ ) e F (` t+⌧ )  ✏(` t + ⌧ )e (|F |+✏)(` t+⌧ ) (22)
for all `2 [t ⌧, t], where we set ✏ = |✏| 1 for brevity. It follows by setting `= t in [START_REF] Perruquetti | Finite time observers: application to secure communication[END_REF] that

|Z(t)|  ✏⌧ h e (|A|+✏)⌧ + e (|F |+✏)⌧ i (23)
for all t ⌧ . From this inequality and Lemma 1 (applied with M = E(⌧ ) and N = Z(t)) we deduce that (8) from Assumption 2 ensures that for all t ⌧ , the matrix

E(⌧ ) + Z(t) is invertible. Also, Lemma 1 implies that |(E(⌧ ) + Z(t)) 1 E 1 (⌧ )|  ✏J(⌧ ) for all t ⌧
, where J is the constant defined in [START_REF] Levant | Robust exact di↵erentiation via sliding mode technique[END_REF]. Then, omitting the argument ⌧ of E to keep our notation simple, [START_REF] Pearson | Nonlinear system identification with limited time data[END_REF] gives

x(t) = [E + Z(t)] 1 R t t ⌧ ⇥ e (t ⌧ `)A f (`) e (t ⌧ `)F g(`) ⇤ d+ [E + Z(t)] 1 R t t ⌧ ⇥ 1 (`, t ⌧ ) e (t ⌧ `)A f (`) + e (t ⌧ `)F 2 (`, t ⌧ ) g(`) ⇤ d`. (24)
We deduce that x(t) = x(t) x(t) satisfies

|x(t)|  |E 1 [E + Z(t)] 1 | R t t ⌧ e (t ⌧ `)A f (`) e (t ⌧ `)F g(`) d+ | [E + Z(t)] 1 | R t t ⌧ ⇥ 1 (`, t ⌧ ) e (t ⌧ `)A |f (`)| + e (t ⌧ `)F 2 (`, t ⌧ ) |g(`)| ⇤ d ✏J R t t ⌧ e (t ⌧ `)A f (`) e (t ⌧ `)F g(`) d+ (✏J + |E 1 |) R t t ⌧ ⇥ 1 (`, t ⌧ ) e (t ⌧ `)A |f (`)| + e (t ⌧ `)F 2 (`, t ⌧ ) |g(`)|
⇤ dẁ here we omitted the dependency of J on ⌧ for brevity. Using [START_REF] Perruquetti | Finite time observers: application to secure communication[END_REF], and setting ✏] = ✏(✏J

+ |E 1 |), we obtain |x(t)|  ✏J R t t ⌧ e (t ⌧ `)A f (`) e (t ⌧ `)F g(`) d+ ✏] R t t ⌧ ⇥ (` t + ⌧ )e (|A|+✏)(` t+⌧ ) |f (`)| + (` t + ⌧ )e (|F |+✏)(` t+⌧ ) |g(`)| ⇤ d ✏J R t t ⌧ ⇥ e ⌧ |A| |f (`)| + e ⌧ |F | |f (`) Ly(`)| ⇤ d+ ✏] ⌧ R t t ⌧ ⇥ e (|A|+✏)⌧ |f (`)| + e (|F |+✏)⌧ |f (`) Ly(`)| ⇤ d ✏J R t t ⌧ ⇥ e ⌧ |A| + e ⌧ |F | |f (`)| + e ⌧ |F | |L||y(`)| ⇤ d+ ✏] ⌧ R t t ⌧ ⇥ e (|A|+✏)⌧ + e (|F |+✏)⌧ |f (`)| + e (|F |+✏)⌧ |L||y(`)| ⇤ d`. (25)
This concludes the first part of the proof, by our choice [START_REF] Khalil | High-gain observers in nonlinear feedback control[END_REF] of the constant c(⌧ ). We now use the preceding result to cover the case where ✏ depends on both t and x to complete the proof of the theorem. Fix any specific solution of ( 1), which we denote by x \ (t). Then we consider the system

⇢ Ẋ(t) = [A + ✏(t, x \ (t))]X(t) + f (t, Y (t), u(t)) Y (t) = CX(t) . (26)
For the system (26), ✏(t, x \ (t)) depends only on t and not on X. Moreover, it is bounded by ✏. Thus, from our previous proof, we deduce that for all solutions X of ( 26) and for all t ⌧ , we have

| X(t) X(t)|  c(⌧ ) Z t t ⌧ |f (`, Y (`), u(`))| + |L||Y (`)| d`( 27
)
with c(⌧ ) as defined in the first part of the proof and with

X(t) = E 1 (⌧ ) Z t t ⌧ h⇣ e (t ⌧ `)A e (t ⌧ `)F ⌘ f (`, Y (`), u(`)) + e (t ⌧ `)F LY (`) i d`. (28)
Since x \ (t) is a solution of (1), it follows that (26) holds with X = x \ . From [START_REF] Shen | Global finite time observers for Lipschitz nonlinear systems[END_REF], it follows that

|x \ (t) x \ (t)|  c(⌧ ) Z t t ⌧ [|f (`, Cx \ (`), u(`))| + |L||Cx \ (`)|] d`( 29
)
with

x\ (t) = E 1 (⌧ ) Z t t ⌧ h⇣ e (t ⌧ `)A e (t ⌧ `)F ⌘ f (`, Cx \ (`), u(`)) + e (t ⌧ `)F LCx \ (`) i d`. (30)
Since x \ was an arbitrary solution of (1), this concludes the proof of Theorem 1.

Stabilization of Systems with Temporary Loss of Measurements

Assumptions and Statement of Result

Theorem 1 relies on the assumption that the output is available for all t 0 and the fact that the matrix A in (1) does not depend on t. In this section, we relax these assumptions. We assume that the output is only available on some specific intervals of time, and we will consider systems whose linear approximation at the origin is time-varying even when no disturbance is acting and y and u are set to zero. Under these assumptions, complemented by a stabilizability assumption of ISS type and a mild restriction on f , we construct a globally exponentially stabilizing dynamic output feedback with an observer for the system

ẋ(t) = [M(t) + ⌘(t, x(t))]x(t) + B(t)u(t) (31)
with x valued in R n for any n 2 N and u valued in R p for any p 2 N. The function u represents a control. We assume that there are constants P > 0 and ✓ 2 (0, P ) such that the R p valued output

y(t) = Cx(t) (32)
is only available when t is in the set

S P,✓ = [ i2N [iP, iP + ✓]. (33)
We assume:

Assumption 3. The functions ⌘, M, and B are locally Lipschitz and bounded, and M is not the zero function. ⇤

For any locally Lipschitz bounded function Z : [0, 1) ! R n⇥p , we can therefore fix a constant k > 0 and choose s matrices A j 2 R n⇥n for j = 1, 2, . . . , s such that the function

D(t, j, x) = M(t) Z(t)C A j + ⌘(t, x) , (34)
possesses the following property: For each i 2 N, there is a j 2 {1, ..., s} (depending on i) such that for all t 2 [iP, iP + ✓] and all x 2 R n , the inequality

|D(t, j, x)|  k (35)
is satisfied. In terms of the matrix C from (31), we also assume that the A i 's from Assumption 3 satisfy: Assumption 4. For all i 2 {1, ..., s} the pair (A i , C) is observable. ⇤ Assumption 4 ensures that there are matrices L i and values ⌧ i 2 (0, ✓) for i = 1 to s such that each of the matrices

E i (⌧ i ) = e Ai⌧i e Fi⌧i where F i = A i + L i C (36)
is invertible for i = 1, ..., s. This follows from the analytic continuity argument from Section 3.1.

For later use, we introduce the matrices

H i = E 1 i (⌧ i ) (37)
and upper bounds for the matrices A i , F i , L i and H i of the form

|A i |  A , |F i |  F , |L i |  L , and |H i |  H ( 38
)
for all i 2 {1, ..., s}.

Assumption 5. There are a locally Lipschitz function u s and a constant µ > 0 such that

|u s (t, x)|  µ|x| ( 39
)
for all t 0 and x 2 R n , and constants  2 (0, 1) and > 0, such that for any piecewise continuous function , all solutions of the system

ẋ(t) = [M(t) + ⌘(t, x(t))]x(t) + B(t)u s (t, x(t) + (t)) (40) satisfy |x(t)|  |x(s)| + sup w2[s,t] | (w)| ( 41
)
for all t 0 and s 2 [t P ✓, t P ]. ⇤ Assumption 5 can be viewed as a generalized Hurwitzness type condition, because in the specific case where u s (t, x) = Kx for some matrix K such that M(t) + ⌘(t, x(t)) + B(t)K is constant and Hurwitz, condition (41) can be checked using variation of parameters. In terms of the constants

ḡ = k ⇣ e (A+k)⌧ + e (F +k)⌧ ⌘ ⇥ J + (kJ + H)⌧ ⇤ (42)
and (47)

J = ⌧ h e (A+k)⌧ + e (F +k)⌧ i H 2 1 k⌧ h e (A+k)⌧ + e (F +k)⌧ i H and ⌧ = max{⌧ i : 1  i  s}, (43)
Thus, the following theorem can be interpreted as a result for switched systems:

Theorem 2. Let the system (31) satisfy Assumptions 3 to 6 and consider the continuous-discrete system

ẋ(t) = M(t)x(t) + B(t)u s (t, x(t)) , t 2 [t i , t i+1 ) x(t i ) = H (ti) (⌧ i ) R ti ti ⌧i h E(t i ⌧ i `, (t i )) Z(`)y(`) + B(`)u s (`, x(`)) + e (ti ⌧i `)F (t i ) L (ti) y(`) i d`( 48
)
where E(m, k) = e mA k e mF k , H (ti) is the matrix defined in (37) and t i = iP + ✓ for integers i 2 N. Then the system (31) in closed loop with the dynamic output feedback

u s (t, x(t)) (49)
is such that the dynamics for (x, x) are globally exponentially stable to 0, where x(t) = x(t) x(t).

Remark 2. Since we chose values ⌧ i 2 (0, ✓), (48) only requires y(t) values at times t 2 S P,✓ . Theorem 2 implies that the dynamic output feedback (49) asymptotically stabilizes the closed loop system to 0. Assumption 3 is inspired by the technique of [START_REF] Yan | Stability analysis of switched time-delay systems[END_REF] and [START_REF] Yan | Robust stabilization of parameter varying time delay systems by switched controllers[END_REF], which makes it possible to stabilize time-varying systems with auxiliary switched systems whose subsystems are time-invariant systems that are a↵ected by disturbances. Notice that if |⌘| 1 = 0, and if M and Z are continuous and periodic, then one can have an arbitrarily small constant k by choosing a su ciently large number of matrices A i , which can facilitate satisfying Assumption 6.

Remark 3. Assumption 5 is a stabilizability assumption. Constructing a feedback u s such that (41) is satisfied can be challenging. However, see [START_REF] Maliso↵ | Constructions of Strict Lyapunov Functions[END_REF] and [START_REF] Mazenc | Asymptotic stabilization for feedforward systems with delayed feedback[END_REF] for techniques for constructing feedback controls for nonlinear time-varying systems. The main di↵erence between (31) and the system studied in [START_REF] Mazenc | Design of continuous-discrete observers for time-varying nonlinear systems[END_REF] is that M depends on t, but [START_REF] Mazenc | Design of continuous-discrete observers for time-varying nonlinear systems[END_REF] does not provide exponentially stabilizing feedbacks and so does not apply to the problems we study here.

Remark 4.

If |⌘| 1 = 0 and k can be chosen arbitrarily small, then by choosing both ⌧ and k su ciently small, the matrix S is Schur stable since  2 (0, 1). In this case, there exists a ⌘ > 0 such that the corresponding matrix S is Schur stable if |⌘| 1  ⌘.

Remark 5. The size of the constant k depends on ✓, and our condition (44) is more easily satisfied when k is su ciently small. Thus, since an arbitrarily small constant ✓ can be selected, it may be useful to choose a ✓ value resulting in a small constant k > 0. Remark 6. Theorem 2 can be extended to the case where C in the output depends on t. For the sake of simplicity, we do not investigate this case. Notice that Z(t) can always be chosen identically equal to zero. However, to decrease the conservatism of the approach, it is worth introducing the function Z that can be freely chosen.

Proof of Theorem 2

Let i 0 be an integer. Bearing in mind (47), we deduce that the system (31) in closed loop with u s (t, x(t)) can be rewritten as

ẋ(t) = [A (ti) + D(t, (t i ), x(t))]x(t) + Z(t)Cx(t) + B(t)u s (t, x(t)) ( 50
)
for all t 2 [iP, (i + 1)P ) and i 0, and our linear growth condition (39) on u s ensures that the closed loop system is forward complete. Assumptions 3 to 6 imply that for each choice of i, Assumptions 1-2 are satisfied by

⇢ ẋ(t) = [A (ti) + D [ (t, (t i ), x(t))]x(t) + Z(t)y(t) + B(t)u(t) y(t) = Cx(t) , (51)
where 9) for suitable choices of f and other functions from ( 9), we deduce from Theorem 1 that

u(t) = u s (t,
|x(t i )|  ḡ Z ti ti ⌧i ⇥ |B(`)u s (`, x(`))| + |Z| 1 + |L (ti) | |y(`)| ⇤ d`(52)
with ḡ defined in (42) and x = x x. As an immediate consequence of the bounds (38), we obtain

|x(t i )|  ḡ R ti ti ⌧i ⇥ |B| 1 µ|x(`)| + (L + |Z| 1 )|C||x(`)| ⇤ d` R ti ti ⌧i [p 1 |x(`))| + p2 |x(`)|] d`(53)
with p1 and p2 defined in (45).

On the other hand, ( 31) and (48) imply that, for all t 2 [t i , t i+1 ) and integers i 0, we have

ẋ(t) = M(t)x(t) ⌘(t, x(t))x(t) (54)
when we use the control (49). This gives

| ẋ(t)|  |M| 1 |x(t)| + |⌘| 1 |x(t)|, which we can integrate to obtain |x(t)|  ⇣ |x(t i )| + R t ti |⌘| 1 |x(`)|d`⌘ + R t ti |M| 1 |x(`)|d`. (55)
Hence, Gronwall's inequality gives

|x(t)|  ⇣ |x(t i )| + R t ti |⌘| 1 |x(`)|d`⌘ e |M|1(t ti)  e |M|1(t ti) |x(t i )| + |⌘| 1 R t ti |x(`)|e |M|1(t `)e |M|1(` ti) d e |M|1(t ti) |x(t i )| + |⌘| 1 e P |M|1 R t ti |x(`)|e |M|1(t `)d`. (56)
We can combine the last inequality in (56) with (53) to obtain

|x(t)|  e |M|1(t ti) R ti ti ⌧ [p 1 |x(`))| + p2 |x(`)|] d`+ |⌘| 1 e P |M|1 R t ti e |M|1(t `)|x(`)|df
or all t 2 [t i , t i+1 ), and so also

|x(t)|  p1 e |M|1(t ti) ⌧ sup w2[ti ⌧,ti] |x(w)| + h p2 e |M|1(t ti) ⌧ + |⌘| 1 e P |M|1 e |M|1 P 1 |M|1 i sup w2[ti ⌧,t] |x(w)| ,
where we used the fact that

Z t ti e |M|1(t `)|x(`)|d` 1 |M| 1 ⇣ e |M|1(t ti) 1 ⌘ sup w2[ti,t] |x(w)|  1 |M| 1 ⇣ e |M|1P 1 ⌘ sup w2[ti ⌧,
for all t P + ⌧ . The system (60) is of the form (6) from Lemma 3 (with (z 1 , z 2 ) = (|x(t)|, |x(t)|), N = P + ⌧ , and M = S), save for the fact that the matrix S from ( 46) is Schur stable but is not necessarily a positive Schur stable matrix (since some of the entries of S can be zero). However, we can majorize all of the entries of S by positive values, in such a way that the new positive matrix that we obtain is a positive Schur stable matrix. This follows from the continuity of the eigenvalues of a matrix as functions of the entries of the matrix, by increasing the entries of S by adding small enough positive constants to the entries. Then it follows from Lemma 3 that (x(t), x(t)) converges exponentially to the origin, which proves the theorem.

Illustrations

Our classes of dynamics ( 1) and ( 31) cover a broad class of systems that were beyond the scope of existing finite time observer approaches. In this section, we illustrate our theorems using the controlled Mathieu equation q(t)

+ (R 1 + R 2 cos(t))q(t) + u(t) = 0 (61)
for positive constants R 1 and R 2 , which arises in the study of vibrations of an elliptic membrane. The Mathieu equation has also been studied in [START_REF] Pearson | Nonlinear system identification with limited time data[END_REF] to illustrate parameter identification for a certain family of nonlinear and time varying systems using data over a limited time interval. See also [START_REF] Insperger | Stability chart for the delayed Mathieu equation[END_REF] for the study of the uncontrolled Mathieu equation corresponding to cases where u = 0. The work [START_REF] Insperger | Stability chart for the delayed Mathieu equation[END_REF] studied domains of stability and instability including Hopf bifurcations along the boundaries of the domains of stability.

Illustration of Theorem 1

The controlled Mathieu equation (61) can be written as

⇢ ẋ(t) = [A + ✏(t, x(t))]x(t) + f (t, y(t), u(t)) y(t) = Cx(t) (62)
with

x =  q q , A =  0 1 R 1 0 , C = e > 1 , ✏(t) =  0 0 R 2 cos(t) 0
, and

f (t, y, u) = e 2 u (63)
where e i is the ith standard basis vector for i = 1 and 2. We consider the case where only an upper bound on R 2 is known. Choosing

L =  0 2R 1 (64)
we can satisfy Assumptions 1-2 for many choices of the parameters R 1 , R 2 , and ⌧ . For instance, if we choose R 1 = 1, any R 2 2 (0, 0.024], and ⌧ = 1, then the matrix E(⌧ ) = e A⌧ e F ⌧ with the choice F = A + LC has determinant det(E(⌧ )) = 0.33254 and so is invertible, and

|✏| 1 ⌧ h e (|A|+|✏|1)⌧ + e (|F |+|✏|1)⌧ i |E 1 (⌧ )| = 0.998418 < 1, (65)
which ensures that Assumptions 1-2 are satisfied for any piecewise continuous locally bounded choice of u.

Then the finite time observer is provided by ( 9) from Theorem 1.

To illustrate our result, Figure 1 shows a MATLAB simulation of our finite time observer (9), using an integration algorithm of the model (62) with a semi-implicit integration step of 0.001, and with the initial conditions x 1 (0) = x 2 (0) = 2 and x1 (0) = x2 (0) = 2, where u(t) = sin(2t), and using the preceding choices of A, f , and L. Since our simulation shows good tracking performance of the estimator x2 for the state component x 2 of (62), it helps to illustrate our general theory in the special case of the system (62).

Illustration of Theorem 2

We rewrite the controlled Mathieu equation (61) as ⇢ ẋ(t) = [M(t) + ⌘(t, x(t))]x(t) + B(t)u(t)

y(t) = Cx(t) (66)
with

x =  q q , M(t) =  0 1 R 1 R 2 cos(t) 0 , C = e > 1 , ⌘(t, x) = 0, and B(t) = e 2 . ( 67
)
We choose P = 2⇡. We consider the case where R 1 = 1 and R 2 = 0.024, and we apply Theorem 2 with s = 1, and we use ⌧ to denote the constant ⌧ 1 from Theorem 2. The preceding choices of ⌘, M, and B ensure that Assumption 3 is satisfied. Let

A 1 =  0 1 R 1 R 2 0 and Z(t) =  0 R 2 (1 cos(t)) . ( 68
)
Then Assumption 4 is satisfied, and our choice of Z gives D(t, j, x) = M(t) Z(t)C A j + ⌘(t, x) = 0. Consequently, our requirements on the functions D are satisfied with k = 0. We choose

L 1 =  0 2R 1 and F 1 = A 1 + L 1 C =  0 1 R 1 R 2 0 . ( 69
)
Then the matrix

E 1 (⌧ ) = e A1⌧ e F1⌧ is E 1 (⌧ ) =  cos (! 1 ⌧ ) 1 !1 sin (! 1 ⌧ ) ! 1 sin (! 1 ⌧ ) cos (! 1 ⌧ )  cosh (! 2 ⌧ ) 1 !2 sinh (! 2 ⌧ ) ! 2 sinh (! 2 ⌧ ) cosh (! 2 ⌧ ) (70)
where

! 1 = p R 1 + R 2 and ! 2 = p R 1 R 2 .
This can be checked by noting that the terms in the di↵erence in (70) are the fundamental matrix solutions of Ṁ = A 1 M and Ṁ = F 1 M , respectively. Thus

E 1 (⌧ ) =  cos (! 1 ⌧ ) cosh (! 2 ⌧ ) 1 !1 sin (! 1 ⌧ ) + 1 !2 sinh (! 2 ⌧ ) ! 1 sin (! 1 ⌧ ) + ! 2 sinh (! 2 ⌧ ) cos (! 1 ⌧ ) cosh (! 2 ⌧ ) . ( 71
)
Then

det E 1 (⌧ ) = 2  1 cos (! 1 ⌧ ) cosh (! 2 ⌧ ) R 2 ! 1 ! 2 sin (! 1 ⌧ ) sinh (! 2 ⌧ ) .
One can easily check that det E 1 (0) = 0 and

d d⌧ det E 1 (⌧ ) = 2R 1 ⇣ 1 p R1+R2 sin( p R 1 + R 2 ⌧ ) cosh( p R 1 R 2 ⌧ ) 1 p R1 R2 cos( p R 1 + R 2 ⌧ ) sinh( p R 1 R 2 ⌧ ) ⌘ > 0 (72) for all ⌧ 2 ⇣ 0, ⇡ 2 p 1.024 ⌘ . ( 73
)
This follows by noting that (72) is equivalent to positivity of 73) holds, which follows because S(0) = 0 and S 0 > 0 on this interval. It follows that E 1 (⌧ ) is invertible when ⌧ 2 (0, ⇡/4]. To check Assumption 5, we select the feedback u s (t,

S(⌧ ) = 1 !1 tan(! 1 ⌧ ) 1 !2 tanh(! 2 ⌧ ) when (
x) = [1 R 1 R 2 cos(t)]x 1 + 2x 2 .
Then the system which corresponds to (40) is ẋ(t) = Gx(t) e 2 a (t), with

G =  0 1 1 2 and a (t) = (1 R 1 R 2 cos(t)) 1 (t) + 2 2 (t) . (74)
One can check that for all ` 0, we have

e G`= e ` 1 + `` `1 ` , (75)
and therefore also for any choices of ⌧ 2 (0, ✓). One can also compute ḡ = p1 = p2 = 0. We deduce that the matrix

S =  0 0 1+4⇡ e 2⇡ (79)
is Schur stable, so Assumptions 3-6 are satisfied for any ⌧ 2 (0, ⇡/4]. Thus, we can apply Theorem 2.

To illustrate our result, Figure 2 shows a MATLAB simulation of the closed loop system with the controller (49) and the continuous-discrete observer (48) using an integration algorithm of the model (66) with a semiimplicit integration of step 0.001 with initial conditions x 1 (0) = x 2 (0) = 2 and x1 (0) = x2 (0) = 2, with ✓ = ⇡/4 and ⌧ = 0.5. We again show the estimate x2 tracking the state component x 2 , and the closed loop system performance and the control values. Since our simulation shows good stabilization and tracking performance, it helps to illustrate our general theory in the special case of the system (66).

Conclusions

We provided new constructions of observers and output feedback controls for time-varying nonlinear systems with intermittent output observations and disturbances. Our feedback control result proved exponentially stable convergence of the closed loop system to the desired equilibrium. This is valuable, because it is common in engineering to encounter systems for which there are periods during which no output measurements are available for use in the control. The presence of the disturbances and nonlinearities makes the observer design much more challenging than standard observer design problems. Our main strategy to overcome these challenges combined finite time observers and a switched systems approach. Our applications of our theorems to the Mathieu equation exhibited the good performance of our methods in simulations. Many extensions can be expected for systems with delay, and for local stability of broader classes of systems.

Appendices A. Results for Matrices

This appendix proves two lemmas, the first of which will be used to prove the second one. The second lemma in this appendix will be used in our proof of Lemma 2. The following lemma includes the inequality |e N I|  e |N | 1 for any square matrix as the special case where M = 0:

Lemma 4. The inequality |e M +N e M |  ⇣ 1 e |N | ⌘ e |M |+|N | (A.1)
is satisfied for any matrices M 2 R n⇥n and N 2 R n⇥n .

Proof. Let U 1 = M and U 2 = N . Then for each pair of integers (i, j) with i 0 and j 1, there is a function j,i : {1, ..., j} ! {1, 2} such that j,1 (k) = 1 for all k 2 {1, ..., j} and such that

(M + N ) j = 2 j P i=1 U j,i (1) ...U j,i (j) (A.2)
which gives

(M + N ) j M j = 2 j P i=2 U j,i (1)
...U j,i (j) . (A.

3)

It follows that

|(M + N ) j M j |  2 j X i=2 |U j,i (1) |...|U j,i (j) | = (|M | + |N |) j |M | j . (A.4)
Observing that Proof. We start the proof with some definitions. Set S(i) = (↵ 1 + . . . + ↵ i )G for all i 1, and set ⌦ 0 = I and ⌦ i = e ↵i(G+⌫i) .e ↵i 1 (G+⌫i 1 ) ....e ↵1(G+⌫1) for i 1, (A.8)

e M +N e M = +1 P j=1 (M +N ) j M j j!  +1 P j=1 |(M+N) j M j | j! (A.
and set ⇠ 0 = 0 and ⇠ i = ⌦ i e S(i) for all i 2 {1, ..., m}. Then the left side of (A. with ⇤ r (t, s) = r (t, s) e Since > 0 is arbitrary, and since we can make r as large as we wish (and therefore p ⇤ as small a positive value as we wish), we deduce that from (A.36) that (t, s) e A(t s)  e (t s)(|A|+|E|1) 1 e (t s)|E|1 (A.37) for all t 2 (T 0 , T ] and s 2 [T 0 , t). Since the constants T and T 0 are arbitrary, we can conclude because 1 e s  s for all s 0.

Figure 1 :

 1 Figure 1: Simulations of finite-time observer (9) for (62): Component x 2 and its estimate x2

  choice P = 2⇡ gives (41) with the choice  = max `2[P,P +✓] e `(1 + 2`) = e P (1 + 2P ) = 1 + 4⇡ e 2⇡ . (77) In terms of the notation from Assumptions 3-6, one has and one can choose k = 0, µ = 2.024, A = |M| 1 = 1.024, |B| 1 = |C| = 1, L = 2, and F = 1, which give k⌧ h e (A+k)⌧ + e (F +k)⌧ i H = 0 < 1 (78)

some j 2 1 P

 21 7) is |⇠ m |. Consider any i 2 {1, ..., m}. Elementary calculations give⇠ i = e ↵i(G+⌫i) ⌦ i 1 e S(i) = ⇥ e ↵i(G+⌫i) e ↵iG ⇤ ⌦ i 1 + e ↵iG ⇠ i 1 (A.9)and therefore also|⇠ i |  e ↵i(G+⌫i) e ↵iG |⌦ i 1 | + e ↵ig |⇠ i 1 |  1 e ⌫ e g+⌫ |⌦ i 1 | + e g |⇠ i 1 | (A.10)where the last inequality is a consequence of Lemma 4 and the fact that ↵ i 2 [0, 1] for all i. Since for all i 2 {0, ..., m}, we have |⌦ i |  e i(g+⌫) , we obtain|⇠ m |  1 e ⌫ e m + e g |⇠ m 1 | |⇠ m 1 |  1 e ⌫ e (m1) + e g |⇠ m 2 | . . . |⇠ 1 |  1 e ⌫ e (A.11) with = g + ⌫. A simple induction argument then gives |⇠ m |  1 e ⌫ {1, 2, . . . , m 1}, then (A.11) gives |⇠ j+1 |  e ḡ (1 e ⌫ ) j k=0 e kḡ e (j k) + (1 e ⌫ )e (j+1) = (1 e ⌫ ) j P k=1 e kḡ e (j k+1) + (1 e ⌫ )e (j+1) = 1 e ⌫ j P k=0 e kg e (j+1 k) , (A.14)which proves the inductive step. Then the geometric sum formula implies that the inequality (A.7) holds.

  {1, ..., s} is any integer such that |D(t, j, x)|  k for all t 2 [iP, iP + ✓] and all x 2 R n , it follows from our choices of k and the A i 's that for all t 2 S P,✓ and x 2 R n , we have |D(t, (t), x)|  k.

	Assumption 6. The inequality	k⌧	h	e (A+k)⌧ + e (F +k)⌧	i	H < 1	(44)
	is satisfied. Moreover, with the choices
		p1 = ḡ|B| 1 µ and p2 = ḡ |B| 1 µ + (L + |Z| 1 )|C| ,	(45)
	the matrix	S =	"	p1 e |M|1P ⌧ p2 e |M|1P ⌧ + |⌘| 1 e P M1 e |M|1 P 1 |M|1	#	(46)
								
	is Schur stable.							⇤
	Defining the switching signal : [0, +1) ! {1, ..., s} by (t) = j for all t 2 [iP, (i + 1)P ) and integers 0, where j 2 Notice that the system (31) can be rewritten as i
		ẋ(t) = [A (t) + D(t, (t), x(t))]x(t) + Z(t)Cx(t) + B(t)u(t) .
	our final assumption is:

  x(t)), and where D [ is any locally Lipschitz function that equals D on [iP, iP + ✓] ⇥ { (t i )} ⇥ R n and that is bounded by k. (We introduce the function D [ because k is not necessarily a global bound on D but provides the required bound for t 2 [iP, iP + ✓].) Therefore, since the xi formula from (48) agrees with the x(t) formula from (

  Let m 1 be an integer. Let G 2 R n⇥n and ⌫ i 2 R n⇥n for all i 2 {1, ..., m} be matrices such that |G|  g and max i |⌫ i |  ⌫, where g > 0 and ⌫ > 0 are two constants. Let ↵ i 2 [0, 1] be a constant for each i 2 {1, ..., m}. Then e ↵m(G+⌫m).e ↵m 1 (G+⌫m 1 ) ....e ↵1(G+⌫1) e ↵ 4 G  e m(g+⌫) 1 e m⌫ (A.7) is satisfied with the choice ↵ 4 = ↵ 1 + . . . + ↵ m .

									5)
	we deduce that	e M +N	e M 	+1 P j=1	(|M |+|N |) j |M | j j!	= e |M |+|N |	e |M |	.	(A.6)
	This allows us to conclude.
	Lemma 5.

  A(t s) . Since t 2 (T 0 , T ] and s 2 [T 0 , t), there exists l 2 N such that t 2 [t l , t l+1 ) and there exists p 2 N such thats 2 [t p , t p+1 ) with p  l. If p = l, then |⇤ r (t, s)| = | r (t, s) e A(t s) | = |e (t s)(A+E(t l ))e A(t s) | .On the other hand, if p < l, then we obtain⇤ r (t, s) = r (t, t l ) r (t l , tl1 ).... r (t p+1 , s) e A(t s) = e (t t l )(A+E(t l )) e p⇤(A+E(t l 1 )) ....e p⇤(A+E(tp+1)) e (tp+1 s)(A+E(tp)) e A(t s) = e t t l p⇤ p⇤(A+E(t l )) e p⇤(A+E(t l 1 )) ....e p⇤(A+E(tp+1)) e Since (t t l )/p ⇤ , (t p+1 s)/p ⇤ , and p ⇤ are all contained in [0, 1], it follows from Lemma 5 with the choices m = ` p + 1, G = p ⇤ A and ⌫ i = p ⇤ E(t i ) that |⇤ r (t, s)|  e (` p+1)p⇤(|A|+|E|1) 1 e (` p+1)p⇤|E|1

								(A.32)
	From Lemma 4, it follows that
	|⇤ r (t, s)| 	⇣	1 e (t s)|E|1	⌘	e (t s)(|A|+|E|1)	.	(A.33)
								(A.34)
						t p+1 s p⇤	p⇤(A+E(tp))	e A(t s)	.
	 e 2p⇤(|A|+|E|1)	e (t s)(|A|+|E|1) 1 e 2p⇤(|A|+|E|1)	e (t s)|E|1	,	(A.35)

where we used the fact that (` p 1)p ⇤  t s. This inequality and (A.

[START_REF] Wang | High-gain observers with limited gain power for systems with observability canonical form[END_REF]

) give (t, s) e A(t s)  + e 2p⇤(|A|+|E|1) e (t s)(|A|+|E|1) 1 e 2p⇤(|A|+|E|1) e (t s)|E|1 . (A.36)

Mazenc and Ahmed were supported by the PHC Bosphore 2016 France-Turkey Project under project numbers 35634QM (France) and EEEAG-115E820 (T ÜB İTAK-The Scientific and Technological Research Council of Turkey). Maliso↵ was supported by US National Science Foundation Grants 1408295 1711299. Mazenc is with Inria, the Laboratoire des Signaux et Systèmes (L2S,

B. Proof of Lemma 3

Set z = (z 1 , z 2 ) > . Let V be a positive eigenvector associated with an eigenvalue 2 [0, 1) of M > (which exist by the Perron-Frobenius Theorem, since M > is a positive Schur stable matrix). Choose an integer k 1 such that k < 1/2, and set R = M k , p = k , and k (t) = (sup s2[t kN ,t] z 1 (s), sup s2[t kN ,t] z 2 (s)) > for all t kN . Then [START_REF] Freidovich | Lyapunov-based switching control of nonlinear systems using high-gain observers[END_REF] can be rewritten as z(t)  M 1 (t) (where inequalities of matrices are taken entry-wise), and we can prove (by induction on k) that z(t)  R k (t) for all t kN . We also have k (t)  2 sup s2[t kN ,t] (z 1 (s), z 2 (s)) > for all t kN , which follows because (z 1 (s 1 ), z 2 (s 2 ))  (z 1 (s 1 ), z 2 (s 1 )) + (z 1 (s 2 ), z 2 (s 2 )) for all s 1 and s 2 in [t kN , t] and all t kN . Since V > R = pV > , we conclude that

for all t kN . Since 2p 2 [0, 1), we can apply [17, Lemma 1] (with the choice w(`) = V > z(`+ kN )) to find a constant c 2 > 0 (that only depends on k, N , and p) such that V > z(t)  sup{V > z(`) : 0  ` kN }e c2(t kN ) for all t kN . The lemma now follows because all components of V are positive.

C. Lemma on Comparison of Flow Maps

This appendix provides a lemma that we later use to prove Lemma 2. Consider the system

with X valued in R n and where N : [0, +1) ! R n⇥n is a continuous function. Let ✏ : R ! R n⇥n be a piecewise continuous function that is bounded everywhere by some constant ✏ 0. Consider

Let and ✏ denote the fundamental solutions of the systems (A. [START_REF] Mazenc | Asymptotic stabilization for feedforward systems with delayed feedback[END_REF]) and (A.17), respectively. Lemma 6. Let T 0 and T T 0 be two real numbers. Let > 0 be a real number. There exists a constant ✏ > 0 such that if ✏  ✏, then for all t 2 [T 0 , T ] and s 2 [T 0 , t], the inequality

holds.

Proof. Observe for later use that the continuity of and 1 implies that there is a constant 0 such that for all a 2 [T 0 , T ] and b 2 [T 0 , a], we have

Through simple calculations, we obtain

where we also used [START_REF] Sontag | Deterministic Finite-Dimensional Systems[END_REF]

, where V 2 R n is any vector satisfying |V | = 1. Then the chain rule and (A.20) give

Using (A. [START_REF] Menhold | Finite time convergent observers for linear time-varying systems[END_REF]) and the upper bound of ✏, we deduce that for all t 2 [T 0 , T ] and s 2 [T 0 , t], we have

D. Proofs of Lemmas 1 and 2 D.1. Proof of Lemma 1

To prove that the matrix M + N is invertible, let us proceed by contradiction. We suppose that it is not invertible. Then there is a nonzero vector V 2 R n such that V > (M + N ) = 0, so invertibility of M gives

, and so also |V |  |V | mn. Since V 6 = 0, we conclude that 1  mn, which contradicts (2). We deduce that M + N is invertible. To prove the inequality (3), we first set R = (M + N ) 1 holds. Recalling that E is locally Lipschitz, let K E 0 be a Lipschitz constant for E(t) on [T 0 , T ]. Let r > T T 0 be a positive integer to be selected later. We define an increasing sequence of real numbers t i by t 0 = T 0 and t i+1 t i = p ⇤ where p ⇤ = (T T 0 )/r. Let r be a switching sequence defined by r (t) = t i when t 2 [t i , t i+1 ) for all integers i 0. Then

Consider any constant > 0. Let r denote the fundamental solution of the time-varying linear system ⇠(t) = [A + E( r (t))] ⇠(t). By letting r ! +1 in (A.29) and using Lemma 6 (with the choices N (t) = A+E(t) and ✏(t) = E( r (t)) E(t)), we deduce that there is an integer r c > 0 such that when r 2 N is such that r r c , then | (a


Abdurrahman Yas

Bora Uc
email: bora.ucar@ens-lyon.fr

Ümit V C ¸atalyürek

SINA: A Scalable Iterative Network Aligner

Given two graphs, network alignment asks for a potentially partial mapping between the vertices of the two graphs. This arises in many applications where data from different sources need to be integrated. Recent graph aligners use the global structure of input graphs and additional information given for the edges and vertices. We present SINA, an efficient, shared memory parallel implementation of such an aligner. Our experimental evaluations on a 32-core shared memory machine showed that SINA scales well for aligning large real-world graphs: SINA can achieve up to 28.5× speedup, and can reduce the total execution time of a graph alignment problem with 2M vertices and 100M edges from 4.5 hours to under 10 minutes. To the best of our knowledge, SINA is the first parallel aligner that uses global structure and vertex and edge attributes to handle large graphs.

I. INTRODUCTION

Fusion of data from different sources can be modeled as merging graphs. In such a case, given two non-isomorphic and comparable graphs, the first step is to identify a graph alignment, where a potentially partial mapping of vertices between two graphs is computed. Graph alignment is a well studied area. However, an important portion of related works are proposed for structurally very similar networks, such as protein protein interaction graphs [START_REF] Klau | A new graph-based method for pairwise global network alignment[END_REF], [START_REF] Singh | Pairwise global alignment of protein interaction networks by matching neighborhood topology[END_REF], [START_REF] Kuchaiev | Topological network alignment uncovers biological function and phylogeny[END_REF], [START_REF] Elmsallati | Global alignment of proteinprotein interaction networks: A survey[END_REF]. Furthermore, the majority of the existing work rely on computationally heavy optimization algorithms or pairwise similarity computations for all vertex pairs which cannot produce results for large graphs in a reasonable time. Hence, they are limited to graphs having less than tens of thousands vertices.

Most of the existing graph alignment techniques are based on the global topology of underlying graphs. More recent ones exploit additional information about vertex and edges for improved alignment, when such information is available. One recent approach is implemented in GSANA [START_REF] Yas | An iterative global structure-assisted labeled network aligner[END_REF] which uses global global structure of the graphs to reduce the problem size. GSANA has been shown to have higher recall than the state of the art algorithms while being orders of magnitudes faster. However, for larger graphs (a few million vertices and edges), its execution time is still in the orders of hours. In this work, we investigate an efficient and scalable parallelization of GSANA without sacrificing its recall. Our aim is to enable high recall alignment of very large graphs on shared memory systems. We implement the findings in a data-parallel, architecture (resource)-aware scalable graph alignment framework called SINA.

GSANA is an iterative algorithm, where each iteration has four pipelined steps. In the first step, shortest path distances to a set of special vertices (called anchor vertices) are computed. Anchor vertices are those for which a mapping is known. If such vertices are available in the given graph, then GSANA uses them; if not, it decides which vertices to use as anchor vertices. Then some of these anchor vertices are singled out and deemed more important. In the second step, GSANA places vertices into a 2D plane (unit circle to be more precise) using the shortest path distances to the important anchor vertices, with the idea that similar vertices will be placed closely. The 2D plane is subdivided into buckets using quadtrees in the third step. In the fourth step, the similarities of vertices that lie in the same or neighboring buckets are computed by combining a set of similarity measures, and a promising set of pairs is mapped greedily. The next iteration then starts with an enlarged set of anchors, where the additional anchors are chosen with the computed mapping. In SINA we parallelize all costly components of GSANA.

The contributions of this paper are as follows:

• We develop parallel algorithms to implement SINA, making the process scalable, and solve the attributed network alignment problem, so that large scale graphs can be handled in a reasonable time. • We propose vertex layout techniques for large scale graphs. These techniques are aimed at decreasing the total execution time by increasing the data access locality. • We optimize GSANA's sequential implementation.

SINA's sequential execution times are about 33 times faster than those reported in [START_REF] Yas | An iterative global structure-assisted labeled network aligner[END_REF]. Experimental results show that our proposed framework, SINA, reaches up to 28× speedup. To the best of our knowledge, SINA is the first parallel aligner that uses global structure and vertex and edge attributes. Therefore we could not compare SINA with other systems. We note that as SINA retains the recall of GSANA, it has higher recall than other state of the art methods.

The organization of the paper is as follows. We introduce the notation and give the formal problem definition in Section II. Then we discuss the parallelization approach in Section III. Section IV presents experimental results, which is followed by a summary of related work in Section V. We conclude the paper in Section VI with a summary and a list of future work.

Symbol Description

V Vertex set E Edge set TV Vertex type set TE Edge type set t[x] Type of the vertex x ∈ V or the edge x ∈ E AV Vertex attribute set AE Edge attribute set a[x] Attribute of the vertex x ∈ V or the edge x ∈ E Ni[u] Neighbor list of vertex u in graph Gi δ(u, v) Distance between u, v ∈ V σ(u, v) Similarity score for u ∈ V1 and v ∈ V2 µ[u] Mapping of u ∈ V1 in V2 S = S1 ∪ S2 Anchor (seed) set where S2 = {v : µ[u] = v, v ∈ S1}

II. PROBLEM DEFINITION AND BACKGROUND

A graph G = (V, E, T V , T E , A V , A E ) consists of a set V of vertices, a set E of edges, two sets T V and T E for vertex and edge types, and two sets A V and A E for vertex and edge attributes, where type and attribute sets can be empty. For instance in a Facebook graph, human, page, group can be used as vertex types, name, location can be defined as vertex attributes, and connection types and connection years can be defined as edge types and attributes respectively. An edge e is referred as e = (u, v) ∈ E, where u, v ∈ V . The neighbor list of a vertex u ∈ V is defined as

N [u] = {v ∈ V : (u, v) ∈ E}.
When discussing two graphs, we will use subscripts 1 and 2 to differentiate them if needed, and ignore those subscript when the intent is clear from the context. For example, N 1 [u] and N 2 [u ] will represent the neighbor lists of vertices u and u in G 1 and G 2 , respectively. Given a vertex x ∈ V or an edge x ∈ E, a[x] represents the set of attributes of x, and t[x] represents the type of x. We also use δ(u, v) to denote the breadth-first search (BFS) distance between vertices u and v. We assume that there are a number of vertices whose mappings are known. These are called anchor vertices, where the initial set of anchor vertices is denoted by S, which can be empty.

Given two different graphs G 1 and G 2 , the similarity score between two vertices u ∈

V 1 and v ∈ V 2 is denoted by σ : V 1 × V 2 → R, where R denotes the set of real numbers. We define µ[u] : V 1 → V 2 as an injective mapping, where µ[u] = v represents mapping of u ∈ V 1 to v ∈ V 2 . If a vertex u ∈ V 1 is not mapped to a vertex of G 2 , we use µ[u] = ⊥,
which is also refereed as nil mapping. Table I displays the notations used in this paper.

Definition 1 (Graph Alignment Problem). Given two graphs

G 1 = (V 1 , E 1 , . . . ) and G 2 = (V 2 , E 2
, . . . ), the graph alignment problem asks for an injective mapping that maximizes

∀u∈V1,µ[v] σ(u, µ[u]) , (1)
with the convention that σ(u, ⊥) = 0, when µ[u] = ⊥.

GSANA, hence SINA, uses quadtrees [START_REF] Finkel | Quad trees a data structure for retrieval on composite keys[END_REF] to partition the 2D plane, hence the vertices of the graphs. Simply put, quadtree is a geometric, tree-based data structure in which each internal node has four children. The points in the plane (in our case the vertices) are stored in the leaf nodes, where a capacity limit is imposed in the number of points in a leaf node. Initially there is a single node, containing the whole plane. Then, the plane is recursively split into four regions with a cross until each of the regions contains at most a predetermined number of points. When a region is split, the node representing that region becomes an internal node, and each of the four sub-regions are associated with a leaf node whose parent is the new internal node. The points that lie in the region are partitioned among the four leaf nodes according to their positions.

III. SINA: PARALLEL ITERATIVE NETWORK ALIGNMENT

The sequential GSANA algorithm is iterative and has four main steps: i) shortest path computations, ii) anchor selection, iii) partitioning, and iv) mapping (see Fig. 1). We parallelize all compute intensive parts which are marked in Fig. 1. We discuss the sequential versions of each of these steps as implemented in GSANA below, and then discuss SINA's parallelization approach in the following subsections. GSANA uses shortest path distances to anchor vertices to reduce the similarity computation requirements. In simple terms, if we know that

v 1 ∈ V 1 is identified with u 1 ∈ V 2 , then v 2 ∈ V 1 is more likely to be mapped to u 2 ∈ V 2 than to u 3 ∈ V 2 , if δ(u 2 , u 1 ) is closer to δ(v 2 , v 1 ) than δ(u 3 , u 1 ).
GSANA assigns unit lengths to edges so that the shortest path distances effectively correspond to the breadth-first search distances. We discuss the parallelization approach of SINA for this first step in Subsection III-A.

When there are a number of different anchor pairs, one has to combine the distances of a vertex pair to the anchor vertices to decide if the two vertices could be similar. GSANA identifies some anchors as better fit for determining similar vertices (these better serving anchor vertices are called vantage anchors in GSANA). Then, GSANA places the vantage anchor vertices on the unit circle in 2D, where an anchor pair is placed diametrically opposite places on the circle. Once the vantage anchors are placed in the circle, the other vertices are placed inside this circle using the distances to the anchors on the circle. Since nearby vertices will be checked for similarity, GSANA divides up this circle into buckets and confines the similarity computations among neighboring buckets. In order to have similar number of vertices per bucket, GSANA inserts the vertices in a quadtree (with an upper bound on the number of vertices per a leaf node of the quadtree) while placing them in the 2D plane. Instead of using one quadtree, SINA uses one quadtree per graph. We needed this deviation from GSANA for better definition of tasks for parallelism (covered in Subsection III-C). We discuss the parallelization approach of SINA for the anchor selection step in Subsection III-B.

Once the vertices are inserted into the buckets corresponding to the leafs of the quadtree, GSANA starts computing similarity between vertices that lie in a common or a neighboring bucket. Vertices in buckets that are not neighbors (even if they are geometrically close by) are not compared, as the alternatives in the neighboring buckets are closer. We discuss the parallelization approach of SINA for the similarity computation in Subsection III-C.

When similarity computation is done, top k most similar vertices are identified for each vertex. Mapping is the fourth step of GSANA, and the goal is to compute potentially a partial mapping between two graphs using identified vertices. This step takes only a tiny fraction of the total execution time (less than 0.5s. per iteration), therefore we do not parallelize it. SINA implements a few optimization techniques, such as data layout and data structure optimizations for increasing practical performance. These are discussed in more detail in Subsection III-D.

A. Parallelizing shortest path computations

SINA needs an efficient breadth first search (BFS) implementation to compute the shortest path distances from each anchor vertex. There are two ways to do this. One of them is to use existing parallel [START_REF] Beamer | Direction-optimizing breadth-first search[END_REF] or vectorized BFS implementations [START_REF] Sarıyüce | Regularizing graph centrality computations[END_REF] per BFS. The second one is to use a sequential BFS per anchor, and parallelize over anchors. Since the number of anchors is usually larger than the potential number of computing threads, yet relatively smaller with respect to the number of vertices (GSANA limits the number of anchors to 2000), this parallelization is more promising. In particular, the sequential tasks on a given graph (a BFS per anchor) have the same complexity and are likely to have similar practical run time. Furthermore, the parallelization overhead is negligible. Because of these reasons, SINA implements the second alternative described in Algorithm 1. As seen in this algorithm, for each newly introduced anchor u ∈ G 1 , a sequential BFS is run starting from u in a parallel for. Once the shortest path distances from the newly introduced anchor vertices of G 1 are computed, a BFS for each newly introduced anchor vertex in G 2 is run in a parallel for. After BFSs, the new set of vantage anchors are determined as in GSANA. We did not parallelize this step of choosing vantage anchors, as its cost is negligible.

B. Placing vertices into 2D and bucketing

SINA needs to place all vertices in a graph into a quadtree. In the quadtree data structure, we keep vertices in the leaf nodes and use internal nodes for routing. Leaf nodes are split into four when the number of vertices in them reaches to a pre-defined size limit, l. In our partitioning, we use two main operations of quadtrees; search and insert. Since the Algorithm 1: Computation of the shortest paths from the anchor vertices with BFS

D1.resize(|S|), D2.resize(|S|)

Memory allocation Computation of shortest paths from new anchors for each u ∈ S1 in parallel do if δ(u, .) is not computed before then δ(u, .) ← BF S(G1, u)

for each u ∈ S1 in parallel do if δ(µ[u], .) is not computed before then (δ(µ[u], .) ← BF S(G2, µ[u])
search operation is read-only, it can be done concurrently. Therefore as the first step, in parallel, SINA computes all vertices' coordinates (p) in the 2D space and instantiates an insert operation, defined in Alg. 2. This algorithm first finds the corresponding leaf node for the given position p. When the leaf node is found, the insertion operation starts. The insertion operation requires some additional effort for two reasons. First, multiple threads may want to insert a vertex into the same queue at the same time. Second, when a leaf node has to be split, children should be created safely by one thread, and all the vertices in the leaf node have to be moved into the newly created children's queues. To overcome the first issue, we use thread-safe queues in leaf nodes, therefore multiple threads can make parallel updates. For the second issue, when a thread observes that a leaf node has to be split, it first locks the leaf node and creates four children. Then, the leaf node is marked as an internal node, and other threads cannot insert any more vertices to that node, or split that node. Finally, the vertices in the queue are moved to the children nodes. Here, we note that more than one thread may reach the split operation at the same time but only one of them creates the children. Then these threads may move the vertices in the queue in parallel.

C. Similarity computation

In the fourth step, GSANA compares all vertices in a bucket of a graph with the vertices of the other graph that appear in the same or neighboring buckets. It checks neighbor buckets to make sure that vertices are close to the border of the buckets are handled appropriately. In SINA, we have two quadtrees. The first one contains the vertices of the first graph, and the second one contains those of the second graph. For each bucket B 1 ∈ QT 1 , we define an associated one in QT 2 by locating the bucket B 2 ∈ QT 2 containing the center of B 1 . Then, the similarity comparisons of vertices in B 1 ∈ QT 1 will be carried out for the vertices in the associated B 2 ∈ QT 2 , and B 2 's neighboring buckets in QT 2 .

We propose two coarse-grain task definitions in SINA. In the all comparison scheme (ALL), the similarity computation between vertices in a bucket of QT 1 and all related buckets in QT 2 are defined as one large coarse-grain task. The second scheme, pairwise comparison (PAIR), is still a coarse task definition but it is finer than the first one. In this scheme, the similarity computation between the vertices in a bucket of QT 1 and one of the related ones in QT 2 is defined as a task. These two definitions are depicted in Figure 2, and explained in more detail below. In ALL, a task is composed of similarity computations between the vertices of a non-empty bucket B 1 ∈ QT 1 , and the vertices of the all the related buckets in QT 2 . To find the related buckets, we first find the bucket B 2 ∈ QT 2 containing the center of B 1 . The related buckets are B 2 and the neighbors of B 2 .

In PAIR a task is composed of similarity computations between the vertices of a non-empty bucket B 1 ∈ QT 1 , and the vertices of a single related bucket B ∈ QT 2 .

Alg. 3 constructs the task lists based on the selected scheme. Then, it executes each task in parallel using COMPSIMALL or COMPSIMPAIR algorithms displayed in Algorithms 4 and 5, respectively. These functions compute the similarity scores for each vertex v ∈ B with every other vertex u ∈ B . GSANA only keeps the top k most similar vertices that are identified for each vertex, and these are stored in a priority queue In COMPSIMALL (Alg. 4) all of the similarity scores for each vertex v ∈ V 1 will be computed by a single thread. Therefore, computed scores can be safely inserted to the priority queue P [v] without any concurrency issues. In COMPSIMPAIR (Alg. 5), the top similarity scores for each vertex v ∈ V 1 might come from multiple threads. Therefore, for each task (i.e., for each bucket pairs that are compared), top k similar vertices first stored in a temporary queue, p. After top k similar ones computed for this task, if the highest similarity score in p is greater than the lowest similarity score in P [v], then we lock the P [v], insert all elements of p to P [v] and finally unlock it. At any time, both p and P [v] only keep just top k most similar vertices. The ALL scheme has two main drawbacks. First, the number of parallel tasks is limited by the number of buckets. Second, due to irregular nature of partitioning and very coarse-grain task composition, this scheme may lead to high load imbalance among tasks. In SINA (Alg. 3), tasks in the task list T are executed in parallel using a dynamic load balancing scheme (using OpenMP or CilkPlus). This will help toi reduce the observed load balance. Furthermore, one can simply sort tasks based on their loads in a non-increasing order, and use this largest-job first heuristic to further reduce the load imbalance, and hence the total execution time.

P [v] for v ∈ V 1 .
The PAIR scheme helps to reduce the load imbalance at the expense of additional synchronization overhead for insertion to the global priority queues. The use of local priority queues helps to reduce the synchronization overheads, but if two or more concurrent tasks have the same source bucket, this overhead may increase. To reduce the probability of such cases, SINA randomly shuffles the task list.

D. Other performance optimization techniques

We applied a few optimizations techniques to better exploit parallelism and improve data locality in SINA. Here, we cover these optimizations under two categories: data layout and data structure.

1) Data Layout: Exploiting data locality is one of the most common optimization techniques to improve the performance of irregular applications. The graphs that are input have already an ordering of the vertices, based on the alphabetical ordering of vertex labels. We call this the natural ordering (NAT). In SINA, we investigate two other layout techniques: "Breadth-First Search (BFS)" based and "Hilbert-Curve Based (HCB)" based. These techniques are presented in Fig. 3. a) BFS Layout: In this layout, shown in Fig. 3(b), the vertices are ordered based on their insertion time to the frontier queue in a BFS, starting from a vertex with the highest degree. This layout gives significant improvement in both similarity and shortest path computation steps over the natural layout. However in this scheme even if the vertices preserve some locality, appearances of the vertices in the quadtree buckets may not preserve BFS locality, which can cause loss of locality during similarity computation.

b) HCB Layout: In this layout, shown in Fig. 3(c), the vertices are ordered based on their Hilbert orders. When all vertices are inserted into the quadtree, SINA sorts the buckets based on a Hilbert order of the buckets. Then, SINA orders all vertices in a bucket according to the bucket's rank. We expect this layout to increase locality during the similarity computation step, where neighboring buckets are likely to be involved with the similarity computations regarding a given bucket in QT 1 .

2) Data-Structure optimizations: The similarity function σ used in GSANA is composed of several metrics. All of these metrics (listed in Table II) measures similarity based on a different property.

Symbol

Description τ : Type similarity α : Anchor similarity ∆ : Relative degree distance [START_REF] Koutra | Big-align: Fast bipartite graph alignment[END_REF] τ V : #Same/#T otal types of adjacent vertices τ E : #Same/#T otal types of adjacent edges C V : Vertex attribute similarity C E : Edge attribute similarity All of the mentioned similarity metrics, except τ and ∆, require accessing two vertices' neighborhood information. For instance, to compute τ V (u, v) a thread has to read

|N [u]|+|N [v]| different elements in T V
to be able to compare types of the neighboring vertices. This type of data access is highly irregular and cache unfriendly. Therefore to improve the cache use, SINA stores "weight" and "type" informations of the neighbors of a vertex in sorted vectors for each vertex.

IV. EXPERIMENTAL EVALUATION

We present several experiments in order to identify the performance trade-offs of the parallel algorithms and optimizations of SINA. Experiments were carried out on a machine that has two, 16-core Intel Xeon E5-2683 2.10GHz processors, 512GB of memory, 1TB disk space, running Ubuntu GNU/Linux with kernel 4.8.0. SINA is implemented in C++ and compiled with GCC 5.4.

A. Dataset

We use real-world graphs obtained from [START_REF] Boldi | WebGraph Datasets: Laboratory for algorithmics[END_REF], [START_REF] Demaine | BigDND: Big dynamic network data[END_REF], [START_REF] Trier | DBLP: Computer science bibliography[END_REF]. We also generated different size DBLP [START_REF] Trier | DBLP: Computer science bibliography[END_REF] graphs. The properties of graphs are listed in Table III, which are described below.

DBLP (2014-2017): As discussed in GSANA [START_REF] Yas | An iterative global structure-assisted labeled network aligner[END_REF], we downloaded two consecutive years of DBLP graphs from 2014 [START_REF] Demaine | BigDND: Big dynamic network data[END_REF] to 2017 [START_REF] Trier | DBLP: Computer science bibliography[END_REF]. The ground-truth is created using key attribute of author elements. Vertices are authors, and two authors have an edge if they have co-authored a paper. As in GSANA, we use publications' cross-ref information to create vertex attributes by splitting a cross-ref by '/' and unionizing initial character of each word as the vertex attribute. Edge attribute between two vertices is the mean of the publication years of co-authored papers between two authors. The other DBLP graphs (DBLP-17-small and DBLP-14-small), listed in Table III, are smaller subgraphs of the original DBLP graph, centered around the highest degree vertex.

Hollywood (2009-2011): Hollywood graphs from years 2009 and 2011 are downloaded from [START_REF] Boldi | WebGraph Datasets: Laboratory for algorithmics[END_REF]. In these graphs, the vertices are actors, and two actors are connected if they have acted in the same movie. We create the ground truth using author names. We give initials of each actor as the vertex attribute. This dataset does not have vertex/edge types nor edge attributes.

B. Workload

In order to investigate load imbalance due to variations of number of vertices in each bucket and variable degrees of the vertices, we plotted bucket and task workload distributions in Fig. 4. Here, the upper bound on the limit of vertices per bucket (i.e., leaf nodes of quadtree) is set to 512. Fig. 4(a) displays a histogram of number of buckets with different maximum number of vertices in a bucket. As expected, the vertex distribution among the buckets is not uniform.

A better estimate for workload for each task is the amount of data that needs to be read, at least once, in order to compute all required similarity scores. Computing similarity score between vertices v and u requires comparing the adjacency lists of those two vertices. Since we have multiple components of similarity scores, we need to multiply this number with a constant. For the sake of simplicity, we present the amount of data that needs to be read, in a single pass as

W (T = B, B ) = ∀v∈B ∀u∈B |N [v]| + |N [u]| . (2)
Using the formula (2) as a metric for workload estimation, in Figure 4(b) and Figure 4(c) we present workload estimations for ALL and PAIR parallelization schemes respectively. In both of these figures x-axis represents tasks in execution order, and y-axis represents size of memory read in MB. We observe from these figures that in PAIR, tasks have better workload since the granularity sizes are smaller. In ALL, some tasks require significantly larger memory reads and hence higher computation times, yielding poor scalability.

C. Effect of Task Granularity and Parallel Programing API

In the next experiment, we compared ALL and PAIR parallelization schemes using two different shared memory programming languages/API: OpenMP and CilkPlus. In order to highlight load imbalance issues, in this experiment we have used our smallest problem (DBLP-17-small and DBLP-14small graphs). Note that the sizes of these graphs are larger than the graphs used in prior studies. Figure 4(d) displays the speedup obtained using these four combinations. When we use static scheduling in ALL, OpenMP achieves only around 14× speedup (not shown), and CilkPlus can only reach to approximately 12× speedup. Using dynamic scheduling in OpenMP increases the speedup to only around 20×. In PAIR scheme and using dynamic scheduling SINA obtains around 26× speedup with OpenMP and around 25× speedup with CilkPlus. This demonstrates that, both of the programming models have competitive performances in PAIR. We observe that even though PAIR incurs additional synchronization overhead, it outperforms ALL, thanks to its finer granularity which yields better load balancing. We again note that in an attempt to reduce the number of locks in PAIR, SINA shuffles the task list before the execution. Since PAIR outperforms ALL, in the remaining experiments we use PAIR as the default scheme.

D. Data Layouts and Speedup

Next we compare three different data layouts obtained using three different vertex ordering on DBLP 14-17 dataset in Fig. 5. Figure 5(a) illustrates the execution time of SINA for this problem for three different data layouts and varying number of threads. In this figure, each grouped stacked-bar in the x-axis represents the number of threads and y-axis represents the execution time of the three parallel components of SINA. Sequential portion of SINA takes less than 10 seconds, hence it is omitted in this figure for simplicity. In the sequential run, more than 75% of the overall execution time comes from similarity computation. When we use 32 threads, this ratio decreases to 60%, because the parallelizations of the other two components do not scale as much as the similarity computation. As expected NAT order has the worst performance, because in this layout SINA portrays a more irregular data access pattern. Both BFS and HCB layouts improve the overall execution time with respect to NAT. BFS and HCB layouts give improvements in different parts of SINA. We observe that in the shortest path computations, BFS layout outperforms HCB layout by nearly 20% while HCB outperforms BFS in the similarity computation by nearly 25%. Since most of the overhead comes from the similarity computation, in the overall case, for 32 threads HCB outperforms BFS by nearly 10% and outperforms NAT by 40%. Figure 6 illustrates execution times and speedup on Hollywood 09-11 dataset. In Figure 6(a), we observe that similarity computation takes more than 85% of the total execution time. As listed in Table III even though the Hollywood graphs have nearly the same number of vertices as the DBLP graphs, the number of edges per vertex is much higher. This is the primary reason as to why similarity computation dominates the total execution time. Figure 6(b) illustrates the overall and the three components' speedup. We observe that in this dataset SINA achieves a speedup of 30× and gets closer to the linear speedup for similarity computation, and achieves 28× speedup for the overall execution. Even though the speedup of the shortest path and partitioning computations are far from the ideal at 14× and 20×, respectively, the overall speedup is not affected significantly, since these two components take a small portion of the total execution time. In Table IV we present speedups of the three parallel components and overall execution time for different datasets. For the similarity computation, SINA obtains 23× to 30× speedup, partitioning's speedup varies from 17× to 19×. One interesting fact that we observe is that for the shortest path computation while the speedup is between 22× to 23× in DBLP graphs, in Hollywood graphs this decreases to 15×. Again, since graphs are denser, frontier queues for BFS's enlarges faster and data access becomes more irregular.

V. RELATED WORK

Existing graph alignment methods are usually classified into four basic groups [START_REF] Conte | Thirty years of graph matching in pattern recognition[END_REF], [START_REF] Elmsallati | Global alignment of proteinprotein interaction networks: A survey[END_REF]: spectral methods [START_REF] Liao | Isorankn: spectral methods for global alignment of multiple protein networks[END_REF], [START_REF] Neyshabur | Netal: a new graph-based method for global alignment of protein-protein interaction networks[END_REF], [START_REF] Patro | Global network alignment using multiscale spectral signatures[END_REF], [START_REF] Singh | Pairwise global alignment of protein interaction networks by matching neighborhood topology[END_REF]; graph structure similarity methods [START_REF] Aladag | Spinal: scalable protein interaction network alignment[END_REF], [START_REF] Kuchaiev | Topological network alignment uncovers biological function and phylogeny[END_REF], [START_REF] Malod-Dognin | L-graal: Lagrangian graphlet-based network aligner[END_REF], [START_REF] Memišević | C-graal: C ommon-neighbors-based global gra ph al ignment of biological networks[END_REF], [START_REF] Milenkovic | Optimal network alignment with graphlet degree vectors[END_REF]; tree search or tabu search methods [START_REF] Chindelevitch | Optimizing a global alignment of protein interaction networks[END_REF], [START_REF] Kpodjedo | Using local similarity measures to efficiently address approximate graph matching[END_REF], [START_REF] Liu | A multiobjective memetic algorithm based on particle swarm optimization[END_REF], [START_REF] Saraph | Magna: maximizing accuracy in global network alignment[END_REF]; and integer linear programming (ILP) methods [START_REF] Bayati | Algorithms for large, sparse network alignment problems[END_REF], [START_REF] El-Kebir | Lagrangian relaxation applied to sparse global network alignment[END_REF], [START_REF] Klau | A new graph-based method for pairwise global network alignment[END_REF]. All of these works have scalability issues. Our algorithms leverage global graph structure and reduce the problem space. Furthermore, with effective parallelization of each compute intensive step of the algorithm, we alleviate most of the scalability issues.

We are aware of two parallel approaches for global alignment. The first approach, termed network similarity decomposition (NSD) [START_REF] Kollias | Network similarity decomposition (nsd): A fast and scalable approach to network alignment[END_REF], decomposes the ranking calculations of IsoRank's similarity matrix [START_REF] Singh | Pairwise global alignment of protein interaction networks by matching neighborhood topology[END_REF] using the singular value decomposition. While NSD shows an order of magnitude improvement over the state of the art sequential algorithms, the largest network tested has less than ten thousand vertices. The second approach is a shared memory parallel algorithm [START_REF] Neyshabur | Netal: a new graph-based method for global alignment of protein-protein interaction networks[END_REF] that is based on the belief propagation (BP) solution for integer program relaxation [START_REF] Bayati | Algorithms for large, sparse network alignment problems[END_REF]. It uses parallel matrix operations for BP iterations and also implements an approximate weighted bipartite matching algorithm. Neither of these approaches can handle additional information (vertex and edge attributes and labels) available in modern graphs.

VI. CONCLUSION

We presented SINA, a scalable iterative graph aligner. SINA is a careful, shared memory parallelization of a recent sequential graph aligner GSANA [START_REF] Yas | An iterative global structure-assisted labeled network aligner[END_REF]. SINA can achieve up to 28.5× speedup on a 32-core machine, reduces the total execution time of a graph alignment problem with 2M vertices and 100M edges from 4.5 hours to under 10 minutes, while retaining the state of the art alignment recall obtained by GSANA.

As a future work, we will investigate how to use both BFS and HCB layouts at the same time. BFS layout is highly beneficial for the shortest path computations, while HCB layout is beneficial for the similarity computations. In order to improve the overall speedup, keeping two copies of the graphs in two different layouts (BFS layout for the shortest path computations, and HCB layout for the similarity computations) could be practical and useful. Another alternative is to use hierarchical ordering that combines HCB and BFS. We will also look at how to make use of existing shortest path distances when computing the shortest paths to the new anchor vertices, in order to reduce the run time both theoretically and practically.

Fig. 1 .

 1 Fig. 1. Overview of the SINA, where boxes marked with "*" are parallelized.

Fig. 2 .

 2 Fig. 2. Task definition and parallel computation schemes.

Algorithm 4 :

 4 COMPSIMALL(B, B , P )for each v ∈ B1 do for each u ∈ B do compute σ(v, u) P [v].insert(u, σ(v, u))Only keeps top k

Algorithm 5 :

 5 COMPSIMPAIR(B, B , P ) for each v ∈ B do p ← ∅ Only keeps top k similar vertices for v for each u ∈ B do compute σ(v, u) p.insert(u, σ(v, u)) if p.top() > P [v].bottom() then lock(P [v]) P [v].insertAll(p) Only keeps top k unlock(P [v])

Fig. 3 .

 3 Fig. 3. Data layouts based on different vertex orderings.

Figure 5 (

 5 b) presents overall relative speedup results for three different layouts. In this figure, blue color represents speedups based on each layout's sequential execution time. Speedup results represented with orange color are computed based on the minimum sequential execution time which is with HCB. When we compare each layout with its sequential execution time, NAT layout gives 25×, BFS gives 24× and HCB layout gives 23× speedup. NAT and BFS gives 14× and 20× speedup when we compare them with the minimum sequential time. Since HCB outperforms other two layout techniques, in the following experiments, we use HCB as the default data layout.

Fig. 4 .

 4 Fig. 4. (a) Histogram of bucket sizes; (b) and (c) task workloads for ALL and PAIR, respectively; (d) overall speedup of different parallelization schemes using OpenMP and CilkPlus.

Fig. 5 .

 5 Fig. 5. Effect of data layouts to execution time and speedup in DBLP 14-17 dataset.

Fig. 6 .

 6 Fig. 6. Execution time and speedup for Hollywood 09-11 dataset using HCB data layout.

  The two similarity computation algorithms update those priority queues. Algorithm 3: PARALLELSIM(QT 1 , QT 2 , scheme) Input : QT1, QT2 are quadtrees of G1, G2. scheme is the comparison scheme Output : P [v]: top k similar vertices of vertex v ∈ V1 in V2. Compute the task list T . T ← ∅ for each non-empty B1 ∈ QT1 do B2 ← QT2.f indLeaf N ode(B1.center) if scheme = ALL then B ← B2 N ∈QT 2 .neigbors(B 2 ) N Execute the task list T in parallel. P [v] ← ∅, for ∀v ∈ V1 for each task T = B, B in parallel do if scheme = ALL then

	T ← T ∪ { B1, B }
	else
	for each B ∈ QT2.neigbors(B2) do
	T ← T ∪ { B1, B }
	COMPSIMALL(B, B , P )
	else
	COMPSIMPAIR(B, B , P )
	return P

TABLE II COMPONENTS

 II OF THE SIMILARITY FUNCTION σ.

TABLE III PROPERTIES

 III OF THE DATASETS. |N [x]| REPRESENTS AVERAGE VERTEX DEGREE, AND |µ| REPRESENT THE SIZE OF GROUND TRUTH MAPPING.

	Data Set	|V |	|E|	|N [x]|	max(|N [x]|)	|N [x]| < 3	|µ|	|S 1 | A V	A E
	DBLP-17-small DBLP-14-small	491,719 366,137	4,089,071 2,542,331	8.31 6.94	2,322 1,782	51,035 46,853	(10%) (13%)	294,531	64
	DBLP-17	1,966,877	9,059,634	4.61	2,322 616,386	(31%)
	DBLP-14	1,464,539	5,906,792	4.03	1,782 491,206	(34%) 1,440,379	64
	Hollywood 2011 1,917,070 114,181,101 Hollywood 2009 1,069,126 56,841,216	59.6 53.2	13,107 119,333 (6.2%) 11,468 72,540 (6.8%)	990,090	64
	64 128 192 256 320 384 448 512 Max. number of vertices in a bucket

TABLE IV SPEEDUP

 IV AND OVERALL EXECUTIONS TIMES FROM 1 TO 32 THREADS.

		Data Set	1	2	4	8	16	32
		Preparation	1.0	1.6	3.2	6.4	12.6 23.2
	DBLP	Partitioning	1.0	1.7	3.2	6.1	11.9 18.8
	14-17	Similarity	1.0	2.0	3.4	7.3	12.4 27.4
	(NAT)	Overall	1.0	1.9	3.4	7.1	12.4 25.8
		Exec. time (s) 6854	3614 2038 978	555	266
		Preparation	1.0	1.5	3.1	6.1	11.6 22.0
	DBLP	Partitioning	1.0	1.7	3.2	5.9	11.3 17.5
	14-17	Similarity	1.0	1.9	2.7	6.9	13.9 26.3
	(BFS)	Overall	1.0	1.9	2.7	6.7	13.2 24.3
		Exec. time (s) 4600	2553 1652 691	349	190
		Preparation	1.0	1.7	3.4	6.6	12.8 23.2
	DBLP	Partitioning	1.0	1.7	3.2	6.3	11.6 19.1
	14-17	Similarity	1.0	1.6	3.5	6.6	12.0 24.2
	(HCB)	Overall	1.0	1.6	3.4	6.6	12.1 23.1
		Exec. time (s) 3851	2354 1121 588	318	167
		Preparation	1.0	1.4	2.7	5.3	9.9	15.5
	Hollywood Partitioning	1.0	1.7	3.3	6.4	11.9 19.1
	09-11	Similarity	1.0	1.9	4.0	7.7	16.1 30.6
	(HCB)	Overall	1.0	1.8	3.9	7.5	15.3 28.5
		Exec. time (s) 16064 8777 4158 2143 1050 564