A Rakotomamonjy
email: alain.rakotomamonjy@insa-rouen.fr

S Chanda

p -norm Multiple Kernel Learning with Low-Rank Kernels

Kernel-based learning algorithms are well-known to poorly scale to large-scale applications. For such large tasks, a common solution is to use low-rank kernel approximation. Several algorithms and theoretical analyses have already been proposed in the literature, for low-rank Support Vector Machine or low-rank Kernel Ridge Regression but not for multiple kernel learning. This work bridges this gap by addressing the problem of scaling p -norm multiple kernel for large learning tasks using low-rank kernel approximations. Our contributions stand on proposing a novel optimization problem, which takes advantage of the low-rank kernel approximations and on introducing a proximal gradient algorithm for solving that optimization problem. We also provide partial theoretical results on the impact of the low-rank approximations over the kernel combination weights. Experimental evidences show that the proposed approach scales better than the SMO-MKL algorithm for tasks involving about several hundred thousands of examples. Experimental comparisons with interior point methods also prove the efficiency of the algorithm we propose.

Introduction

Kernel methods such as Support Vector Machines or Kernel Rigde Regression have now become classical tools for classification or regression problems. In these methods, the choice of the kernel is of primary importance for achieving good performances. Multiple kernel learning (MKL) is a learning framework that transfers the choice of the kernel from the practitioner to the algorithm since the latter provides as outputs a learned linear combination of some base kernels and a decision function [START_REF] Lanckriet | Learning the kernel matrix with semidefinite programming[END_REF][START_REF] Bach | Multiple kernel learning, conic duality, and the SMO algorithm[END_REF].

Since the last few years, the literature on MKL have been flourishing and many efforts have been spent on proposing and analyzing novel regularizers for kernel combinations [START_REF] Szafranski | Composite kernel learning[END_REF][START_REF] Kloft | lp-norm multiple kernel learning[END_REF][START_REF] Hinrichs | Q-mkl: Matrix-induced regularization in multi-kernel learning with applications to neuroimaging[END_REF]. Improving the learning algorithm computational efficiency [START_REF] Rakotomamonjy | SimpleMKL[END_REF][START_REF] Sonnenburg | Large scale multiple kernel learning[END_REF][START_REF] Vishwanathan | Multiple kernel learning and the smo algorithm[END_REF][START_REF] Aflalo | Variable sparsity kernel learning[END_REF][START_REF] Suzuki | SpicyMKL : A fast algorithm for multiple kernel learning with thousands of kernels[END_REF][START_REF] Tian | A multiple kernel framework for inductive semi-supervised svm learning[END_REF] has also been on the focus of many interests. While these works have pushed the boundaries of multiple kernel learning method scalability, we believe that some efforts are still needed for making MKL capable of better handling large-scale learning problems.

One common way for handling large-scale problems in kernel methods is to consider low-rank approximation of the kernel matrix and to take advantage of this approximation for accelerating the learning process [START_REF] Williams | Using the nystrom method to speed up kernel machines[END_REF][START_REF] Drineas | On the nyström method for approximating a gram matrix for improved kernelbased learning[END_REF][START_REF] Kumar | Sampling methods for the nyström method[END_REF]. As an example, for Kernel ridge regression, the Woodbury matrix inversion formula [START_REF] Golub | Matrix computations[END_REF] makes it easy to derive that efficient algorithm. For SVM, Fine et al. [START_REF] Fine | Efficient svm training using low-rank kernel representations[END_REF] have proposed an interior point optimization method that can leverage on the low-rank kernel approximation. A recent work has also considered kernel approximation for multiple fisher discriminant analysis [START_REF] Wang | A novel multiple nystrm-approximating kernel discriminant analysis[END_REF]. However, in this latter work, the kernel weights are not learned according to the data and the problem at hand, but defined according to the quality of the kernel approximation. While kernel matrix low-rank approximations are often computed without any supervision on the labels, some works also proposed to improve the kernel approximation by taking into account distance or similarity constraints over the training examples [START_REF] Kulis | Learning low-rank kernel matrices[END_REF] or even by considering their labels [START_REF] Bach | Predictive low-rank decomposition for kernel methods[END_REF]. In this work, we propose to make a another leap towards the goal of very largescale multiple kernel learning. For this purpose, we develop an algorithm which is able to benefit from low-rank approximation of the base kernels used in a multiple kernel learning framework instead of directly learning a single low-rank kernel as in [START_REF] Kulis | Learning low-rank kernel matrices[END_REF][START_REF] Bach | Predictive low-rank decomposition for kernel methods[END_REF]. In particular, we focus on p -norm MKL because of its theoretical soundness [START_REF] Kloft | On the convergence rate of p-norm multiple kernel learning[END_REF] and because it has been proven to be useful in some practical applications [START_REF] Kloft | lp-norm multiple kernel learning[END_REF]. Within this context of p -norm MKL, we make the following contributions :

• we introduce a novel framework for MKL with low-rank kernels and we adapt, in a non-trivial way, a recently proposed technique for convex optimization to this problem of multiple kernel learning with low-rank kernels.

By exhibiting a trick on projection on affine sets, we are able to drastically reduce the computational complexity of our algorithm. In addition, our approach is parallelizable, a property that comes from the nature of the algorithm itself as it is essentially based on matrix-vector multiplications.

• we provide a theoretical result that bounds the difference between the minimizer of the MKL problem with the exact kernel and and the one with low-rank kernel approximation, of each component of the MKL decision function. This bound gives us intuitions on how the kernel approximations impact the weights in the kernel linear combination.

• We empirically demonstrate that our algorithm for solving the low-rank optimization problem is more efficient than interior point methods. In addition, we show that our MKL with low-rank kernel approach can be significantly faster than state-of-the-art p MKL solvers such as the SMO-MKL algorithm [START_REF] Vishwanathan | Multiple kernel learning and the smo algorithm[END_REF] while providing similar performance accuracies.

Framework

We introduce in this section the multiple kernel learning framework we are interested in. Let us consider a classification learning problem from data {x i , y i } i=1 where x i belongs to some input space X and y i = {+1, -1} denoting the class label of examples {x i }. We denote Y as the diagonal matrix composed from the label vector.

In a multiple kernel learning framework which involves a linear combination of m kernels, one looks for a decision update {d k } according to Equation 35: until convergence is met positive definite kernel K k . One way to learn these functions f k (•) is to solve the following MKL primal problem as first proposed in [START_REF] Zien | Multiclass Multiple Kernel Learning[END_REF] and [START_REF] Rakotomamonjy | SimpleMKL[END_REF] and extended by Kloft et al. [START_REF] Kloft | lp-norm multiple kernel learning[END_REF] :

function of the form f (x) = m k=1 f k (x) + b with each function f k ∈ H K k , H K k being the RKHS associated with
min d k ≥0,f k ,b,ξi≥0 m k=1 f k 2 H K k d k + C i=1 ξ i s.t y i k f k (x i ) + y i b ≥ 1 -ξ i ∀i d 2 p ≤ 1 (1)
where each d k regularizes the squared norm of f k in the objective function. Hence, a smaller d k indicates a smoother f k . The p-norm constraint on the weight vector d controls the amplitudes of {d k } while the choice of p eventually enforces the kernel linear combination to be sparse (p = 1), non-sparse (p > 1) or have equal weights (p = ∞).

Problem ( 1) is usually solved by considering an alternating optimization strategy which consists in optimizing a slave problem with respect to {f k } k with fixed weights {d k } and then in optimizing the weights while the functions f k are fixed [START_REF] Kloft | lp-norm multiple kernel learning[END_REF]. This alternating strategy, depicted in Algorithm 1, has been proved to be efficient owing to several features. The first one is that the slave problem boils down to be an SVM problem where the kernel matrix is a fixed combination of kernels min α

1 2 α Y k d k K k Yα -1 α st y α = 0 0 ≤ α i ≤ C (2)
and thus off-the-shelf efficient SVM solvers can be used. The other point is that the optimization with respect to a weight d k admits a closed-form solution :

d k = f k 2 p+1 H K k ( m k=1 f k 2p p+1 H K k ) 1 p ∀k (3)
where

f k (•) = d k i=1 α i y i K k (x i , •).
Although this strategy can be made further efficient, by for instance interleaving the SVM and the weight optimization solvers [START_REF] Kloft | lp-norm multiple kernel learning[END_REF], it still suffers because of the size of kernel matrix, when large-scale learning problems are in play. The SMO-MKL algorithm [START_REF] Vishwanathan | Multiple kernel learning and the smo algorithm[END_REF] departs from this alternating optimization strategy and instead solves the dual of a modified version of Problem [START_REF] Aflalo | Variable sparsity kernel learning[END_REF]. A very nice feature of this latter algorithm is that it can handle large-scale problems since kernel matrix entries can be computed on the fly.

One way to deal with such large-scale training datasets in kernel-based learning algorithms is to resort to low-rank kernel approximation of the kernel matrix. Several methods are available for performing low-rank approximation of a semi-definite positive matrix including Nystrom methods [START_REF] Williams | Using the nystrom method to speed up kernel machines[END_REF][START_REF] Kumar | Sampling methods for the nyström method[END_REF] or Incomplete Choleski decomposition [START_REF] Golub | Matrix computations[END_REF]. In this study, we will consider Nystrom method that we briefly remind in what follows. Suppose we want to approximate a positive definite Gram matrix K. A rank-R Nystrom approximation K of K is obtained by randomly sampling N n training examples among all available ones and by defining K = CW † R C = VV , where C is the matrix formed by the columns of the matrix K related to the N n samples, W is the Gram matrix of those examples, W † R being the pseudo-inverse of the best rank-R approximation of W with R ≤ N n and V a matrix of size × R so that V = CLD -1 2 with D the diagonal matrix of size R × R composed of the largest R eigenvalues of W and L the matrix of the corresponding eigenvectors. Computational complexity of the Nystrom method is dominated by the pseudo-inverse computation which is O(N 3 n ). The impact of the kernel approximation on the kernel method performance is of course of interest, since using lowrank approximation boils down in trading some computational efforts against eventually some losses of generalization performances. For some families of kernel methods, these impacts have been formally analyzed [START_REF] Cortes | On the impact of kernel approximation on learning accuracy[END_REF][START_REF] Bach | Sharp analysis of low-rank kernel matrix approximations[END_REF]. In this paper, we have addressed algorithmic issues and have studied the impact of the kernel approximation on the minimizer of problem (2) on which the weight d k depends. The theoretical study of the low-rank approximation impact on the performances is left for future works. However, our experimental results show that when the rank of the approximation becomes sufficiently large, the loss in generalization performances can become small compared to the use of the full kernel matrix.

Our contribution stands on developing an algorithm for large-scale multiple kernel learning by taking advantage of the low-rank decomposition of the base kernels. The algorithm we propose in the sequel follows the classical trend in MKL which consists in optimizing alternatively the inner SVM problem with fixed kernel weights and then in optimizing these weights while keeping the SVM optimization variables fixed. Hence, our main objective is to solve efficiently problem [START_REF] Bach | Sharp analysis of low-rank kernel matrix approximations[END_REF] given that we are provided with some Nystrom low-rank approximations {V k } m k=1 or all the {C k } and {W k } of all base kernels {K k } m k=1 . For the sake of clarity and simplicity, we suppose that all {V k } are of the same dimensions.

Algorithms

There exists two ways for solving the SVM problem (2) given low-rank kernel approximations. Indeed, the sum of kernel can either be approximated as a sum of low-rank kernels or as a low-rank approximation of k d k K k , that needs to be computed at each MKL iteration. This section presents the most efficient one when used within a MKL framework, which is based on approximating each kernel once and then in considering the sum of these low-rank kernels.

The approach is based on the idea that low-rank approximation of each single kernel can be computed beforehand and the sum of kernel K d is then approximated as the weighted sum of low-rank approximation of each single kernel.

Formally, we consider the following approximation :

K d = m k=1 d k K k ≈ m k=1 d k C k W † k C k = m k=1 d k V k V k
From this kernel approximation, Problem (2) boils down to be

min α 1 2 α Y k d k V k V k Yα -1 α st y α = 0 0 ≤ α i ≤ C (4) and each f k 2 H K k
can be obtained as :

f k 2 H K k = d k α YV k V k Yα (5)
Let us highlight that Problem (4), because of the sum of kernels, does not fit into the framework introduced by Fine et al. [START_REF] Fine | Efficient svm training using low-rank kernel representations[END_REF] for solving a low-rank kernel SVM. Indeed, because of the linear combination of kernels, the low-rank property may be lost as the rank of

k d k V k V k is upper bounded by k rank(V k V k ).
Hence, the gain achieved through the rank-R approximation of each V k may be spoiled ( for instance, one can easily construct a rank-1 kernel

approximation of the each V k V k which sum k V k V k is of rank m.
) Note that while SMO-MKL can be fed with the set of {V k V k } kernels, it is not able to take computational advantage of these low-rank approximations since each approximate kernel will still be considered as a full Gram matrix. This can be easily experimentally verified. Our objective is thus to derive an algorithm that is able to leverage from the low-rank approximations {V k } of all base kernel matrices.

For this purpose, we propose an equivalent formulation of Problem (4) by introducing some auxiliary variables and additional equality constraints :

min α,{γ k } m k=1 1 2 m k=1 γ k γ k -1 α γ k = √ d k V k Yα ∀k ∈ 1, • • • , m y α = 0 0 ≤ α i ≤ C ∀i ∈ 1, • • • , (6)
This optimization problem is a quadratic programming problem with R×m+ variables, R×m+1 equality constraints and 2 inequality constraints, thus when the number of training examples is large, we have to face with a pretty largescale problem. We can resort to classical QP solvers that handle such large-scale problems, for instance interior point methods [START_REF] Potra | Interior-point methods[END_REF]. However, in this work, we depart from this classical route by using a proximal gradient algorithm for solving problem [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF]. In the next paragraphs, we first describe the algorithm before giving its motivations. Indeed, as they are related to computational complexities, they will be clearer after the algorithm exposition.

The proximal method that we are considering here is the generalized forward-backward splitting algorithm [START_REF] Raguet | Generalized forward-backward splitting[END_REF] which is tailored for minimizing the sum of an arbitrarily large number of convex functions as in

min z F 0 (z) + M i=1 F i (z)
where one of these functions (here F 0 ) is smooth with gradient Lipschitz of constant L while the others here the {F i } M i=1 are possibly non-smooth but whose proximal operators can be simply computed [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF].

The Generalized Forward-Backward (GFB) splitting algorithm [START_REF] Raguet | Generalized forward-backward splitting[END_REF] follows the same principle as the Forward-Backward splitting [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF] in the sense that, at a given iteration, with an iterate z n , it first takes an explicit gradient step and then makes a implicit step where proximal operators are computed. The main feature of the Generalized Forward-Backward (GFB) splitting algorithm is that, since we have several proximal operators related to the constraints or regularizers to compute, these computations are applied in parallel on some auxiliary variables {Z i }. These variables are finally averaged so as to yield the next iterate z n+1 . The following (simplified) proposition makes the GFB algorithm explicit :

Proposition 1. [22, Th 2.1] Let {F i } M i=0 be M + 1 convex lower semi-continuous functions of R d such
that a standard domain qualification holds and such that the set of minimizers of

M i=0 F i (z) is not empty. Set {Z i } M i=1 ∈ R d , z 0 = 1 M M i=1 Z i and build at each iteration n ≥ 0 Z i ← Z i + λ n prox M ζnFi (2z n -Z i -ζ n ∇F 0 (z n )) -z n for all i ∈ 1, • • • , M , and
z n+1 ← 1 M M i=1 Z i with ζ n ∈ (0, 2 L
) and λ n ∈ (0, 1], then every sequence {z n } generated by this algorithm weakly converges towards a minimizer of

F 0 + M i=1 F i .
Note that this algorithm is exactly the forward-backward splitting algorithm when M = 1 and in the same way, it is robust to noise on the proximal operators and on the gradient computation under weak conditions on the noise.

By defining the vector z = [γ 1 ; γ 2 ; • • • ; γ m ; α], with z ∈ R R×m+ we can see that this GFB algorithm can be straightforwardly applied to our problem by choosing F 0 (z) as the objective function [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF]. The constraints in that problem can form our functions F

1 2 k γ k γ k -1 α of problem
prox I A •=0 (v) = v -A (AA ) -1 Av
Computing the term A (AA ) -1 A involved in the proximal operator, at a cost of the order of O(m 3 R 3 + m 2 R 2 ), is the main burden of the algorithm, although it is computed only once. Applying the proximal operator to a vector v involves about O(m 2 R 2 ) multiplications. Hence, because of the matrix inversion, the number of equality constraints is critical as it introduces a cubic dependency with respect to the number of kernels. In order to break down this complexity, we take advantage of one property of the GFB algorithm. Indeed, since this algorithm can handle any number of constraints and by noticing that the following simple but essential equality holds

I Az=0 = j=1 I A C j ,: z=0
where the {C j } forms a partition of the row index of A, one can consider the constraints in problem [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF] as the intersection of the m + 2 sets defined by : {γ k = √ d k V k Yα} m k=1 , α y = 0 and 0 α C1. Basically, we have splitted the matrix A in submatrices A k which applies only on some components of z. Under this novel perspective, Problem (6) still fits into the framework of the generalized forward-backward splitting [START_REF] Raguet | Generalized forward-backward splitting[END_REF]. However, we now have m + 1 affine sets and a box constraint on which we have to project the vector z. The m affine sets related to the linear

constraints {γ k = √ d k V k Yα} m k=1
, are formally defined as

I A k z=0 with, for k ∈ 1, • • • , m A k = 0 • • • I R • • • 0 - √ d k V k Y (7)
with I R being an identity matrix of size R × R. The constraint related to α y = 0 is defined as

A m+1 = 0 • • • 0 • • • 0 y since the optimization variables are [γ 1 ; γ 2 ; • • • ; γ m ; α].
Because of its peculiar structure, prox 6) if its objective function is uniformly convex. This can be easily achieved through the addition of a term η 2 α α with η > 0 at the expense of slightly perturbing the original problem.

I A k •=0 (
Remark 2. In the above algorithm, we have decided to split the original linear equality constraints in m affine sets as they are naturally related to a given low-rank kernel approximation. However, we could have chosen different number of sets. For instance, the extreme case is the case where a proximal operator is associated to a single linear equality constraint making the computation of the related

(A k A k ) -1 in O(1 + ).
The number of proximal operators involved in the for-loop of the algorithm is now R • m but each application of the proximal operator only costs O(1 + ). The drawback of this approach is that the memory complexity drastically increases since we have to store all the copies Algorithm 2 Efficient adaptation of the Generalized Forward-Backward for SVM with sum of low-rank kernels for i = 1 → m do 12:

1: Input : {V k } k matrices so that K k ≈ V k V k . 2: set ζ n = 1, λ = 1 3: initialize auxiliary variables {Z i } m+2 i=1 to random values where each Z i is a vector of size R × m + 4: set z 0 = 1 M +2 M +2 i=1 Z i . 5: set [γ 0 ; α 0 ] = z 0 6: precompute A k and A k (A k A k ) -1 for k = 1, • • • , m + 1 7: n ← 0 8: repeat 9: ∇F 0 = [γ n ; -1]
z n = 2z n -Z i -ζ n ∇F 0 13: Z i ← Z i + λ(prox I A i •=0 (z n ) -z n ) 14:
end for 15:

% Compute the forward-backward step related to the constraint α y = 0 16:

z n = 2z n -Z m+1 -ζ n ∇F 0 17: Z m+1 ← Z m+1 + λ(prox I A m+1 •=0 (z ) -z) 18:
% Compute the forward-backward step related to the box constraint 19:

z n = 2z n -Z m+2 -ζ n ∇F 0 20: Z m+2 ← P b L ,b U (z n ) 21:
% Define the next iterate of z as the average of the auxiliary variables 22:

z n+1 = 1 m+2 m+2 i=1 Z i 23: [γ n+1 ; α n+1 ] = z n+1 24:
n ← n + 1 25: until convergence is met {Z i }, which would be in (R•m+ )×R•m. In addition, while the computational complexity of each proximal operator is small, we have no guarantee regarding the number of iterations needed by the global algorithm to converge towards a solution. Hence a trade-off has be made. We have not thoroughly analyzed this issue in this work and have left it for future works. This possibility of splitting the linear constraints and then in dealing with each constraint separately is one of the main reason that makes us prefer the Generalized Forward-Backward algorithm over a standard interior-point approach. Indeed, for instance, in the interior-point approach proposed by Mehrotra [START_REF] Mehrotra | On the implementation of a primal-dual interior point method[END_REF], at each iteration, one needs to solve several linear systems of the size of the number of linear constraints (in our case R × m). While a careful implementation, using for instance a Choleski factorization and a proper use of the diagonal Hessian, helps in reducing computational complexity, the cubic dependency in the number of kernels can hardly be avoided.

Let us emphasize that from our experimental study of interior point and proximal methods, when an accurate solution of problem ( 6) is needed, IP methods may do a better job as they tends to provide (slightly) lower objective values. However, following the lines of Bottou and Bousquet [START_REF] Bottou | The trade-offs of large scale learning[END_REF], in large-scale machine learning settings, finding an approximate solution of the optimization problem is sufficient, and in this situation, we advocate the use of proximal gradient methods as they are far more efficient.

Analysis

Our objective in this section is to derive some theoretical understandings of the kernel approximation impact on the SVM problem (2) and thus on the multiple kernel learning. Similar analyses have already been carried out by Fine et al. [START_REF] Fine | Efficient svm training using low-rank kernel representations[END_REF] and Cortes et al. [START_REF] Cortes | On the impact of kernel approximation on learning accuracy[END_REF]. The former work proposed a bound on the variation of the SVM optimization problem objective value when the kernel is replaced by an approximated one. The latter one instead, analyzed how the decision function f (x) varies with respect to the Frobenius norm of kernel difference. In our work, we focus our interest on the norm difference between the minimizer of the exact SVM problem (2) and the minimizer of the approximated one (4). We justify our interest for the vector α as it plays a major role in the computation of the weights d k (see Equation ( 3)). While we agree that bounding α -α does not tell how |d k -d k | varies, it provides a partial result that gives us some intuitions about the key components in the kernel approximation that help in reducing the norm difference.

Proposition 2. Suppose that K d = k d k K k is positive definite. Denote K k as a low-rank approximation of K k .
Denote α as the solution of problem (2) and α the solution of the same problem but with kernel

K d = k d k K k .
Then, we have,

α -α 2 ≤ 2s λ min (8)
where λ min is the smallest eigenvalue of K d and s is so that

s = max α:y α=0,0≤αi≤C 1 2 α 2 2 k d k (K k -K k ) F
The complete proof is given in the supplementary material. It relies on bounds on variation of minimizers of perturbed strictly positive definite quadratic programming problems [START_REF] Hx Phu | Some properties of boundedly perturbed strictly convex quadratic functions[END_REF]. This proposition gives us the intuition that the impact of the low-rank kernel approximation is controlled by one term independent of the approximation (λ min )

and by the Frobenius norm of the kernel difference.

Numerical experiments

The objectives of these experiments are three-fold. At first, we want to provide empirical evidences that in largescale situations the low-rank MKL method we propose is faster than the state of the art algorithm SMO-MKL while yielding to similar performances as long as the kernel approximations are accurate enough. Secondly, we provide support that proximal gradient algorithms are more efficient than interior points methods for solving problem [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF].

Finally, experiments depicting the gain in efficiency brought by our non-trivial application of the GFB algorithm compared to its naive implementation have also been reported.

General settings

Computations have been carried out on a 16-core Intel Xeon E5530@2.4 GHz machine with 24 GB of memory. All the codes have been written in Matlab and have been run on a single core of the above machine. For the SMO-MKL, we have used the C implementation of the authors. We have considered datasets from the UCI and Libsvm repository and their statistics are summarized in Table 1.

For each of these datasets, we have randomly split 70% -30% the examples in training and test sets. 50 Gaussian kernels, with bandwidth σ ranging logarithmically from 10 0 to 10 2 , have been considered. Since we are interested in relative performances, we have arbitrarily set C = 1000, which is a value that yields good accuracies. For SMO-MKL, we used the default stopping criterion and all its default parameter. Stopping criterion of our MKL algorithm is based on the maximal variation of the weights between two iterations, which should be lower than 10 -4 . For our approach, the inner low-rank SVM algorithm is stopped when the relative variation of its objective value is lower than 10 -4 or when 3000 iterations have been reached.

Comparing with SMO-MKL

For supporting our claims, we have evaluated the computational speed (including the time for computing low-rank kernels) of our algorithm as well as the accuracy of the resulting decision function with respect to the rank of the kernel approximations. We show that for large-scale datasets our low-rank multiple kernel algorithm is more efficient than SMO-MKL while achieving similar or slightly worse performances. Note that the setting we are interested in is the large-scale number of examples/few kernels ones and in such a situation SMO-MKL is more efficient than the spectral projected gradient MKL [START_REF] Jain | Spg-gmkl: generalized multiple kernel learning with a million kernels[END_REF].

Averaged performances (accuracy and running time) over 5 runs and increasing quality of approximation are reported in Figure 1 and Figure 2 for p = 2 and in Figure 3 and Figure 4 for sparse multiple kernel learning with p = 1.1.

Note that we have chosen this latter value since SMO-MKL is not able to handle the exact sparse case p = 1. For the smallest datasets, running time comparison is in favour of the SMO-MKL, and we explain this fact through the advantage brought by the C implementation and the memory cache handling. As soon as, the number of examples is larger than 10 5 , our algorithm can be significantly faster (up to an order of magnitude in some cases) while achieving similar performances. Similar behaviors are observed for p = 1.1, and this can be expected as p does not influence the problem (6) we are solving.

Scalability with respects to the number of kernels is reported in Figure 5. Again, running time is in favour of our low-rank approach.

Proximal gradient vs interior-point

For all experiments described above, we were not looking for a precise solution of the optimization problem but as advocated by Bottou et al. [START_REF] Bottou | The trade-offs of large scale learning[END_REF], we were searching for approximate solutions of the problem [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF].

The next experiments provide some empirical evidences on our claim regarding efficiency of the gradient proximal algorithm with respects to interior point (IP) methods for solving problem [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF], at least for the precision we looked for.

The interior point method we have considered is similar to the one proposed by Fine et al. [START_REF] Fine | Efficient svm training using low-rank kernel representations[END_REF]. We have implemented a Mehrotra predictor-corrector algorithm [START_REF] Mehrotra | On the implementation of a primal-dual interior point method[END_REF] but capable of handling a larger number of equality constraints. A careful attention has also been brought to the fact that the Hessian matrix of problem ( 6) is diagonal. Despite this, we emphasize that two linear systems involving a dense R × R positive definite matrix has to be solved at each iteration making the method computationally very expensive.

Remind that the stopping criterion of our low-rank SVM problem is based on relative variation of the objective value that should be lower than 10 -4 . For the sake of fair comparison, we have early stopped the IP methods when the difference between their objective value and the one of the proximal method is less than 0.1%. Note that at this precision the infinite norm of the difference of the two minimizers is in the order of 10 -6 . The experimental set-up is the same as for the large-scale experiments except that the rank R has been set to 500 so as to emphasize the difference in running times and that we have considered only 10 kernels and we have arbitrarily set d k = 1 m . We have compared these two methods on their efficiency for solving problem [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF]. Figure 6 depicts the running time of both algorithms as the number of training examples increases. We note that GFB is always more efficient with a gain factor varying from 1.1 to 8.

Splitting equality constraints

One of our contribution is to propose a non-trivial use of the GFB algorithm for solving our learning problem [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF]. This contribution is based on the idea that splitting the constraints in several blocks lead to reduce computational complexity. In what follows, we provide empirical supports of this fact.

The running time of two GFB algorithms, one which considers all linear constraints of problem (6) at once and one which splits the constraints (as described in Algorithm 2) have been reported in Figure 7. Again, we note the drastic gain in efficiency brought by our specific GFB approach compared to a naive GFB. We gain about 1 to 3 order of magnitude in efficiency.

Conclusion

The paper addresses the problem of large-scale MKL learning. In a large-scale setting, kernel methods usually resort to low-rank approximation of the original Gram matrices. In this work, we make the hypothesis that our base kernels for MKL are low-rank approximations of the original ones. We propose a novel MKL optimization problem which takes advantage of such approximations as well as a gradient proximal algorithm for solving the the related problem. Partial theoretical results dealing with the impact of the kernel approximation on the kernel combination weights are also provided as well as empirical evidences on the better scalability of our approaches compared to stateof-the-art MKL solvers.

As for future work, we plan to apply our work to real-world large-scale problems related to computer visions and audio signal classifications.

Let us first introduce the theorem which bounds the norm difference between the minimizer of a strictly convex quadratic programming problem and its perturbed version.

where K d = k d k K k , with K k a low-rank approximation of K k . Note that the above problem ( 9) is a perturbed version of Problem (2) with the following perturbation

p(α) = 1 2 α Y(K d -K d )Yα
From the definition of s, in our specific case, we have

∞ > s ≥ sup α:y α=0,C≥αi≥0 1 2 α Y(K d -K d )Yα
Let us now show that such s exists and that it can be upper-bounded by a term depending on K d -K d F .

First note that the perturbation function p(•) is continuous and that the constraints is subset of a closed subset of R . Hence there exists s so that

≤ Yαα Y F K d -K d F = α 2 2 K d -K d F
Hence, we can choose s as s = max α:y α=0,C≥αi≥0

1 2 α 2 2 ( K d -K d F )
Then from Theorem 1 and the definition of s, we get the desired inequality given in Equation ( 8).

Algorithm 1

 1 Alternate optimization algorithm for p MKL 1: set k=1, initialize the kernel weights {d k } k 2: repeat 3: optimize an SVM with kernel k d k K k 4:

  1 and F 2 by choosing F 1 (z) = I b L ≤z≤b U the indicator function on the box constraints delimited at each coordinate by the coordinates of b L and b U (in our case, we have b L = 0 and b U = C) and F 2 (z) = I Az=0 , the indicator function on an affine set defined by a matrix A. In our case, A is the matrix of size (R × m + 1) × (R × m + ) defining the linear equality constraints as given in Problem (6). The proximal operators of F 1 (•) and F 2 (•) are simple to derive and in particular, we have :

Remark 1 .

 1 z) leaves unchanged all components but γ k and α. The computationalcomplexity for computing A k (A k A k ) -1 A k is now of the compelling order of O(R 3 + R 2 ).And since we have to compute this proximal operator m times, the complexity is now only linear in m. A detailed instantiation of our nontrivial application of the GFB algorithm to our particular problem is given in Algorithm 2. In this algorithm, P b L ,b U (•) defines the projection on the box constraints defined by the lowest coordinate b L and highest coordinate b U . Note that in practice we have set the step size ζ n to 1 for all iterations since the Lipschitz constant of F 0 is 1. According to Raguet et al. [22, Th 2.1], the sequence {z n } can be made strongly convergent to the minimizer of Problem (

10 :%

 10 Compute the forward-backward steps related to the constraints {γ k = √ d k V k Yα} m k=1 11:

Theorem 1 . 2 α 2 α

 122 [20, Corollary 4.2] Let α be the solution of strictly convex quadratic programming problem (P ) min α∈C 1 Qα + c α with Q a positive definite matrix R × , c a vector of R and C a non-empty convex set. Let α be any global infimizer the perturbed version of this problem, where the perturbation stands in adding a function p(α) to the objective function, then the following inequality holds α -α 2 ≤ 2s λ min with λ min being the smallest eigenvalue of the matrix Q and s is so that sup α∈C |p(α)| ≤ s < ∞. Consider now α as the solution of Problem (2) and α the solution of problem min α:y α=0,0≤αi≤C 1 YK d Yα -1 α

  sup α:y α=0,0≤αi≤C |p(α)| ≤ max 0≤αi≤C |p(α)| = s < ∞Now, the following inequalities provide us with an upper-bound on p(α).2|p(α)| = α Y(K d -K d )Yα = tr(Yαα Y(K d -K d ))

Table 1 :

 1 Dataset statistics

	Data sets	# examples Dimension	comments
	svmguide 1	7089	4
	mushrooms	8124	112
	magic	19020	10
	pokerhand	25010	10	class 0 against all
	codrna	59535	8
	seismic	78823	50	class 3 against all
	adult	48842	14
	ijcnn1	141691	22

Figure 5: Running time with respect to the number of kernels of our low-rank MKL and the SMO-MKL. Rank of approximations have been fixed to 100 so that both algorithms provide similar performance accuracies. [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF]. GFB merged constraints and GFB respectively refer to a naive implementation generalized forward-backward algorithm which handles all constraints at once and a GFB as described in Algorithm 2. The number of kernels is m = 20 and each kernel has been approximated using rank-500 matrices. All of these results in 5001 linear constraints to deal with.


J Neggers

• J P M Hoefnagels

• F Hild

• S Roux

• M G D Geers

Direct stress-strain measurements from bulged membranes using topography image correlation

Keywords: digital image correlation, surface profilometry, thin film, membrane, full-field measurement, strain, stress

This paper discusses an experimental method to characterize thin films as they are encountered in micro-electronic devices. The method enables the measurement of the stress and strain of pressure deflected bulged membranes without using a priori defined bulge equations. An enrichment to the Global Digital Image Correlation method is detailed to capture the membrane strain and curvature while robustly dealing with acquisition noise. The accuracy of the method is analyzed and compared to the standard bulge test method. The method is applied to a proof of principle experiment to investigate its applicability and accuracy. Additionally, it is shown for two experimental cases that the method provides accurate results, although the bulge equations do not hold.

Introduction

Thin films are key in miniaturizing components in consumer electronic devices. It is known that the properties of most thin films are size dependent, since their intrinsic microstructural length scale is of comparable size to the design length scale [START_REF] Nix | Mechanical properties of thin films[END_REF][START_REF] Vinci | Mechanical behavior of thin films[END_REF]. Moreover, the processing techniques, and the adjacent materials, also influence the mechanical response [START_REF] Arzt | Size effects in materials due to microstructural and dimensional constraints: a comparative review[END_REF][START_REF] Gruber | Size effects on yield strength and strain hardening for ultra-thin Cu films with and without passivation: A study by synchrotron and bulge test techniques[END_REF]. It is therefore important to perform an in-situ characterization of these films, i.e. with the selected manufacturing processes as used in the device.

The bulge-test is a well known technique for testing thin film properties, especially at small scales [START_REF] Nix | Mechanical properties of thin films[END_REF]. In a bulge-test, a membrane of specific shape (rectangular, circular) is deflected under pressure. The deflection and pressure can be converted into stress and strain using shape specific equations. Generally, bulge equations are derived using energy minimization methods [START_REF] Beams | Mechanical properties of thin films of gold and silver[END_REF][START_REF] Tabata | Mechanical property measurements of thin films using load-deflection of composite rectangular membranes[END_REF][START_REF] Timoshenko | Theory of plates and shells[END_REF] resorting to approximate membrane displacement fields. In these derivations several assumptions are used to make the system solvable. Although these assumptions may be clearly specified, they are often only approximately fulfilled and hence the resulting error/uncertainty on the estimates of the mechanical properties of these films is difficult to evaluate.

For micro-electronic systems the miniaturized rectangular bulge equations developed by Vlassak et al. [START_REF] Vlassak | A new bulge test technique for the determination of Young's modulus and Poisson's ratio of thin films[END_REF][START_REF] Xiang | Plane-strain bulge test for thin films[END_REF] are popular. The reason is that, for these systems, it is possible to manufacture the rectangular membranes with high precision by back etching a Si wafer using KOH (potassium hydroxide), thereby improving the accuracy significantly. Vlassak et al. have shown that for rectangular membranes (width 2a along the x direction, length 2b along the y axis, and thickness t), with sufficiently large width-to-length aspect ratio b/a > 6, the center part of the membrane closely approximates a plane-strain state and deflects cylindrically. The bulge equations then obtained to compute the stress and strain from the pressure and deflection are,

κ xx = 2δ a 2 + δ 2 , (1)
σ xx = P tκ xx , (2)
ε xx = 1 aκ xx sin -1 (aκ xx ) -1, (3)
where, P is the pressure and δ the deflection at the center (or apex) of the membrane and κ xx its curvature along direction x. Equation ( 2) is also known as the hoop stress equation, and is derived from static equilibrium assuming negligible flexural moment. Equation ( 3) is purely geometrical and based on the assumption that the membrane deflects to a circular profile. It has been shown that this method is accurate for films with large width-to-thickness ratios a/t > 1000 (depending on the level of strain) because these thin films approximate a free rotating hinge at the boundary [START_REF] Neggers | A Global Digital Image Correlation Enhanced Full-Field Bulge Test[END_REF]. Thicker films conform poorly to the circular deflection profile that connects the two edges. In particular close to their edge, flexural moments are the largest and hence the evaluation of the curvature from Equation (1) may lead to significant errors. Alternatively, as long as the film edges are removed from the analysis, the hoop stress equation, Equation (2) may still be valid.

More importantly, all bulge equation methods are limited to homogeneous, unstructured membranes, with specific shapes. Nonetheless, in current microelectronic applications more complex thin film shapes are often found [START_REF] Madou | Fundamentals of Microfabrication[END_REF][START_REF] Rogers | Materials and mechanics for stretchable electronics[END_REF], even with multiple layered and structured heterogeneous thickness materials. Since a bulgeequation based method can only provide a global (or averaged) material response, it would not apply to such thin films.

In this paper, a new method is proposed for measuring thin films that are loaded with a pressure medium similar to the bulge test method. In this new method, full-field Optical Confocal Profilometry is combined with dedicated Global Digital Image Correlation (GDIC). The proposed method enables the full-field strain and curvature fields to be measured for any type of deformed sample measured with a profilometric full-field technique.

Using the fact that the pressure is uniform under the sample, the full-field membrane stresses are derived from the pressure and curvature if the bending stresses are negligible and the local stress state (e.g., planestrain, biaxial loading) is known (Section 2.3). The proposed method utilizes the wealth of full-field data that allows for noise attenuation, resulting in improved accuracy in cases where the bulge equations apply, as discussed in Section 3. Moreover, the proposed method is direct, i.e. measuring the local stress and strain without assuming any global deformed shape. The method is therefore more widely applicable in terms of e.g. membrane thickness, membrane shapes, and structured membranes, as demonstrated in Section 4.

Methods

The proposed method has been developed for thin films as they are used in micro-electronics applications. For these films, the small strain regime is usually the most relevant for service conditions. This regime also implies small curvatures, which puts high demands on the used experimental technique. Therefore, a microscopic surface profilometry technique (such as confocal profilometry or phase shift interferometry) is chosen as a measurement method. This type of system measures the "surface profile" of the sample, which is the height of the sample at every pixel location. The resulting data are images (i.e. 2D matrices), where the pixel information is not a gray value but height. As a result, in this study, GDIC is applied to identify the 3D displacement vector field, from (2D) height-images, and therefore called Quasi-3D. Such a Quasi-3D GDIC method has been used before on Atomic Force Microscopy data [START_REF] Han | Measuring nanoscale stress intensity factors with an atomic force microscope[END_REF]. Additionally, the proposed method aims to capture the curvature field, which is the second derivative of the position field. It is well known that noise has a detrimental effect on derivatives. It is therefore imperative to deal with the acquisition noise robustly, for which particular choices are made in the GDIC procedure.

Global Digital Image Correlation

Digital Image Correlation (DIC) in general consists of measuring the displacement field between two images obtained from two increments of loading in an experiment [START_REF] Chu | Applications of digital-image-correlation techniques to experimental mechanics[END_REF][START_REF] Bruck | Digital image correlation using Newton-Raphson method of partial differential correction[END_REF][START_REF] Hild | Digital Image Correlation: from Displacement Measurement to Identification of Elastic Properties -a Review[END_REF]. The first image is usually a picture of the reference configuration f , and the second image a picture of the deformed configuration g. If the pictures capture a (speckle) pattern, which is deformed with the underlying sample, then the gray level conservation for a pixel at location x reads

η(x) = f (x) -g(x + u(x)e x + v(x)e y ), (4)
where u and v are the in-plane components of the displacement respectively in x and y direction, and η is the acquisition noise. In the case of profilometric data, the pixel values contain height information, which may evolve between two measurement increments. Thus the gray level conservation is rewritten as "topography" conservation,

η(x) = f (x) -g(x + u(x)e x + v(x)e y ) -w(x) , (5)
where w denotes the out-of-plane displacement, and η now represents noise, but also height distortions not captured in w. Note that the topography is in data structure similar to a 2D image, thus the position vector remains 2D,

x = xe x + ye y , (6)
in contrast with the 3D displacement vector,

u(x) = u(x)e x + v(x)e y + w(x)e z . (7)
To write the problem in terms of a limited number of unknowns, the displacement field is obtained through the minimization of a functional over of subspace of functions ϕ,

u(x) ≈ u * (x, λ) = 3n i=1 λ i ϕ i (x), ( 8
)
where λ is a column vector with the degrees of freedom for a chosen set of basis functions ϕ, detailed later on. The displacement field u * is then obtained by minimizing the squared residual with respect to λ,

λ = argmin λ (Φ) = argmin λ Ω η(x, λ) 2 dx, ( 9
)
where η is now dependent on the degrees of freedom through the substitution of Equation ( 8) into Equation [START_REF] Beams | Mechanical properties of thin films of gold and silver[END_REF].

For this purpose, an iterative procedure is performed, where λ consists of a current value λ k , and an unknown iterative update δλ

λ k+1 = λ k + δλ. ( 10
)
where k denotes the current iteration. The update is found by solving the linearized form of Equation ( 9),

∀ j 3n i=1 Ω (ϕ j •G)(G•ϕ i ) dx δλ i = Ω ϕ j •G (f -g) dx, ( 11
)
where G is the image gradient, defined below, and,

g = g(x + u * (x, λ k )e x + v * (x, λ k )e y ) -w * (x, λ k ),
is the corrected version of g(x) using the current estimate of the displacement field. Equation ( 11) is usually written in matrix form,

M δλ = b, (12)
where the tangent matrix M and the right-hand member b are conceptually similar to established GDIC [START_REF] Hild | Digital Image Correlation: from Displacement Measurement to Identification of Elastic Properties -a Review[END_REF][START_REF] Hild | Comparison of Local and Global Approaches to Digital Image Correlation[END_REF]. The main difference is the vector G, which is the image gradient enriched with an out-of-plane component

G = ∂f ∂x e x + ∂f ∂y e y -e z , (13)
and the three-dimensional nature of the basis functions

ϕ = ϕ x e x + ϕ y e y + ϕ z e z . (14)
To build g(x), the image g(x) needs to be estimated at locations x + u 2 (x), for which a cubic spline interpolation scheme is used, where u 2 denotes the in-plane displacement field. This interpolation is not exact, and is a source of error [START_REF] Schreier | Systematic errors in digital image correlation caused by intensity interpolation[END_REF][START_REF] Schreier | Systematic Errors in Digital Image Correlation Due to Undermatched Subset Shape Functions[END_REF]. Interestingly, the out-of-plane component does not enter the interpolation scheme, and does not suffer from interpolation errors in a direct way.

The iterative system is solved until convergence is reached. The chosen criterion is such that the norm of the degree of freedom increments, δλ, becomes sufficiently small, i.e. δλ < 10 -6 .

The implementation used in this paper is considered a "global" DIC method because the system is minimized globally, instead of minimizing each zone of interest independently. However, the out-of-plane enrichment can also be applied to "local" DIC methods by following exactly the same procedure.

Basis Functions

In a general DIC problem, a displacement vector is to be determined for each pixel in an image. This makes the problem ill-posed, even when not considering image noise. The solution is to reduce the number of unknowns to less than the number of pixels, for instance as was done in Equation ( 8). If there is uncorrelated (i.e. white) noise, then using more pixels per unknown will further reduce the noise sensitivity, thereby improving the accuracy. However, using too few unknowns restricts the kinematics of the solution, which decreases the accuracy. For that reason, it is important to use a basis that adequately captures the full kinematics of clamped bending membranes, inside the Region Of Interest (ROI), with a minimum number of degrees of freedom. Inspired by continuum solutions for loaded plates, a 2D polynomial basis is used. The C ∞ continuity property of polynomials makes these functions well suited for calculating curvature fields, noted that polynomials of at least second order are required to be sensitive to curvature. The basis functions are vector functions where one component of the field is expressed as the product of two terms of a polynomial series,

ϕ = χ 0 ζ 0 ϕ = χ 0 ζ 1 ϕ = χ 0 ζ 2 ϕ = χ 1 ζ 1 ϕ = χ 1 ζ 4 ϕ = χ 2 ζ 3
ϕ x = ϕ y = ϕ z = χ a ζ b , (15)
where χ and ζ are normalized coordinates for x and y such that -1 ≤ χ, ζ ≤ 1. The complete basis is formed by using all combinations of a and b up to a certain order p, see Figure 1 for some examples. Each basis function is used three times, each with its own degree-of-freedom (for each component of the displacement vector). For example for p = 2,

u * = λ 1 χ 0 ζ 0 + λ 2 χ 0 ζ 1 + λ 3 χ 1 ζ 0 + λ 4 χ 0 ζ 2 + λ 5 χ 1 ζ 1 + λ 6 χ 2 ζ 0 , v * = λ 7 χ 0 ζ 0 + λ 8 χ 0 ζ 1 + λ 9 χ 1 ζ 0 + λ 10 χ 0 ζ 2 + λ 11 χ 1 ζ 1 + λ 12 χ 2 ζ 0 , w * = λ 13 χ 0 ζ 0 + λ 14 χ 0 ζ 1 + λ 15 χ 1 ζ 0 + λ 16 χ 0 ζ 2 + λ 17 χ 1 ζ 1 + λ 18 χ 2 ζ 0 .
This particular decomposition into separate degrees of freedom for each component of the displacement vector is one of the possible ways to define the basis functions.

In particular, this will give 3 basis functions for p = 0, 9 basis functions for p = 1, 18 basis functions for p = 2, 30 basis functions for p = 3, etc. Because of the fact that the matrix elements M ij involve the scalar product of pairs of functions with a "weight" that is the tensor G ⊗ G, there is no gain in choosing an orthogonal polynomial basis. The same space of function will be generated by the above polynomials and say Legendre polynomials of the same order. Hence the resulting displacement, and convergence rate is independent of the particular choice provided the generated space is the same.

Curvature, Stress and Strain

Let us revisit the assumptions in the bulge equations discussed in the introduction. Equation ( 1) is a purely geometrical formulation to obtain the curvature from the deflection, assuming the membrane deflects according to a circular profile. For thin membranes and deflections larger than the membrane thickness this assumption holds, attributing to the established accuracy of the bulge test method. However, with the GDIC method discussed in this paper the curvature tensor κ(x) at each pixel location can be computed without resorting to the circularity assumption (see also [START_REF] Neggers | A Global Digital Image Correlation Enhanced Full-Field Bulge Test[END_REF]), by first defining

κ(x) = ∇n(x), (16)
where n(x) is the normal vector, which in turn is the gradient of the position field z(x) (corrected for rotations)

n(x) = ∇z(x) ||∇z(x)|| . ( 17
)
The position field is measured as the topography g(x), however, it also includes the additional pattern. The membrane topography (i.e. the required position field) is more accurately obtained by by applying the measured displacement field to the initial (flat) membrane position. Because the displacement fields are forced to be smooth, by the polynomial basis, noise in the curvature is here controlled by the basis order, as will be shown further down. Equation ( 3) is a similar geometrical formulation as Equation [START_REF] Nix | Mechanical properties of thin films[END_REF], which expresses the strain as function of the curvature assuming a circular deflection profile. For thicker or inhomogeneous membranes the circular deflection profile assumptions is violated, especially at the membrane boundary and heterogeneous points. However, the GDIC method recovers the displacement field, from which the strain field is easily obtained without assumptions, by e.g. computing the Green-Lagrange strain tensor E in the membrane plane.

Equation ( 2) is a simple balance equation evaluated without recourse to the small displacement hypothesis. It is valid as long as the bending moment can be neglected, but the actual size of the membrane does not come into play. Its validity depends on the flexural moment within the region over which the curvature is evaluated. However, the flexural moment is proportional to the change in curvature, which is measured and thus can be assessed.

At the apex of bulged square membranes, the pressure is balanced with membrane stresses in two directions. If the curvatures in both directions are equal then the principal stresses are also equal if the material is isotropic,

σ xx = σ yy = P 2t κ . ( 18
)
If a part of the membrane deforms axisymmetrically, the principal stresses of that membrane part are related to the meridional and circumferential curvatures [START_REF] Hill | A theory of the plastic bulging of a metal diaphragm by lateral pressure[END_REF][START_REF] Hsu | Identification of response functions from axisymmetric membrane inflation tests: implications for biomechanics[END_REF],

σ m = P 2t κ c , (19)
σ m = P t κ c 1 - κ m 2κ c . ( 20
)
where c and m denote the meridional and circumferential directions respectively. For kinematically more complex cases, closed-form expressions cannot be obtained, and numerical schemes are required, such as Finite Element Method Updating [START_REF] Avril | Overview of identification methods of mechanical parameters based on full-field measurements[END_REF].

To conclude, the assumption of neglecting flexural stiffness remains in the proposed methodology, yet only for relating stress to curvature. The assumption is appropriate if a restrained region is analyzed, such that the change in curvature is insignificant over the analyzed region.

Experimental Procedure

Two types of samples are used in the present paper, a square (1 × 1 mm 2 ) membrane and a rectangular (1 × 6 mm 2 ) membrane (Figure 2). Both samples are manufactured by deposition a 100 nm thick Si 3 N 4 layer on a monocrystalline 200 µm thick Si wafer. A freestanding membrane is created by back KOH etching the wafer up to the Si 3 N 4 layer. The membrane dimensions are chosen such that the bulge equations are known to perform well [START_REF] Vlassak | A new bulge test technique for the determination of Young's modulus and Poisson's ratio of thin films[END_REF], enabling a comparison with the GDIC method presented in this paper.

These samples are processed with high precision micromanufacturing techniques, resulting in a surface roughness that is smaller than measurable with the used optical profilometer. Since a height correlation is to be

(a) (b) (c) (d) (e) (f) (g)
Fig. 3 The bulge test apparatus with, (a) the sample, (b) the sample holder, (c) the 270 ml chamber, (d) the 30 ml chamber, (e) the linear actuator connected to the piston, (f) the three pressure sensors, (g) the heating stage.

performed on these samples, a height pattern is required. This pattern is created by coating the membrane surface with a colloidal suspension of Ag particles (80-500 nm diameter) dispersed in ethanol. After evaporation of the suspending fluid, the deposited particles (or rather clusters of particles as they tend to aggregate) showed enough physical adhesion to withstand the bulge test without detaching from nor moving on the surface. Note that the basis functions in the GDIC method have support over the entire ROI, allowing the GDIC method to work efficiently with relatively sparse patterns. The sparse pattern allowed the evaluation of the unloaded freestanding film with and without the pattern. Application of the pattern did not cause a measurable deformation, confirming the insignificant influence on the mechanics.

To test the membranes a custom made low profile bulge test apparatus is specifically designed to allow for optical access to the sample surface (Figure 3). The apparatus consists of two ethanol filled stainless steel chambers of 270 ml and 30 ml respectively (Figure 3cd), which are connected with a valve. The sample is mounted on the smaller chamber with an interchangeable sample holder (Figure 3b). Additionally, the apparatus is equipped with three pressure sensors with respective pressure ranges [0-5], [0-15], [0-50] bar (Figure 3f). The pressure sensors are high performance millivolt output transducers from GE Sensing, with an accuracy of 0.04% of their rated maximum. The two large range sensors are mounted on the smaller cham-ber while the most sensitive sensor is mounted on the larger chamber. The pressure difference is applied with a piston radius of 1.6 mm, connected to the smallest chamber, driven by a linear actuator with a closed-loop DC motor (Figure 3e). The actuator steps in 50 nm increments with a total stroke of 25 mm, and a maximum applied force of 70 N. The two chamber design allows experiments in two pressure ranges, i.e. by closing the valve to the second chamber the most sensitive pressure sensor is naturally protected from overloading, while at the same time, the effect of the piston motion on the pressure is amplified. Finally, a heating stage (Figure 3g) serves to keep the entire system at a constant temperature, slightly above the room temperature (∼ 28 • C), to eliminate pressure changes due to thermal expansion.

The bulge test apparatus fits under the Sensofar Optical Profilometer, which is used to measure the surface roughness of each pressure increment. The experimental procedure is controlled via custom NI LabVIEW code that follows a user-defined pressure loading curve in a closed loop, with intermittent pressure hold periods for the profile measurements. The height profile is captured by a Sensofar Plµ2300 Optical Confocal Microscope using a CCD camera (definition: 557 × 557 pixels) with a Nikon EPI 50x objective lens, resulting in a square field of view of 184 × 184 µm 2 . The objective lens has a numerical aperture of 0.80, which together with the monochromatic blue light (λ = 470 nm) results in an in-plane resolution of 358 nm (i.e. slightly larger than the resolution of one pixel of the CCD sensor (332 nm)). For the patterned membranes discussed in this paper, the out-of-plane resolution has been measured to be ∼20 nm from the RMS of the difference between two height profiles of the same area on the sample. The difference between in-plane and out-of-plane resolution will also be apparent in the accuracy analysis of Section 3.

Virtual Experiment

In this paper, proof of principle experiments on square and rectangular membranes are performed. To make a thorough evaluation of the accuracy of the method, virtual experiments are used. A virtual experiment is a synthetic procedure that produces similar data as expected from the real experiment, yet, with the additional knowledge of the reference.

To create realistic virtual height profiles, the virtual experimental procedure consists of deforming a real experimental height profile f (Figure 4b) with an FE displacement field (Figure 4c) to obtain a virtual deformed height profile g (Figure 4d). Different increments in the Fig. 4 The reference profile f (x) and the deformed profile g(x) are measured with the experimental setup (a,b,e). Additionally, the reference profile is virtually deformed with an FE displacement field (c) to obtain a virtually deformed profile (d), which approximates the real experiment.

FE simulation are used to obtain the incremental pressure height profiles similar to a real experiment.

Two different FE geometries are used, corresponding to the rectangular (Figure 2a) and to the square (Figure 2b) membrane. In the FE simulations, a dense mesh is chosen to match the pixel discretization level of the height images. This results in a mesh of approximately 50,000 3D 4-node bilinear Mindlin shell elements, which is densest at the field of view (FOV) and opens towards the boundaries. Note that, for illustrative reasons, the grid shown in (Figure 4c) is not the actual used mesh, which would be too dense to show any meaningful detail.

At this stage it is possible to determine the order of the basis functions required to accurately capture the kinematics in the virtual experiment. Figure 5 shows the mean absolute displacement error between the FE displacement field u fem and displacement fields approximated with increasing polynomial order u * . An error measure is introduced as the average of the lengths Fig. 5 The mean absolute error as a function of polynomial order obtained by comparing the FE displacement field with approximated displacement fields with increasing polynomial order.

of all difference vectors between the FE displacement field and the approximated displacement fields, evaluated over all pixels k,

E u = 1 n n k=1 ( u fem (x k ) -u * (x k ) ) , (21)
where the norm is defined as the Euclidean norm of, in this case, the difference vector. As expected, the mean absolute error decays for increasing order, and reaches an error level of less than one nanometer for sets of basis functions with p ≥ 4. Since this error is smaller than the resolution of the measurement system, it is chosen to truncate the series of basis functions at p = 4.

3 Accuracy Analysis

Error Fields

The use of the virtual experiment allows for direct comparisons between the results obtained with the proposed GDIC method and the FE reference solution.

Similarly to Equation ( 21), a pixel-wise displacement error is computed as the vector difference between reference and the displacements estimated by GDIC. The different components and magnitude of the displacement error fields are shown in Figure 6 for the correlated ROI. More importantly, Figure 6a shows that the maximum error in displacement is approximately 10 nm while the average absolute error is much smaller. Even though there is no noise present in this virtual experiment, the result is based on a real experimental image f with the same pixel discretization and pattern as the real experiment. Note the difference in the in-plane accuracy versus the out-of-plane accuracy (Figure 6b-d).

The error in w is approximately 20 times smaller than in u and v. This coincides with the observed accuracy of the confocal microscopy set-up itself.

Error vs. number of DOFs

The error sources in DIC are in general threefold:

1. acquisition noise or pattern changes; 2. subpixel interpolation; 3. displacement approximation.

To minimize the impact of the second error sources, a local cubic spline subpixel interpolation method is used. More interestingly, the first and third error source are dependent on the chosen basis. Using too few (or incorrect) basis functions will introduce systematic errors because the kinematics of the experiment cannot be captured. Conversely, giving the system more freedom, allows the solution to be sensitive to noise, resulting in statistical errors.

In Figures 7ab, image correlations are performed on virtual experiments for an increasing number of degrees of freedom, for various levels of (white) noise added to the deformed profile g. The virtual experiment method allows the application of artificial acquisition noise by adding generated noise fields to the images. The various noise levels are quantified as the standard deviation of the white noise, relative to the image gray level dynamics. It is clearly shown that for the noise-free situation, increasing the number of degrees of freedom only improves the result, while when noise is included, an optimum accuracy is found using a conservative number of degrees of freedom (i.e. 10 to 20 DOFs). These results demonstrate that choosing the right kinematic basis for the GDIC procedure enhances the accuracy. The preferred basis is the one that captures the full kinematics with the minimum number of degrees of freedom. The applicability of a narrow basis, which is ideal in terms of accuracy and robustness, is limited to experiments within the range of this narrow basis. In contrast, the polynomial basis, discussed in this paper, is applicable to many types of experiments, provided they have slowly varying displacement fields, making the proposed method relatively general. In a real experiment, no reference case exists. Therefore, the recommended procedure to select a proper set of basis functions, is to perform the correlation for increasingly rich sets of basis functions. The residual (Figure 7c) will decrease for richer bases (if converged), yet the decrease will flatten at some level of richness. At this point, adding more degrees of freedom to the system tends to make the correlation less accurate due to noise sensitivity. Another method for evaluating the chosen basis, is by examining the residual field. A correlation with too few basis functions will result in a residual field exhibiting long wavelength modulations, while the residual field for an appropriate basis will only contain the acquisition noise, and thus will resemble white noise.

Demonstration experiments

Proof of Principle Experiments

To investigate the accuracy of the method on real cases, four experiments are performed: two experiments deal with rectangular membranes, and two with square ones. The samples are chosen for allowing optimal accuracy when using the bulge equations, to make a fair comparison with the discussed GDIC method, see Section 2.4.

The membranes are pressurized with the bulge test apparatus to a pressure of 1 bar in 100 increments (above atmosphere), and then unloaded to 0 bar in 20 increments. One increment takes approximately 90 seconds including the pressure change, and scanning of the confocal microscope. On each of the recorded images the polynomial GDIC procedure is performed with all basis functions up to p = 2 for the square membranes (i.e. 18 DOFs). For the rectangular membrane the DOFs that operate in y-direction are removed leaving 15 DOFs.

Figure 8 shows the stress and strain results after processing the topographies with both the GDIC method and the bulge equations. The square membrane response is more stiff, which is expected since it is in a more confined biaxial strain state when compared to the plane strain state at which the rectangular membrane is deformed. Consequently, the obtained moduli are the plane strain modulus and the biaxial modulus that are directly related to the Young's modulus (E ≈ 234 GPa) and the Poisson's ratio, (ν ≈ 0.16) and are as expected for these thin Si 3 N 4 films. For the square membranes the stress-strain curves in both x and y-directions are shown, using κ xx and κ yy respectively. The curves overlap showing that the membranes are not anisotropic, at least not within the membrane plane. More importantly, the figure shows that both the bulge equations and the GDIC method give the same results, confirming the validity of the method.

Beyond the bulge equations

The previous results show that the method is comparable to the accuracy obtained by using the bulge equations. However, the goal of this method is to go be- Fig. 9 GDIC method and the bulge equations applied to virtual experimental data with the same material properties, but with various membrane thicknesses. The GDIC results do not deviate from the reference elastic response up to a thickness of 20 µm.

yond the validity regime of the bulge equations. Therefore, the GDIC method is applied to data of thicker membranes which are created by the virtual experiment framework.

A number of virtual experiments are performed (see Section 2.5). All virtual experiments use the same linear elastic material properties with Young's modulus of 235 GPa. Figure 9 shows the stress strain response taken from the center of membranes for all virtual experiments. As expected the bulge equations start deviating from the correct response for thicknesses above 2 µm, for this 1 mm wide membrane. In contrast, applying the GDIC method to the same data shows correct stress strain response up to thicknesses of 20 µm. The GDIC method considers only the membrane inside the ROI, and thus is much less restricted by the flexural stiffness assumption.

For even thicker membranes the stress across the thickness is far from uniform, and no longer satisfies the assumptions made in Equation [START_REF] Vinci | Mechanical behavior of thin films[END_REF]. It may be possible to obtain a corrected stress from the curvature for thicker membranes, but this requires a modification that is outside the scope of this paper.

Inhomogeneous Films

The previous results use only the stress and strain at the center of the membrane, but the GDIC method gives full-field stress strain results. To show that the GDIC method captures features in the membranes, the virtual experiment framework is again applied to generate experimental data. In this case a 24 µm wide strip along the y-axis is modeled using a Young's modulus of E str. = 235 GPa, while the remaining membrane material is modeled three times with different moduli, i.e. E mem. = [130, 170, 235] GPa.

The generated topography images are processed with the GDIC method twice, once masking the structure, and once masking everything but the structure. For all three virtual experiments the stress-strain curves are collected in one point on the structure and in two points on the membrane (See Figure 10). This shows that it is possible to obtain the correct stress strain results for two different materials from a single full-field measurement.

With the proposed method it is always possible to obtain curvature fields and strain fields. However, obtaining the stress fields relies on applying either Equation (2), [START_REF] Schreier | Systematic errors in digital image correlation caused by intensity interpolation[END_REF], or [START_REF] Schreier | Systematic Errors in Digital Image Correlation Due to Undermatched Subset Shape Functions[END_REF], which are only valid for particular cases (e.g. plane-strain, biaxial, or axisymmetric states).

Discussion and Conclusions

The global DIC method has been extended to images containing height profiles, relaxing the brightness conservation principle and enabling for quantitative measurement of in-plane and out-of-plane displacements. This extension is useful for any topographic measurement, in particular optical confocal profilometry, but also e.g. Scanning Probe Microscopes like Atomic Force Microscopy. Furthermore, the out-of-plane quantity is not limited to topography, it can be any measured field, e.g., temperature [START_REF] Maynadier | One-shot measurement of thermal and kinematic fields: Infra-red image correlation (iric)[END_REF].

In the proposed GDIC method the image pattern and the quantified out-of-plane displacement field are both taken from the same topographical data. In order to avoid interference between them, a separation in length scales is required between the length scale of the out-of-plane displacement fluctuations and the pattern length scale. For the presented case of bulged membranes this separation of length scales is achieved by applying low order polynomial basis functions with support over the entire region of interest and by controlling the surface marking with a deposition of small clusters of nanoparticles. Considering that the membrane mechanics are the goal of the analysis, it is important that the added pattern does not mechanically influence the membrane. For the discussed cases it was shown that the pattern did not significantly influence the mechanics.

The introduced polynomial basis is particularly suited for the bulging of membranes, and not limited to square or rectangular membranes. The polynomial basis functions have a wide support and capture the membrane kinematics with only a few degrees of freedom. This makes the method highly robust with respect to noise, but also with respect to the initial guess. Moreover, the measured displacement fields are C ∞ continuous, making them ideal for the calculation of curvature fields.

Through the use of virtual experiments, where the input (FE) displacement field is known, the reference solution can be compared with the measured (via GDIC) displacement fields, through which the accuracy of the method was analyzed. The accuracy depends on the level of the acquisition noise and the richness of the kinematics that is controlled through the number of basis functions. Too restricted kinematics will give inaccurate results because the kinematics of the experiment is not adequately captured; too rich kinematics will make the solution sensitive to noise also leading to inaccurate results.

In a real experiment no reference displacement field exists to asses if an optimal basis is applied. To identify this optimum, it is suggested to perform DIC calcula-tions for increasing number of degrees of freedom until the residual field stabilizes. The residual field should always decrease for increasing number of degrees of freedom, but when the residual decrease becomes small, it suggests that modes that cannot be measured are progressively included thereby inducing a degradation of the conditioning of the tangent GDIC problem, and hence an increased noise sensitivity.

The proposed method recovers the displacement fields, from which strain fields are trivially computed, regardless of the shape or thickness of the membrane. Moreover, the proposed method recovers curvature fields, and it is shown that for particular cases (i.e. planestrain, biaxial, or axisymmetric strain-states) the stress fields can be recovered. Consequently, the bulge membrane shape can be any shape, and as long as the interrogated area fulfills one of these particular strain-states, then it is possible to obtain the stress in that respective area. For cases with more complex kinematics, a coupling with numerical method is to be made in order to obtain the stress.

As a proof-of-principle, bulge test experiments are performed for 100 nm thick Silicon Nitride films. Both the GDIC method and the plane-strain bulge equations are applied to the data, and give comparable results. This shows that the GDIC method at least matches the accuracy of the bulge equations for this particular case. However, the applicability regime of the GDIC method is much less restrictive.

The method applies the same assumption (negligible flexural moments) as the plane-strain bulge equations to obtain stresses. However, by considering only a small portion of the membrane, the assumption is less restrictive. Therefore, application of the method to inhomogeneous membranes with less restrictions on the membrane thickness has shown to accurately recover the stress-strain response of the various materials in the field of view. Consequently, the proposed method is less restricted in terms of membrane thickness, membrane homogeneity and membrane shapes, thereby, alleviating these shortcomings of the original bulge test.

Fig. 1

 1 Fig. 1 Example 2D basis-functions ϕ.

Fig. 2

 2 Fig.2The two sample types used in the experiments, the Si 3 N 4 membrane is 100 nm thick in both types, (a) a drawing of a rectangular membrane (1 × 6 mm 2 ), (b) a photograph of a square membrane (1 × 1 mm 2 ).

  z-component of the displacement error field

Fig. 6

 6 Fig.6Error in the estimated displacement fields for the ROI for a full 4 th order set of shape functions without additional noise. The Euclidean norm is used as the vector norm shown in subfigure (a). Note that the accuracy of the out-of-plane displacement is approximately 20 times better than the accuracy of the in-plane displacement, which coincides with the observed accuracy of the confocal microscopy set-up.

Fig. 7

 7 Fig. 7 Effect of increasing the number of degrees of freedom on (b) the error in displacement, (c) the error in curvature, (a) the residual.

Fig. 8

 8 Fig.8Measured stress strain response for two square and two rectangular 100 nm thick membranes. The GDIC method and the bulge equations perform equally well for such thin membranes.

Fig. 10

 10 Fig. 10 GDIC applied to three heterogeneous membranes, where only the Young's modulus of the membrane material is varied, and the center strip of material has a constant modulus of 235 GPa. Note that the method can capture full-field strain (a) and curvature (b) fields. The red dots in (a,b) show the location where the stress-strain depicted in (c) is taken.

Eindhoven University of Technology, Department of Mechanical Engineering, 5612 AZ Eindhoven, Netherlands Email: j.p.m.hoefnagels@tue.

nl • 2 LMT Cachan, ENS Cachan / CNRS / UPMC / PRES UniverSud Paris 61 avenue du Président Wilson, 94235 Cachan Cedex, France

Acknowledgment

This work has been supported by the Dutch Technology Foundation (STW) and the Dutch Organization for Scientific Research (NWO).