Sébastien Bougleux

Benoit Gaüzère

David B Blumenthal

Luc Brun

Fast Linear Sum Assignment with Error-Correction and no Cost Constraints

come

Introduction

Finding correspondences between structured or unstructured data is a fundamental problem in data processing and analysis, in particular in computer vision, pattern recognition and machine learning. In various application scenarios, data is represented by attributed graphs, and the correspondence problem hence translates to the task to compute those graphs G, from a given collection, that "correspond best" to a given graph H.

One way to formalize correspondence between graphs is to say that two attributed graphs G and H correspond well, if and only if they are close to each other w. r. t. some distance measure. A very sensitive and therefore widely used 1 distance measure is the graph edit distance (GED) [START_REF] Bunke | Inexact graph matching for structural pattern recognition[END_REF][START_REF] Sanfeliu | A distance measure between attributed relational graphs for pattern recognition[END_REF][START_REF] Riesen | Structural Pattern Recognition with Graph Edit Distance. Advances in Computer Vision and Pattern Recognition[END_REF]. It is defined as the minimum cost of an edit path between G and H. An edit path between G and H is a sequence of elementary edit operations that transform G into H, where the elementary edit operations are removal or insertion of an edge or an isolated node and substitution of a node's or an edge's attribute. Each elementary edit operation comes with an associated non-negative edit cost.

Since it is N P -hard to exactly compute GED [START_REF] Zeng | Comparing stars: On approximating graph edit distance[END_REF], algorithms that compute lower and upper bounds are important. Many existing algorithms [START_REF] Riesen | Approximate graph edit distance computation by means of bipartite graph matching[END_REF][START_REF] Zeng | Comparing stars: On approximating graph edit distance[END_REF][START_REF] Riesen | Structural Pattern Recognition with Graph Edit Distance. Advances in Computer Vision and Pattern Recognition[END_REF][START_REF] Gaüzère | Approximate graph edit distance guided by bipartite matching of bags of walks[END_REF][START_REF] Zheng | Efficient graph similarity search over large graph databases[END_REF][START_REF] Carletti | Approximate graph edit distance computation combining bipartite matching and exact neighborhood substructure distance[END_REF][START_REF] Benjamin | Correcting and speeding-up bounds for non-uniform graph edit distance[END_REF][START_REF] Blumenthal | Improved lower bounds for graph edit distance[END_REF][START_REF] Takami | Accurate and fast computation of approximate graph edit distance based on graph relabeling[END_REF][START_REF] Daller | Approximate graph edit distance by several local searches in parallel[END_REF] proceed as follows: In a first step, they lossily transform the problem of computing the graph edit distance between G and H into an instance C of the linear sum assignment problem with error-correction (LSAPE), where C is a real-valued matrix. LSAPE is similar to the linear sum assignment problem for bipartite graphs (LSAP), but also allows node removals and insertions in addition to node substitutions. In a second step, a LSAPE solver is called to compute an optimal solution X for C. Subsequently, X is interpreted as an edit path between G and H and hence gives an upper bound for GED. By using the optimality of X , some approaches also allow to derive a lower bound of the GED from the cost of X under C [START_REF] Zeng | Comparing stars: On approximating graph edit distance[END_REF][START_REF] Zheng | Efficient graph similarity search over large graph databases[END_REF][START_REF] Benjamin | Correcting and speeding-up bounds for non-uniform graph edit distance[END_REF][START_REF] Blumenthal | Improved lower bounds for graph edit distance[END_REF]. In addition to that, LSAPE arises as a subproblem in several approaches that refine the upper bound by local or hybrid greedy search strategies [START_REF] Riesen | Improving graph edit distance approximation using various search strategies[END_REF][START_REF] Riesen | Improved graph edit distance approximation with simulated annealing[END_REF][START_REF] Ferrer | A first step towards exact graph edit distance using bipartite graph matching[END_REF]. This is also the case for approaches based on relaxations of quadratic programming formulations of GED and conditional gradient descent [START_REF] Bougleux | Graph edit distance as a quadratic assignment problem[END_REF][START_REF] Bougleux | A hungarian algorithm for errorcorrecting graph matching[END_REF][START_REF] Daller | Approximate graph edit distance by several local searches in parallel[END_REF].

The runtime performance of the approaches mentioned in the previous paragraph crucially depend on the performance of the LSAPE solver they employ. Existing LSAPE solvers fall in two categories. Solvers of the first kind reduce LSAPE to LSAP [START_REF] Riesen | Bipartite graph matching for computing the edit distance of graphs[END_REF][START_REF] Riesen | Approximate graph edit distance computation by means of bipartite graph matching[END_REF][START_REF] Serratosa | Fast computation of bipartite graph matching[END_REF][START_REF] Serratosa | Speeding up fast bipartite graph matching through a new cost matrix[END_REF][START_REF] Serratosa | Computation of graph edit distance: Reasoning about optimality and speed-up[END_REF]: Using different strategies, they first transform an instance C of LSAPE into an instance C of LSAP and then use standard approaches available for LSAP such as the Hungarian Algorithm for computing an optimal solution X for C. Finally, X is transformed into an optimal solution X for the LSAPE instance C. Algorithms of the second kind directly solve LSAPE by adapting existing approaches for LSAP [START_REF] Jones | Revisiting Volgenant-Jonker for approximating graph edit distance[END_REF][START_REF] Jones | Optimising the Volgenant-Jonker algorithm for approximating graph edit distance[END_REF][START_REF] Bougleux | A hungarian algorithm for errorcorrecting graph matching[END_REF]. Furthermore, LSAPE solvers can be separated in general methods [START_REF] Riesen | Bipartite graph matching for computing the edit distance of graphs[END_REF][START_REF] Riesen | Approximate graph edit distance computation by means of bipartite graph matching[END_REF][START_REF] Jones | Revisiting Volgenant-Jonker for approximating graph edit distance[END_REF][START_REF] Jones | Optimising the Volgenant-Jonker algorithm for approximating graph edit distance[END_REF][START_REF] Bougleux | A hungarian algorithm for errorcorrecting graph matching[END_REF] and cost-constrained methods that are only applicable to those instances of LSAPE that respect the triangle inequality [START_REF] Serratosa | Fast computation of bipartite graph matching[END_REF][START_REF] Serratosa | Speeding up fast bipartite graph matching through a new cost matrix[END_REF][START_REF] Serratosa | Computation of graph edit distance: Reasoning about optimality and speed-up[END_REF].

The fastest currently available LSAPE solver is the algorithm FBP which reduces an instance C ∈ R (m+1)×(n+1) of LSAPE to an instance C ∈ R n×m of LSAP. However, FBP requires C to respect the triangle inequality. In this paper, we propose the algorithm FLWC, a fast solver for LSAPE without cost constraints. Like FBP, FLWC reduces C to a (n × m)-sized instance of LSAP. Unlike FBP, FLWC does not require C to respect the triangle inequality. FLWC builds upon the insight that node substitutions whose costs do not respect the triangle inequality can be factorized into removals and insertions. A thorough experimental evaluation shows that, for instances that do respect the triangle inequality, FLWC is as fast as FBP, while, for unconstrained instances, FLWC clearly outperforms all competitors.

The remainder of the paper is organized as follows: In Section 2, we formally introduce LSAPE and LSAP. In Section 3, we discuss related work. In Section 4, we present FLWC. In Section 5, we empirically evaluate the performance of FLWC w. r. t. the performance of existing competitors.

Preliminaries

Linear Sum Assignment Problem

The linear sum assignment problem (LSAP) is defined on complete bipartite graphs K n,m = (U, V, U × V), where U = {u i } i∈I and V = {v j } j∈J are sets of nodes, and I = {1, . . . , n} and J = {1, . . . , m} are index sets. In the following, we will assume w. l. o. g. that n ≤ m. This can easily be enforced by exchanging the roles of U and V, if necessary. An edge set M ⊂ U × V is called a maximum matching for K n,m if and only if all nodes in U are incident with exactly one edge in M and all nodes in V are incident with at most one edge in M. A maximum matching M can be encoded with a binary matrix X = (x i,j ) ∈ {0, 1} n×m by setting x i,j = 1 just in case (u i , v j ) ∈ M. In the following, we will always represent matchings with their induced binary matrices. The set Π n,m of all maximum matchings for K n,m is hence given as follows:

Π

n,m = X ∈ {0, 1} n×m : ∀j ∈ J, n i=1 x i,j ≤ 1 ∀i ∈ I, m j=1 x i,j = 1 (1)
When n = m, maximum matchings become perfect matchings, and Π n,n becomes the set of all n × n permutation matrices.

Given a cost matrix C = (c i,j ) ∈ R n×m for the edges of K n,m , LSAP asks to compute a maximum matching for K n,m which minimizes the induced edge costs.

Problem 1 (LSAP) Given a cost matrix C ∈ R n×m , LSAP consists to find a maximum matching X ∈ Π n,m with X ∈ argmin X∈Πn,m L(X, C), where L(X, C) = n i=1 m j=1 c i,j x i,j . This is a classical combinatorial optimization problem, also known as the minimum cost maximum (or perfect) matching problem in bipartite graphs. It can be solved in polynomial time and space complexities with several algorithms. For instance, if the cost matrix C is balanced, i. e., if n = m, it can be solved in O(n 3 ) time and O(n 2 ) space with the Kuhn-Munkres Hungarian Algorithm [START_REF] Harold W Kuhn | Variants of the hungarian method for the assignment problem[END_REF][START_REF] Munkres | Algorithms for the assignment and transportation problems[END_REF]. For the unbalanced case, it can be solved in O(n 2 m) time [START_REF] Bourgeois | An extension of the Munkres algorithm for the assignment problem to rectangular matrices[END_REF]. See [START_REF] Lawler | Combinatorial Optimization: Networks and Matroids[END_REF][START_REF] Burkard | Assignment Problems[END_REF] for more details.

Linear Sum Assignment Problem with Error-Correction

Given a maximum matching X for K n,m = (U, V, U × V) and a cost matrix C ∈ R n×m , a cell x i,j of X with x i,j = 1 can be interpreted as a substitution of the node u i ∈ U by the node v j ∈ V with associated substitution cost c i,j . LSAP can hence be viewed as the task to substitute all the nodes of U by pairwise distinct nodes of V such that the substitution cost is minimized. Note that, under this interpretation, LSAP does not require all nodes in V to be taken care of, as, if n < m, there are always nodes in V that do not substitute any node in U.

There are scenarios such as approximation of GED, where one is faced with a slightly different matching problem: First, in addition to node substitutions, one also wants to allow node insertions and node removals. Second, one wants to enforce that all the nodes in V are taken care of. The linear sum assignment problem with error-correction (LSAPE) models these settings. To this purpose, the sets U and V are extended to U = U ∪ { } and V = V ∪ { }, where is a dummy node. An edge set M = S ∪ R ∪ I ⊂ U × V is called error-correcting matching for K n,m, = (U , V , U ×V ) if and only if it does not contain the edge ( , ) and all nodes in U and V are incident with exactly one edge in M. The set S ⊂ U × V contains all node substitutions, R ⊂ U × { } contains all removals, and I ⊂ { } × V contains all insertions. The set Π n,m, of all error-correcting matchings for K n,m, is hence given as

Π n,m, = X = X sub x rem x ins 0 ∈ {0, 1} (n+1)×(m+1) : ∀j ∈ J, x j + n i=1 x i,j = 1 ∀i ∈ I, x i + m j=1 x i,j = 1 , (2)
where X sub = (x i,j ) ∈ R n×m encodes node substitutions, x rem = (x i, ) ∈ R n×1 encodes node removals, and x ins = (x ,j ) ∈ R 1×m encodes node insertions. Assume now that C ∈ R (n+1)×(m+1) is a cost matrix of the form

C = C sub c rem c ins 0 , (3)
where C sub = (c i,j ) ∈ R n×m contains the substitution costs, c rem = (c i, ) ∈ R n×1 contains the removal costs, and c ins = (c ,j ) ∈ R 1×m contains the insertion costs. Then LSAPE asks to compute an error-correcting matching for K n,m, such that the sum of the induced substitution, removal, and insertion costs is minimized.

Problem 2 (LSAPE) Given a cost matrix C ∈ R (n+1)×(m+1) of form (3), LSAPE consists to find an error-correcting matching

X ∈ Π n,m, with X ∈ argmin X∈Πn,m, L(X, C), where L(X, C) = n i=1 m j=1 c i,j x i,j + n i=1 c i, x i, + m j=1 c ,j x ,j .
Table 1 summarizes the notations introduced in this section. The standard algorithm extended bipartite matching (EBP) [START_REF] Riesen | Structural Pattern Recognition with Graph Edit Distance. Advances in Computer Vision and Pattern Recognition[END_REF][START_REF] Riesen | Approximate graph edit distance computation by means of bipartite graph matching[END_REF][START_REF] Riesen | Bipartite graph matching for computing the edit distance of graphs[END_REF] solves LSAPE by transforming it to LSAP. Given an instance C of LSAPE, EBP constructs an instance C ∈ R (n+m)×(m+n) of LSAP as

C = C sub ω(1 n×n -I n ) + diag(c rem ) ω(1 m×m -I m ) + diag(c ins ) 0 m×n ,
where ω denotes a very large value and the operator diag maps a vector (v i ) k i=1 to the diagonal matrix (d i,j ) k i,j=1 with d i,i = v i and d i,j = 0 for all i = j. EBP then calls a LSAP solver to compute an optimal maximum matching X ∈ Π n+m,m+n for LSAP. By construction, X is of the form

X = X sub X rem X ins X comp , (4)
where X sub ∈ {0, 1} n×m , X rem ∈ {0, 1} n×n , and X ins ∈ {0, 1} m×m . The southeast quadrant X comp ∈ {0, 1} m×n contains assignments from dummy nodes to dummy nodes which are not needed for encoding node removals or insertions but are anyway computed by LSAP algorithms. X is then transformed into an optimal error-correcting matching X ∈ Π n,m, for LSAPE with L(X , C) = L(X , C) by setting X sub = X sub , x rem = X rem 1 n , and x ins = 1 T m X ins . The time complexity of EBP is dominated by the complexity of solving the LSAP instance C. Therefore, EBP runs in O((n + m) 3 ) time.

Cost-Constrained Reductions to LSAP

The algorithms fast bipartite matching (FBP) [START_REF] Serratosa | Fast computation of bipartite graph matching[END_REF] and square fast bipartite matching (SFBP) [START_REF] Serratosa | Computation of graph edit distance: Reasoning about optimality and speed-up[END_REF][START_REF] Serratosa | Speeding up fast bipartite graph matching through a new cost matrix[END_REF] build upon more compact reductions of LSAPE to LSAP than EBP. However, both FBP and SFBP are applicable only to those instances C of LSAPE which respect the following triangle inequalities:

∀(i, j) ∈ I × J, c i,j ≤ c i, + c ,j (5)
In other terms, FBP and SFBP can be used for instances of LSAPE where substituting a node u i ∈ U by a node v j ∈ V is never more expensive than removing u i and inserting v j . The following proposition is the key-ingredient of both FBP and SFBP. Recall that we have assumed w. l. o. g. that n ≤ m.

Proposition 1 (Cf. [START_REF] Serratosa | Fast computation of bipartite graph matching[END_REF][START_REF] Serratosa | Computation of graph edit distance: Reasoning about optimality and speed-up[END_REF][START_REF] Serratosa | Speeding up fast bipartite graph matching through a new cost matrix[END_REF]) Let X ∈ Π n,m, be an optimal error-correcting matching for an instance C ∈ R (n+1)×(m+1) of LSAPE which satisfies [START_REF] Bougleux | Graph edit distance as a quadratic assignment problem[END_REF]. Then X contains no node removal, i. e., satisfies x rem = 0 n . Note that this implies that X contains exactly m -n node insertions and hence that

x ins = 0 m , if n = m.
Given an instance C of LSAPE that satisfies [START_REF] Bougleux | Graph edit distance as a quadratic assignment problem[END_REF], FBP constructs an instance

C ∈ R n×m of LSAP by setting C = C sub -c rem 1 T n -1 m c ins .
Subsequently, FBP computes an optimal maximum matching X ∈ Π n,m for C. Because of Proposition 1, it then holds that the matrix X ∈ Π n,m, defined by X sub = X , x rem = 0 n , and

x ins = 1 T m -1 T n X is an optimal error-correcting matching for C.
The variant FBP 0 of FBP transforms C into a balanced instance C 0 ∈ R m×m of LSAP, by adding the matrix 0 m-n,n below C defined for FBP. After solving this instance, FBP 0 transforms the resulting optimal maximum matching X for C 0 into an optimal error-correcting matching X ∈ Π n,m, for C, by applying the same transformation rules as FBP on X .

Like FBP 0 , SFBP transforms an instance C into a balanced instance C ∈ R m×m of LSAP. However, C is now defined as follows:

C = C sub 1 m-n c ins
In the next step, SFBP computes an optimal maximum matching X = X sub X ins for C. SFBP then constructs the matrix X ∈ Π n,m, by setting X sub = X sub , x rem = 0 n , and x ins = 1 T m X ins . Again, Proposition 1 ensures that X is indeed an optimal error-correcting matching for LSAPE.

The time complexities of FBP, FBP 0 , and SFBP are dominated by the complexities of solving the LSAP instances C. Therefore, FBP runs in O(n 2 m) time, while FBP 0 and SFBP run in O(m 3 ) time. These are significant improvements over EBP. However, recall that FBP, FBP 0 , and SFBP can be used only if the cost matrix C respects the triangle inequalities [START_REF] Bougleux | Graph edit distance as a quadratic assignment problem[END_REF].

Adaptations of Classical Algorithms LSAPE can also be solved directly by adapting algorithms originally designed for LSAP. An adaptation of the Jonker-Volgenant Algorithm is proposed in [START_REF] Jones | Optimising the Volgenant-Jonker algorithm for approximating graph edit distance[END_REF]. An adaption of the Hungarian Algorithm, denoted HNG in this paper, has been suggested in [START_REF] Bougleux | A hungarian algorithm for errorcorrecting graph matching[END_REF]. Both modifications lead to an overall time complexity of O(n 2 m).

A Compact Reduction from LSAPE to LSAPE without Cost Constraints

The main contribution of this paper is to show that LSAPE without cost constraints can be reduced to an instance of LSAP of size n × m. The reduction translates into the algorithm FLWC (fast solver for LSAPE without cost constraints), which, like FBP, runs in O(n 2 m) time, but, unlike FBP, FBP 0 , and SFBP, does not assume the costs to respect the triangle inequalities [START_REF] Bougleux | Graph edit distance as a quadratic assignment problem[END_REF]. The following Theorem 1 states the reduction principle. It relies on a cost-dependent factorization of substitutions, removals and insertions. In Section 4.1, we discuss special cases and present FLWC. In Section 4.2, we present the proof of Theorem 1.

Theorem 1 (Reduction Principle) Let C ∈ R (n+1)×(m+1) be an instance of LSAPE and let X ∈ Π n,m be an optimal maximum matching for the instance C of LSAP defined by

∀(i, j) ∈ I × J, c i,j = δ C i,j, c i,j + (1 -δ C i,j, )(c i, + c ,j ) -δ n<m c ,j , (6)
where

δ n<m = 1 if n < m 0 else , δ C i,j, = 1 if c i,j ≤ c i, + c ,j 0 else
Furthermore, let the function f C : Π n,m → Π n,m, be defined as follows: For instance, we have

∀(i, j) ∈ I × J, x i,j = f C (X) i,j = δ C i,j, x i,j (7)
∀i ∈ I, x i, = f C (X) i, = 1 - m j=1 δ C i,j, x i,j (8)
∀j ∈ J, x ,j = f C (X) ,j = 1 - n i=1 δ C i,j, x i,j (9)
Then X = f C (X ) is an optimal error-correcting matching for C with cost L(X , C) = L(X , C) + δ n<m m j=1 c ,j . Example 1 Assume that n = 2, m = 3,
C =   i j 1 2 3 1 3 5 1 4 2 8 9 4 4 2 4 0 0   apply (6) =====⇒ C = i j 1 2 3 1 1 1 1 2 4 4 4 u 1 u 2 v 1 v 2 v 3 1 4 (a) X u 1 u 2 v 1 v 2 v 3 3 4 4 0 (b) X u 1 u 2 v 1 v 2 v 3 9 4 4 (c) u 1 u 2 v 1 v 2 v 3 3 4 2 (d)
c 1,1 = δ C 1,1 3 + (1 -δ C 1,1 )(2 + 4) -δ 2<3 2 = 3 -2 = 1 and c 2,2 = δ C 2,2 9 + (1 -δ C 2,2 )(4 + 4) -δ 2<3 4 = 8 -4 = 4.
Fig. 1(a) shows an optimal maximum matching X ∈ Π n,m for C, and Fig. 1(b) shows an optimal error-correcting matching X = f C (X ) ∈ Π n,m, for C. Note that f C factorizes the substitution (u 2 , v 2 ) into the removal (u 2 , ) and the insertion ( , v 2 ), since c 2,2 > c 2, + c ,2 (Fig. 1(c)). On the other hand, the substitution

(u 1 , v 1 ) is not factorized, since c 1,1 ≤ c 1, + c ,1 (Fig. 1(d)). Furthermore, we have L(X , C) = 11, L(X , C) = 5, and δ n<m m j=1 c ,j = 6. Therefore, it holds that L(X , C) = L(X , C) + δ n<m m j=1 c ,j
, as stated by our reduction principle.

Discussion of Special Cases and Presentation of FLWC

Theorem 1 states that a general instance C ∈ R (n+1)×(m+1) of LSAPE, which is not required to respect the triangle inequalities [START_REF] Bougleux | Graph edit distance as a quadratic assignment problem[END_REF], can be reduced to a (n×m)sized instance of LSAP. However, there is still room for improvement if C does respect the triangle inequalities.

Proposition 2 Let C ∈ R (n+1)×(m+1) be an instance of LSAPE that respects the triangle inequalities [START_REF] Bougleux | Graph edit distance as a quadratic assignment problem[END_REF]. Then the reduction principle specified in Theorem 1 can be carried out without knowledge of the removal costs c rem . If n = m holds, too, then knowledge of the insertion costs c ins is not necessary, either.

Proof 1 The proposition immediately follows from the facts that δ C i,j, = 1 holds for all (i, j) ∈ I × J if C respects [START_REF] Bougleux | Graph edit distance as a quadratic assignment problem[END_REF], and that

δ n<m = 0 if n = m.2.
Proposition 2 is useful, because in many application scenarios for LSAPE, C is not given but has to be computed. Furthermore, C is often known a priori to respect the triangle inequalities, for instance, because its entries contain the distances between elements in a metric space. In such settings, the overall performance of an algorithm that calls a LSAPE solver as a subroutine can improve significantly if only parts of the cost matrix C have to be computed. Table 2 summarizes which parts of C have to be known for our reduction from LSAPE to LSAP. The case n > m can straightforwardly be obtained from the case n < m by transposing C. Recall that both FBP and FBP 0 always require the entire cost matrix C to be known. SFBP requires the same parts of C as our approach, but reduces LSAPE to a larger instance of LSAP (max{n, m} × max{n, m} vs. n × m).

n = m C sub C sub , cins, crem n < m C sub , cins C sub , cins, crem n > m C sub , crem C sub , cins, crem
Algorithm 1 shows the algorithm FLWC, which turns our reduction from LSAPE to LSAP into a method for computing an optimal error-correcting matching. FLWC only uses those parts of C which are really required by the reduction (cf. Table 2).

Algorithm 1: FLWC input: an instance C ∈ R (n+1)×(m+1) of LSAPE output: an optimal error-correcting matching X ∈ Πn,m, for C 1 if n > m then 2 C ← C T ; (n, m) ← (m, n); 3 initialize C ∈ R n×m ; 4 for i ∈ {1, . . . , n} do 5 for j ∈ {1, . . . , m} do 6 if C respects triangle inequalities then 7 if n = m then 8 ci,j ← ci,j; 9 else 10 ci,j ← ci,j -c ,j ; 11 else 12 ci,j ← δ C i,j, ci,j + (1 -δ C i,j, )(ci, + c ,j ) -δn<m c ,j ;
13 call LSAP solver to compute opt. max. matching X ∈ Πn,m for C;

14 X ← f C (X ); 15 if C was transposed in lines 1-2 then 16 X ← X T ; 17 return X ;
If the adaption of the Hungarian Algorithm to unbalanced instances of LSAP is used in line 13, FLWC runs in O(min{n, m} 2 max{n, m}) time and O(nm) space. Table 3 compares FLWC's time and space complexities to the complexities of existing competitors. Note that our reduction principle can also be used for the fast computation of a suboptimal solution for LSAPE. To this end, it suffices to replace the optimal LSAP solver in line 13 of Algorithm 1 by a suboptimal one such as one of the the greedy heuristics suggested in [START_REF] Riesen | Approximation of graph edit distance in quadratic time[END_REF].

Proving the Correctness of the Reduction Principle

The first step towards the proof is the following Proposition 3, which constitutes a relation between error-correcting matchings and maximum matchings.

Proposition 3 Let X ∈ Π n,m, be an error-correcting matching. Furthermore, let Z X = (U, V, {(u i , v j ) ∈ U × V : x i, x ,j = 1}) be the bipartite graph between U and V whose edges encode all combinations of node removals and insertions, let Z be the set of of maximum matchings for Z X , and let the set Y X be defined as follows:

Y X = {X sub + Z : Z ∈ Z } (10)
Then Y ∈ Π n,m holds for each Y ∈ Y X .

Proof 2 Let I sub = {i ∈ I : m j=1 x i,j = 1} and J sub = {j ∈ J : n i=1 x i,j = 1} be the set of indices of those nodes of U and V that are substituted by X. Furthermore, let s = |I sub | be the number of substitutions encoded by X, and let Z be a maximum matching for Z X . We observe that Z X can be viewed as the complete bipartite graph between the nodes U rem = {u i : i ∈ I \ I sub } and V ins = {v j : j ∈ J \ J sub } that are removed and inserted by X, and that we have

|U rem | = n -s ≤ m -s = |V ins |.
These observations imply that we have j∈J z i,j = 0 for each i ∈ I sub and j∈J z i,j = 1 for each i ∈ I \I sub . Similarly, we have i∈I z i,j = 0 for each j ∈ J sub and i∈I z i,j ≤ 1 for each j ∈ J \ J sub . This gives us j∈J x i,j + z i,j = 1 for each i ∈ I and i∈I x i,j + z i,j ≤ 1 for each j ∈ J, which implies

X sub + Z ∈ Π n,m .2
Example 2 Consider the error-correcting matching X shown in Fig 2(a). Since X removes u 2 and inserts v 2 and v 3 , Z X contains the edges (u 2 , v 2 ) and (u 2 , v 3 ) (Fig 2(b)). There are exactly two maximum matchings for Z X , namely 2(c)) and Y 2 = {(u 1 , v 1 ), (u 2 , v 3 )} (Fig. 2(d)).

Z 1 = u 1 u 2 v 1 v 2 v 3 (a) X u 1 u 2 v 1 v 2 v 3 (b) Z X u 1 u 2 v 1 v 2 v 3 (c) Y 1 u 1 u 2 v 1 v 2 v 3 (d) Y 2
{(u 2 , v 2 )} and Z 2 = {(u 2 , v 3 )}. This implies that Y X = {Y 1 , Y 2 } with Y 1 = {(u 1 , v 1 ), (u 2 , v 2 )} (Fig.
We now introduce the notion of a minimally-sized error-correcting matching. To this purpose, we call two error-correcting matchings X, X ∈ Π n,m, equivalent w. r. t. an instance C of LSAPE (in symbols: X ∼ C X ) if and only if, for all (i, j) ∈ I × J, x i,j = x i,j or c i,j = c i, + c ,j and x i,j = x i, x ,j or x i,j = x i, x ,j . By definition of ∼ C , the cost L is invariant on the equivalence classes induced by ∼ C . Definition 1 (Minimally-Sized Error-Correcting Matching) Let C ∈ R (n+1)×(m+1) be an instance of LSAPE and X ∈ Π n,m, be an error-correcting matching. Then X is called minimally-sized if and only if | supp(X)| < | supp(X )| holds for all X ∈ [X] ∼ C , where supp(X) is the support of X.

In other words, X is minimally-sized just in case it always favours substitution over removal plus insertion, if the costs are the same. By construction, each equivalence class [X] ∼ C contains exactly one minimally-sized errorcorrecting matching. In particular, there is always a minimally-sized optimal error-correcting matching.

The next step is to characterize a subset Π n,m, (C) ⊆ Π n,m, which contains all optimal minimally-sized error-correcting matchings for a given instance C of LSAPE.

Proposition 4 Let C ∈ R (n+1)×(m+1) be an instance of LSAPE, and let the set Π n,m, (C) ⊆ Π n,m, of cost-dependent error-correcting matchings be defined as follows:

Π n,m, (C) = { X ∈ Π n,m, : ∀(i, j) ∈ I × J, (1 -δ C i,j, )x i,j = 0, ∀(i, j) ∈ I × J, δ C i,j, x i, x ,j = 0}
Then Π n,m, (C) contains all minimally-sized optimal error-correcting matchings for the instance C of LSAPE.

Proof 3 Assume that there is a minimally-sized optimal error-correcting matching X ∈ Π n,m, \ Π n,m, (C). Then there is a pair (i, j) ∈ I × J such that (1 -δ C i,j, )x i,j = 1 or δ C i,j, x i, x ,j = 1. Assume that we are in the first case, i. e., that δ C i,j, = 0 and x i,j = 1. This implies c i,j > c i, + c ,j . Now consider the error-correcting matching X , which, instead of substituting u i by v j , removes u i and inserts v j (x i, x ,j = 1). Since δ C i,j, = 0, X is cheaper than X . This contradicts X 's optimality. If we are in the second case, we have x i, x ,j = 1 and c i,j ≤ c i, + c ,j . Since X is minimally-sized, we can strengthen the last inequality to c i,j < c i, + c ,j . Consider the error-correcting matching X , which, instead of removing u i and inserting v j , substitutes u i by v j (x i,j = 1). Again, X is cheaper than X , which is a contradiction to X 's optimality.2

The following Proposition 5 shows that the transformation function f C defined in Theorem 1 indeed maps maximum matchings to error-correcting matchings and that it is surjective on Π n,m, .

Proposition 5 Let f C : Π n,m → Π n,m, be defined as in Theorem 1. Then it holds that img(f C ) ⊆ Π n,m, and that img(f C ) ⊇ Π n,m, (C).

Proof 4 Consider a maximum matching X ∈ Π n,m and let X = f C (X). From ( 7) and (8), we have x i, + m j=1 x i,j = 1 for each i ∈ I. From ( 7) and ( 9), we have x ,j + n i=1 x i,j = 1 for each j ∈ J. This implies X ∈ Π n,m, and thus img(f C ) ⊆ Π n,m, .

For showing img(f C ) ⊇ Π n,m, (C), we fix an error-correcting matching X ∈ Π n,m, (C). From Proposition 3, there is a set Y X of maximum matchings representing X. Consider a matching X ∈ Y X , i. e. X = X sub + Z with Z ∈ Z X We will show that X = f C (X), which proves the proposition. To this end, we first show the following equalities:

∀(i, j) ∈ I × J, δ C i,j, z i,j = 0 ( 11
)
Consider a pair (i, j) ∈ I × J with z i,j = 1. From Z X ∈ Z X , we know that x i, x ,j = 1. As X ∈ Π n,m, (C), this implies δ C i,j, = 0 and hence proves [START_REF] Ferrer | A first step towards exact graph edit distance using bipartite graph matching[END_REF]. Now let (i, j) ∈ I × J. It holds that f C (X) i,j = δ C i,j, x i,j + δ C i,j, z i,j = δ C i,j, x i,j = x i,j , where the first equality follows from the definitions of f C and X, the second equality follows from [START_REF] Ferrer | A first step towards exact graph edit distance using bipartite graph matching[END_REF], and the third equality follows from X ∈ Π n,m, (C). We have hence shown that

X sub = f C (X) sub . Next, we show that f C (X) ins = x ins . Let j ∈ J. We have f C (X) ,j = 1 - n i=1 δ C i,j, x i,j - n i=1 δ C i,j, z i,j = 1 - n i=1 δ C i,j, x i,j = 1 - n i=1
x i,j = x ,j , as required. Again, the first equality follows from the definitions of f C and X, the second equality follows from [START_REF] Ferrer | A first step towards exact graph edit distance using bipartite graph matching[END_REF], and the third equality follows from X ∈ Π n,m, (C). The last equality follows from the fact that X is an error-correcting matching. The argument for showing that f C (X) rem = x rem is analogous.2 Now we prove the correctness of our reduction principle.

Proof of Theorem 1 1 Let X ∈ Π n,m be a maximum matching. Its cost w. r. t. C is given by:

L(X, C) = n i=1 m j=1 (δ C i,j, c i,j + (1 -δ C i,j, )(c i, + c ,j ) -δ n<m c ,j )x i,j = n i=1 m j=1 c i,j δ C i,j, x i,j + (c i, + c ,j )(1 -δ C i,j, )x i,j -δ n<m m j=1 c ,j n i=1 x i,j + m j=1 c ,j - m j=1 c ,j = n i=1 m j=1 c i,j (δ C i,j, x i,j ) =f C (X)i,j + n i=1 c i,   m j=1 (1 -δ C i,j, )x i,j   =Ai + m j=1 c ,j δ n<m + n i=1 (1 -δ C i,j, -δ n<m )x i,j =Bj -δ n<m m j=1 c ,j
Since X is a maximum matching for Π n,m , we know that

A i = 1- m j=1 δ C i,j, x i,j = f C (X) i,
for each i ∈ I. We now distinguish the cases δ n<m = 1 and δ n<m = 0. In the first case, we immediately have

B j = 1 - n i=1 δ C i,j, x i,j = f C (X)
,j for each j ∈ J. In the second case, we have n = m and B j = f C (X) ,j holds, too, since, for balanced instances, a maximum matching contains an edge (u i , v j ) for each j ∈ J. We have thus shown the following equality:

∀X ∈ Π n,m , L(X, C) = L(f C (X), C) -δ n<m m j=1 c ,j (12)
Now let X be a minimally-sized optimal error-correcting matching for C, X be an optimal maximum matching for C, and X = f C (X ). From ( 12), we know that L(X , C) = L(X , C) + δ n<m m j=1 c ,j . Therefore, the theorem follows if we can show that L(X , C) = L(X , C). The ≥-direction of the desired equality follows from X 's optimality and the fact that, from Proposition 5, we know that X ∈ Π n,m, . For showing the ≤-direction, we pick an X ∈ Π n,m with f C (X ) = X . Such an X exists because, from Proposition 4, we know that X ∈ Π n,m, (C), and, from Proposition 5, we have img(f C ) ⊇ Π n,m, (C). Assume that L(X , C) > L(X , C). Then [START_REF] Gaüzère | Approximate graph edit distance guided by bipartite matching of bags of walks[END_REF] implies that L(X , C) > L(X , C), which contradicts X 's optimality.2

Experimental Validation

In our experiments, we compared FLWC to all existing competitors mentioned in Table 3. We evaluated the runtime of the different methods (Section 5.1) and their performances when used within approximate approaches for the computation of GED (Section 5.2). All procedures use the same C++ implementation of the Hungarian Algorithm as LSAP solver, which is based on the version presented in [START_REF] Lawler | Combinatorial Optimization: Networks and Matroids[END_REF][START_REF] Burkard | Assignment Problems[END_REF]. HNG is also based on this version, so that all procedures are comparable. 1

Time Comparison

To illustrate the differences between the methods, we recorded their execution time on three types of instances, some of which do and some of which do not respect the triangle inequalities: Machol-Wien instances, flat instances, and random instances. Experiments on further types yielded similar results. Fig . 3 shows the results for Machol-Wien instances and flat instances, Fig. 4 displays the results for random instances. For each instance type, the tested methods' execution time is reported for several values of n and m. In the first and third columns of the figures, n is fixed and m is varied, while in the second and fourth columns, m is fixed and n is varied. In any case, we have n ≤ m. The general methods FLWC, EBP, and HNG were tested on all instances, the cost-constrained algorithms FBP, FBP 0 , and SFBP only on those that respect the triangle inequalities. Since FBP 0 was slower than FBP across all instances, our plots do not contain curves for FBP 0 .

Machol-Wien Instances [START_REF] Machol | A 'hard' assignment problem[END_REF] Machol-Wien instances are defined by c i,j = i * j ∀(i, j) and thus satisfy the triangle inequalities. They are known to be hard to optimize for classical LSAP solvers such as the Hungarian Algorithm. This is also the case for the LSAPE solver HNG , but not for the other methods. This is explained by the fact that all methods except for HNG transform an instance C of LSAPE into an instance C of LSAP, which reduces the difficulty of the problem. The first row of Fig. 3 shows the results of our experiments on Machol-Wien instances. We observe that FLWC and FBP perform very similarly, provide the best results in all cases, and are more stable than the other methods. EBP and SFBP are more time consuming, in particular when n and m increases.

Flat Instances

In a first series of experiments, we considered flat instances of the form C = α1 n×m , with α = 10. These instances satisfy the triangle inequalities. Moreover, they are easy to solve, which implies that a large part of the execution time is spent on the initialization step of the Hungarian Algorithm and the cost transformations. The results for instances of this kind are displayed in the second row of Fig. 3. As before, FLWC and FBP are more stable and efficient than the other methods. In a second series of experiments, we varied the construction of C in order to enforce that k% of the nodes contained in the smaller set be removed. This can be achieved by randomly selecting m -n nodes in the larger set and setting their insertion cost to 0. Subsequently, k% of the remaining elements are selected and their insertion cost is set to (α/2) -1. Similarly, the removal costs of k% of the nodes contained in the smaller set are set to (α/2) -1. The resulting instances do not satisfy the triangle inequalities, and thus FBP, FBP 0 , and SFBP are not tested (they do not compute an optimal solution). The results are shown in the third and the fourth row of Fig. 3. Both FLWC and HNG clearly outperform EBP, and FLWC is slightly more stable than HNG .

Random instances

We also carried out experiments on random instances similar to the ones presented in [START_REF] Serratosa | Computation of graph edit distance: Reasoning about optimality and speed-up[END_REF][START_REF] Serratosa | Speeding up fast bipartite graph matching through a new cost matrix[END_REF]. These instances are very similar to the ones that occur in the context of approximation of GED. The cost c i,j encodes the cost of substituting a node v i and its set of incident edges in a graph G by a node v j and its set of incident edges in a graph H. an integer value in {1, . . . , l max }, where l max = √ n * m/10. Similarly, edges are labeled with a binary value. Then c i,j is defined as the cost of substituting the labels (0 if they are the same, or the constant c s else) plus the cost of an optimal error-correcting matching between the sets of incident edges. The cost of substituting two edges is defined as before, and the cost of removing or inserting an edge is defined by the constant c eri . Similarly, removal and insertion costs are defined as the cost of removing the node plus the cost of removing its incident edges, which is given by d i c eri + c nri , where d i is the node degree and c nri is its removal or insertion cost. The three following set of parameters (c s , c nri , c eri ) are considered: (40,[START_REF] Riesen | Structural Pattern Recognition with Graph Edit Distance. Advances in Computer Vision and Pattern Recognition[END_REF][START_REF] Riesen | Structural Pattern Recognition with Graph Edit Distance. Advances in Computer Vision and Pattern Recognition[END_REF], [START_REF] Riesen | Structural Pattern Recognition with Graph Edit Distance. Advances in Computer Vision and Pattern Recognition[END_REF]40,0), and (40, 2, 2). The first parameter setting (first row of Fig. 4) satisfies the triangle inequality. The second (second row of Fig. 4) and the third parameter setting (third row of Fig. 4) do not satisfy the triangle inequalities. However, for the second setting, none of the computed optimal error-correcting contains removals and insertions, and so all methods were tested. On the contrary, optimal error-correcting matchings for the third setting indeed contain up to 20% of removals, which is why we did not carry out tests for the cost-constrained methods FBP, FBP 0 , and SFBP. We again observe that our method FLWC is globally the most stable algorithm and obtains the best results on average.

Effect on Approximation of GED

As mentioned in the introduction, many methods approximate GED by lossily transforming the problem of its computation into an instance C of LSAPE [START_REF] Riesen | Structural Pattern Recognition with Graph Edit Distance. Advances in Computer Vision and Pattern Recognition[END_REF][START_REF] Gaüzère | Approximate graph edit distance guided by bipartite matching of bags of walks[END_REF][START_REF] Benjamin | Correcting and speeding-up bounds for non-uniform graph edit distance[END_REF][START_REF] Riesen | Approximate graph edit distance computation by means of bipartite graph matching[END_REF][START_REF] Zeng | Comparing stars: On approximating graph edit distance[END_REF][START_REF] Zheng | Efficient graph similarity search over large graph databases[END_REF][START_REF] Carletti | Approximate graph edit distance computation combining bipartite matching and exact neighborhood substructure distance[END_REF]. Lossy transformations from GED to LSAPE can be used for computing both upper and lower bounds for GED. Given an error-correcting matching for C computed by any LSAPE algorithm, an upper bound for GED is derived by computing the distance associated to the edit path induced by the considered matching. Since this edit path may be suboptimal, the computed distance might be an overestimation of the exact GED, and thus constitutes an upper bound. With this paradigm, each error-correcting matching for C yields a valid upper bound, although the ones induced by optimal matchings are usually tighter. Note that, if the LSAPE instance C constructed by a lossy transformation has only one optimal solution, each optimal LSAPE solver yields the same upper bound for GED. If C has more than one optimal solution, there might be differences in accuracy, since different LSAPE solvers might pick different optimal solutions depending on the organization of the input data. However, these differences in accuracy are arbitrary and hence disappear once the upper bounds are averaged across enough pairs of input graphs. In other words, when it comes to the approximation of GED, the only relevant property of an exact LSAPE solver is its runtime behaviour.

This observation is empirically confirmed by the experiments reported in Table 4, which shows how different methods for solving LSAPE affect the performance of the algorithm BP suggested in [START_REF] Riesen | Approximate graph edit distance computation by means of bipartite graph matching[END_REF], which computes an upper bound for GED. Given two graphs G 1 and G 2 on n 1 and n 2 nodes, respectively, BP constructs a LSAPE instance C ∈ R (n1+1)×(n2+1) . For computing the cell c i,j for (i, j) ∈ {1, . . . , n 1 } × {1, . . . , n 2 }, BP has to solve another LSAPE instance of size (deg G1 (u i ) + 1) × (deg G2 (v j ) + 1), where deg G1 (u i ) is the degree of node u i in G 1 and deg G2 (v j ) is the degree of node v j in G 2 . So altogether, BP has to solve 1 + n 1 n 2 instances of LSAPE. The tests were carried out on the datasets Acyclic and MAO from the ICPR GED contest [START_REF] Abu-Aisheh | Graph edit distance contest 2016: Results and future challenges[END_REF]. In order to ensure that also the cost-constrained LSAPE solvers compute optimal solutions, we defined metric edit costs by setting the cost of substituting nodes and edges to 1 and the cost of deleting and inserting nodes and edges to 3. We see that, as expected, there are no significant differences between the different solvers w. r. t. the tightness of the produced upper bounds. In terms of runtime, our solver FLWC performs best, followed by FBP and HNG .

In contrast to the situation for upper bounds, lossy transformations from GED to LSAPE which aim at the computation of lower bounds [START_REF] Benjamin | Correcting and speeding-up bounds for non-uniform graph edit distance[END_REF][START_REF] Blumenthal | Improved lower bounds for graph edit distance[END_REF] as well as methods based on conditional gradient descent [START_REF] Bougleux | Graph edit distance as a quadratic assignment problem[END_REF][START_REF] Bougleux | A hungarian algorithm for errorcorrecting graph matching[END_REF] where LSAPE occurs as a subproblem crucially depend on the optimality of the computed error-correcting matching. Therefore, these methods cannot use the existing fast LSAPE solvers FBP, FBP 0 , and SFBP, unless the triangle inequalities are known to be satisfied. For instance, lower bounds for GED are obtained from the cost L(X , C) of an optimal error-correcting matching for the LSAPE instance C. If X is not optimal, L(X , C) is not in general a valid lower bound. Table 5 shows how different methods for solving LSAPE affect the performance of the algorithm BRANCH suggested in [START_REF] Benjamin | Correcting and speeding-up bounds for non-uniform graph edit distance[END_REF][START_REF] Blumenthal | Improved lower bounds for graph edit distance[END_REF], which computes a lower bound for GED. Like BP, given to graphs G 1 and G 2 on n and n 2 nodes, respectively, BRANCH has to solve 1 + n 1 n 2 instances of LSAPE. One of them is of size (n 1 + 1) × (n 2 + 1), the other ones are of size (deg G1 (u i ) + 1) × (deg G2 (v j ) + 1). The tests were again carried out on the datasets Acyclic and MAO, but this time, edit costs that do not satisfy the triangle inequality were used (cf. [START_REF] Abu-Aisheh | Graph edit distance contest 2016: Results and future challenges[END_REF] for more details on the edit costs, especially equation ( 1) and setting 3 in Table 2). First of all, we can see that the classic approach EBP requires more computational time than other optimized approaches, and that FLWC, FBP, and HNG are the fastest methods. Second, we observe that the cost-constrained methods FBP, FBP 0 and SFBP are not able to deal with costs which do not satisfy the triangle inequality: Very often, they compute invalid lower bounds that exceed the exact GED, which we computed using a binary linear programming approach [START_REF] Lerouge | Exact graph edit distance computation using a binary linear program[END_REF]. This is explained by the fact that, if the triangle inequality is not satisfied, optimal error-correcting matchings might well include both insertions and removals. However, those error-correcting matchings are not considered by cost-constrained methods for LSAPE.

Figure 5 shows a case where the cost-constrained methods FBP, FBP 0 and SFBP compute an invalid lower bound. Considering the two graphs G 1 and G 2 extracted from Acyclic dataset, cost-constrained methods provide a matching ϕ FBP that favours node substitution over removal plus insertion even if removal plus insertion is cheaper. For instance, ϕ FBP matches the nodes 7 and 8 of G 1 to the nodes 9 and 5 of G 2 , although it is cheaper to first remove 7 and 8 and then insert 9 and 5, as done by the matching ϕ FLWC computed by FLWC. In conclusion, cost-constrained methods do not guarantee valid lower bounds for general edit costs, while our algorithm FLWC does.

Conclusions

In this paper, we presented FLWC, a new efficient method for solving the linear sum assignment problem with error correction. The technical backbone of our method is a cost-dependent factorization of substitutions into removals and insertions. FLWC runs in O(min{n, m} 2 max{n, m}) time and O(nm) space. Its complexities are hence the same as the ones of the most efficient state of the art methods FBP and HNG .

The advantage of our method FLWC over FBP is that, unlike FBP, FLWC remains valid for any configuration of the cost matrix and does not require the triangle inequality to holds. Therefore, our method allows to efficiently retrieve a lower bound for the graph edit distance irrespectively whether or not the edit costs respect the triangle inequality. This is especially important in settings where the edit costs are deduced from calculus and where one does not have an a priori guarantee that the triangle inequality will be satisfied. Situations of this kind occur, for instance, in some quadratic approximations of the graph edit distance problem [START_REF] Bougleux | Graph edit distance as a quadratic assignment problem[END_REF][START_REF] Bougleux | A hungarian algorithm for errorcorrecting graph matching[END_REF].

One advantage of FLWC over HNG is that, although both algorithms have the same time complexity, FLWC is slightly faster in practice. In particular, FLWC is more stable than HNG , in the sense that it also allows to quickly solve difficult LSAPE instances with whom HNG struggles. A second advantage is that FLWC is much easier to implement: While implementing HNG 's adaptation of the Hungarian Algorithm to LSAPE requires a thorough knowledge of matching theory, FLWC can be implemented by slightly transforming the cost matrix and then calling a solver for LSAP. Since LSAP is a very famous combinatorial optimization problem, libraries are available for all major programming languages.

For future work, we are planning to integrate the factorization of the assignment matrix, which is the core of our paper, in a modeling of the graph edit distance problem as a quadratic assignment problem.

  and consider the instance C of LSAPE and the induced instance C of LSAP:

Figure 1 :

 1 Figure 1: Illustration of the reduction principle stated in Theorem 1

Figure 2 :

 2 Figure 2: Illustration of Proposition 3

13

 13

Figure 3 :

 3 Figure 3: Time comparison for Machol-Wien instances (first row) and flat instances (all other rows).

  Graphs are constructed locally by assigning node degrees randomly from 1 to d max = (3/10)|V G |, where |V G | is the number of nodes in the graph G. Nodes are labeled randomly with

Figure 4 :

 4 Figure 4: Time comparison for random instances of LSAPE (see text).

Figure 5 :

 5 Figure 5: Example of lower bound computation where based methods do not allow to compute a valid lower bound.

Table 1 :

 1 Frequently used notations

	syntax	semantic
	C ∈ R (n+1)×(m+1) instance of LSAPE
	C ∈ R n×m	instance of LSAP
	X ∈ Πn,m,	error-correcting matching
	X ∈ Πn,m	maximum matching
	X ∈ Πn,m,	optimal error-correcting matching
	X ∈ Πn,m	optimal maximum matching
	L(X, C)	cost of error-correcting matching X w. r. t. C
	L(X, C)	cost of maximum matching X w. r. t. C
	3 State of the Art for LSAPE
	3.1 Unconstrained Reduction to LSAP

Table 2 :

 2 Required parts of C for our reduction from LSAPE to LSAP triangle ineqs. hold triangle ineqs. do not hold

Table 3 :

 3 Time and space complexities of existing algorithms for LSAPE under the assumptions that reductions to LSAP use the Hungarian Algorithm for solving LSAP O(min{n, m} 2 max{n, m}) O(nm) FLWC [our approach] O(min{n, m} 2 max{n, m}) O(nm)

	method	time	space
	cost-constrained methods
	FBP [27]	O(min{n, m} 2 max{n, m}) O(nm)
	FBP0 [27]	O(max{n, m} 3 )	O(max{n, m} 2 )
	SFBP [29]	O(max{n, m} 3 )	O(max{n, m} 2 )
	general methods
	EBP [21]	O((n + m) 3 )	O((n + m) 2 )
	HNG [4]

Table 4 :

 4 Effect of LSAPE methods on algorithm BP suggested in[START_REF] Riesen | Approximate graph edit distance computation by means of bipartite graph matching[END_REF] that computes an upper bound for GED on datasets with metric costs.

	method time in µs average UB time in µs average UB
		Acyclic		MAO
	cost-constrained methods
	FBP	1.61	38.86	6.62	107.97
	FBP0	2.04	38.89	8.10	107.85
	SFBP	2.14	39.28	11.10	107.48
	general methods
	EBP	4.84	39.14	23.30	107.48
	HNG	1.87	38.97	7.46	107.59
	FLWC	1.53	38.86	6.05	108.04

Table 5 :

 5 Effect of LSAPE methods on algorithm BRANCH suggested in[START_REF] Benjamin | Correcting and speeding-up bounds for non-uniform graph edit distance[END_REF][START_REF] Blumenthal | Improved lower bounds for graph edit distance[END_REF] that computes a lower bound for GED on datasets with non metric costs. method time in µs # invalid LB time in µs # invalid LB

		Acyclic		MAO
	cost-constrained methods
	FBP	1.28	4	8.30	416
	FBP0	1.64	4	10.80	416
	SFBP	1.91	4	12.81	416
	general methods
	EBP	5.34	0	21.10	0
	HNG	1.54	0	7.59	0
	FLWC	1.23	0	8.50	0

Procedures are written both in C++ and MATLAB (or GNU Octave). The source code is available at bougleux.users.greyc.fr/lsape/.


Younes Faydi
email: younes.faydi@ensam.eu

Loic Brancheriau

Guillaume Pot

Robert Collet

Prediction of Oak Wood Mechanical Properties Based on the Statistical Exploitation of Vibrational Response

Keywords: Wood mechanical properties, Vibration test, Partial least squares, Stepwise regression, Oak

In the European Union, timber is used in structural applications and must be graded with a Conformité Européene (CE) mark. To achieve standard, machine strength grading is used. A common technology for these machines is based on using the vibrational response of each wood board to estimate the timber modulus of elasticity and modulus of rupture. The first Eigen frequency is usually used to predict these mechanical properties. However, in heterogeneous wood species such as oak, this parameter is less correlated with mechanical properties. The current study proposes two new methods based on an extended exploitation of the vibrational response that predicts oak wood mechanical properties. The first method was based on the mechanical parameters deduced from several Eigen frequencies that were chosen with regards to a stepwise regression. The second method was based on the full vibrational spectrum and used a partial least squares method. The first method slightly improved the prediction of the modulus of elasticity compared with the first Eigen frequency in edgewise transversal vibration. Both methods significantly improved the prediction of the modulus of rupture.

INTRODUCTION

The characterization of timber mechanical properties is needed for structural applications of wood. However, wood heterogeneity requires the individual characterization of each board. European standards have defined different grades based on three characteristics (EN 338 (2016);EN 384 (2016)): density, modulus of elasticity (MOE), and modulus of rupture (MOR). The reference for these properties is based on a four point bending test according to the EN 408 (2012) standard. However, in the frame of machine strength grading described in EN 14081 (2016), these properties must be determined nondestructively. Different mechanical grading methods and machines are commercially available and based on technologies including X-rays, vibrational methods, or stress ratings [START_REF] Hanhijarvi | Potential of Strength Grading of Timber with Combined Measurement Techniques Report of the Combigrade project -Phase 1[END_REF][START_REF] Baillères | Strength and stiffness assessment technologies for improving grading effectiveness of radiata pine wood[END_REF]. Mechanical properties are typically predicted from the board's density or dynamic measurement of the MOE, and sometimes with the help of the detection of singularities like knots or fiber direction.

The board mass and dimensions allow for the measurement of the global board density. MOE can also be estimated by dynamic methods [START_REF] Hanhijarvi | Potential of Strength Grading of Timber with Combined Measurement Techniques Report of the Combigrade project -Phase 1[END_REF][START_REF] Baillères | Strength and stiffness assessment technologies for improving grading effectiveness of radiata pine wood[END_REF]. The most challenging characteristic to predict is the MOR. The efficiency of prediction is typically evaluated with the coefficient of determination (R²) between indicating properties and destructively obtained variables. However, it is also essential to estimate the mean error of prediction to evaluate each model. The coefficient of determination between the measured MOR and MOR indicating property obtained by non-destructive techniques can reach R² = 0.68 for Spruce and R² = 0.58 for Douglas fir [START_REF] Viguier | Improving strength grading of lumber by grain angle measurement and mechanical modeling[END_REF]. However, there is little data available in the literature concerning the prediction of hardwood mechanical properties and particularly oak wood. It could be mentioned, nonetheless, that the coefficient of determination between static MOE and MOR for chestnut was R² = 0.36 in [START_REF] Vega | Proposal for visual grading criteria of structural timber of sweet chestnut from Spain[END_REF], which is lower than spruce with R² = 0.71 in [START_REF] Viguier | Modelling mechanical properties of spruce and Douglas fir timber by means of X-Ray and grain angle measurements for strength grading purpose[END_REF]. Grading oak is a major difficulty for the large commercialization of structural applications [START_REF] Collet | The processing of small lowgrade French oaks into solid cross laminated panels[END_REF]. Thus, more effort is needed to improve its characterization and non-destructive grading.

This paper describes the use of vibrational analysis to predict MOE and MOR. The first Eigen frequency allows the prediction of the MOE [START_REF] Hanhijarvi | Potential of Strength Grading of Timber with Combined Measurement Techniques Report of the Combigrade project -Phase 1[END_REF][START_REF] Baillères | Strength and stiffness assessment technologies for improving grading effectiveness of radiata pine wood[END_REF]. Guindos andGuaita (2013, 2014) showed that MOR depends on knot characteristics including shape, position, etc. Detecting these characteristics is essential for MOR prediction. Roohania et al. (2015) observed the influence of the position and orientation of a local heterogeneity on the resonance frequencies by comparing clear beams and beams that were drilled. The behavior of the free flexural vibration was different depending on the position and orientation of the heterogeneity. When the drilling axis lies in the vibration bending plane, the weakening of frequency was maximal at the location of an antinode of vibration. However, the frequency offset was maximal in the vibration node when the drilling axis was orthogonal to the bending plane. This particular behavior could be used for grading lumber because it takes into account the presence of defects by the estimation of the static MOE and the MOR. For the same purpose, [START_REF] Sobue | Identification of defect position in a wooden beam from the power spectrum of longitudinal vibration[END_REF] detected the position of defects by looking at the resonance frequency shifts in longitudinal vibrating. Indeed, the shift was large when the defect coincided with a node of vibration, while no shift was observed when the defect coincided with an antinode of vibration. The proposed method worked only to identify a maximum of two dominant defects. Another study from [START_REF] Yang | Application of modal analysis by transfer function to nondestructive testing of wood I: Determination of localized defects in wood by the shape of the flexural vibration wave[END_REF] detected some defects such as knots by comparing clear and defected boards. The criterion used was the variation of the mode shape of transversal vibration waves. The theoretical and experimental mode shapes coincided with clear beams and differed near a defect. None of the previous studies went beyond heterogeneity detection to predict MOR from vibrational signals.

The current study proposed two analysis techniques based only on vibrational measurements to predict MOE and MOR. The first method was based on a stepwise regression and considered several output parameters from the vibration test. Stepwise regression detects the most representative parameters that explain most of the variance in MOE and MOR. The second technique was related to frequency spectrum standardization in frequency and in amplitude. In the second method, the partial least squares (PLS) statistical method was used, such that each frequency of the spectrum was a predictive variable. [START_REF] Brancheriau | Use of the partial least squares method with acoustic vibration spectra as a new grading technique for structural timber[END_REF] conducted a similar analysis for larch species, but they did not standardize the frequency vectors. Thus, their analysis was dependent on the board dimensions. The current study standardized the amplitude and the frequency, which resulted in the removal of the percussion impact effect and board dimensions effect. These two analysis techniques were compared to the more classical use of the dynamic MOE, based on the first Eigen frequency in longitudinal vibration as an indicating property of static MOE and MOR.

EXPERIMENTAL

Sampling

A total of 164 boards of French oak wood (Quercus petraea and Quercus robur) were tested in normal laboratory conditions of 20 °C and 65% humidity as noted in the EN 408 (2012) standard. The boards were selected according to the NF B52 001-1 (2011) visual grading standard to ensure an almost equal distribution between the three grades of D30, D24, and D18 (Table 1). The waste category was over-represented because the goal of machine grading is to reduce the underestimation of visual grading. Four different board dimensions were chosen (Table 2). The number of boards was different in each section because the wood board manufacturing process reduced the final number of exploitable boards.

Apparatus

Four points bending

A four point bending test was performed according to EN 408 (2012) using the distance between edge supports equals to 18 times the board's height, and a distance between the central loading heads equals to 6 times the board's height. These distances were modified according to each cross-section. The longitudinal position of the boards was chosen so that the supposed critical board's zone was between loading heads. The tension edge was selected at random. Global MOE and MOR were determined in accordance with the EN 408 (2012) standard and constituted the observations for statistical methods.

Vibrational measurements

Vibration tests were carried out using the BING device developed by CIRAD [START_REF] Brancheriau | Bing: Beam identification by nondestructive grading (Cirad)[END_REF]. The technique consisted of generating an impulse with a hammer in beam extremity and placing a microphone as a receiver on the other side. The beam was placed on two elastic supports to ensure free vibrations. These supports have been positioned at the theoretical node points at the first mode in bending. Three types of vibrational tests were done. An impact along the sample involved longitudinal vibrations representative of a compressive test (Fig. 1); bending or transversal vibrations were generated by edgewise and flatwise impacts (Fig. 1). The longitudinal vibration test represents industrial measurements of the dynamic MOE using the first Eigen frequency. The analog output signal was converted to a digital one with a sampling frequency of 39 kHz, acquisition of 1678 ms, and 12 bits resolution. A discrete Fourier transformation converted the temporal response to a frequency spectrum.

Method Based on Eigen Frequencies

Global output parameters

The relationship between Eigen frequencies and MOE or shear modulus depend on the vibration direction and the mechanical model applied. Several output parameters were computed from these different relationships and used later as predictive variables.

Equations 1, 2 and 3 were demonstrated in [START_REF] Brancheriau | Natural vibration analysis of clear wooden beams: Theoretical review[END_REF], who made a review and theoretical comparison of the different models.

In longitudinal vibrations, the MOE is computed thanks to the first Eigen frequency only, following Eq. 1,

𝑀𝑂𝐸 𝐶 = 4𝜌𝐿²𝑓 𝑚𝑜𝑑𝑒 1 (1)
where 𝑀𝑂𝐸 is the MOE in longitudinal vibrations (compression), 𝜌 is the density, L is the length, and 𝑓 is the frequency of the first mode in longitudinal vibration.

In transversal vibrations (edgewise and flatwise), four Eigen frequencies were chosen. In this case, two common theories were used to determine the MOE: Bernoulli (1748) and [START_REF] Timoshenko | On the correction for shear of the differential equation for transverse vibrations of prismatic bars[END_REF]. Bernoulli theory does not consider beam shearing. To compute the MOE according to Bernoulli theory, equation 2 was used,

𝑀𝑂𝐸𝐵 𝑖 = 4𝜋²𝜌𝑆𝐿 4 𝑓 2 𝑚𝑜𝑑𝑒 𝑖 𝐼𝑃 𝑖
(2) where 𝑀𝑂𝐸𝐵 is the MOE following Bernoulli theory for the i th mode in transversal vibrations; 𝜌 is the density; L is the length; 𝑓 is the frequency of the i th mode in transversal vibration; 𝑆 is the cross section; 𝐼 is the second moment of area; 𝑃 is a scalar relative to the i th mode (i th solution of cos √𝑃 cosh √𝑃 = 1). Applying Timoshenko bending theory, [START_REF] Bordonné | Module Dynamique et Frottement Intérieur dans le Bois, Mesures Sur poutres Flottantes en Vibrations Naturelles; Dynamic Modulus and Internal Friction in Wood, Measurements on Floating Beams in Natural Vibrations[END_REF] and more recently [START_REF] Brancheriau | Natural vibration analysis of clear wooden beams: Theoretical review[END_REF] gave the following solution of the equation of motion for a beam in transversal vibration,

MOET 𝜌 - MOET 𝐾G 𝑥 𝑖 = 𝑦 𝑖 (3)
where is the MOE following Timoshenko theory, 𝜌 is the density, 𝐾 is the shear factor (𝐾 = for a rectangular cross-section), G is the shear modulus, 𝑥 and 𝑦 are parameters that depend on the Eigen frequency [START_REF] Brancheriau | Natural vibration analysis of clear wooden beams: Theoretical review[END_REF] for detailed expression). Thanks to a linear regression between 𝑥 and 𝑦 parameters for the four selected Eigen frequencies, the MOE and shear modulus were found. The coefficient of determination of these regressions were also computed to be used as predictive variables ( , in Table 3).

Fig. 1. Illustration of vibrational measurements

As suggested by Roohania et al. (2015), a ratio between edgewise and flatwise vibrational MOE was computed for each mode. The inharmonicity factor, proposed by [START_REF] Aramaki | Sound quality assessment of wood for xylophone bars[END_REF], was computed. It represented the ratios between the second, third, and fourth Eigen frequencies divided by the first Eigen frequency. To take into account the viscoelastic behavior, the loss factor (tan δ) was determined for the first Eigen frequency [START_REF] Brancheriau | Internal friction measurement of tropical species by various acoustic methods[END_REF]. The loss factor could have been calculated for each Eigen frequency, but the first was the only one that was considered because of its high energy and low computational uncertainty. Finally, it was possible to extract 27 variables, as presented in Table 3. This list of parameters organized the input data for the stepwise regression method.

Stepwise regression

The stepwise regression [START_REF] Ing | A stepwise regression method and consistent model selection for high dimensional sparse linear models[END_REF] aimed to fit a model based on the vibrational parameters of Table 3 to predict the static MOE and MOR. The stepwise regression was conducted with R software [START_REF] Venables | An Introduction to R: Notes on R: A Programming Environment for Data Analysis and Graphics (Version 3.3[END_REF] with the Akaike information criterion method described by [START_REF] Venables | Modern Applied Statistics with S-Plus[END_REF].

Method Based on the Full Vibrational Spectrum

This analysis considers standardizing amplitudes of the whole spectrum for each vibration test as predictive variables. Thus, the spectrum was standardized in its frequencies and amplitudes. Note that the present method only used the edgewise vibrational measurements for reasons that are discussed later.

Standardization of spectrum frequencies

The determination of the fundamental frequency f 0 allowed standardized frequencies by computing the ratio of each frequency f to the fundamental frequency. The standardization made the analysis independent from the board's dimensions. Thus, each board had its own frequency scale (f/f 0 ). By using a linear interpolation, only one frequency scale was constituted for all measured spectra.

Standardization of spectrum amplitudes

Explicative variables were composed of amplitudes which were standardized to remove impact influence. The signal energy was defined as follows (Eq. 4), (4) where S(f) is the magnitude of the Fourier transform at the frequency f. In a linear system, the responses of the system under two different impacts were as follows, (5) where 𝑆 is the Fourier transform of the output signal, 𝐸 is the Fourier transform of the input signal (impact), and is the impulse response of the linear system. Assuming there is a linear relationship between 𝐸 1 and 𝐸 2 (α a constant): 𝐸 = 𝐸 , the output spectrum is then equal to: 𝑆 = 𝑆 and the signal energy: = . To remove the influence of impact, an appropriate standardization is to divide each spectrum amplitude by the square root of the energy:

(6)
Because this expression is non-dimensional, the MOE calculated from the first Eigen frequency extracted from edgewise tests was included in the final amplitude standardization

The explicative and explained (MOE and MOR) variables were, in turn, in the same unit.

Partial least squares

The partial least squares (PLS) were applied on treated spectral signals to predict MOE and MOR. The PLS was based on multiple linear regressions (MLR) represented by Eq. 7. MLR looks at a linear relationship between explicative variables (elements of the spectrum) to explain a single dependent goal variable (MOE or MOR) according to Eq. 7, which is an example for MOR, (7) where is the i th observation of goal variable, n is the number of observations, m is the number of explicative variables, 𝑥 is the j th standardized amplitude of spectrum in the i th observation, is the j th element of regression coefficient and is the i th residual element.

The goal was to avoid the collinearities between the variables by compressing the information and maximizing the explained variance of the dependent variable (MOE or MOR). Equation 7was then converted to Eq. 8. The details of this method are available in [START_REF] Naes | A User-Friendly Guide to Multivariate Calibration and Classification[END_REF]. Equation 8is an example of two latent variables that are selected in the case of MOR.

(8)

Energy = 𝑆(𝑓) ² 𝑓 𝑆 1 = 𝐸 1 . and 𝑆 2 = 𝐸 2 . 𝑆 2 Energy 2 = 𝑆 1 Energy 1 ( 𝑆 Energy MOE B 1 EW ). 𝑀𝑂𝑅 1 𝑀𝑂𝑅 2 ⋮ ⋮ 𝑀𝑂𝑅 𝑛 = 1 1 𝑥 1,2 𝑥 2,2 ⋮ ⋮ 1 ⋮ ⋮ 𝑥 𝑛,2 … … 𝑥 1,𝑚 -1 … … 𝑥 2,𝑚 -1 𝑥 1,𝑚 𝑥 2,𝑚 … … ⋮ … … ⋮ … … 𝑥 𝑛,𝑚 -1 ⋮ ⋮ 𝑥 𝑛,𝑚 1 2 ⋮ ⋮ 𝑚 + 1 2 ⋮ ⋮ 𝑛 𝑀𝑂𝑅 1 𝑀𝑂𝑅 2 ⋮ ⋮ 𝑀𝑂𝑅 𝑛 = 𝑡 1,1 𝑡 2,1 𝑡 1,2 𝑡 2,2 ⋮ ⋮ 𝑡 𝑛,1 ⋮ ⋮ 𝑡 𝑛,2 𝛽 1 𝛽 2 + 𝜖 1 𝜖 2 ⋮ ⋮ 𝜖 𝑛
where is the i th observation of goal variable, 𝑡 is the j th element of the i th latent variable which is a linear combination of initial explicative variables, 𝛽 is the j th new coefficient element of regression and 𝜖 is the i th residual element.

The number of latent variables was selected in order to minimize the error between the predicted and actual values of the goal variable (cross validation with 10 random segments). The population was separated into ten random folds. Calibrations and validations were done ten times with different samples. Calibration was completed with 9 folds, while validation was completed with the last one. The principle is detailed in the publication of [START_REF] Naes | A User-Friendly Guide to Multivariate Calibration and Classification[END_REF]. After cross validation, the root mean square error of cross validation (RMSECV) was calculated (Eq. 9). The selected number of latent variables correspond to the lowest RMSECV. The associated coefficient of determination of cross validation (

) was also determined.

(

) 9
where ̂ is the i th estimate value of at the j th fold of cross validation (note that the same equation was used for MOE), and N j is the number of samples used in validation (j th fold).

The final model was computed using all the observations. The root mean square error of calibration (RMSEC) was defined on the total number of observations (n) as follows, (10) where k is the number of latent variables and n is the total number of observations. The computation was performed with R software with the PLS package [START_REF] Mevik | The PLS package: Principal component and partial least squares regression in R[END_REF].

RESULTS AND DISCUSSION

Comparison of Destructives and Nondestructive Tests

Table 4 shows the descriptive statistics of density, MOE, and MOR. Destructive tests showed a large representation of high MOR values (MOR > 30 MPa), which constitutes 88% of the total number of boards. However, boards with a low MOR were not negligible. Indeed, 12% of boards were under 30 MPa. Despite the large range between the maximum and minimum values of density (275 kg.m -3 ), its dispersion, which is represented by the coefficient of variation (COV), was low (COV = 7.9%) with a mean value of 711 kg.m -3 . Close values were observed in the PhD of [START_REF] Viguier | Mechanical Classification of Structural Wood[END_REF] where the mean of density equals to 725 kg.m -3 and the COV = 9.3%. The dispersion of MOE was more pronounced than the density with a COV = 20%, a range of 10669 MPa, and a mean value of 10301 MPa. The values observed in the PhD of [START_REF] Viguier | Mechanical Classification of Structural Wood[END_REF] where the mean of MOE equals to 10140 MPa and the COV = 24.8% are close to those obtained in this study.

The same observation was made for MOR with a large dispersion (COV = 39%), a range of 96 MPa and a mean value of 56 MPa. The mean value of MOR observed in the PhD of Viguier (2016) were equals to 40.7 MPa and the COV% = 37.8% which are lower comparing to this study.

Coefficients of determination values (with interval of confidence p < 0.05) and their degree of significance between density, MOE, and MOR are shown in Fig. 2. The correlation between static MOE and density was low (R² = 0.09) with a large interval of confidence. Thus, predicting MOE from the density appears to be inappropriate. A similar observation was made between the density and static MOR (R² = 0.01). However, the r-squared value between the static MOE and static MOR was equal to 0.40. As a result, it was interesting to determine the MOE to predict the MOR by using nondestructive techniques. Typically, grading machines based on vibration measurements use the dynamic MOE obtained from the first Eigen frequency in longitudinal vibration to determine MOE and MOR. In the current study, MOE and MOR were predicted from the dynamic MOE calculated by using transversal (edgewise vibrational MOE) and longitudinal (longitudinal vibrational MOE) vibration tests. A cross-validation was performed to estimate the prediction error of the statistical modeling. The models were stable in both cases of longitudinal and edgewise-transversal vibrations. Indeed, the RMSECV was close to the RMSEC. In the longitudinal vibration model, the prediction of MOE leads to a RMSEC of 1100 MPa and a RMSECV of 1035 MPa and the prediction of MOR leads to RMSEC of 19.1 MPa and a RMSECV of 19.2 MPa. The same behavior is observed in edgewise transversal vibrations model where the prediction of MOE leads to a RMSEC of 896 MPa and a RMSECV of 835 MPa and the prediction of MOR leads to RMSEC of 18.7 MPa and a RMSECV of 18.9 MPa. For both MOE and MOR prediction, the errors were higher for the longitudinal vibration than for the edgewise-transversal vibration.

Either in the longitudinal or transversal vibration, the vibrational method was a good way to predict the MOE for oak species ( > 0.75). In contrast, the MOR was not estimated accurately by the vibration tests ( < 0.27). Thus, it is interesting to show how far spectrum analyses improve the prediction of MOE and MOR, first by using parameters based on Eigen frequencies, then the full spectrum.

Method based on Eigen Frequencies

This analysis used a stepwise regression method. The goal consisted of finding the best global output parameters that explain the variance of observations, i.e., MOE and MOR. Each model contained significant parameters with p-values that were under 0.001.

Stepwise regression for MOE

The selected regression parameters to predict MOE were and Ratio 2 which is the ratio / . The determination coefficient of the regression was = 0.86 (F-statistic = 496, p-value < 2.2 × 10 -16 ) with a residual error of RMSEC = 812 MPa. The result of the cross-validation procedure was RMSECV = 817 MPa ( = 0.86). This model was found to be stable because the RMSECV was very close to the RMSEC. Figure 3a shows the predicted values according to static MOE values.

The first selected parameter was the edgewise dynamic modulus, which corresponded to the elastic property measured in the same configuration than the static test (edgewise flexural loading). The second selected parameter was the second ratio between edgewise and flatwise dynamic moduli. This ratio was affected by the presence and positions of natural defects in the boards (Roohania et al. 2015). Due to the collinearity of the explicative variables, the first ratio ( / ) or the coefficient of determination between Eigen frequencies ( ) could be used too (in addition to parameter). By doing so, for the first ratio the model results were = 0.86, RMSEC = 825 MPa, RMSECV = 830 MPa; and for the model results were = 0.85, RMSEC = 831 MPa, RMSECV = 833 MPa.

Stepwise regression for MOR

The selected parameters to predict MOR were and . The regression coefficient was equal to = 0.35 (F-statistic = 42.8, p-value = 1.23 × 10 -15 ), and the residual error was RMSEC = 18.0 MPa. The associated cross-validation gave a RMSECV value of 18.1 MPa ( = 0.32). This model was stable because the RMSECV was close to the RMSEC. Figure 3b shows the regression scatterplot for the MOR values.

As well as the static MOE model, the first parameter selected was an edgewise elastic property (

). The second parameter was the coefficient of determination associated with the Bordonné's solution of the Timoshenko equation [START_REF] Brancheriau | Natural vibration analysis of clear wooden beams: Theoretical review[END_REF]. This coefficient was affected by the values of the Eigen frequencies. If the material was perfectly homogeneous, would be equal to one. Otherwise, the presence and position of natural defects shifted the Eigen frequencies, lowering the value [START_REF] Roohnia | Orientation and position effects of a local heterogeneity on flexural vibration frequencies in wooden beams[END_REF]. The dynamic MOE corresponding to the first and second Eigen frequencies were correlated ( , R-squared value of 0.95). Thus, replacing with in the multiple regression with lead to similar results: = 0.34, RMSEC = 18.1 MPa, and RMSECV value of 18.2 MPa.

Method Based on the Full Vibrational Spectrum

Because multiple linear regressions showed that edgewise transversal vibration parameters were mainly selected to predict MOE and MOR, only the edgewise transversal vibrations were used in this section. The number of predictive variables (standardized frequencies) was reduced according to the bilateral R² threshold between the considered frequency and the dependent variable (MOE or MOR) in order to ensure stability to the model. We kept the frequencies associated with a bilateral R² greater than 0.025 (several R² thresholds were tested).

PLS regression for MOE

In regards to the bilateral R² criterion on the frequency selection, 597 frequencies were chosen for the MOE prediction (Fig. 4a). Figure 4a shows that the relevant frequencies kept for MOE prediction were below 10 f/f 0 . This corresponded to a maximum frequency range of 1250 Hz because the first Eigen frequency was 125 Hz maximum (average of 91 Hz). The bilateral R² value did not exceed 0.40. The selected frequencies were not only focused on the harmonic frequencies of the flexural vibrations but considers the whole spectrum. At low frequencies the vibration energy was at a maximum, but the wavelengths were large. Consequently it had little sensibility to punctual defects (high standardized frequencies were more sensible to punctual defects). However, the presence of natural defects might cause mode conversions from flexural to torsional and compressional vibrations, even in low frequencies. Four latent variables were needed to minimize the root mean square error of cross validation (RMSECV = 1458 MPa, = 0.54). The partial least squares regression was significant with a value of 0.71 (F-statistic = 97.8, p-value < 2 × 10 -16 ) and a RMSEC = 1175 MPa. However, the difference between the RMSECV and RMSEC highlighted the instability of the model (ratio of 1.24). The realistic error and determination coefficient of this model were those obtained by cross-validation. Figure 4b shows the reference MOE values according to predicted values of the partial least squares model.

PLS regression for MOR

The process of variable selection by bilateral R² was conducted to 290 selected frequencies (Fig. 5a). The selected frequencies were not all identical to those associated to the MOE. However, they remain 5 times lower than the first Eigen frequency, thus mostly below 625 Hz. For this model, thirteen observations were considered as outliers and were removed. The outlier detection was based on the distance of each observation to the centroid of the population computed in the basis defined by the selected scores (Euclidian distance applied in an orthogonal basis, which was equal to the Mahalonobis distance in this case).

To predict MOR, two latent variables minimized the error of cross-validation with a RMSECV value of 17.6 MPa (R² CV = 0.33). The partial least squares regression was significant with a value of 0.44 (F-statistic = 58.1, p-value < 2 × 10 -16 ) and a RMSEC value of 16.25 MPa. RMSECV was found to be close to RMSEC with a ratio of 1.08. Figure 5b shows the reference MOR values according to the predicted values of the partial least squares model.

MOE and MOR predictive variables were obtained by multiplying the amplitudes of the frequency spectrum by the dynamic modulus. From a statistical point of view, there was only an interaction between the variables of amplitude and the dynamic modulus. The dynamic modulus was linked to the Eigen frequency, the density, and the dimensions. The amplitude of the vibrations depended on the internal viscosity of the material, the presence of natural defects, the wood stiffness, the positions of the supports, and the sensor. The interaction between these variables provided complementary information on the mechanical behavior of the wooden beams. EW . It is interesting to notice that the longitudinal vibration is the more common method used in industry for practical reasons. Finally, the PLS method based on the full vibrational spectrum did not improve the prediction of MOE regarding RMSECV and

Comparison of Models

(1458 MPa, 0.54, respectively). For the MOR prediction, the PLS method based on the full vibrational spectrum conducted to the lowest RMSECV value of 17.6 MPa, which was similar to the one obtained when the static MOE was used as a predictive criterion. The corresponding value of 0.33 was greater than those of the others models. Indeed, with the stepwise regression method based on Eigen frequencies, equals 0.32 (RMSECV of 18.1 MPa). Using the first Eigen frequency in transversal vibration (MOE B1 EW ), equals 0.26 (RMSECV of 18.9 MPa). The first Eigen frequency in longitudinal vibration (MOEC) leads to of 0.23 (RMSECV of 19.2 MPa). As a result, the PLS method based on the full vibrational spectrum was better suited for MOR prediction. Indeed, the full spectrum contains the local information about wood defect (magnitudes taken by the frequency spectrum) which allows a better determination of MOR. But, the MOE is mainly dependent on the Eigen frequencies values, and this information was lost when the frequency scale was standardized this is why the stepwise regression method based on Eigen frequencies or the first edgewise transversal vibration (MOEB1EW) were better for the MOE prediction.

Because Since the PLS method was based on the full vibrational spectrum in edgewise transversal vibration, it is possible to build a single grading machine that can predict MOR and MOE with the best values in a single measure. The obtained errors would be significantly lower than the current grading machine based on longitudinal vibration. For MOE prediction the would increase from 0.75 to 0.86; and for MOR prediction the would increase from 0.22 to 0.33. However, the PLS model is quite unstable due to the model sensibility to noise in the signals. To overcome this problem, specific filtering procedures and/or a signals average with multiple acquisitions can be tested. A larger experimental sample can also increase the robustness of the model. CONCLUSIONS 1. By using transversal vibrations and statistical methods, it is possible to significantly improve MOE and MOR prediction with respect to the usual model based on longitudinal vibrations.

2. The stepwise regression method based on Eigen frequencies allowed an improvement of the MOE prediction by choosing specific global vibrational parameters. The main parameter is the dynamic MOE calculated from the first resonance frequency in the edgewise bending configuration. Concerning the MOR estimation, this method had a better efficiency than the usual models based only on the dynamic MOEs (bending or compression).

3. The PLS method based on the full vibrational spectrum did not improve the MOE prediction. However, the prediction of the MOR was very efficient and led to the lowest prediction error.

4. A combination of the two proposed statistical methods is possible in order to obtain the best prediction of MOE and MOR from a single edgewise vibrational measurement.

Fig. 2 .

 2 Fig. 2. Scatterplots of density, static MOE, static MOR, and vibrational MOE. R-squared values are in the upper right panel, with the significance in upper case (*** for p-value <0.001, nothing for 0.1 < p-value < 1). The confidence intervals at the 0.05 level are shown in brackets.

Fig. 3 .

 3 Fig. 3. (a) Static MOE vs. predicted MOE. (b) Experimental MOR vs. predicted MOR. In both panels, the solid line is the first bisector.

Fig. 4 .

 4 Fig. 4. (a) Bilateral coefficient of determination between standardized frequencies and MOE. The hollow circles are the selected frequencies (R² ≥ 0.025). (b) Relationship between static and predicted MOE. Partial least squares regression on the full vibrational spectrum. The solid line is the first bisector.

Fig. 5 .

 5 Fig. 5. (a) Bilateral coefficient of determination between standardized frequencies and MOR. The hollow circles are the selected frequencies (R² ≥ 0.025). (b) Relationship between static and

Table 1 .

 1 Board Grading based on Visual Method (NF B52 001-1 (2011))

	Grades	Number of Boards
	1 (D30)	30
	2 (D24)	37
	3 (D18)	25
	Waste	72

Table 2 .

 2 Board Numbers for Each Section

	Number of Boards	Dimensions (mm)
	51	2200 x 76 x 23
	52	2000 x 97 x 24
	38	2500 x 105 x 23
	23	3000 x 170 x 25

Table 3 .

 3 List of Global Output Vibrational Variables

	Parameters	Reference	Number of Parameters
	Density	ρ	1
	Longitudinal vibrational MOE in compression	MOE C	1
	Transversal vibrational MOE (flatwise and edgewise) relative to i th Eigen frequency based on Bernoulli theory Transversal vibrational MOE (flatwise and edgewise) based on Timoshenko theory	MOEB i MOEB i EW , FW MOET EW , MOET FW	8 2
	Shear moduli relative to edgewise vibration	G	EW	1
	Frequency ratios between first Eigen frequencies and i th Eigen frequency, with 2<i<4 (flatwise and edgewise)	Inha i Inha i EW , FW	6
	Ratio between flatwise and edgewise vibrational MOE relative to i th Eigen frequency	Ratio i	4
	Coefficient of determination between Eigen frequencies (flatwise and edgewise)		,	2
	loss factor (tan δ)	D C , D EW	2

Table 4 .

 4 Statistics of Density, MOE, and MOR Obtained by Destructive Tests

	Basic statistics	Min	1 st quartile Mean Std. deviation Median 3 rd quartile Max COV (%)
	Density (kg.m	-3 )	582	682	711	56	712	744	857	7.9
	MOE (MPa)	5385	8780	10301	2159	10114	12000	16054	20
	MOR (MPa)	13	38	56	22	55	75	109	39

Table 6 .

 6 Summary of Prediction Methods

			MOE		MOR
	Model	RMSECV (MPa)		RMSECV (MPa)
	Static MOE			17.2	0.39
	MOE B1 EW	835	0.85	18.9	0.26
	MOE	C	1035	0.75	19.2	0.23
	Stepwise regression based on Eigen frequencies	817	0.86	18.1	0.32
	PLS based on the full vibrational spectrum	1458	0.54	17.6	0.33

  Table 6 shows the results of MOE and MOR predictions for the different models in cross-validation. Concerning the MOE predictions, the stepwise regression method based on Eigen frequencies presented the lowest RMSECV, which equals 817 MPa. The associated equals to 0.86, which demonstrates that this method explains a large part of the variance of MOE. However, a model based only on the first Eigen frequency in edgewise transversal vibration (MOE B1 EW ) presented similar values. This result was expected because MOE B1 EW was the main variable used in the stepwise regression model. The model based on the first Eigen frequency in the longitudinal vibration (MOE C ) presents a higher RMSECV value of 1035 MPa. The reason why the models based on

	MOE B1	EW were the most efficient must be because the static MOE was measured in a
	bending test like MOE B1

ACKNOWLEDGMENTS

This study was funded by the regional council of Bourgogne Franche-Comté, Carnot ARTS institute, French ministry of ecology, sustainable development and energy and ministry of housing and territorial equality, the CODIFAB, and "France Bois Forêt" interbranch organization.