B Mohammadi
email: bijan.mohammadi@um2.fr

Principal angles between subspaces and reduced order modelling accuracy in optimization

Keywords: robust optimization, reduced order modelling, principal angles, uncertainty quantification Mathematics Subject Classification (2000) 49N30, 49M29, 65Y05, 65D30, 14N20, 15A18, 76B75, 76D55, 90C52, 93B51

The paper considers robust parametric optimization problems using multi-point formulations and addresses the issue of the approximation of the gradient of the functional by reduced order models. The question of interest is the impact of such approximations on the search subspace in the multi-point optimization problem. The mathematical concept used to evaluate these approximations is the principal angles between subspaces and practical ways to evaluate these are provided. An additional indicator is provided when a descent minimization algorithm is used. The approach appears also to be an interesting tool for uncertainty quantification of the design in the presence of models of increasing complexity. The application of these concepts is illustrated in the design of the shape of an aircraft robust over a range of transverse winds.

Introduction

Reduced order models are widely used in both simulation and design. These can be built, for instance, by assimilation of high fidelity data or experimental results by a parametric model. Once the model built it can be used in an optimization procedure providing cheap estimation of the functional and its sensitivity with respect to the model parameters. It can also be locally used in trust region approaches where the reduced order model needs to be dynamically enriched. Reduced order modelling receives different denominations following the field of research it is issued from: learning in neural networks, fitting with least squares in experimental plans or higher degrees polynomials, reduced-order modeling with proper orthogonal decomposition, etc. Exhaustive referencing is therefore out of scope and we indicate a few review texts such as [START_REF] Schilders | Model order reduction: Theory, research aspects and applica[END_REF][2][3].

The performance of a system designed for given functioning conditions often seriously degrades when these conditions are modified. Typical situations of interest involve a few (typically one or two) parameters describing the functioning of the system. Multi-point optimization is widely used to address robustness issues in engineering and we showed [START_REF] Gallard | Robust parametric design by multi-point optimization[END_REF][START_REF] Mohammadi | Reduced sampling and incomplete sensitivity for low-complexity robust parametric optimization[END_REF] how to use such formulation to address the robustness issue when a few functioning parameters are involved. In particular, we showed how to introduce in this formulation what we would like for the outcome of the design through a 'target-based' weighting in the functional. We also discussed optimal sampling of the functioning parameters intervals and quantitative confidence levels on the quality of our search direction through Gram-Schmidt orthonormalization of the multi-point gradients. This naturally brought into light the question of the sampling size which needs to be monitored to maintain the calculation complexity low. The aim was to avoid the worst-case theoretical limitation indicating that the sampling size should be larger by one the size of the control space [START_REF] Li | Robust Airfoil Optimization to Achieve Consistent Drag Reduction Over a Mach Range[END_REF][START_REF] Gallard | Robust parametric design by multi-point optimization[END_REF][START_REF] Desideri | Multiple-gradient descent algorithm for multiobjective optimization[END_REF]. We showed that this is only necessary if all the associated gradients (i.e. evaluated at the sampling points) are linearly independent which is never the case in optimizations involving a state equation. Hence, we showed, again through Gram-Schmidt procedure, that large dimensional parametric optimization problems can be treated with very small sampling of the functioning parameters range with marginal losses on the gradient information.

However, one often uses approximation in the definition of the gradients. And this is the starting point of this paper. One would like to analyze the impact of these using reduced order models on the search space defined by the gradients of the functional at the different sampling points of the functioning parameters. The mathematical tool we use for this analysis is the principal angles between subspaces. Reduced order modelling also includes situations where lower accuracy and modelling complexity are accounted for in the linearization step than for the direct simulation chain used for the definition of the functional. This is often the situation where an approximate simulation chain is considered for linearization, droping or approximating some of the ingredients in the direct chain. Of course, if the full gradient (i.e. based on the complete direct model) is accessible this discussion in pointless. But, calculation complexity and also use of black box simulation tools make that the linearization often takes place for a variant of the simulation loop.

After discussing algorithmic evaluations of the principal angles between the exact and approximate search spaces, we will compare the evolution of these through optimization iterations for a complex shape optimization problem for two scenarios of reduced order modelling. Even if nonintuitive, it will be shown why the modelling errors will have less impact on the design when using descent minimization algorithms. This approach appears to be a powerful tool to evaluate the pertinence of an increase in the modelling complexity and address related uncertainty quantification issues.

Robust parametric optimization

We are interested in a class of optimization problems where the cost function involves a parameter α not considered as a design parameter:

min x∈O ad j(x, α), α ∈ I ⊂ IR n , O ad ⊂ IR N , n N. ( 1
)
where x is the actual design vector belonging to O ad the optimization admissible domain. Usually, the functioning parameters α are just a few. Typical situations of interest are where n = 1 or n = 2.

In [START_REF] Mohammadi | Reduced sampling and incomplete sensitivity for low-complexity robust parametric optimization[END_REF] we showed how to use multi-point optimization to address such optimization problem. The aim is to remove during the optimization the dependency in α. This is done minimizing a functional J(x) encapsulating this dependency:

J = µ, such that σ ≤ σ 0 , where µ = 1 ω α l ∈Im ω l j(x, α l ), ω = m l=1 ω l , (2)
with I m a given sampling of I whose choice has been discussed in [START_REF] Mohammadi | Reduced sampling and incomplete sensitivity for low-complexity robust parametric optimization[END_REF] as well as the weights ω l which account for the kind of performance we eventually want for the design: constant performance over the functioning parameters ranges.

For the sake of simplicity and without loss of generality one can consider here both the sampling and the weighting to be uniform.

As in First-Order Second Moment (FOSM) methods [START_REF] Melchers | Structural Reliability Analysis and Prediction[END_REF], σ is to monitor the regularity of the final performance which should be as regular as possible with respect of α:

σ = 1 2 α l ∈Im < ∇ α j(x, α l ), ∇ α j(x, α l ) >, (3)
where <, > is the Euclidean scalar product. σ 0 is the level of variability for the initial design or of a reference configuration. As n is small, ∇ α j is evaluated by finite differences on I m . Such minimization problems have brought new interest to descent methods and this not only because of their lower computational complexity than the gradient free methods. Indeed, we will see that gradients are useful to see what should actually be the search space in a context of robust multi-point design. More precisely, the definition of the suitable search space is only possible after identification of the maximal free subspace in the space generated by the gradients of the functional at the sampling points I m . Indeed, as mentioned in the introduction, previous theory suggested that the size of the sampling should be larger than the dimension of the control space (i.e. m = N + 1) while we noted [START_REF] Gallard | Robust parametric design by multi-point optimization[END_REF][START_REF] Mohammadi | Reduced sampling and incomplete sensitivity for low-complexity robust parametric optimization[END_REF] that the dimension of the maximal free search subspace is usually much smaller than N .

A multi-point descent algorithm

We consider the following iterative descent algorithm for our constrained minimization problem involving a direct simulation chain linking the parameters (x, α) to the state U solution of a state equation F (U (q(x, α)) = 0 and to a functional j is:

                                       Given x 0 , 0 < ρ, 0 < η 1, I m , p max , 0 < ε 1, optimization iterations p = 1, ..., p max                              -m parallel state equation solutions F (U (q(x p ), α l )) = 0, α l ∈ I m , -m parallel evaluations of j(x p , α l ), α l ∈ I m , -m parallel solutions of the adjoint state V equation: V t F U (U (q(x p ), α l )) = j t U , α l ∈ I m , -m parallel evaluations of ∇ x j(x p , α l ) = j x + (V t F x ) t , α l ∈ I m , -define d the descent direction, d = ∇ x µ-< ∇ x µ, ∇ x σ > ∇ x σ + η∇ x σ, -control parameter variation, x p+1 = x p -ρd, Stop if d ≤ ε,
where a = a/ a is the normalized vector a and

∇ x µ = 1 ω α l ∈Im ω l ∇ x j(x, α l ), (4)
∇ x σ = α l ∈Im < ∇ α j(x, α l ), ∇ αx j(x, α l ) > . (5)
∇ xα j(x, α l ) is obtained from ∇ x j(x, α l ) by finite differences on I m , component by component, following what has been done for ∇ α j. The definition of the descent direction permits to make sure that both µ and σ decrease for small descent steps. Indeed, a first order development in x gives:

σ(x p+1 ) -σ(x p ) = ∇ x σ ∇ x σ.(x p+1 -x p ) = -ρη ∇ x σ ≤ 0,
and we have

µ(x p+1 ) -µ(x p ) = ∇ x µ ∇ x µ.(x p+1 -x p ) = -ρ ∇ x µ (1 -ζ 2 + ηζ),
where ζ =< ∇ x µ, ∇ x σ >. Therefore, µ is also decreasing as 1 -ζ 2 + ηζ ≥ 0 for |ζ| ≤ 1 and 0 < η 1 as chosen in the algorithm. Despite the natural presence of parallelism in this algorithm in the m independent evaluations of the state, functional and its gradient, computational complexity remains an issue. Possible solutions for the reduction of the sampling size have been presented in [START_REF] Mohammadi | Reduced sampling and incomplete sensitivity for low-complexity robust parametric optimization[END_REF] together with the use of incomplete sensitivity concept in the evaluation of the gradients. This latter permits to avoid the solution of the m adjoint equations in the algorithm. This is particularly suitable when using black-box state equation solvers not providing the adjoint of the state variables. We will mention these issues in section 7.

Beyond individual gradient accuracy (i.e. at each of the sampling point), what is important in these multi-point problems is the global search space defined by the ensemble of the gradient vectors. This means that one might tolerate higher error levels in each of the gradient defined at the different sampling point than for a single-point optimization situation as what is important is for the global search space to remain nearly unchanged. An interesting mathematical concept which permits to measure the deviation between two subspaces is the principal angles between subspaces.

Angles between subspaces

We use the mathematical concept of 'principal angles' between subspaces in the Euclidean space IR N ) initially introduced by C. Jordan [4]. If the maximum principle angle between the two subspaces is small, the two are nearly linearly dependent. Geometrically, this is the angle between two hyperplanes embedded in a higher dimensional space.

Let us briefly recall the concept of principal angles and how to practically compute them [5, 8]. For simplicity, suppose A and B are two subspaces of dimension k of IR N , N ≥ 2k, although this is not a prerequisite to define the principal angles. The k principal angles {θ i , i = 1, ..., k} are recursively defined as:

cos(θ i ) = < a i , b i > a i b i = max{ < a, b > a b : a ⊥ a m , b ⊥ b m ; m = 1, ..., i -1},
where a j ∈ A and b j ∈ B.

The principal angles θ i are between 0 and π/2. This is an important point and will be used later to take advantage of the positivity of the cosine of the angles. The procedure finds unit vectors a 1 ∈ A and b 1 ∈ B minimizing the angle θ 1 between them. It then takes the orthogonal complement of a 1 in A and b 1 in B and iterates. This procedure is not useful in practice as computationally inadequate. We would like to be able to find the angles θ i from the inner products < a i , b j > of the elements of two bases of A and B [9]. This would be interesting in our multi-point optimization context where we can exhibit an orthonormal basis of the global search space for the multi-point optimization problem using Gram-Schmidt orthonormalization. Now, let {a i , i = 1, ..., k} and {b i , i = 1, ..., k} be two arbitrary orthonormal bases for A and B.

Orthonormal bases are easy to obtain through the Gram-Schmidt orthonormalization procedure. Consider M being the matrix of the projection operator P r A of B onto A defined by:

P r A (b i ) = k j=1 < b i , a j > a j , M = (< b i , a j >) i,j .
The principal angles can be linked to this operator [9] through:

M = GΣH t ,
where G and H are orthogonal matrices and Σ = diag(cos(θ i )).

As G and H are orthogonal matrices, this is a Singular Vector Decomposition (SVD) of M . G and H are unknown at this point. But, we will show that we do not need them to get the θ i . Otherwise, the approach will be again computationally useless.

We recall that the columns of G are the left-singular vectors of M and eigenvectors of M M t and the columns of H are the right-singular vectors of M and eigenvectors of M t M . And most important that cos 2 (θ i ) are the eigenvalues of P r t A P r A which writes in matrix form as:

M t M = (GΣH t ) t (GΣH t ) = HΣ 2 H t with Σ 2 = diag(cos 2 (θ i )).
Therefore, to find the principal angles between subspaces A and B, knowing an orthonormal basis in each subspace, one should calculate M and find the eigenvalues of M t M and take the square root of them. This last operation is valid as the angles are between 0 and π/2, and their cosine therefore always positive.

We presented the approach for subspaces of the same dimension k, but it is not necessary for the two subspaces to be of the same size in order to find the angles between them. We need N ≥ 2k to be able to exhibit two orthogonal subspaces. If N < 2k, some principal angles necessarily vanish and for N = k they all vanish. This procedure is still valid if the subspaces have different dimensions. The projection operator can be defined as well as its transpose and the eigenvalues of M t M are real as this is a symmetric square matrix.

In our optimization applications we always proceed first with a reduction in size of the search space using a sampling reduction size algorithm [START_REF] Mohammadi | Reduced sampling and incomplete sensitivity for low-complexity robust parametric optimization[END_REF]. This makes the calculation of the whole set of eigenvalues feasible in terms of calculation complexity. However, if this is out of reach, one can evaluate the bounds on the angles to see the global pertinence of our reduced order models and gradient approximations. This can be done without an exact calculation of the all eigenvalues. It is sufficient to use the Gershgorin circle theorem to find these bounds as every eigenvalue of M t M lies within at least one of the Gershgorin discs D((M t M ) ii , R i ) centred on (M t M ) ii and with radius R i = j =i |(M t M ) ij |. And because M t M is symmetric, the eigenvalues being real, we only consider the intersection of the discs with the x-axis. Alternatively, the largest and smallest principal angles can be found using iterative power and inverse power methods applied to M t M .

One should however be aware that these bounds might not be sufficiently sharp to discriminate between two reduced order models and decide, for instance, which one is the more adequate for sensitivity analysis. Figure 1 shows a typical sketch. It represents principal angles calculated between a first subspace generated by the exact gradients of a transport model for a 10 points sampling of one of the functioning parameters of the model and two subspaces generated by the sensitivities derived from two approximations of this model for the same sampling. Details of the models can be found in [START_REF] Mohammadi | Value at Risk for confidence level quantifications in robust engineering optimization, optimal Control: Applications and Methods[END_REF]. But the modelling problem is not of a main concern here. What is important is that if one only considers the first and last principal angles, model M 2 is found being a better approximation to be used in a linearization procedure. However, with the whole spectrum in hand the picture is quite different and M 1 appears be more suitable if one intends to use this reduced order model for sensitivity analysis.

Search space dimension

The first angle is larger for model M1, but it represents a better approximation overall.

M1 M2 (deg)

Search space dimension

Fig. 1 Principal angles between the subspaces generated by an exact gradient calculation and the linearization of two reduced order models. This permits to quantify the pertinence of an approximation from the whole spectrum. Model M 1 is found to be a better approximation even with a first principal angle slightly larger than with M 2 .

Principal angles between multi-point search spaces are interesting to measure the pertinence of sensitivity definitions based on reduced order models. Indeed, the design will be unaffected by a reduction in the model's complexity if the search subspaces, generated by the gradients at the sampling points of the functioning parameter interval and their approximations, remain the same. In the sequel we discuss the application of these ideas in our multi-point robust optimization context.

Definition of the subspaces

Let us define subspaces A and B in the context of our multi-point optimization problems defined in section 2. In our problem the size N is very large and is given by the dimension of the control space. But, the subspaces A and B are of dimension k and we have 2k

N . For the sake of simplicity we consider that both subspaces are of the same dimension. These are free subspaces and we make sure by Gram-Schmidt orthonormalization that dim

(A) = dim(B) = k.
Let us denote by  an approximation of the functional j. The orthonormalization is applied to a set of exact G(α l ) and incomplete g(α l ) gradient vectors with respect to the control variable x evaluated at α l ∈ I m the sampling of size m:

G(α l ) = ∇ x j(x, α l ), g(α l ) = ∇ x (x, α l ).
The choice of m is such that k < m N . The issue of minimal sampling still providing an enough large search space has been discussed in [START_REF] Mohammadi | Reduced sampling and incomplete sensitivity for low-complexity robust parametric optimization[END_REF].

So we have A = Span(G(α l ), α l ∈ I m ) and B = Span(g(α l ), α l ∈ I m ). Now, the orthonormal bases {a i , i = 1, ..., k} and {b i , i = 1, ..., k} are defined by the Gram-Schmidt orthonormalization procedure applied to {G(α l ), α l ∈ I m } and {g(α l ), α l ∈ I m }. Hence, we will have

A = Span(a i , i = 1, ..., k) and B = Span(b i , i = 1, ..., k) but with k < m.
The size k depends on the accuracy requested in the Gram-Schmidt procedure we recall briefly (in practice modified Gram-Schmidt is applied for numerical stability). This algorithm stops when a given tolerance is achieved and k will be the rank at which the rest after successive projections is found smaller than this tolerance:

               Given I m , T OL, a 1 = G(α 1 ), iterations k = 1, ... a k = G(α k ) - k-1 j=1 <G(α k ),aj > <aj ,aj > a j , until a k < T OL. (6)
This gives us the basis {a i , i = 1, ..., k}. The same procedure is applied to g s to get the basis {b i , i = 1, ..., k}. Again, the theory is valid even if the two bases do not have the same number of vectors for the same accuracy T OL (i.e. A and B have different dimension).

Evolution of the principal angles during optimization

We would like to use the principal angles concept to evaluate gradients approximation by reduced order modelling. These can be evaluated versus exact gradients when available or against each other. This former permits to estimate, for instance, how much deviation is caused in the search space by different levels of approximation in an incremental modelling process.

To implement these we modify the algorithm in section (3) introducing a few extra steps in order to dynamically evaluate during the optimization the principal angles between the subspaces generated by the set of exact and incomplete gradients over the functioning parameters sampling I m and described in section 5. We think this procedure is necessary because unlike in the single-point optimization, the approximation here needs to be valid over a range of the functioning parameters and one would like also to make sure that the approximation holds over the optimization history for the different control parameter values visited.

The algorithm given in section 3 is modified as:

                                               Given x 0 , ρ, η, I m , p max , ε, T OL, optimization iterations p = 1, ..., p max                                        -m parallel state equation solutions F (U (q(x p ), α l )) = 0, α l ∈ I m , -m parallel evaluations of j(x p , α l ), α l ∈ I m , -m parallel evaluations of G l (α l ), α l ∈ I m , -m parallel evaluations of g l (α l ), α l ∈ I m , -define the maximal free subspaces A ⊂ Span(G(α l ), α l ∈ I m )
and B ⊂ Span(g(α l ), α l ∈ I m ), -evaluate the principal angles between A and B, -define d the descent direction as in 3, -control parameter variation,

x p+1 = x p -ρd, Stop if d ≤ ε,
As A and B are the 'approximate' largest free subspaces (in the sense given in section 5), they generate the respective gradient-based search spaces and we have dim(A) = dim(Span(G(α l ), α l ∈ I m )) and dim(B) = dim(Span(g(α l ), α l ∈ I m )) assuming that what falls beyond the threshold T OL is related to numerical artefacts and should not be considered as potential search directions. Indeed, we give in section 8 an example of how to numerically chose T OL. It is shown that its value should not be too small in order to avoid adding to the maximal free search subspace additional directions related to numerical noises and not containing any pertinent information. Avoiding these extra dimensions also improves the convergence of the optimization algorithm for the solution of our multi-point minimization problem. Indeed, extra gradients would make an admissible descent direction harder to find.

Sensitivity analysis

In the algorithms given in sections 3 and 6, at each iteration of minimization we need to provide ∇ x j(x p , α l ) for different values of α l . This is computer intensive even if it is fully parallel. Our aim is to see how to reduce the complexity of this optimization problem using reduced order modelling. Beyond the computational complexity issue, the necessity for such alternatives also comes from the fact that it is not always possible to proceed with the linearization of the direct simulation chain used for the definition of j. Indeed, consider a general simulation loop linking the independent parameters (α, x) to a functional j. As described in section 2, α is not a design parameter, but one would like the design to be robust over a range of α:

j : (α, x) → (α, q(x)) → U (α, q(x)) → j(α, x, q(x), U (α, q(x))), (7)
where q and U are dependent variables. Suppose q is cheap to compute and U expensive. For instance, q represents geometrical quantities and U state related variables related to the physics of the problem and solution of some costly governing equation. The gradient of j with respect to x is:

∇ x j(x, α) = j ,x + (j ,q + j ,U U ,q )q ,x , (8)
where j ,x , j ,q and j ,U are easy to access and are usually provided by the user as external modules in an industrial simulation platform. U q and q x are, on the other hand, difficult to access. The major part of the cost of this evaluation is due to U ,q and one usually uses an adjoint method to make this cost independent of the size of the control space as recalled in algorithm 3. However, today's industrial simulation platforms are more and more based on black-boxes and commercial packages not enabling the user for a direct access to the source of the codes. Differentiating the simulation codes, by automatic differentiation [START_REF] Griewank | Computational derivatives[END_REF] for instance, is therefore off the table. In the same way, it is quite unconceivable to develop in-house adjoint solvers when the cost function calculation relies on commercial packages. The only feasible gradient calculation approach with black-boxes is with finite difference approximation which has a cost proportional to the size of the control parameter space. Despite this limitation, this technique remains the most widely used.

Here, we consider it only as an alternative when an adjoint solution is not available and it is only used to validate sensitivity analysis based on reduced order modelling together with the ingredients presented above. Suppose we have q and Ũ two reduced order models for q and U expressed at iteration p of optimization as Φq(x p ) = q(x p ) and Ψ Ũ (α, q(x p )) = U (α, q(x p )). Φ and Ψ are scaling or transfer functions which will be frozen during linearization. The model involves control parameters values x p at iteration p of optimization. This indicates that the reduced order model is enriched during the iterations assimilating the new high-fidelity data by (7). Simple forms of the transfer functions are Φ = q/q and Ψ = U/ Ũ . We consider the approximate simulation loop  where only black-box and high complexity terms have been approximated:

 : (α, x p ) → (α, Φq(x p )) → Ψ Ũ (α, q(x p )) → j α, x p , Φq(x p ), Ψ Ũ (α, q(x p )) . (9)
We need both q and Ũ as q is now often part of a black-box package. Actually, considering a shape optimization situation for instance, the q-code which includes the CAD manipulation tools (CAD: Computer Aided Design) is even more difficult to access than the U -code. An example of Φ linking 3D surface meshes and the corresponding CAD-based geometry definitions is given in the problem of section 8. Φ -1 describes the path used to go from the CAD shape definition to the mesh and for a given mesh M, Φ(M) is defined in order to minimize the distance between M and Φ -1 ΦM meshes. An example of these relations is shown in figure 2. The notation Φ -1 is formal and this operator is not built as the inverse of Φ (same for the operator Ψ and Ψ -1 below). We consider the following linearization of (9) at iteration p of the optimization algorithms in 3 or 6:

∇ x (x p , α) = j ,x + (j ,q + j ,U (Ψ Ũ,q (Ψ -1 U )))Φq ,x , (10)
where Ũ,q (Ψ -1 U ) indicates the linearization of the reduced order model around the high-fidelity solution U restricted to the domain of definition of Ũ . Again, the notation Ψ -1 is formal and this transfer operator is not built as the inverse of Ψ (see the multi-grid context in section 7.1.2.

Examples of reduced order models

The ingredients of the paper will be tested on two examples of low complexity models in the context of shape optimization. The first approximate modelling concerns Hadamard incomplete sensitivity formulation and the second is based on grid coarsening for sensitivity analysis. Below we describe both briefly using the notations above and refer to mentioned references for details of each approach. These two approximations represent two major classes of complexity reduction based either on a model simplification on the continuous level or use of lower accuracy in the discretization for the sensitivity analysis.

Model reduction on the continuous level

One example of gradient approximation useful in our shape optimization problems is Hadamard incomplete sensitivity which addresses the following context [10, 15]:

the cost function j and control x have the same domain of definition D(x) (e.g. a shape and an aerodynamic coefficient defined over it), j is a product of functions such as j(x) = f (x, q(x)) g(U (q(x))).

If these requirements hold, we can use an incomplete evaluation of this gradient, neglecting the sensitivity with respect to the state, leading to the approximation ∇ x  = j ,x + j ,U q ,x = ∇ x f g. This is very interesting as ∇ x f can be analytically calculated in some situations. If not we will again need an approximate model q as described above leading to ∇ x  = j ,x + j ,U (Φq ,x ). An example of such implementation is given in the problem of section 8 for the linearization of the normal to the shape with respect to the shape parameterization necessary to get the gradient of the aerodynamic coefficients, such as the lift and drag, with respect to the shape. Finally, the approach is also interesting because the quantities involved are all locally defined on the domain of definition D(x) of x and do not involve the full domain of definition of the state variable U or even q.

One can go one step further from incomplete sensitivities introducing simple algebraic relations in Ũ [10]. Examples of such algebraic relations are, for instance, the Newton relation for the pressure distribution over the shape or wall functions describing the tangent velocity and temperature distribution along the direction normal to the wall.

Let us illustrate our purpose with the former relation. Designing a shape with minimum drag with respect to a parameter x (scalar for clarity) involves an integral on the shape of p(q(x), α)u ∞ (α) • n(x):

C d (x, α) = 1 2ρ ∞ u ∞ (α) 2 shape(x) p(q(x), α)(u ∞ (α).n(q(x))dγ, (11)
where superscript ∞ indicates inflow conditions. The Newton or the cosine-square law for the pressure is:

p(x, α) = p tot (α)(u ∞ (α) • n(x)) 2
where p tot (α) is the total pressure function of the inflow conditions α (e.g. inflow Mach number and fluid density). We therefore have p

(x, α)u ∞ (α).n(x) = p tot (α)(u ∞ (α)•n(x)) 3 . Its derivative with respect to x is (pu ∞ (α)•n) ,x = (pu ∞ (α))•n ,x +p ,x (u ∞ (α)• n) = 3p tot (α)u ∞ (α)(u ∞ (α) • n) 2 n ,x .
The first term in the sum is what we called above incomplete sensitivity and it worths (pu

∞ (α))n ,x = p tot (α)u ∞ (α)(u ∞ (α).n) 2 n ,x .
We see that if the pressure is defined by the cosine-square law, the exact and incomplete derivatives only differ by a factor of 3 and have the same sign. Now what happens when the pressure is through the solution of the Euler equations for instance. The expression above can be rewritten as pu ∞ (α

).n = p|u ∞ (α)| cos( u ∞ (α) |u ∞ (α)| .n).
The incomplete gradient is therefore p(u ∞ (α).n

) x = -p|u ∞ (α)| sin( u ∞ (α) |u ∞ (α)| .n) = 0 when n is aligned with u ∞ (α).
Otherwise, the incomplete sensitivity fails (e.g. in the area near the leading edge for an airfoil at no incidence). On the other hand, the model tells us that the pressure sensitivity with respect to shape variations vanishes if those are along the normal to the shape such that n ,x = 0. But this is compatible with the pressure boundary condition p ,n = 0 one uses with the Euler equations as in the problem of section 8. We therefore expect the incomplete sensitivity to be a good approximation of the gradient if the shape deformation parameterization is along the normal to the shape. Finally, if Ũ is the cosine-square law, Ψ is defined as a scaling term ratio of the pressure by the Euler equations over this expression.

Complexity reduction on the discrete level

Another widely used strategy for gradient approximation is using different space and time discretizations for the high-fidelity calculation with (7) and its approximation (9). This is particularly adapted to multigrid implementations [START_REF] Hackbush | Multi-Grid Methods and Applications[END_REF]6,7] where Ψ = I h H is the transfer (interpolation) operator from the coarse H to the fine level h. The two meshes can be separated by n intermediate levels and in this case

Ψ = I h h 1 U h1 h2 ...I hn-1 hn I hn H . Ψ -1 = R H
h denotes the restriction operator from the fine h to the coarse level H. And in the presence of n intermediate levels

Ψ -1 = R H hn R hn hn-1 ...R h2 h1 R h1 h .
Hence, principal angles permit to evaluate the impact on the multi-point search space of the different transfer operators between grids. These operators concern, for instance, reconstruction by high order interpolations or different smoothing/filtering strategies to remove the high frequency features of the solution. The angles will tell how the introduction of a new operator impacts the search space defined by the multi-point gradients in a context of robust design.

Evaluation of an incremental or adaptive modelling procedures

Beyond mesh refinement, the most interesting application of this approach is the evaluation of the pertinence of an incremental or adaptive modelling procedure through the evolution of the associated search spaces generated by the gradients at the sampling of the functioning parameter. Indeed, the interest of an increase in modelling complexity, and therefore cost, should be evaluated through its true impact on the search space which is the truly important entity in a design process.

Defining the descent direction

We need to define d the descent direction in both algorithms given in 3 and 6. Either with an exact or incomplete evaluation of the gradients, one remarks that one should not use all the available vectors G(α l ) or g(α l ) for α l ∈ I m , but only those 'clearly' linearly independent (at given tolerance T OL) and generating subspaces A or B which are of dimension k such as k < m. This is because what falls below this tolerance threshold does not bring pertinent information on the problem and is rather due to numerical artifacts. With the approximate models, the descent direction d given in section 3 is evaluated with  instead of j giving d:

d = ∇ x µ()-< ∇ x µ(), ∇ x σ() > ∇ x σ() + η∇ x σ(). (12)
To make sure we only consider the relevant part in the descent direction we proceed with a projection on subspace A or B when working respectively with the exact or incomplete gradients. Recalling that {a i , i = 1, ..., k} and {b i , i = 1, ..., k} are orthonormal bases of A and B, the descent directions then becomes:

P r A (d) = i=1,...,k < d, a i > a i and P r B ( d) = i=1,...,k < d, b i > b i .
In a context of optimization, these are eventually the most important quantities. Therefore, instead of the principal angles between A and B we might consider the reduced order modelling pertinent if the deviation between these descent directions remains small. This is the case if the following holds through the optimization iterations:

1-< P r A (d) P r A (d) , P r B ( d) P r B ( d) > 1. (13)
This estimator avoids the calculation of the eigenvalues of M t M described in section 4. We show in section 8 an example of the evolution of this indicator during optimization iterations for a shape design problem.

Full aircraft shape optimization

In this section, we show how the pertinence of gradient approximations can be evaluated using the principal angles between the subspaces generated by the gradients of different approximate models and the exact gradient.

The numerical example concerns a shape optimization problem for a full aircraft in transonic cruise condition. For such configurations, an efficient single-point optimization is already a difficult task as it involves several software elements. It also requires coherent geometrical manipulation of the shape parameterization and the surface and volume meshes.

Several sources of variability exist for these problems where large deviations appear between the prescribed shape and the shape during the flight, which we suppose stationary. These may be due to a change in the weight of the aircraft or due to variability in the flight conditions. One example of the latter is when the aircraft cruises against transverse winds which are very common. Usually aircraft are designed for a range of angle of incidence. These designs are usually realized with the sideslip angle set to zero.

We previously used this problem [START_REF] Mohammadi | Reduced sampling and incomplete sensitivity for low-complexity robust parametric optimization[END_REF] to illustrate robustness and optimal sampling issues in parametric optimization. And also how to provide quantitative confidence levels on the search direction compatible with a search space defined by parametric gradient evaluations. Here, we would like to go one step further and evaluate the impact of gradient approximations using principal angles between subspaces.

Let us briefly link our shape optimization procedure to the direct dependency chain (7) linking independent variables (α, x) to geometry and state dependent variables (q, U ) and to the cost function j and constraints C i , i = 1, ..., n c .

Shape parameterization and other independent variables

In our dependency chain x denotes a CAD-free parameterization [START_REF] Mohammadi | Shape Optimization in Fluid Mechanics[END_REF] which does not require a priori local regularity assumptions on the shape as it is implicitly the case in CAD-based shape definitions. More precisely, x represents shape deformations along the normal to the triangular faces of the surface mesh as shown in figure 2. This search space has a dimension N ∼ 5000. This parameterization receives different denominations and belongs to the same class as node-based or free-form shape definitions. In all these approaches the regularity of the deformation needs to be controlled [START_REF] Firl | Regularization of shape optimization problems using fe-based parametrization[END_REF]10]. This is because, unlike with a CAD definition, the shape ∂Ω of an object Ω and a gradient-based deformation do not belong to the same function space in terms of regularity and, actually, the second is always less regular. Let us illustrate this on a simple example with

J(x) = Ax -b 2 taking x ∈ H 1 (∂Ω), Ax and b in L 2 (∂Ω). The gradient J x = 2A T (Ax -b)
belongs to H -1 (∂Ω). Therefore, any variation along J x will have less regularity than x: δx = -ρJ x = -ρ(2(Ax -b)A) ∈ H -1 (∂Ω). We therefore need to project (or filter or smooth) into H 1 (∂Ω). But, why it is not the case with a CAD parameterization? Suppose the shape is described in a finite dimensional parameter space, as for instance with a polynomial definition of a surface (this is like a CAD parameterization). When we consider as control parameters the coefficients of the polynomial, changes in those do not change the regularity as the new shape will always belong to the same polynomial space. Sobolev inclusions give the key for the choice of the regularity operator with the CAD-free parameterization [10]. In our case, because we are using a piecewise linear discretization, a second-order elliptic system is sufficient. This operator will also be used below to enforces some of the geometrical constraints.

The parameter α is the sideslip angle inducing fully 3D effects on the flow around the plane making the consideration of a full aircraft necessary during the design. However, because the airplane geometry is symmetric spanwise, it is not necessary to consider a symmetric range for the transverse wind. We consider the sideslip angle α in I = [0, 10 o ]. In our multi-point analysis, we consider a uniform sampling I m of [0, 10 o ] with m = 30 points. The other flow conditions are given by a Mach number of 0.8 and an inflow incidence of zero degrees. Together with the sideslip angle these parameters fully describe a 3D inviscid flow around the aircraft.

Figure 2 shows the initial shape and a view of the initial triangular surface mesh. This correspondence is given by Φ -1 (described in section 7) while Φ permits to construct the corresponding optimized shape from this mesh and the shape deformation over it.

Geometric and state dependent variables

Let q(x) denotes the auxiliary unstructured mesh related geometrical quantities. When the shape is modified, this change must be propagated through the mesh keeping it admissible and we need to recalculate all related geometrical quantities. Admissible and positive mesh defirmation is achieved by a 3D torsional spring analogy method [START_REF] Farhat | A three-dimensional torsional spring analogy method for unstructured dynamic meshes[END_REF].

U (q(x), α) denotes flow variables depending also on the extra parameter α not part of those involved in the definition of the shape. More precisely, U (q(x), α) = (ρ, ρu, ρE) t represents the conservative flow variables solution of the Euler equations in conservation form. Our flow solver is based on a finite volume Galerkin method on unstructured tetrahedral meshes. The Roe [START_REF] Roe | Approximate Riemann solvers, parameters vectors and difference schemes[END_REF] flux is used for the approximation of the advection operator together with MUSCL reconstruction with Van Albada limiters [START_REF] Van Albada | Flux vector splitting and Runge Kutta methods for the Euler equations[END_REF] in presence of shocks. We target steady solutions and use time marching to reach these. The time integration procedure is explicit and is based on a low-storage Runge-Kutta scheme (RK4). The details of this implementation are available in [10]. Of course, these are not the only choices possible and the literature on numerical methods for compressible flows is huge.

Cost function and constraints

We consider a classical aerodynamic shape optimization which aims at minimizing the drag coefficient j = C d given by formula [START_REF] Griewank | Computational derivatives[END_REF]. This minimization is usually performed under equality constraints on the lift and volume of the aircraft. Let us denote these by C 1 = C l -C 0 l = 0 and C 2 = V -V 0 = 0 where superscript 0 indicates the corresponding values for the initial shape. The volume of an object Ω (here the aircraft) is expressed through the boundary integral formula:

V = Ω 1 = Ω 1 3 ∇.(X) = ∂Ω X.n,
where X = (x 1 , x 2 , x 3 ) t is the local coordinate. In addition to these, we consider some extra purely geometric constraints. Minimum surface regularity is prescribed for the wings leading edges. This is enforced through the CAD-free parameterization regularity operator. Other examples are the wings by-section thicknesses and the wings plan-form which are requested to remain unchanged. The latter does not introduce any difficulty but the by-section thickness is less obvious to enforce a posteriori. We define a by-section (figure 3) definition of the shape where the number of sections n s is free (here n s = 100) and depends on the complexity of the geometry. Each node in the parameterization is associated to a section Σ i , and for each section, we define the maximum thickness ∆ i . This last operation requires a projective analysis based on the projection of the upper-surface nodes over the lower surface for each section. This constaint can therefore be expressed as:

C 3 = ns i=1 |∆ i -∆ 0 i |.

Gradient of the functional and constraints

As we saw in algorithms given in sections 3 and 6, to built the multi-point descent direction, we need m evaluations of the state and functional and constraints sensitivities with respect to the shape. All the sensitivities are computed by automatic differentiation in reverse mode using tapenade [START_REF] Hascoet | Tapenade 2.1 user's guide[END_REF].

Of course, the constraints can be accounted for by introducing a penalty term in the cost function: J = j + i=1,3 a i C i , a i ∈ IR + . But this should be avoided when possible. Also, one classical technique to recover the lift during optimization is to change the flow incidence taking advantage of the linear relationship between the incidence and lift away from stall conditions. Suppose we do not want to use either penalty or such an approximation. An alternative would be to consider a locally admissible gradient orthogonal to Π = Span(∇ x C i , i = 1, 2, 3) with dim(Π) ≤ 3. Let us denote by π an orthonormal basis of this subspace obtained by the Gram-Schmidt procedure applied to the gradients of the constraints. The admissible gradient will be used in algorithms 3 and 6 and is given by: ∇

x j(x, α) ← ∇ x j-< ∇ x j, π > π,
where <, > indicates the scalar product over subspace Π. We need, therefore, to provide ∇ x C l receive a same treatment and we take advantage of the capability for multi-right-hand-side adjoint calculation of tapenade in reverse mode to acces these two gradients without necessitating the solution of two separate adjoint problems. Both also enter the validity domain of incomplete sensitivities as presented in section 8.5. Our direct Euler code uses time marching to the steady solution with local time steps. The corresponding reverse (or adjoint) code can be used easily to produce the incomplete sensitivity by setting to one the number of reverse time iterations for the steady direct state variables. This is interesting as it permits to implement this approximation in existing adjoint based optimization platform with no extra coding. Another optimization of the reverse mode code comes from the fact that, our situations of interest being stationary in time, there is no interest in storing the forward states for backward integration. [START_REF] Christianson | Reverse accumulation and implicit functions[END_REF]10].

∇ x C i , i = 1, 2, 3 which require ∇ x C l , ∇ x V and ∇ x ∆ i . ∇ x j = ∇ x C d and
Figure 4 shows the contours of a multi-point descent direction based on sensitivities evaluated on a fine and a coarse meshes using the state evaluated on the fine mesh and restricted on the coarse mesh. Below we compare the impact on the multi-point search space of gradient approximations using either the incomplete sensitivity concept or a two-level evaluation of the complete sensitivity.

Exact vs. incomplete sensitivities

Let us start with the evaluation of the impact on the multi-point search space for our shape optimization problem of a model reduction using the Hadamard incomplete sensitivity concept described in section 8.5. We think this approximation is valid because the cost function and the constraints are either purely geometric or enter the validity domain of Hadamard incomplete sensitivity as they can be written as a boundary integral over the shape for a product of geometry by state functions. One interest of incomplete sensitivity as presented in section is that there is no need for field variable linearization as everything is defined on the shape. This locality issue is very important in parallel computing where data locality is aimed as much as possible to reduce the communications.

To evaluate the pertinence of this gradient approximation in robust parametric optimization we follow the algorithm given in section 6. In particular, at each iteration of optimization, we identify the size of the free subspaces generated by the exact and incomplete sensitivities. This is achieved using algorithm (6) and is done for a given tolerance level T OL.

Figure 5 shows the dimension of the maximal free search subspace found by the incomplete Gram-Schmidt orthonormalization algorithm (6) for the original shape for different values of the tolerance level T OL. Small values of this parameter lead to large dimension search spaces. But these do not necessarily contain pertinent physical information and are partly related to numerical artifacts and noises. One considers T OL = 10 -4 which gives a dimension of the search subspace less than 10. This choice appears also a posteriori pertinent from figure 7 where one sees that the dimensions of the maximal free search subspaces decrease during optimization to 4. This analysis is enforced by figure 6 where one follows dim(Span(∇ x j(x, α k ), α k ∈ I 30 )) with T OL = 10 -4 , T OL = 10 -6 and T OL = 10 -7 during optimization iterations where the descent directions have been built in the maximal free search subspaces obtained for T OL = 10 -4 . We see that in all cases the dimension of the maximal free search subspace remains below 25 making 30 a safe choice for the sampling. Again, we consider that the dimensions above 10 are related to numerical artifacts and that these search directions should be avoided.

Figure 7 also shows the impact of gradient approximations on the dimension of the maximal free search subspaces during optimization. The evolution of the principal angles between the subspaces generated by the exact and the incomplete sensitivities shows a maximum angular deviation of around 30 degrees but with always some collinear directions between the subspaces. Most important the pertinence of the gradient approximation in the context of optimization is established following indicator [START_REF] Mohammadi | Value at Risk for confidence level quantifications in robust engineering optimization, optimal Control: Applications and Methods[END_REF] showing a deviation of about 15 degrees. This suggests, and this is surprising and non-intuitive, that descent methods will be less sensitive to modelling inacurracy. Convergence histories for the first and second momentum of the functional over the 10 degrees sideslip angle range demonstrate that the minimization problem is properly addressed and that the choice of the descent direction d by ( 12) is both pertinent and robust validating the analysis presented in section 3.

Sensitivity evaluation on coarser meshes

The next complexity reduction approach is, as described in section 7.1.2, when using different discretizations of meshes for the state and sensitivity evaluations. One could simply use a direct solver on a fine mesh h, restrict the state U h to a coarser mesh using operator

U H = Ψ -1 U h = R H
h U h , use the adjoint solver to access the gradient G H on this coarse mesh and interpolate back to the initial mesh with G h = Ψ G H = I h H G H . Practical meshing and geometry handling issues are not of importance here and any method, based either on structured or unstructured meshes, can be considered.

Figure 8 shows the histories of Gram-Schmidt orthonormalizations during optimization to find the size of the maximal free search subspaces with the exact gradient on the two meshes. The principal angles between these subspaces during optimization show a larger deviation than when using the incomplete sensitivity evaluation on the same fine mesh in figure 7. This somehow recovers what is expected from the incomplete sensitivity theory: it is better to perform an incomplete evaluation of the sensitivity on a fine mesh avoiding any state approximation than an exact evaluation of the sensitivity on a coarse mesh, even with a state restricted from a fine mesh [10]. This can be seen as an a posteriori validation of the dominant importance of the geometric part in these sensitivities over the state dependency. This analysis is also confirmed by the histories of the the mean and variance during an optimization based on the gradients on the coarse mesh, which appear less regular than in figure 7. Comparison of the evolution of indicator [START_REF] Mohammadi | Value at Risk for confidence level quantifications in robust engineering optimization, optimal Control: Applications and Methods[END_REF] in the two optimizations also confirms sometimes a larger deviation between the exact and approximate [START_REF] Mohammadi | Reduced sampling and incomplete sensitivity for low-complexity robust parametric optimization[END_REF] descent directions in the present optimization. The descent direction is however eventually better for the second approach.

This evolution can be explained a posteriori as in both optimizations, the principal angles between the subspaces based on the exact and approximate gradients decrease during optimization (see the middle pictures in figures 7 and 8). But, this is even more visible for the second case where we eventually also have more common directions between the two maximal free search subspaces shown by vanishing principal angles. This means that the second approximate models eventually better coincides with the exact model and therefore explains why the optimization performances have been recovered (lower-left picture in figure 8).

→ 

-1 φ  ← φ

Initial shape

Initial mesh

Shape deformation

Reconstructed optimal shape

Fig. 2 Initial shape and a view of the initial triangular surface mesh. This correspondence is given by Φ -1 while Φ permits to construct the corresponding optimized shape from this mesh and the shape deformation over it.

Concluding remarks

The impact of different gradient approximations by reduced order models has been evaluated in the context of robust shape design by multi-point optimization. In particular, the deviations between the multi-point search spaces based on the exact evaluation of the gradients and these approximations have been analyzed using the concept of principal angles between subspaces. It has been shown how this can be seen as a singular value decomposition problem which can be reduced in our situation to an eigenvalue evaluation problem for a symmetric matrix of small size. The reduction in size has been made through the extraction of the pertinent information in the multi-point vector search space. This approach can therefore be seen as a tool to validate reduced order models in the context of robust optimization. Indeed, what is important is to make sure that model reduction does not introduce any major perturbation in the free search subspace. Finally, Fig. 3 An example of the spanwise slices of an aircraft used to enforce the by-section maximum thickness constraint.

Multi-point descent direction based on the complete sensitivity -fine mesh

Multi-point descent direction based on the incomplete sensitivityfine mesh

Multi-point descent direction based on the complete sensitivity -coarse mesh & interpolated on the fine.

Fig. 4 Iso-contours of a multi-point descent direction based on the complete and incomplete sensitivities on the fine and coarse levels. The sensitivity on the coarse level has been evaluated using the state restricted from the fine to the coarse level with operator Ψ -1 = R H h and the sensitivity transfered back to the fine level using interpolation operator Ψ = I h H .

it has been shown that the situation can be more favorable when using a descent minimization algorithm as the reduced order model gradients need then to verify a less constraining criteria.

Beyond the examples presented, the approach provides a valuable tool to a posteriori measure the 'distance' between models. In particular, this could be used with adaptive models. Indeed, in practice the high-fidelity model is not available or too expensive and one would like to be able to identify the 'right' level of complexity for the model to be used during design. This concerns both the continuous model and its discretizations. In that sense, we showed examples of the impact of grid refinements on the deviation between the search spaces generated by the sensitivities for the different meshes. Hence, the approach permits to quantify the pertinence of an increase in the modelling complexity through the deviation between the search spaces. Free search subspace dimension Fig. 7 Upper: histories of Gram-Schmidt orthonormalizations during optimization to find the size of the maximal free search subspaces with the exact and incomplete gradients. Middle: principal angles between these subspaces during optimization (a curve per iteration). Lower-left: optimization histories of the mean and variance of the functional over the 10 degrees sideslip angle range. Lower-right: angles between the two descent directions given by indicator [START_REF] Mohammadi | Value at Risk for confidence level quantifications in robust engineering optimization, optimal Control: Applications and Methods[END_REF].

Mean / Mean0

Evolution of principal angles during 20 iterations of optimization between subspaces with the exact sensitivity on a coarse and fine mesh.

Variance / Variance0

Histories of Gram-Schmidt orthonormalizations (with the exact sensitivities on a coarse (o) and a fine (+) grid) during optimization giving a size of the free subspace always below 8 for TOL=1E-4.

Deviation between subspaces decreases during optimization

Fig. 8 Upper: histories of Gram-Schmidt orthonormalizations during optimization to find the size of the maximal free search subspaces with the exact gradient on these meshes. Middle: principal angles between these subspaces during optimization (a curve per iteration). Lower-left: optimization histories of the mean and variance of the functional over the 10 degrees sideslip angle range. Lower-right: angles between the two descent directions given by indicator [START_REF] Mohammadi | Value at Risk for confidence level quantifications in robust engineering optimization, optimal Control: Applications and Methods[END_REF]

Fig. 5 Fig. 6

 56 Fig.5Dimension of the maximal free search subspace versus the tolerance value T OL in the incomplete Gram-Schmidt orthonormalization algorithm (6) evaluated for the original shape. This evaluation helps choosing T OL in order to avoid non pertinent search directions due to numerical artifacts. The curves are for two different sortings of the gradients during orthonormalization.

  tions, Springer Math in Industry series Vol.13, Berlin, (2008). 2. Obinata, G. Andersonn B. Model reduction for control system design, Springer, Berlin (2000). 3. Qu, Z. Model Order Reduction Techniques with Applications in Finite Element Analysis, Springer, Berlin (2004). 4. Jordan, C. Essay on geometry in n dimensions,Bull. Soc. Math. France, 3, 103:174 (1875). 5. Gluck, H. Warner, F. Great circle fibrations of the three-sphere,Duke Math. J., 50,107-132 (1983). 6. Morano, E. Mavriplis, D. Venkatakrishnan, V. Coarsening strategies for unstructured multigrid techniques with application to anisotropic problems, SIAM J. Sci. Comput., 20-2, 393:415 (1998).

  . 7. Lallemand, M.-H. Steve, H. Dervieux, A. Unstructured multigridding by volume agglomeration: Current status, Computers and Fluids, 21, 397:433 (1992). 8. Jiang, S. Angles between Euclidean subspaces,Geometricae Dedicata, 63(2), 113-121 (1996). 9. Shonkwiler, C. Poincare duality angles for Riemannian manifolds with boundary, PhD. thesis, Univ. Pennsylvania (2009). 10. Mohammadi B., Pironneau O. Applied Shape Optimization for Fluids (2nd Edition), Oxford Univ. Press, London

Acknowledgements The author would like to thank F. Gallard, M. Meaux and M. Montagnac from Airbus, EADS-IW and Cerfacs for their feedbacks.


Michal Koperski
email: michal.koperski@inria.fr

Piotr Bilinski
email: piotr.bilinski@inria.fr

Francois Bremond
email: francois.bremond@inria.fr

3D TRAJECTORIES FOR ACTION RECOGNITION

Keywords:

Recent development in affordable depth sensors opens new possibilities in action recognition problem. Depth information improves skeleton detection, therefore many authors focused on analyzing pose for action recognition. But still skeleton detection is not robust and fail in more challenging scenarios, where sensor is placed outside of optimal working range and serious occlusions occur. In this paper we investigate state-of-the-art methods designed for RGB videos, which have proved their performance. Then we extend current state-of-the-art algorithms to benefit from depth information without need of skeleton detection. In this paper we propose two novel video descriptors. First combines motion and 3D information. Second improves performance on actions with low movement rate. We validate our approach on challenging MSR DailyActivty3D dataset.

INTRODUCTION

Human action recognition has been an active research topic for many years. It has also found applications such as: video surveillance, video data indexing, patient monitoring, humancomputer interface. Apart from many contributions in last years, action recognition is still a challenging task. The major problem is to find model which would be discriminative, but still flexible enough to handle intraclass variation, because same actions can be done in different way.

Recent studies have shown that methods based on local space-time approach (such as trajectories) and Bag-of-Words, reached high accuracy rate. These methods model motion by detecting points of interest on each frame and then compute trajectories by tracking them in time space. However this methods fail to recognize similar actions, as they ignore spatial relationship between features. This problem has been addressed by Bilinski et al. [START_REF] Bilinski | Relative dense tracklets for human action recognition[END_REF]. They use head as reference point and compute relative position of trajectories according to head position. Even though this method improves action recognition accuracy, it still has problems using 2D information to distinguish actions which are performed in depth plane.

Recently, thanks to rapid development in cost effective depth sensors it is feasible to capture real-time depth information. Compared to conventional RGB cameras, the depth camera has several advantages: depth images are insensitive to changes in lighting conditions, moreover depth information simplifies task of object segmentation, which eases human skeleton detection. Recently many authors [START_REF] Wu | Mining actionlet ensemble for action recognition with depth cameras[END_REF], [START_REF] Oreifej | Hon4d: Histogram of oriented 4d normals for activity recognition from depth sequences[END_REF] have focused on action recognition using pose detection. Unfortunately skeleton detection is not robust when occlusions occur and it can lead to inaccurate or even missing detentions. What is more, direct employing state-of-the-art local interest points detectors on depth map (instead of skeleton detector) is not feasible due to high noise of depth map comparing to RGB sequences [START_REF] Oreifej | Hon4d: Histogram of oriented 4d normals for activity recognition from depth sequences[END_REF]. Recently many dedicated depth map descriptors have been proposed [START_REF] Oreifej | Hon4d: Histogram of oriented 4d normals for activity recognition from depth sequences[END_REF], [START_REF] Hadfield | Hollywood 3d: Recognizing actions in 3d natural scenes[END_REF], [START_REF] Li | Action recognition based on a bag of 3d points[END_REF], but they are still limited in terms of high order search space and sensitivity to depth map noise or missing measurements.

To address the above limitations of methods based on RGB sequences and depth sequences we propose a novel method which can benefit from both data sources. As detection of interest points is easier on RGB sequences we propose to extend Wang et al. [START_REF] Wang | Action Recognition by Dense Trajectories[END_REF] "Dense trajectories" method by adding depth information to each trajectory point. Such approach improves detection accuracy of actions mainly performed in depth plane. To model spatial relationship we employ [START_REF] Bilinski | Relative dense tracklets for human action recognition[END_REF] "Relative trajectories". We propose to use head as center of dynamic coordinate system. To improve discriminative power of such descriptor the positions of head and trajectories is given in 3D (x, y, z) space.

Descriptors based on trajectory detection require a certain amount of movement in processed video, because all computed features depend directly on move. In case where given action needs only little movement or action is occluded, proposed descriptors fail. To overcome this issue we propose a novel descriptor which combines features computed on RGB sequences and depth map. We use SURF key point detector [START_REF] Bay | Speeded-up robust features (surf)[END_REF] on RGB sequences and on each detected SURF point we compute Local Depth Pattern based on depth map.

We evaluate our approach on challenging MSR DailyAc-tivity3D dataset. The experiments shows that our approach improves action recognition performance without need of skeleton detection.

The contribution of this paper are summarized as follows:

• We propose a new way to combine motion and depth information.

• We improve discriminative power of dense trajectories by adding depth information.

• We use relative trajectories to encode spatial information of features.

• We propose a novel descriptor for actions with low movement rate.

• Our approach does not need skeleton detection. Which is a very important feature as in real life video surveillance environment, reliable skeleton detection is not feasible due to camera viewpoint angle or high depth map noise.

OUR APPROACH

To extract trajectories information we employ (similarly to [START_REF] Wang | Action Recognition by Dense Trajectories[END_REF]) dense trajectories sampled on RGB sequences. For each frame we sample feature points with step W pixels. Such extracted points are tracked using optical flow and median filter kernel. We limit length of trajectory to L frames to avoid drifting problem. If trajectory exceeds length L we remove it from tracking process. To assure dense coverage, in situation when we detect that there are no tracked points in W ×W area, we sample a new point from this area and add it to tracking process.

Trajectory Shape Descriptor (TSD)

Local motion patterns can be encoded by the shape of trajectory. We describe shape of the trajectory of length L as sequence S = (∆P t , ∆P (t+1) , ..., ∆P t+L-1 ), where ∆P t is a displacement vector ∆P t = (x t+1 -x t , y t+1 -y t ). Vector S is normalized by sum of the magnitudes of the displacement vectors. Thus for each video v we obtain descriptor set:

Ω v = {S 1 , S 2 , ..., S Nv } (1)
where N v is number of trajectories detected in video v.

3D Trajectory Shape Descriptor (3DTSD)

To improve discriminative power of Trajectory Shape Descriptor we add depth information z to each trajectory point p = (x, y). The z value is computed as mean of area of size N around point p. Such approach allows to estimate depth value of missing measurement points and to filter out noise which is an issue in depth map. It is a very important step as displacement vector cannot be computed for points which do not have measured depth. In addition, trajectories are very often detected at the edges, where depth missing values are most likely to present. Before computing displacements, trajectories which contain at least one point with missing depth value are being removed. For remaining trajectories 3D Trajectory Shape Descriptor is computed in the analogical way to Trajectory Shape Descriptor. Thus for each video v we obtain descriptor set:

Ψ v = {S 1 , S 2 , ..., S Nv } (2)
where N v is number of trajectories which remained in the video v.

Relative Trajectory Descriptor (RTD)

Both Trajectory Shape Descriptor and 3D Trajectory Shape Descriptor encodes only displacement information ignoring important spatial position of trajectory. A common way to resolve this issue is to use either spatio-temporal grids or multiscale pyramids. But those methods provide only coarse information. To resolve this problem similarly to [START_REF] Bilinski | Relative dense tracklets for human action recognition[END_REF] we employ relative trajectories. We use detected head as a center of dynamic coordinate system. Such descriptor encodes both shape characteristics and spatial position, this approach helps to distinguish similar trajectories detected at different positions. Given trajectory t i = [(x j , y j ), ..., (x j+L , y j+L )] where L is the length of trajectory and head trajectory t h = [(x ′ j , y ′ j ), ..., (x ′ j+L , y ′ j+L )]. We can define relative trajectory as: R = [(x j , y j , ..., x j+L , y j+L ) -(x ′ j , y ′ j , ..., x ′ j+L , y ′ j+L )].

(3) Thus for each video v we obtain descriptor set:

Φ v = {R 1 , R 2 , ..., R Nv } (4)
where N v is a number of trajectories detected in video v. RTD descriptor may be combined with other non-relative descriptors. Such solution would allow action recognition even if head detection is missing.

3D Relative Trajectory Descriptor (3DRTD)

To improve discriminative power of Relative Trajectory Descriptor we add depth information to each trajectory point and each head trajectory point. We do it in same way as in 2.2. Thus for each video v we obtain descriptor set:

Π v = {RD 1 , RD 2 , ..., RD Nv } (5)
where N v is number of trajectories extracted from the video v.

Trajectory Appearance Descriptors

Descriptors proposed in sections: 2.1, 2.2, 2.3, 2.4 provide only trajectory shape information. To empower the motion information of dense trajectories, similarly to [START_REF] Wang | Action Recognition by Dense Trajectories[END_REF] we compute HOG, HOF, MBH (motion boundary histogram) in spatiotemporal volume around each trajectory point. The volume size is N × N pixels and L frames long.

SURF Key Points Appearance Descriptor

Descriptors proposed in previous sections rely on trajectory detection. This implies requirement that there must be certain amount of movement in the processed video to compute enough discriminating features. In case where given action needs only little movement or action is occluded proposed descriptors fail.

To overcome this issue we propose another descriptor which will catch appearance features even when there are no trajectories detected. To achieve this goal we propose to use SURF key point detector [START_REF] Bay | Speeded-up robust features (surf)[END_REF]. On each sampled frame SURF key points are detected inside bounding box of detected person. For each detected key point SURF descriptor is computed based on RGB appearance.

Local Depth Pattern Descriptor

To improve discriminating power similarly to [START_REF] Zhao | Combing rgb and depth map features for human activity recognition[END_REF] we compute Local Depth Pattern in the neighborhood of detected SURF points. At given frame f for each detected SURF point p we divide space around p into N x × N y cells. Each cell has size of S x × S y pixels. Then we compute average depth value for each cell. To create the feature vector we compute difference between each cell pair.

Action recognition

We use standard bag-of-words approach. First we construct codebook for each descriptor using k-means algorithm. We empirically set the size of codebook to 4000 words. Descriptors are then assigned to the closest word from codebook in means of Euclidean distance. For each video histogram of occurrences of codebook words is computed. Which is then normalized by L1 norm. For classification we use a non-linear SVM with χ 2 -kernel.

As mentioned before trajectory based descriptors fail when there is no enough motion, due to either characteristics of action or occlusions. To distinguish actions with high and low movement and benefit from both trajectories descriptor and key points descriptor we use hierarchical framework (see figure 1). First we compute descriptor which combines Trajectory Shape Descriptor (TSD) and 3D Trajectory Shape Descriptor (3DTSD) and train classifier which distinguish actions with high movement from actions with low movement rate. Then for actions with low movement rate we combine SURF Key Points Descriptor with Local Depth Pattern Descriptor and then train classifier using described bag-ofwords approach. For high movement rate actions we re-use descriptor from top hierarchy level and train classifier.

EXPERIMENTS

In this section we evaluate the performance of our approach on challenging MSRDailyActivity3D dataset. We compare the results to state-of-the-art methods.

MSRDailyActivity3D

In this section we evaluate our approach on MSRDailyAc-tivity3D [START_REF] Wu | Mining actionlet ensemble for action recognition with depth cameras[END_REF] dataset. The data set consists of 16 actions performed by 10 subjects. Each action is performed in standing and sitting position which brings additional intraclass variation. We use Leave One Subject Out setup, where in each split one person is selected for testing and training is performed on the remaining subjects. In this case we have 10 splits. For computing 2D dense trajectories, HOG, HOF and MBH features we use LEAR's implementation 1 . The following parameters values has been: W = 5, L = 15, N = 32, as they gave the best experimental results.

We have selected the actions such: write on paper, use laptop, sit still as low movement actions, due to low number of detected trajectories (see figure 3). Play game, play guitar are also considered as low movement actions as most trajectories are detected on head or due to whole body move which is not discriminative for those actions. In play game action subject plays on pad involving small fingers movement, the same case is with play guitar action. The remaining actions such as drink, eat, read book, call phone, use vacuum, cheer up, toss paper, lie down on sofa, walk, stand up, sit down are considered as high movement actions. The top level classier (see figure 1) achieved 90% of accuracy. In table 2 we provide comparison of different descriptors combination performance for high movement actions. In table 1 we provide performance comparison for low movement actions. The final results are in table 3. Note that methods printed in italic require skeleton detection which is not the case in our approach. It is a very important feature as in real life video surveillance environment reliable skeleton detection is not feasible due to camera viewpoint angle or high depth map noise (as people often stay outside of optimal sensor range).

Fig. 1 .

 1 Fig. 1. Classification framework

Fig. 2 .Fig. 3 .

 23 Fig. 2. Sample frames from MSR DailyActivity3D dataset

Table 1 .

 1 Low movement action classification -descriptors comparison for DailyActivity3D datset.

	Method	Accuracy
	Trajectory Shape Desc (TSD)	0.78
	3D Trajectory Shape Desc (3DTSD)	0.74
	TSD + 3DTSD	0.85
	TSD + Relative Trajectory Descriptor (RTD)	0.83
	3DTSD + 3D Relative Trajectory Descriptor (3DRTD)	0.83
	HOG	0.68
	HOF	0.80
	MBH	0.84

Table 2 .

 2 High movement actions classification -descriptors comparison for DailyActivity3D datset.

	Method	Accuracy
	Dynamic Temporal Warping [9]	0.54
	HON4D [3]	0.80
	Actionlet Ensemble [2]	0.85
	TSD (without hierarchical classifier)	0.58
	3DTSD (without hierarchical classifier)	0.55
	TSD + RTD (without hierarchical classifier)	0.65
	TSD + 3DTSD (without hierarchical classifier)	0.63
	Our Approach (with hierarchical classifier)	0.72

Table 3 .

 3 Recognition Accuracy Comparison for DailyActiv-ity3D datset. Methods in italic require full skeleton detection.

http://lear.inrialpes.fr/people/wang/dense_ trajectories (Second version)

Acknowledgments

This work is supported by Dem@care and SafEE projects. Method Accuracy SURF 0.57 Local Depth Pattern (LDP) 0.56 SURF + LDP 0.60