H Ould-Baba
email: hanany.ould-baba@ut-capitole.fr

J Antoni

V Robin

Multidimensional Independent Component Analysis with Higher-order cumulant matrices for vector sources with possibly differing dimensions ✩,✩✩

Keywords: Multidimensional Independent Component Analysis, Independent Component Analysis, Higher-order cumulant, Independent Subspace Analysis, Joint Block-Diagonalization

The paper addresses the separation of multidimensional sources, with possibly different dimensions, by means of higher-order cumulant matrices. First, it is rigorously proved, in a general setting, that contracted cumulant matrices of any order are all block-diagonalizable in the same basis. Second, a family of joint blockdiagonalization algorithms is proposed that separate multidimensional sources by combining contracted cumulant matrices of arbitrary orders. Third, a specific solution is given to determine the source dimensions when they are unknown but all different. The performances of the proposed algorithms are compared between them and with algorithms of the literature based on orders 3 and 6.

Introduction

In its original formulation independent component analysis (ICA) assumes the mutual independence of the sources to be separated (e.g. see Comon [START_REF] Comon | Independent component analysis, a new concept?[END_REF]). Unfortunately, there are many instances in the real-world where this assumption is not fulfilled, which precludes the recovery of sources that are (partly) dependent. This has brought in the concept of multidimensional independent component analysis (MICA), where multidimensional rather than scalar sources are considered Cardoso [START_REF] Cardoso | Multidimensional independent component analysis[END_REF]. Formally, the MICA model reads

x = As = A 1 A 2 • • • A p      s 1 s 2 . . . s p      = x 1 + x 2 + • • • + x p (1)
where

                      
x ∈ R m is the vetor of observations, A ∈ R m×n is the unknwon mixing matrix, s ∈ R n is the unknown vector of sources, s i ∈ R ni denotes the i-th source, A i ∈ R m×ni denotes the i-th block of matrix A, x i = A i s i ∈ R m denotes the i-th independent component, n = p i=1 n i is the dimension of the source vector, p ∈ N denotes the number of multidimensional sources.

(

) 2
That is to say MICA consists in decomposing an arbitrary random vector, x, into its independent components (ICs), x i , i = 1, ..., p where a vector v is an IC of a random vector u if there exists an invertible matrix F and an decomposition u = F v w where v and w are independent random vectors. Up to trivial indeterminacies (see subsection 2.1), this amounts to estimating the unknown mixing matrix A (or its inverse, the separation matrix B) and sources s in model ( 1) from the sole observations x. In other words, when x accepts a decomposition into one-dimensional (scalar ) sources, the problem boils down to classical ICA [START_REF] Comon | Independent component analysis, a new concept?[END_REF]. More generally, when the decomposition involves k-dimensional sources, this is referred to as k-MICA, which is the topic of this paper. So far MICA has attracted much less attention than ICA, although it is also known under (or related to) the following names in the literature: independent subspace analysis (ISA) [START_REF] Hyvärinen | Emergence of phase-and shift-invariant features by decomposition of natural images into independent feature subspaces[END_REF], independent feature subspace analysis (IFSA) [START_REF] Kim | Membership scoring via independent feature subspace analysis for grouping coexpressed genes[END_REF], subspace ICA [START_REF] Sharma | Subspace independent component analysis using vector kurtosis[END_REF], group ICA [START_REF] Theis | Blind signal separation into groups of dependent signals using joint block diagonalization[END_REF].

The application fields where MICA finds interest are indeed numerous. Examples comprise electrocardiography (ECG) [START_REF] Cardoso | Multidimensional independent component analysis[END_REF][START_REF] Theis | Blind signal separation into groups of dependent signals using joint block diagonalization[END_REF][START_REF] Sharma | Subspace independent component analysis using vector kurtosis[END_REF][START_REF] Chawla | Detection of indeterminacies in corrected ecg signals using parameterized multidimensional independent component analysis[END_REF], f-MRI electroencephalography (EEG) [START_REF] Ma | Independent subspace analysis with prior information for fmri data[END_REF], topography [START_REF] Hyvärinen | Topographic independent component analysis[END_REF], texture classification [START_REF] Santos | Using independent subspace analysis for selecting filters used in texture processing[END_REF], action recognition [START_REF] Le | Learning hierarchical invariant spatio-temporal features for action recognition with independent subspace analysis[END_REF], face recognition [START_REF] Li | View-based clustering of object appearances based on independent subspace analysis[END_REF], motion segmentation [START_REF] Fan | Motion segmentation based on independent subspace analysis[END_REF], audio source separation [START_REF] Casey | Separation of mixed audio sources by independent subspace analysis[END_REF], gene decoding [START_REF] Kim | Membership scoring via independent feature subspace analysis for grouping coexpressed genes[END_REF][START_REF] Kim | Tree-dependent components of gene expression data for clustering[END_REF], metabolomic data analysis [START_REF] Gutch | An isa algorithm with unknown group sizes identifies meaningful clusters in metabolomics data[END_REF].

Moreover, although most previous works have concentrated on separating ICs with identical dimensions (k-MICA with k = dim(s i ), ∀i) [START_REF] Hyvärinen | Topographic independent component analysis[END_REF][START_REF] Theis | Multidimensional independent component analysis using characteristic functions[END_REF][START_REF] Theis | Uniqueness of complex and multidimensional independent component analysis[END_REF][START_REF] Theis | Blind signal separation into groups of dependent signals using joint block diagonalization[END_REF][START_REF] Hyvärinen | FastISA: A fast fixed-point algorithm for independent subspace analysis[END_REF][START_REF] Sharma | Subspace independent component analysis using vector kurtosis[END_REF] (in which case the problem is known as k-MICA with k = dim(s i ) = n i ), this situation is by no means general. This paper considers the general case where p multidimensional ICs with arbitrary and unknown dimensions, n i = dim(s i ) = dim(x i ) (not necessarily identical), are to be separated. Several approaches have been considered in the literature to solve the general MICA problem. An early solution consists of regrouping a posteriori the components separated by traditional ICA according to their residual dependence. Although the idea was globally covered in Ref. [START_REF] Cardoso | High-order contrasts for independent component analysis[END_REF], no proof of separability and unicity was provided. The first proof of separability of the k-MICA problem was given by Theis [51] and is based on a theorem of prob ability introduced in Ref. [START_REF] Gutch | Second order subspace analysis and simple decompositions[END_REF], in the sixteens, that generalizes the Darmois-Skitovitch theorem [START_REF] Févotte | Orthonormal approximate joint block-diagonalization[END_REF]48] to the multidimensional case (square matrices). The complete proof of unicity was later given by Theis [50] and Gutch and Theis [START_REF] Hyvärinen | Emergence of phase-and shift-invariant features by decomposition of natural images into independent feature subspaces[END_REF] in the general case of independent subspaces with possibly different dimensions. It is there remarked that the initial model introduced in Ref. [START_REF] Cardoso | High-order contrasts for independent component analysis[END_REF] is not complete and requires an additional assumption for being separable. Indeed, for any observation x, a decomposition of x into independent subspaces with some internal dependence in each subspace is given by x itself. The additional condition required to avoid this trivial decomposition has been coined "irreducibility of multidimensional sources". The condition is revisited in subsection 8 of this paper (see Definition 8) Note that the proof of unicity of the MICA separation has recently been completed in Ref. [START_REF] Hyvärinen | FastISA: A fast fixed-point algorithm for independent subspace analysis[END_REF]. Besides, some joint block-diagonalization (JBD) algorithms have been proposed that potentially provide solutions to MICA by reducing a set of matrices to diagonal blocks with different sizes [50, [START_REF] Nion | A tensor framework for nonunitary joint block diagonalization[END_REF][START_REF] Sharma | Subspace independent component analysis using vector kurtosis[END_REF]. They make use of second-order statistics [49, [START_REF] Kim | Membership scoring via independent feature subspace analysis for grouping coexpressed genes[END_REF][START_REF] Mac Cullagh | Tensor Methods in Statistics[END_REF], of the functional form of the probability density function or its generating function [52,49], or of dispersion matrices [START_REF] Theis | Towards a general independent subspace analysis[END_REF]. Two such algorithms have been applied to the MICA problem in Ref. [49]. The first one, which is a direct extension of SOBI [START_REF] Albera | Asymptotic performance of contrast-based blind source separation algorithms[END_REF], jointly blockdiagonalizes a set of correlation matrices taken at different time lags and has been coined Multidimensional SOBI (MSOBI). The second one jointly block-diagonalizes a set of Hessian matrices of the first (or second) generating function and has been coined Multidimensional Hessian ICA (MHICA). However, solutions to MICA based on higher-order cumulants are still relative ly seldom although they offer several advantages. Indeed, Theis in [START_REF] Theis | Towards a general independent subspace analysis[END_REF] has proposed to simply solve the MICA problem by jointly block-diagonalizing a set of fourth-order cumulant matrices in the spirit of the JADE algorithm [START_REF] Cardoso | Blind beamforming for non gaussian signals[END_REF], which he coined SJADE (for Subspace-JADE). Yet, the following questions have so far remained unanswered:

-are fourth-order cumulant matrices (as initially defined in [START_REF] Cardoso | Blind beamforming for non gaussian signals[END_REF] and used in Theis [START_REF] Theis | Towards a general independent subspace analysis[END_REF]) block-diagonalizable in the basis spanned by A? -more generally, are cumulant matrices of arbitrary order r (to be defined later on in the paper) blockdiagonalizable in the basis spanned by A?

In other words, can the MICA problem be formally formulated as a joint block-diagonalization problem on arbitrary order r? On the second-order, the answer is positive and has been given in [START_REF] Theis | Blind signal separation into groups of dependent signals using joint block diagonalization[END_REF][START_REF] Gutch | Second order subspace analysis and simple decompositions[END_REF][START_REF] Lahat | Optimal performance of second-order multidimensional ica[END_REF]. On higher-orders, the answer is yet not trivial and is the subject of the present paper.

The remainder of the paper is organized as follows. First, the usual assumptions that sustain the MICA problem are reviewed in section 2. Next, the definitions of cumulant matrices and contracted cumulant matrices are introduced in section 3 and 4, respectively. Section 5 contains the main result of the paper, which proves that contracted cumulant matrices of any order are all block-diagonalizable in the same basis. This result is then exploited in section 6 to generalize the JADE algorithm to the multidimensional case with arbitrary orders, thus leading to a family of algorithms coined SJADE r . The advantage of combining various orders is illustrated in section 7. Section 8 addresses the important question as how to determine the dimensions of the vector sources when they are unknown. The proposed algorithms are finally compared in section 9 by means of numerical experiments. Simple examples are given to demonstrate the advantage of combining different orders in order to improve the separation. All proofs are collected in the Appendices.

Generalities about MICA

From the onset, some generalities about ICA are to be reminded in the context of multidimensional sources.

Uniqueness and indeterminacies

The indeterminacies underlying MICA are obvious generalizations of those of ICA. They are reminded here for the sake of completeness. It suffices to note that, for any invertible matrices

D i of size n i × n i , i ∈ [[1, p]],
and for any permutations σ ∈ S p , model (1) still holds since

A σ(1) D σ(1) • • • A σ(p) D σ(p)    D -1 σ(1) s σ(1)
. . .

D -1 σ(p) s σ(p)    = x σ(1) + • • • + x σ(p) = x. (3)
Therefore, the multidimensional ICs x i are determined up to an arbitrary permutation (that applies on components with identical dimension), whereas sources are determined up to an arbitrary permutation and invertible matrices. Taking such indeterminacies into account, a matrix B will be recognized as a separation matrix (i.e. a solution to the MICA problem) if BA = PD, where P is a permutation matrix and

D = bdiag(D 1 , D 2 , . . . , D p ) a block-diagonal matrix.
Besides, the following definition will play an important role later on.

Definition 1 ( [START_REF] Theis | Towards a general independent subspace analysis[END_REF][START_REF] Gutch | Uniqueness of linear factorizations into independent subspaces[END_REF]). A random vector, x, of dimension n is said reducible if it can be expressed as

x = A y 1 y 2 , ( 4
)
where A is an invertible matrix and y 1 and y 2 are two independent random vectors of dimensions k = 0 and nk, respectively. A random vector that is not reducible is said irreducible.

Assumptions

One is now in position to announce the assumptions that will guarantee the uniqueness of the separation model:

1. matrix A is of full-rank n, 2. no more than one independent (scalar ) source can have a nil cumulant of order r, c r (s i ) = 0, when working on order r, 3. independent sources are all irreducible.

Whitening

Whitening is a classical preprocessing used in blind source separation. It is used here to simplify the theoretical developments. Without loss of generality, let us assume that the source vectors are centered, E(s) = 0, and standardized, R s = E(ss ′ ) = I n . Thus R x = E(xx ′ ) = AA ′ . Since the covariance matrix R x is symmetric and semi-positive definite, it is well-known that it admits an eigenvalue decomposition R x = UDU ′ where D is an n × n diagonal matrix, with n the number of non-zero eigenvalues, and U m×n an m × n matrix satisfying U ′ U = I n . Following the usual practice, let us define the whitening matrix W n×m = D -1 2 U ′ and the whitened components x = (Wx) n×1 ; thus

R x = E( x x ′ ) = WR x W ′ = D -1 2 U ′ (UDU ′ )UD -1 2 = I n . (5)
Since R x = AA ′ , it follows that WAA ′ W ′ = I n ; in other words, A = WA is an n × n orthogonal matrix. Upon pre-multiplication with the whitening matrix, model (1) then becomes x = As. Therefore, it will be assumed from now on and without loss of generality that matrix A is orthogonal (with m = n). Solving the MICA problem then amounts to finding, using higher-order statistics, the orthogonal matrix A.

Notations

This last subsection introduces some notations that will be used in the remaining of the paper. Let us first note that matrix A can be partitioned into n × n i sub-matrices A i or as a collection of raw vectors ℓ j :

A = A 1 A 2 . . . A p =      ℓ 1 ℓ 2 . . . ℓ n      . (6)
Alternatively, matrix A can also be partitioned into raws as

A i =      ℓ i 1 ℓ i 2 . . . ℓ i n      (7)
where, ∀k ∈ [ [1, n]], ℓ i k with dimension 1 × n i is the k-th raw of sub-matrix A i . Similarly, each raw of matrix A can be expressed as a function of the raws of matrices A i , as follows

ℓ k = ℓ 1 k ℓ 2 k . . . ℓ p k . ( 8
)
One is now ready to introduce specific statistical tools -cumulant matrices, their contracted forms, and the related algebra -dedicated to the purpose of this paper.

Cumulant matrices

The main tools used in this paper are cumulants. In the scalar case it is a well-known fact that cumulants can be deduced from moments by means of the Leonov and Shiryaev formula (see e.g. Mac Cullagh [START_REF] Mac Cullagh | Tensor Methods in Statistics[END_REF] et Albera and Comon [START_REF] Albera | Asymptotic performance of contrast-based blind source separation algorithms[END_REF]). However, it is often convenient to rearrange cumulants into matrices (see e.g. Albera et al. [START_REF] Albera | Blind identification of over-complete mixtures of sources (biome)[END_REF]), for instance in order to exploit their algebraic properties (such as matrix redundancies Albera et al. [START_REF] Albera | Sixth order blind identification of underdetermined mixtures (birth) of sources[END_REF][START_REF] Albera | Icar: a tool for blind source separation using fourth-order statistics only[END_REF][START_REF] Albera | Blind identification of over-complete mixtures of sources (biome)[END_REF]) or to diagonalize them Kollo [START_REF] Kollo | Multivariate skewness and kurtosis measures with an application in ica[END_REF]. In the multidimensional case, the notion of a cumulant matrix needs to be carefully defined for it will play a key role in the remaining of the paper.

Definition 2 (cumulant matrix). Let x be a random vector with dimension n × 1. Then, the cumulant matrix of order r is defined as

c r (x) = n i1,i2,...,ir=1 cum(x i1 , x i2 , • • • , x ir )E i1,i2,...,ir (9)
with

E i1,i2,...,ir =        e i1 ⊗ e i3 ⊗ • • • ⊗ e i 2k-1 e i2 ⊗ e i4 ⊗ • • • ⊗ e i 2k ′ if r = 2k e i1 ⊗ e i3 ⊗ • • • ⊗ e i 2k+1 e i2 ⊗ e i4 ⊗ • • • ⊗ e i 2k ′ if r = 2k + 1, (10)
where (e i ) i∈[ [1,n]] denotes the canonical basis of R n and ⊗ the Kronecker product.

Closed-form expressions of cumulant matrices in terms of lower-order moments were given in Ref. [START_REF] Kollo | Multivariate skewness and kurtosis measures with an application in ica[END_REF] up to orders 3 and 4, and for the first time in Ref. [START_REF] Ould-Baba | Concise formulae for the cumulant matrices of a random vector[END_REF] up to order 6. They are reproduced in Appendix (Appendix B) for the sake of completeness. Not only are such expressions easy to handle due to their compactness, but they also lead to faster numerical computation with matrix-based languages such as Matlab (interested readers are invited to consult [START_REF] Ould-Baba | Concise formulae for the cumulant matrices of a random vector[END_REF] for the systematic derivation of these formulae). It is noteworthy that Definition 2 produces n k × n k square and symmetric cumulant matrices of even orders (r = 2k) which, as compared to other possible definitions, will turn out advantageous; for instance, on order r = 4, c 4 (x) is a n 2 × n 2 matrix known as the "quadricovariance" [START_REF] Albera | Icar: a tool for blind source separation using fourth-order statistics only[END_REF]; similarly, on order r = 6, c 6 (x) is a n 3 × n 3 matrix known as the "hexacova riance" [START_REF] Albera | Sixth order blind identification of underdetermined mixtures (birth) of sources[END_REF]. For odd orders r = 2k + 1, the cumulant matrices have dimensions n k+1 × n k . Several useful properties can be deduced from Definition 2, a few of which are reminded here below.

Proposition 1. (i) If z is a Gaussian random vector, then

∀r > 2, c r (z) = 0. (11)
(ii) If x and y are two independent vectors with identical dimensions, then

∀r ∈ N c r (x + y) = c r (x) + c r (y). (12)
(iii) If x = As, where A is a matrix with dimension m × n and x and s are two random vectors with dimensions m × 1 and n × 1, respectively, then

c r (x) =              (A ⊗ • • • ⊗ A) k fois c r (s) (A ⊗ • • • ⊗ A) ′ k fois = ⊗ k A c r (s) ⊗ k A ′ if r = 2k (A ⊗ • • • ⊗ A) k+1 fois c r (s) (A ⊗ • • • ⊗ A) ′ k fois = ⊗ k+1 A c r (s) ⊗ k A ′ if r = 2k + 1 . ( 13
) (iv) If x = As = A 1 • • • A p    s 1 . . . s p    ∈ R m
, where s i ∈ R ni are independent random vectors,

A i ∈ R m×ni , n = n 1 + • • • + n p and A ∈ R m×n , then for any integer r, c r (x) =    A k bdiag (c r (s 1 ), . . . , c r (s p )) A ′ k if r = 2k A k+1 bdiag (c r (s 1 ), . . . , c r (s p )) A ′ k if r = 2k + 1 ( 14
)
where matrices A k are returned by

A k = ⊗ k A 1 ⊗ k A 2 • • • ⊗ k A p ∈ R m k × p j=1 n k j . ( 15
)
Proof. See Appendix A. Properties (i) to (iv) will be used to establish the main results of this paper.

Contracted cumulant matrices

The higher-order-statistic solution to MICA is based on an extension of the former cumulant matrix which is referred herein as the "contracted" cumulant matrix. This generalizes a notion initially introduced on the fourth-order by Cardoso and Soulomiac in Ref. [START_REF] Cardoso | Blind beamforming for non gaussian signals[END_REF] in the JADE algorithm. The definition is first given hereafter on the fourth-order, and then on arbitrary order.

Definition 3 (contracted matrix cumulant on order 4). Let

M = (m i1,i2 ) (i1,i2)∈[[1,n]] 2 be an n × n matrix and x an n × 1 random vector. The contracted cumulant matrix of order four of dimension n × n, denoted as Q x (4) [M], has generic term Q x (4) [M] i3,i4 = n i1,i2=1 cum(x i1 , x i2 , x i3 , x i4 )m i1,i2 . (16)
The following lemma shows how the contracted cumulant matrix Q x (4) can be directly expressed as the star product MacRae [START_REF] Macrae | Matrix derivatives with an application to an adaptive linear decision problem[END_REF] (see Appendix Appendix C) of matrix M with the cumulant matrix c 4 . This will be needed in Definition 4.

Lemma 2. Let M = (m i1,i2 ) (i1,i2)∈[[1,n]] 2 be a matrix with dimension n × n; then Q x (4) [M] = M * c 4 (x) = n i1,i2=1 m i1,i2 [c 4 (x)] i1,i2 (17)
where

[c 4 (x)] i1,i2 is an n × n matrix corresponding to block (i 1 , i 2 ) of the cumulant matrix c 4 (x).
Proof. It suffices to realize that element

(i 3 , i 4 ) of matrix [c 4 (x)] i1,i2 is [c 4 (x)] i1,i2 i3,i4 = cum(x i1 , x i2 , x i3 , x i4 ). (18)
An alternative proof is to apply properties (P 1) to (P 3) (see Appendix C.2, C.3 and C.4) of the star product * to the expression of c 4 (x) given in Appendix B; it follows that

M * c 4 (x) = E(x ′ Mxxx ′ ) -ΣMΣ -ΣM ′ Σ -tr(ΣM)Σ (19)
which is the explicit expression of Q x (4) [M] (see [START_REF] Cardoso | High-order contrasts for independent component analysis[END_REF]).

The generalization of the definition to any order r is obtained as follows.

Definition 4 (contracted cumulant matrix of order r). Given q ∈ N, any matrix matrix M (q) expands as

M (q) = n i1,i2,...,iq=1 m i1,i2,••• ,iq E i1,i2,...,iq (20)
where matrices E i1,i2,...,iq have been introduced in Eq. ( 10) and m i1,i2,••• ,iq are arbitrary constants. Then, the n × n contracted cumulant matrix of order r is defined as

Q x (r) [M (r-2) ] = M (r-2) * c r (x) = n i1,i2,...,ir-2=1 m i1,i2,••• ,ir-2 c r (x) i1,i2,••• ,ir-2 (21)
where c r (x)

i1,i2,••• ,ir-2
denotes the matrix indexed by

(i 1 , i 2 , • • • , i r-2 ) in the partition of c r (x) into ma- trices of dimensions n × n; that is, element (i r-1 , i r ) of c r (x) i1,i2,••• ,ir-2 is c r (x) i1,i2,••• ,ir-2 ir-1,ir = cum(x i1 x i2 • • • x ir ). (22)
It results immediately that element

(i r-1 , i r ) of the contracted cumulant matrix Q x (r) [M (r-2) ] reads Q x (r) [M (r-2) ] ir-1,ir = n i1,i2,...,ir-2=1 cum(x i1 x i2 • • • x ir )m i1,i2,••• ,ir-2 , (23)
which corresponds to the natural generalization of ( 16).

r-th order MICA by means of joint block diagonalization

Equipped with the tools introduced hitherto, ones is now in a position to prove the main result of this paper concerning the block-diagonalization of the contracted cumulant matrices. The first step towards this perspective is to note that the cumulant matrix of the observations in the MICA model (1) reads

c r (x) =    A k C 2k (s)A ′ k if r = 2k A k+1 C 2k+1 (s)A ′ k if r = 2k + 1 (24) with C r (s) = bdiag c r (s 1 ), c r (s 2 ), • • • , c r (s p ) (25)
where use has been made of property (iv) of proposition 1 and of the statistical independence of the sources. The structure of matrices A k in Eq. ( 24) is given in the following lemma:

Lemma 3 (structure of A k ). Given any integer k 1, matrix A k = ⊗ k A 1 ⊗ k A 2 • • • ⊗ k A p of dimension n k × p i=1 n k i is
returned by the block matrix

A k =                           AD 11•••11 AD 11•••12 . . . AD 11•••1n AD 11•••21 AD 11•••22 . . . AD 11•••2n . . . AD 1n•••nn AD 2n•••nn . . . AD nn•••nn                           (26)
where, for any

(i 1 , i 2 , • • • , i k-1 ) ∈ [[1, n]] k-1 , the n × p i=1 n k i matrix D i1i2•••i k-1 reads D i1i2•••i k-1 = bdiag ℓ 1 i1 ⊗ • • • ⊗ ℓ 1 i k-1 ⊗ I n1 , ℓ 2 i1 ⊗ • • • ⊗ ℓ 2 i k-1 ⊗ I n2 , • • • , ℓ p i1 ⊗ • • • ⊗ ℓ p i k-1 ⊗ I np ( 27
)
with ℓ i k the k-th raw of matrix A i (see subsection 2.4). Proof. See Appendix D.
It directly results from lemma 3 and Eq. ( 24) that the cumulant matrix is made of n × n blocks

c r (x) i1,i2,••• ,ir-2 =    AD i1i3•••ir-3 C 2k (s)D ′ i2i4•••ir-2 A ′ if r = 2k AD i1i3•••ir-2 C 2k+1 (s)D ′ i2i4•••ir-3 A ′ if r = 2k + 1. (28)
One is now ready to state the main result of the paper.

Theorem 4 (joint block-diagonalization of contracted cumulant matrices).

Let s = (s ′ 1 , s ′ 2 , • • • , s ′ p )
′ be a random vector where the s i 's are mutually independent with dimensions dim(s i ) = n i , n = p i=1 n i . Let also A be an n × n orthogonal matrix such that x = As. Then, for any integer r > 2, the contracted cumulant matrix reads

Q x (r) [M (r-2) ] = A bdiag ∆ (r) s1 , ∆ (r) s2 , • • • , ∆ (r) sp A ′ (29)
with

(∀i ∈ [[1, p]]) ∆ (r) si =        Q si (2k) ⊗ k-1 A ′ i M (2k-2) ⊗ k-1 A i if r = 2k Q si (2k+1) ⊗ k A ′ i M (2k-1) ⊗ k-1 A i if r = 2k + 1 (30)
Proof. See Appendix E.

Theorem 4 clearly answers the questions raised in the introduction: contracted cumulant matrices are block-diagonalizable in the basis spanned by matrix A at any order r. In addition, it also returns the general structure of the diagonal blocks in terms of the cumulant matrices ∆ (r) si of the sources.

A particular consequence of Theorem 4 is given in the following corollary.

Corolary 5 (The ICA case). In the case of scalar sources, p = n and n i = 1 for any i ∈ [ [1, p]], the contracted cumulant matrix reads

Q x (r) [M (r-2) ] = A diag δ (r) s1 , δ (r) s2 , • • • , δ (r) sp A ′ (31)
where

(∀i ∈ [[1, n]]) δ (r) si =      ⊗ k-1 a ′ i M (2k-2) ⊗ k-1 a i k (2k) i if r = 2k ⊗ k a ′ i M(2k -1) ⊗ k-1 a i k (2k+1) i if r = 2k + 1 ( 32
)
with a i the i-th column of A and k

(r) i = c r (s i ) = cum(s i , • • • , s i ) the r-order cumulant of source s i .
In particular, when all sources are scalar and 2k = 4 (fourth-order statistics), corollary 5 returns the classical ICA result originally proved in [8, page 5],

Q x (4) M (2) = A diag k (4) 1 a ′ 1 M (2) a 1 , k (4) 2 a ′ 2 M (2) a 2 , • • • , k (4) p a ′ p M (2) a p A ′ (33)
with k (4) i the kurtosis of source s i . It is emphasized here that corollary 5 extends classical fourth-order ICA to any order r and theorem 4 generalizes it to any dimension. For instance, the block diagonal forms involved on orders r ∈ {3; 4; 5; 6} (which are to be used in the experimental section of the paper) are readily found as:

Q x (3) [M (1) ] = A bdiag Q s1 (3) A ′ 1 M (1) , • • • , Q sp (3) A ′ p M (1) A ′ (34)
Q x (4) [M (2) ] = A bdiag Q s1 (4) A ′ 1 M (2) A 1 , • • • , Q sp (4) A ′ p M (2) A p A ′ (35)
Q x (5) [M (3) ] = A bdiag Q s1 (5) (A ′ 1 ⊗ A ′ 1 ) M (3) A 1 , • • • , Q sp (5) (A ′ p ⊗ A ′ p ) M (3) A p A ′ (36)
Q x (6) [M (4) ] = A bdiag Q s1 (6) (A ′ 1 ⊗ A ′ 1 ) M (4) (A 1 ⊗ A 1 ) , • • • , Q sp (6) (A ′ p ⊗ A ′ p ) M (4) (A p ⊗ A p ) A ′ (
37) with M (1) , M (2) , M (3) and M (4) of dimensions n × 1, n × n, n 2 × n, and n 2 × n 2 , respectively.

SJADE r : a family of MICA algorithms

Theorem 4 makes it possible to propose a generalization of the JADE algorithm originally introduced by Cardoso and Soulomiac [START_REF] Cardoso | Blind beamforming for non gaussian signals[END_REF] to any dimension and to any order. Indeed, since it has been proven that contracted cumulant matrices Q -2) ] are block-diagonalizable in the basis spanned by matrix A for any order r, it suggests that joint block-diagonalization will effectively return the unknown mixing matrix in model [START_REF] Albera | Asymptotic performance of contrast-based blind source separation algorithms[END_REF]. Two families of JBD algorithms have actually been proposed in the literature:

x (r) [M (r
(a) algorithms that estimate an orthogonal block diagonalizer ( [START_REF] Févotte | Orthonormal approximate joint block-diagonalization[END_REF][START_REF] Theis | Towards a general independent subspace analysis[END_REF][START_REF] Maehara | Algorithm for error-controlled simultaneous block-diagonalization of matrices[END_REF]); (b) algorithms that avoid the whitening step by estimating an non-orthogonal block diagonalizer ( [START_REF] Nion | A tensor framework for nonunitary joint block diagonalization[END_REF][START_REF] Ghennioui | Gradient-based joint block diagonalization algorithms: Application to blind separation of fir convolutive mixtures[END_REF][START_REF] Ghennioui | A nonunitary joint block diagonalization algorithm for blind separation of convolutive mixtures of sources[END_REF]).

Both approaches are applicable to the results of this paper. However, for the sake of consistency with the assumption of section 2 (matrix A is orthogonal as a result of pre-whitening), the first strategy only will be considered from now on.

As for the set of cumulant matrices to be jointly block-diagonalized, a natural choice is to consider the contracting matrices

M (r-2) = E i1,i2,••• ,ir-2 which return Q x (r) [E i1,i2,••• ,ir-2 ] = [c r (x)] i1,i2,••• ,ir-2 .
The synopsis of the SJADE r algorithm that solves the MICA problem on order r is as follows.

Algorithm 1 (SJADE r ).

-Whiten the observations, -Construct of a set of contracted cumulant matrices: for instance, the set of n r-2 matrices,

M (r) = Q x (r) E i1,i2,••• ,ir-2 ; 1 i 1 , i 2 , • • • , i r-2 n , (38)
estimated from the data according the formulae worked out in the previous sections where empirical moment matrices are substituted for their theoretial versions [START_REF] Ould-Baba | Concise formulae for the cumulant matrices of a random vector[END_REF]. -Joint block-diagonalization (JBD): estimation of the orthogonal matrix A ′ which jointly blockdiagonalizes M (r) . -Separation of sources:

s = A ′ x.
It is noteworthy that the particular case on the fourth-order coincides with the so-called SJADE algorithm introduced in Ref. [START_REF] Theis | Towards a general independent subspace analysis[END_REF]. It should also be noted at this point that the set

Q x (r) [E i1,i2,••• ,ir-2 ] ; 1 i 1 , i 2 , • • • , i r-2 n
of contracted matrices is highly redundant due to the symmetries of the cumulant matrices c r (x). Specifically, it is seen from expression ( 23) that for any permutation σ ∈ S r-2 ,

Q x (r) [σ • M (r-2) ] ir-1,ir = n i1,i2,...,ir-2=1 cum(x i1 x i2 • • • x ir ) =cum(xi σ(1) •••xi σ(r-2) xi r-1 xi r ) m i σ(1) ,i σ(2) ,••• ,i σ(r-2) = Q x (r) [M (r-2) ] ir-1,ir . (39)
As a consequence, the set of matrices to be jointly block-diagonalized is indeed

M (r) = Q x (r) E i1,i2,••• ,ir-2 ; 1 i 1 i 2 • • • i r-2 n , (40)
with cardinality

card(M (r) ) = n + r -3 r -2 = n(n + 1) • • • (n + r -3) (r -2)! . (41)
For instance, on order 4, instead of block-diagonalizing the set of n 2 matrices (38), a wiser choice is to consider the set of n(n+1)

2 matrices M (4) = Q x (4) [E i,j ], 1 i j n (42)
where

E i,j =    e i e ′ i if i = j 1 √ 2 e i e ′ j + e j e ′ i if i < j . (43)
This basis was first proposed in Ref. [START_REF] Cardoso | High-order contrasts for independent component analysis[END_REF]. Similarly, on order 5, instead of block-diagonalizing the set of n 3 (38) matrices, it suffices to consider the set of n(n+1)(n+2) 6 matrices

M (5) = Q x (5) [E i,j,k ], 1 i j k n (44)
where

E i,j,k =                      e i ⊗ e i e ′ i if i = j = k 1 √ 3 (e k ⊗ e i e ′ i + e i ⊗ e k e ′ i + e i ⊗ e i e ′ k ) if i = j < k 1 √ 3 (e i ⊗ e k e ′ k + e k ⊗ e i e ′ k + e k ⊗ e k e ′ i ) if i < j = k 1 √ 6 e i ⊗ e j e ′ k + e i ⊗ e k e ′ j + e j ⊗ e i e ′ k + e j ⊗ e k e ′ i + e k ⊗ e i e ′ j + e k ⊗ e j e ′ i if i < j < k . ( 45
)
7. On the advantage of combining several orders Theorem 4 generalizes the JADE algorithm to any order and to any dimension. It also justifies other generalizations, such as the combination of orders. Indeed, since it has been proved that contracted cumulant matrices are block diagonalizable in the same basis whatever their orders, this makes possible to solve the MICA problem by considering the union of sets,

M = q l=1 M (r l ) , (46)
of contracted cumulant matrices at several orders r 1 , . . . , r q . Let us refer to the corresponding algorithm as SJADE r1,...,rq . The motivation beyond such an approach is to allow the separation of sources characterized by higher-order statistics of different orders. By way of an example, let us consider the case of independent scalar sources

s i (i = 1, • • • , 4) such that k (4) 1 = k (4) 2 = k (6) 3 = k
(6) 4 = 0, k (6)
1 = 0, k (6)
2 = 0, k (4)
3 = 0, and k

(4) 4 = 0, where k (4) i = cum(s i , s i , s i , s i ) et k (6) i = cum(s i , s i , s i , s i , s i , s i ).
According to corollary 5, the contracted cumulant matrices of orders 4 and 6 of the observations x = As, A ∈ R 4×4 , s = (s 1 , s 2 , s 3 , s 4 ) ′ , have expressions

Q x (4) [M (2) ] = A diag 0, 0, a ′ 3 M (2) a 3 k (4) 3 , a ′ 4 M (2) a 4 k (4) 4
A ′ (47)

Q x (6) [M (4) ] = A diag (a 1 ⊗ a 1 ) ′ M (4) (a 1 ⊗ a 1 )k (6) 1 , (a 2 ⊗ a 2 ) ′ M (4) (a 2 ⊗ a 2 )k (6) 2 , 0, 0 A ′ , (48)
respectively. It is clear that in such a scenario matrix A cannot be estimated neither from order 4 nor from order 6 alone because of the existence of some zero cumulants in each case (see assumption (H 2 ) of section 2.1). Nevertheless, resorting to the combination 4) ] of orders 4 and 6, one has

Q (4,6) = Q x (4) [M (2) ] + Q x (6) [M (
Q (4,6) = A diag (a 1 ⊗ a 1 ) ′ M (4) (a 1 ⊗ a 1 )k (6) 1 , (a 2 ⊗ a 2 ) ′ M (4) (a 2 ⊗ a 2 )k (6) 2 , a ′ 3 M (2) a 3 k (4) 3 , a ′ 4 M (2) a 4 k (4) 4
A ′ (49) which shows that matrix A can be estimated by block-diagonalizing contracted cumulant matrices of the type Q [START_REF] Albera | Blind identification of over-complete mixtures of sources (biome)[END_REF][START_REF] Cardoso | Multidimensional independent component analysis[END_REF] . This is illustrated with the numerical values k

Two set of contracted cumulant matrices, M (4) and M (6) , are then constructed according to Eq. ( 40), and the joint diagonalization algorithm (FFdiag) of Ref. Ziehe et al. [START_REF] Ziehe | A fast algorithm for joint diagonalization with non-orthogonal transformations and its application to blind source separation[END_REF] is used to return three estimates of the inverse of matrix A:

-B 4 , from the joint diagonalization of set M (4) ; -B 6 , from the joint diagonalization of set M (6) ; -B 4,6 , from the joint diagonalization of the union of sets M (4,6) = M (4) ∪ M (6) .

The separation results are assessed by the product B • A, • = (4), ( 6), [START_REF] Albera | Blind identification of over-complete mixtures of sources (biome)[END_REF][START_REF] Cardoso | Multidimensional independent component analysis[END_REF]:

B 4 A =    
-0.5109 -0.6400 -0.0000 0.0061 0.5266 0.6311 -0.0000 0.0007 0.0420 0.2293 0.0001 0.0000 -0.0479 -0.2329 0.0039 0.0000

    (51)
B 6 A =    
0.0027 -0.0000 -0.0631 -0.1051 -0.0000 0.0000 -0.3507 -0.0808 -0.0001 0.0000 0.6606 0.2647 -0.0000 0.0043 0.3823 0.0952

    (52) B 4,6 A =    
0.0000 -0.1792 -0.0000 -0.0000 0.0000 0.0000 -0.1813 0.0000 -0.1435 -0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 -0.0801

    . ( 53
)
As expected, orders 4 and 6 alone are unable to separate the 4 sources whereas their combination returns a perfect separation (up to the ICA undeterminacies).

Towards automatic sizing of sources

The previous section has proposed a family of algorithms apt to separate sources of different dimensions. This section now addresses the issue of determining the respective source dimensions from the observations, an objective of considerable importance in practice. An original solution is devised based on the joint diagonalization of the contracted cumulant matrices introduced previously. This is addressed here on order 4 for the sake of conciseness, yet it is easily generalizable to other orders. It proceeds from the concept of irreducibility of the sources.

Irreducibility of vector sources

It has been shown in Ref. [START_REF] Theis | Towards a general independent subspace analysis[END_REF][START_REF] Gutch | Independent subspace analysis is unique, given irreducibility[END_REF][START_REF] Gutch | Uniqueness of linear factorizations into independent subspaces[END_REF] that the irreducibility of sources is a fundamental assumption of MICA that guaranties the unicity of the solution. However, testing for the irreducibilty of vector sources from Definition 1 is not obvious. A simple but sufficient condition of irreducibility is proposed hereafter that makes use of the rank of the cumulant matrix. Lemma 6. Let y ∈ R d be a random vector satisfying the condition

r(c 4 (y)) = d(d + 1) 2 . ( 54
)
Then y is irreducible.

Proof. Appendix F Remark 1. Lemma 6 can be easily extended to other orders; for instance, on order 6, it reads

r(c 6 (y)) = d(d + 1)(d + 2) 6 . (55)
It must be highlighted that conditions (54) and (55) are sufficient, but not necessary. In particular, one may be fulfilled whilst the other one is not.

This leads to the following proposition.

Proposition 7. Suppose that vector sources s 1 , s 2 , • • • , s p are irreducible in the sense of lemma 6 and let V be the modal matrix in the eigenvalue decomposition VΛV ′ of the cumulant matrix c 4 (x). Then, matrix

R = n k=1 V ′ V k,k (56)
where

V ′ V k,k stands for the k-th diagonal block of dimension n × n in matrix V ′ V has factorization R = A      n1+1 2 I n1 0 • • • 0 0 n2+1 2 I n2 . . . 0 . . . . . . . . . . . . 0 0 • • • np+1 2 I np      A ′ . ( 57
)
Proof. Appendix F Proposition 7 means that the sought mixing matrix A that jointly block-diagonalizes the set M (4) also diagonalizes matrix R. The eigenvalues of R are then returned by the p positive quantities n1+1

2 , • • • , np+1 2
with respective multiplicities n 1 , • • • , n p . Two particular cases of proposition 7 are the following: for ICA, R = I n , whereas for k-MICA, R = k+1 2 I n . In both cases, matrix R alone does not allow the recovery of the unknown permutation of independent subspaces in general.

The implications of Proposition 7 are twofold. First, it allows the determination of unknown source dimensions in the mixture by inspecting the eigenvalues of matrix R. Second, it allows the grouping of source components by ordering the eigenvalues. This is now illustrated on a numerical example.

A numerical example

The following mixture is considered

x = A     s 1 s 2 s 3 s 4     , (58)
where sources s i , i = 1, • • • , 4 are independent, verify condition (54), and have dimensions dim(s 1 ) = dim(s 2 ) = 1, dim(s 3 ) = 2, and dim(s 4 ) = 3. The corresponding cumulant matrices are

c 4 (s 1 ) = 1.2331; c 4 (s 2 ) = -0.5013; c 4 (s 3 ) =     0.
              .
Note that the condition of irreducibility of the sources given in lemma 6 is perfectly satisfied since r(c 4 (s 1 )) = r(c 4 (s 2 )) = 1, r(c 4 (s 3 )) = 3, and r(c 4 (s 4 )) = 6. Beside, the mixing matrix A is the 7 × 7 orthogonal matrix,

A =           -0.
          .
These data are then used to compute the cumulant matrix c 4 (x) by means of Eq. 24, which is from now on considered as the only available observation. Note that r = r(c 4 (x)) = r(c 4 (s 1 )) + r(c 4 (s 2 )) + r(c 4 (s 3 )) + r(c 4 (s 4 )) = 1 + 1 + 3 + 6 = 11. The joint diagonalization algorithm of Ref. [START_REF] Cardoso | Blind beamforming for non gaussian signals[END_REF] is then applied to the contracted cumulant matrices M (4) and returns

E =           0.
         
. The next step is to compute the eigenvalue decomposition VΛV ′ of c 4 (x) with r = 11 eigen-elements and to construct matrix R from Eq. ( 56):

R =          
1.5813 0.0631 0.0909 0.4541 0.0386 -0.0503 -0.0253 0.0631 1.5314 0.0127 -0.0752 0.4556 0.0315 -0.0479 0.0909 0.0127 1.6059 -0.1254 -0.0334 0.0993 0.2499 0.4541 -0.0752 -0.1254 1.502 -0.0178 0.0223 0.0542 0.0386 0.4556 -0.0334 -0.0178 1.4471 0.1756 0.0038 -0.0503 0.0315 0.0993 0.0223 0.1756 1.5247 0.0546 -0.0253 -0.0479 0.2499 0.0542 0.0038 0.0546 1.8077

          .
Finally, using the above entries of E and R,

E ′ RE =           2 -0 -0 0 -0 0 0 -0 2 0 -0 0 -0 -0 -0 0 1.5 -0 0 -0 0 0 -0 -0 1 -0 0 -0 -0 0 0 -0 1 -0 0 0 -0 -0 0 -0 1.5 -0 0 -0 0 -0 0 -0 2           .
As indicated by proposition 7, matrix E ′ RE is made of 3 diagonal blocks with n 1 = n 2 = 1, n 3 = 2, and n 4 = 3. This returns the a priori unknown dimensions of the vector sources. Eventually, the permutation that rearrange the eigenvalues in decreasing order is

P =          
0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0

          such that P ′ E ′ A =          
-0.8289 0.4456 0.3382 0 0 -0 -0 -0.5545 -0.5748 -0.6017 -0 -0 0 0 -0.0737 -0.6863 0.7236 0 -0 -0 -0 0 0 0 -0.7239 0.6899 -0 -0 -0 -0 -0 0.6899 0.7239 0 0 0

0 0 -0 -0 -0 1 -0 0 -0 0 0 -1 -0          
, which correctly groups together the components of the independent vector sources. This proves that B = P ′ E ′ is a solution to the MICA problem up to the usual indeterminations.

Estimation in the presence of noise

The previous example has been designed to satisfy exactly the condition of irreducibility of the sources,

r(c 4 (x)) = ni(ni+1) 2
. This assumes the data are free of noise. In order to address the more realistic situation were noise is present, let us add to matrix c 4 (x) a small perturbation ǫc 4 (z), where z is a random vector independent of s and ǫ a small positive parameter:

c 4 = c 4 (x) + ǫc 4 (z). ( 59
)
In such a situation, one will generally have r(c 4 ) = n(n+1)

2 = ni(ni+1) 2
. After computing Ê by joint diagonalization, it remains to find the effective rank r of c 4 (x) as would be observed in the absence of noise. A natural idea is to try all values from r = n (the ICA case with n mutually independent sources) to r = n(n-1) 2 + 1 (2 independent sources with dimensions (n -1) and 1) and select that value which makes matrix Ê′ R r Ê the closest to a diagonal matrix. This can be tested by means of Amari's index [START_REF] Amari | A new learning algorithm for blind signal separation[END_REF] which takes values within 0 and 1: the closer it is to 0, the more diagonal the matrix.

The following example generates z from a standardized Gaussian and sets ǫ = 10 -2 ; the rank of the cumulant matrix is r(c 4 ) = r(c 4 (x) + ǫc 4 (z)) = n(n+1) 2 = 28. The corresponding Amari's indices are {0.038, 0.048, 0.022, 0.016, 0.001, 0.025, 0.035, 0.035, 0.040, 0.044, 0.048, 0.046, 0.040, 0.039, 0.033, 0.032} for r varying from n = 7 to n(n-1) 2 + 1 = 22. The minimum is found in the fifth entry, which corresponds to r = n + 4 = 11 and, accordingly,

Ê′ R 11 Ê =          
2.0000 0.0000 0.0098 -0.0027 0.0015 -0.0045 -0.0001 0.0000 2.0000 0.0054 -0.0016 -0.0003 0.0016 -0.0000 0.0098 0.0054 1.4996 -0.0012 0.002 0.0003 -0.0106 -0.0027 -0.0016 -0.0012 1.0002 -0.0013 -0.0015 -0.0010 0.0015 -0.0003 0.002 -0.0013 1.0002 -0.0004 0.0029 -0.0045 0.0016 0.0003 -0.0015 -0.0004 1.4998 0.0096 -0.0001 -0.0000 -0.0106 -0.001 0.0029 0.0096 2.0002

          .
It is seen that matrix Ê′ R 11 Ê is quite close to a diagonal matrix. After rearranging the diagonal elements in decreasing order, one finds the corresponding permutation matrix P, which finally returns the separation matrix

B = P ′ Ê′ =          
0.5301 0.4373 -0.2347 0.4274 0.4145 0.0581 -0.3388 -0.4454 0.5069 0.0356 -0.5085 0.4736 0.2288 -0.0902 0.2472 0.1609 0.5172 0.1539 0.1997 0.2564 0.7185 0.033 -0.2975 0.2746 0.0534 -0.083 0.8063 -0.4185 0.2815 0.1595 0.6997 -0.2666 -0.1303 -0.3953 -0.4019 0.4357 -0.5368 -0.144 -0.4975 0.4974 -0.0512 0.0606 -0.4351 -0.3518 0.3008 0.4622 0.5412 -0.2637 -0.1433

          .
The quality of the separation is assessed by the matrix product,

BA =          
-0.5544 -0.5749 -0.6018 0.0000 0.0002 -0.0003 0.0003 -0.0742 -0.6861 0.7237 0.0001 -0.0001 -0.0001 0.0006 -0.829 0.4459 0.3377 0.0001 0.0005 -0.0001 -0.0005 0.0005 -0.0002 -0.0001 0.6897 0.7241 -0.0001 -0.0004 0.0003 -0.0001 0.0000 -0.7241 0.6897 0.0005 -0.0029 -0.0002 0.0008 -0.0000 -0.0018 0.0023 0.0003 1.0000 0.0003 0.0002 0.0001 -0.0004 0.0003 -1.0000 0.0003

         
, which clearly singles out four independent subspaces of dimensions 3, 2, 1, and 1.

The synopsis of the algorithm is as follows.

Algorithm 2 (A simple Fourth-order Multidimensional ICA).

-Whiten the observations, -Construct of a set of contracted cumulant matrices: for instance, the set of n(n+1) 2 matrices,

M (4) = Q x (4) [E i1,i2 ] = E(x ′ E i1,i2 xxx ′ ) -2(E i1,i2 ) -tr(E i1,i2 )I n ; 1 i 1 i 2 n (60)
-Joint diagonalization (JD): estimation of the orthogonal matrix E ′ by JD of M (4) .

-Compute the eigenvalue decomposition

VΛV ′ of c 4 (x), where V ∈ R n 2 × n(n-1) 2 +1
.

-For r = n to r = n(n-1)

2 + 1 compute V r = v 1 • • • v r (where v i is the i-th column of V), compute V r V ′
r and R r (by means of Eq. 56), compute the matrices D r = E ′ R r E and the Amari-index of D r , noted error(r), -End.

-Select r 0 that minimizes error(r) and the corresponding D r0 .

-Estimate the sources dimensions.

-Compute the permutation matrix P by the eigenvalue decomposition of diag(diag(D r0 )).

-Separation of sources: s = P ′ E ′ x.

Experimental analyzes

This section now investigates the performance of the proposed SJADE r algorithm on two numerical examples.

Performance index for MICA

In order to measure the quality of the MICA separation, it is proposed to generalize the performance index originally introduced in Moreau [START_REF] Moreau | A generalization of joint-diagonalization criteria for source separation[END_REF]. The idea is to assess the proximity of G = BA to the product PD of a permutation matrix and a block-diagonal matrix. Assume without loss of generality that the estimated sources are sorted according to their dimensions, i.e.

n 1 ≤ n 2 ≤ • • • ≤ n p . Next, let us partition the n × n matrix G into blocks G ij of dimension n i × n j , with (i, j) ∈ [[1, p]] 2 . Then, define the p × p matrix G as G ij = 1 n i n j ni,nj k,l=1 |(G ij ) kl |. (61)
Having matrix G block-diagonal (for blocks of size

n 1 ≤ n 2 ≤ • • • ≤ n p ) is equivalent to having matrix G diagonal.
Therefore, the performance index I(G) for MICA is defined as

I(G) = 1 2p(p -1)   p i   p j G ij max l G il -1   + p j p i G ij max l G lj -1   . (62)
A perfect separation corresponds to I(G) = 0 (or -∞ on a logarithmic scale). Note that for ∀i, n i = 1, the proposed performance index reduces to the one introduced in [START_REF] Moreau | A generalization of joint-diagonalization criteria for source separation[END_REF] for scalar sources.

Comparisons 9.2.1. First simulation

This first example considers the case of three synthetic sources, of which two are vectors of dimension 2:

s 1 =
exp(2(0.007t + 0.5 -floor(0.007t + 0.5)) -1) 2(0.007t + 0.5 -floor(0.007t + 0.5)) -1

s 2 = exp(cos(0.3t)) cos(0.3t) (63)
s 3 = square(0.2t)
with t = 0, 1, • • • , 2999 and "floor" the operator that rounded down to the nearest whole number. Square mixing matrix A i , i = 1, • • • , 100, were randomly generated by sampling a uniform distribution in the interval [0, 1], so as to produce the observations x i = A i s. No noise was added. The JBD algorithm of Ref. [START_REF] Theis | Towards a general independent subspace analysis[END_REF] was used with threshold θ = 0.06. Performances of SJADE 3 , SJADE [START_REF] Theis | Towards a general independent subspace analysis[END_REF], SJADE 5 , SJADE 6 , MHICA, [START_REF] Theis | Blind signal separation into groups of dependent signals using joint block diagonalization[END_REF] and MSOBI [START_REF] Theis | Blind signal separation into groups of dependent signals using joint block diagonalization[END_REF] are compared in Fig. 1. It is seen that SJADE 6 evidences the best performance on this example.

-

SJADE_3 SJADE_4 SJADE_5 SJADE_6 MSOBI I (G) 50 -45 -40 -35 -30 -25 -20 -15 -10 MHICA
Figure 1: Box-plot of I(G) for SJADE 3 , SJADE [START_REF] Theis | Towards a general independent subspace analysis[END_REF], SJADE 5 , SJADE 6 , MHICA, [START_REF] Theis | Blind signal separation into groups of dependent signals using joint block diagonalization[END_REF] et MSOBI [START_REF] Theis | Blind signal separation into groups of dependent signals using joint block diagonalization[END_REF] (100 Monte-Carlo runs).

Second simulation

The second example is inspired from Ref. [START_REF] Póczos | Independent subspace analysis using k-nearest neighborhood distances[END_REF][START_REF] Theis | Towards a general independent subspace analysis[END_REF], where the random vector s = (s ′ 1 , s ′ 2 , s ′ 3 , s ′ 4 ) ′ , of dimension 8 is composed of 4 sources s 1 , s 2 , s 3 and s 4 , of dimensions 1, 3, 2 and 2, respectively. Source s 1 is computed from formula s 1 (ω) = exp(2(0.007ω + 0.5 -floor(0.007ω + 0.5)) -1),

whereas sources s 2 , s 3 and s 4 take the shape of a trihedron, a "β", and a "µ", respectively. The performance indices of the separations obtained with SJADE 3 , SJADE [START_REF] Theis | Towards a general independent subspace analysis[END_REF], SJADE 5 , SJADE 6 , MHICA, [START_REF] Theis | Blind signal separation into groups of dependent signals using joint block diagonalization[END_REF] and MSOBI [START_REF] Theis | Blind signal separation into groups of dependent signals using joint block diagonalization[END_REF] The separated sources are displayed in Fig. 3 for SJADE 3 to SJADE 6 . The product of the estimated separation matrix B with the mixing matrix A is shown in Fig. 4 for SJADE and SJADE 5 . It is clearly seen that B 4 and B 5 are equal up to a permutation, that is the indeterminacy of MICA.

Conclusion

The present work addressed the problem of multidimensional ICA in the most general setting where sources may have different dimensions. The objective was to propose a solution that does not make use of the strong assumptions of mutually independent scalar sources (ICA) or independent groups of components of equal size (k-MICA). In particular, it has been shown that contracted cumulant matrices are all blockdiagonalizable in the same basis. This opens the door to a family of algorithms coined SJADE r1,...,rq that jointly block-diagonalize a set of contracted cumulant matrices at various orders. One advantage is that different orders can be easily combined together according to the nature of the sources in order to improve their separation. Another result of the paper is to provide a solution to infer the dimensions of the sources when these are unknown. The idea proceeds from the concept of irreducibility of vector sources and leads to an algorithm that provides th e source dimensions as well as a permutation matrix that groups together the components of each vector source. The theoretical results of the paper have been illustrated by means of a few numerical examples, which have clearly evidenced the advantage of exploiting various higher orders. multilinearity of cumulants with x j = n i=1 a j,i s i for j ∈ [ [1, m]], where A = (a j,i ) (j,i)∈[ [1,m]]×[ [1,n]] :

c r (x) = m j1,j2,...,j 2k =1 cum(x j1 , • • • , x j 2k )E j1,j2,...,j 2k = m j1,j2,...,j 2k =1 cum n i1=1 a j1,i1 s i1 , • • • , n i 2k =1 a j 2k ,i 2k s i 2k E j1,j2,...,j 2k = m j1,j2,...,j 2k =1   n i1,i2,...,i 2k =1 a j1,i1 • • • a j 2k ,i 2k cum(s i1 , • • • , s i 2k )   E j1,j2,...,j 2k = n i1,i2,...,i 2k =1 cum(s i1 , • • • , s i 2k )   m j1,j2,...,j 2k =1 a j1,i1 • • • a j 2k ,i 2k E j1,j2,...,j 2k   = n i1,i2,...,i 2k =1 cum(s i1 , • • • , s i 2k ) m j1,j2,...,j 2k =1 a j1,i1 • • • a j 2k ,i 2k e j1 e ′ j2 ⊗ • • • ⊗ e j 2k-1 e ′ j 2k . (A.2)
The notation (e j ) j∈[ [1,m]] denotes the canonical basis of R m . Now, denoting

(f i ) i∈[[1,n]] the canonical basis of R n , it comes Af i = m j=1 a j,i e j , for any i ∈ [[1, n]]; thus, c r (x) = n i1,i2,...,i 2k =1 cum(s i1 , • • • , s i 2k )     m j1=1 a j1,i1 e j1     m j2=1 a j2,i2 e j2   ′ ⊗ • • • • • • ⊗   m j 2k-1 =1 a j 2k-1 ,i 2k-1 e j 2k-1     m j 2k =1 a j 2k ,i 2k e j 2k   ′   = n i1,i2,...,i 2k =1 cum(s i1 , • • • , s i 2k ) (Af i1 ) (Af i2 ) ′ ⊗ • • • ⊗ Af i 2k-1 (Af i 2k ) ′ = n i1,i2,...,i 2k =1 cum(s i1 , • • • , s i 2k ) (A ⊗ • • • ⊗ A) (f i1 f ′ i2 ⊗ • • • ⊗ f i 2k-1 f ′ i 2k ) (A ⊗ • • • ⊗ A) ′ = (A ⊗ • • • ⊗ A)   n i1,i2,...,i 2k =1 cum(s i1 , • • • , s i 2k )(f i1 f ′ i2 ⊗ • • • ⊗ f i 2k-1 f ′ i 2k )   (A ⊗ • • • ⊗ A) ′ = (A ⊗ • • • ⊗ A) c r (s) (A ⊗ • • • ⊗ A) ′ . (A.3)
Point (iv) is again proved for r = 2k only, the other case following similarly. One has

x = As = p j=1 A j s j (A.4)
where the independence of the sources s j implies the independence of the components A j s j . Using (ii) and

(iii), it comes

c r (x) = c r (As) = c r   p j=1 A j s j   = p j=1 c r (A j s j ) = p j=1 ⊗ k A j c r (s j ) ⊗ k A ′ j = ⊗ k A 1 • • • ⊗ k A p =A k    c r (s 1 ) • • • 0 . . . . . . . . . 0 • • • c r (s p )       ⊗ k A ′ 1 . . . ⊗ k A ′ p    A ′ k . (A.5)
where K pq , is the pq × pq commutation matrix (see Magnus and Neudecker [START_REF] Magnus | The commutation matrix: Some properties and applications[END_REF]) and the vec operator defined by vecA = vec a 1 . . .

a p =    a 1 . . . a p   .
Appendix C. The star product Definition 5. Let A be a matrix of size p × q and B a matrix of size rp × sq. The star product of A and B, noted A * B, is the matrix of size r × s defined as

A * B = p,q i1,i2=1 a i1,i2 B i1,i2 (C.1)
where B i1,i2 is the block indexed by (i 1 , i 2 ) in the partition of B into matrices of size r × s.

Some useful properties of the star product are listed below. (P3) Let A, B and C be three matrices of size p × q, q × r and r × s, respectively. Then Proof. Properties (P 1), (P 2) et (P 3) are given without proof in Ref. [START_REF] Macrae | Matrix derivatives with an application to an adaptive linear decision problem[END_REF]; a proof is found in [25, pp. 34-37]. Property (P 4) is proved hereafter because it will be used subsequently in the proof of theorem 4 (see Appendix E). Matrix M, of size pr × qs, is partitioned into blocks of size r × s as follows

ABC = B * vecAvec ′ C ′ = B ′ * (C ⊗ I p )K ps (A ⊗ I s ). (C.4) (
M =      M 1,1 M 1,2 • • • M 1,q M 2,1 M 2,2 • • • M 2,q . . . . . . . . . . . . M p,1 M p,2 • • • M p,q      . (C.6) If a = (a 1 , ..., a p ) ′ and b = (b 1 , ..., b q ) ′ , one has (a ′ ⊗ I r ) M (b ⊗ I s ) = a 1 I r a 2 I r • • • a p I r      M 1,1 M 1,2 • • • M 1,q M 2,1 M 2,2 • • • M 2,q . . . . . . . . . . . . M p,1 M p,2 • • • M p,q           b 1 I s b 2 I s . . . b q I s      = p i=1 a i M i,1 p i=1 a i M i,2 • • • p i=1 a i M i,p      b 1 I s b 2 I s . . . b q I s      = b 1 p i=1 a i M i,1 + b 2 p i=1 a i M i,2 + • • • + b q p i=1 a i M i,p = p i=1 q j=1 a i b j M i,j = ab ′ * M. (C.7)
Appendix D. Proof of lemma 3

For the sake of simplicity, the proof is given for k ∈ {2, 3} only, the general proof thus resulting by mathematical induction.

Case k = 2. One has

A 2 = A 1 ⊗ A 1 A 2 ⊗ A 2 . . . A p ⊗ A p =      B 1 B 2 . . . B n      (D.1)
where, for i ∈ [ [1, n]], matrix B i of size n × p j=1 n 2 j has expression

B i = ℓ 1 i ⊗ A 1 ℓ 2 i ⊗ A 2 • • • ℓ p i ⊗ A p , (D.2)
with ℓ j i , of size 1 × n j , is the i-th raw of matrix A j (voir ( 8)). Making use of the properties of the Kronecker product, it comes

ℓ j i ⊗ A j = (1ℓ j i ) ⊗ (A j I nj ) = A j ℓ j i ⊗ I nj . (D.3)
Inserting (D.3) into (D.2), one gets

B i = A 1 ℓ 1 i ⊗ I n1 A 2 ℓ 2 i ⊗ I n2 • • • A p ℓ p i ⊗ I np = A 1 A 2 • • • A p =A      ℓ 1 i ⊗ I n1 0 n1×n 2 2 . . . 0 n1×n 2 p 0 n2×n 2 1 ℓ 2 i ⊗ I n2 . . . 0 n2×n 2 p . . . . . . . . . . . . 0 np×n 2 1 0 np×n 2 2 . . . ℓ p i ⊗ I np      =Di = A bdiag ℓ 1 i ⊗ I n1 , ℓ 2 i ⊗ I n2 , . . . , ℓ p i ⊗ I np . (D.4) Case k = 3. A 3 = A 1 ⊗ A 1 ⊗ A 1 A 2 ⊗ A 2 ⊗ A 2 • • • A p ⊗ A p ⊗ A p =                           B 11 B 12 . . . B 1n B 21 B 22 . . . B 2n . . . B n1 B n2 . . . B nn                           =      B 1 B 2 . . . B n      (D.
B i =      B i1 B i2 . . . B in      = ℓ 1 i ⊗ A 1 ⊗ A 1 ℓ 2 i ⊗ A 2 ⊗ A 2 • • • ℓ p i ⊗ A p ⊗ A p = (A 1 ⊗ A 1 )(ℓ 1 i ⊗ I n 2 1 ) • • • (A p ⊗ A p )(ℓ p i ⊗ I n 2 p ) = A 2 bdiag ℓ 1 i ⊗ I n 2 1 , ℓ 2 i ⊗ I n 2 2 , • • • , ℓ p i ⊗ I n 2 p . (D.6)
Replacing A 2 by its expression given in case k = 2,

B i =      AD 1 AD 2 . . . AD n      bdiag ℓ 1 i ⊗ I n 2 1 , ℓ 2 i ⊗ I n 2 2 , • • • , ℓ p i ⊗ I n 2 p , (D.7)
which implies

B ij = AD j bdiag ℓ 1 i ⊗ I n 2 1 , ℓ 2 i ⊗ I n 2 2 , • • • , ℓ p i ⊗ I n 2 p = A bdiag ℓ 1 j ⊗ I n1 , • • • , ℓ p j ⊗ I np bdiag ℓ 1 i ⊗ I n 2 1 , • • • , ℓ p i ⊗ I n 2 p = A bdiag 1 ⊗ ℓ 1 j ⊗ I n1 , • • • , 1 ⊗ ℓ p j ⊗ I np bdiag ℓ 1 i ⊗ I n1 ⊗ I n1 , • • • , ℓ p i ⊗ I np ⊗ I np = A bdiag ℓ 1 i ⊗ ℓ 1 j ⊗ I n1 , ℓ 2 i ⊗ ℓ 2 j ⊗ I n2 , • • • , ℓ p i ⊗ ℓ p j ⊗ I np =Dij . (D.8)
Appendix E. Proof of theorem 4

Again the proof is given for r = 2k only, the case r = 2k + 1 following similar lines. Let replace A k in c 2k (x) = A k C 2k (s)A ′ k by its expression given in lemma 3. It is then seen that the block of size n × n in matrix c 2k (x) indexed by (i 1 , j 1 , • • • , i k-1 , j k-1 ) has expression

c 2k (x) i1,j1,••• ,i k-1 ,j k-1 = AD i1•••i k-1 C 2k (s)D ′ j1•••j k-1 A ′ . (E.1)
Using definition [START_REF] Hyvärinen | Topographic independent component analysis[END_REF] of the contracted cumulant matrix of order 2k, it comes

Q x (2k) [M] = n i1,j1,••• ,i k-1 ,j k-1 =1 m i1,j1,••• ,i k-1 ,j k-1 c 2k (x) i1,j1,••• ,i k-1 ,j k-1 = n i1,j1,••• ,i k-1 ,j k-1 =1 m i1,j1,••• ,i k-1 ,j k-1 AD i1•••i k-1 C 2k (s)D ′ j1•••j k-1 A ′ = A   n i1,j1,••• ,i k-1 ,j k-1 =1 m i1,j1,••• ,i k-1 ,j k-1 D i1•••i k-1 C 2k (s)D ′ j1•••j k-1   =∆ A ′ . (E.2) It is noted that ∆ is block-diagonal since D i1•••i k-1 , C 2k (s) and D ′ j1•••j k-1 are block-diagonal. Thus ∆ = n i1,j1,••• ,i k-1 ,j k-1 =1 m i1,j1,••• ,i k-1 ,j k-1 D i1•••i k-1 C 2k (s)D ′ j1•••j k-1 = bdiag ∆ (2k) s1 , • • • , ∆ (2k) sp , (E.3)
where, upon invoking ( 25) and [START_REF] Lahat | Optimal performance of second-order multidimensional ica[END_REF],

∆ (2k) si = n i1,j1,••• ,i k-1 ,j k-1 =1 m i1,j1,••• ,i k-1 ,j k-1 ℓ i i1 ⊗ • • • ⊗ ℓ i i k-1 ⊗ I ni c 2k (s i ) ℓ i j1 ⊗ • • • ⊗ ℓ i j k-1 ⊗ I ni ′ .
(E.4) Making use of property (C.5) of the star product,

∆ (2k) si =   n i1,j1,••• ,i k-1 ,j k-1 =1 m i1,j1,••• ,i k-1 ,j k-1 ℓ i i1 ⊗ • • • ⊗ ℓ i i k-1 ′ ℓ i j1 ⊗ • • • ⊗ ℓ i j k-1   * c 2k (s i ) =   n i1,j1,••• ,i k-1 ,j k-1 =1 m i1,j1,••• ,i k-1 ,j k-1 A ′ i e i1 ⊗ • • • ⊗ A ′ i e i k-1 (e ′ j1 A i ⊗ • • • ⊗ e ′ j1 A i )   * c 2k (s i ) =   n i1,j1,••• ,i k-1 ,j k-1 =1 m i1,j1,••• ,i k-1 ,j k-1 (⊗ k-1 A ′ i ) e i1 ⊗ • • • ⊗ e i k-1 (e ′ j1 ⊗ • • • ⊗ e ′ j1 ) (⊗ k-1 A i )   * c 2k (s i ) =   (⊗ k-1 A ′ i )   n i1,j1,••• ,i k-1 ,j k-1 =1 m i1,j1,••• ,i k-1 ,j k-1 E i1,j1,••• ,i k-1 ,j k-1   (⊗ k-1 A i )   * c 2k (s i ) = (⊗ k-1 A ′ i )M (2k-2) (⊗ k-1 A i ) * c 2k (s i ) = Q si (2k) ⊗ k-1 A ′ i M (2k-2) ⊗ k-1 A i . (E.5)
Lemma 9. Let y be a random vector verifying (54) and UΛU ′ be the eigenvalue decomposition of c 4 (y), where Λ = diag(λ 1 , λ 2 , • • • , λ d(d+1)

2

). Then,

UU ′ = 1 2 (K dd + I d 2 ), (F.4)
with K dd the commutation matrix (see [START_REF] Magnus | The commutation matrix: Some properties and applications[END_REF]).

Proof of lemma 9. Since c 4 (y) is a symmetric matrix with rank d(d+1) , and UU ′ KUU ′ = UU ′ . Therefore UU ′ is t he pseud o-inverse of K which is equal to K.

Proof of proposition 6. Assume that all multidimensional sources fulfill condition (54), i.e. r(c 4 (s i )) = ni(ni+1) 2

. It then follows that none of them is Gaussian since r(c 4 (s i )) = 0. In addition, it follows from Eq. ( 24) that

0 • • • Λ p    U ′ D ′ 1 A ′ • • • U ′ D ′ n A ′ . (F.8)
One recognizes the eigenvalue decomposition of c 4 (x) = VΛV ′ with

V =      AD 1 U AD 2 U . . . AD n U      . (F.9)
The next step is to evaluate the product,

V ′ V = n k=1 U ′ D ′ k A ′ A =In D k U = U ′ n k=1 D ′ k D k U.
(F.10)

  3614 and the mixing matrixA = a 1 a 2 a 3 a 4 =

Figure 2 .

 2 (a) displays 3000 samples of the sources and Fig. 2.(b) their observations after mixing with a matrix A ∈ R 8×8 whose entries are sampled in a uniform distribution on [0, 1].

  Projection of mixed sources on 4 planes.

Figure 2 :

 2 Figure 2: MICA, second simulation.

  Separated sources with SJADE 6 .

Figure 3 :

 3 Figure 3: Separated sources (MICA, second simulation)

Figure 4 :

 4 Figure 4: Product of the estimated separation matrix with the mixing matrix (MICA, second simulation)

Proposition 8 .

 8 (P1) Let A and B be two matrices of same dimension. Then A * B = tr(A ′ B). (C.2) (P2) Let A, B and C be three matrices of size p × q, pr × sq and m × n, respectively. Then A * (B ⊗ C) = (A * B) ⊗ C. (C.3)

  P4) Let a and b be two column vectors of dimensions p and q, respectively, and let M be a matrix of size pr × qs. Then (a ′ ⊗ I r ) M (b ⊗ I s ) = ab ′ * M. (C.5)

2 , 2 ) 2 .

 222 it accepts the eigenvalue decomposition UΛU ′ , whereU ∈ R d 2 × d(d+1) 2 is a semi-orthogonal matrix (U ′ U = I d(d+1) matrix made of the d(d+1) 2 non-zero eigenvalues of c 4 (y). let K = 1 2 (K dd + I d 2 ). Since Kc 4 (y) K = c 4 (y), one has KUΛU ′ K = UΛU ′ . Thus, it comes U ′ K KU = I d(d+1) Since K is idempotent ( K K = K), it comes U ′ KU = I d(d+1) 2

1 A 1 2

 11 ′ • • • D ′ n A ′ . (F.5) Let us now introduce U i Λ i U ′ i , i = 1, • • • , p, the eigen-elements of c 4 (s i ) of lemma 9, with U i U ′ i = (K nini + I n 2 i ). Thus U = bdiag U 1 , • • • , U p (F.6) and Λ = bdiag Λ 1 , • • • , Λ p ,

Table 1 :

 1 are displayed in Table1. It is seen that SJADE 5 achieves the best performance in this case. Performance index (MICA, seconde simulation)

	Méthode MSOBI SJADE 3 SJADE 4 SJADE 5 SJADE 6 MHICA
	I(G)[dB]	-4.05	-26.76	-30.10	-32.09	-24.21	-24.95

Appendix A. Proof of proposition 1

Point (i) is a direct consequence of the property that ∀r > 2, cum(z 1 , z 2 , • • • , z r ) = 0, for a Gaussian vector. Point (ii) is a direct consequence of the property (see e.g. [35, p. 280]) that

for independent random vectors x and y. Point (iii) is proved here for r = 2k only, the case r = 2k + 1 following similar lines. It suffices to use the

Appendix B. Cumulant matrices

Let x be a real-valued random vector of dimension p × 1, with assumed zero-mean for simplicity (E [x] = 0). Let us further assume that E [|x| n ] < ∞ for a given integer n ≥ 6 so as to guarantee the existence of moments and cumulants up to order six. Next let us define the four matrices:

which depend only on dimension p. Thus, the cumulant matrices of orders 2 to 6 read

respectively, where

Proof of lemma 6. Let us give a proof by contraposition. Let assume y is reducible, then according to definition 1, there exists an invertible matrix A = A 1 A 2 ∈ R d×d and independent vectors y 1 ∈ R d1 and

This implies from ( 14) that

holds true since r(A) = r(A 1 ) + r(A 2 ) (see [START_REF] Tian | Some rank equalities and inequalities for kronecker products of matrices[END_REF]). Thus, the majoration of rank (r c 4 (y i ) di(di+1) 2

, ∀i ∈ {1, 2} ). This means that r c 4 (y) < d(d+1)

when y is reducible. Therefore, if r(c 4 (y)) = d(d+1)

2

, then y is irreducible. Substituting Eq. ( 27) for D k , one has

Besides, from lemma 9,

so that,

) and, finally upon using property (P4) of the star product (see (Appendix C)),


Andreas Fischer

Kieu Van Cuong

Muriel Visani
email: muriel.visani@univ-lr.fr

Ching Y Suen
email: suen@encs.concordia.ca

Generation of Learning Samples for Historical Handwriting Recognition Using Image Degradation

Keywords: Historical Documents, Handwriting Recognition, Learning Sample Generation, Synthetic Images, Image Degradation, Hidden Markov Models

Historical documents pose challenging problems for training handwriting recognition systems. Besides the high variability of character shapes inherent to all handwriting, the image quality can also differ greatly, for instance due to faded ink, ink bleed-through, wrinkled and stained parchment. Especially when only few learning samples are available, it is difficult to incorporate this variability in the morphological character models. In this paper, we investigate the use of image degradation to generate synthetic learning samples for historical handwriting recognition. With respect to three image degradation models, we report significant improvements in accuracy for recognition with hidden Markov models on the medieval Saint Gall and Parzival data sets.

INTRODUCTION

The interest in handwriting recognition for historical documents has been growing rapidly in recent years [START_REF] Antonacopoulos | Special issue on the analysis of historical documents[END_REF]. Automatic reading of historical manuscripts would provide access to the contents of a vast amount of digitized historical documents stored at libraries worldwide. However, current handwriting recognition systems are still far from being perfect. Typical problems include the high variability in character shapes, the large number of word classes, and the inability to segment touching and broken characters prior to recognition [START_REF] Koerich | Large vocabulary off-line handwriting recognition: A survey[END_REF]. For historical documents, an additional problem are differences in the image quality stemming for instance from faded ink, ink bleed-through, wrinkled and stained parchment. While this problem has been mainly addressed in the context of layout analysis [START_REF] Likforman-Sulem | Text line segmentation of historical documents: A survey[END_REF], it has also an effect on learning-based recognition systems. Indeed, the morphological character models have to incorporate different degradations at learning stage in order to generalize well to unseen test images.

Image degradation is a promising approach to distort available learning samples in order to enlarge the training set and incorporate more variability in the image quality. Diverse image degradation models have been proposed in document analysis [START_REF] Baird | The State of the Art of Document Image Degradation Modeling[END_REF], for example barcode printing defects [START_REF] Loce | Halftone Banding due to Vibrations in A Xerographic Image Bar Printer[END_REF], ink bleed-through [START_REF] Moghaddam | Low Quality Document Image Modeling and Enhancement[END_REF], and "hard pencil noise" [START_REF] Jian Zhai | A Line Drawings Degradation Model for Performance Characterization[END_REF]. A typical application scenario is the evaluation of document analysis systems under different degrees of distortion, for instance to assess symbol recognition and spotting [START_REF] Delalandre | Generation of Synthetic Documents for Performance Evaluation of Symbol Recognition & Spotting Systems[END_REF], music score staff removal [START_REF] Visani | The ICDAR/GREC 2013 Music Scores Competition on Staff Removal[END_REF], and handwriting recognition [START_REF] Mori | Generating New Samples from Handwritten Numerals Based on Point Correspondence[END_REF].

In our application scenario, image degradation is used to generate learning samples. This approach has been pursued in [START_REF] Varga | Effects of Training Set Expansion in Handwriting Recognition Using Synthetic Data[END_REF][START_REF] Varga | Generation of Synthetic Training Data for an HMM-based Handwriting Recognition System[END_REF] for modern handwriting. Several perturbation models were applied to the handwriting images to distort the shape of the characters and the baseline of the text lines. Promising improvements of the recognition accuracy are reported with synthetic learning samples. However, aiming at modern handwritings, this degradation model does not provide the means to generate typical noise observed in digitized historical manuscripts.

In this paper, we investigate the case of historical handwriting based on three degradation models, namely Kanungo noise [START_REF] Kanungo | Global and Local Document Degradation Models[END_REF], a character degradation model [START_REF] Kieu | A Character Degradation Model for Grayscale Ancient Document Images[END_REF], and a geometric distortion model [START_REF] Liang | Geometric Rectification of Camera-Captured Document Images[END_REF]. They are applied to binary text line images to generate additional learning samples for hidden Markov model based recognition. On the medieval Saint Gall [START_REF] Fischer | Transcription alignment of latin manuscripts using hidden Markov models[END_REF] and Parzival [START_REF] Fischer | Lexicon-free handwritten word spotting using character HMMs[END_REF] data sets, we demonstrate that the recognition accuracy can be significantly improved with the synthetic learning samples.

The remainder of this paper is organized as follows. First, the data sets are introduced in Section 2. Next, the image degradation models are detailed in Section 3. Section 4 provides information about the recognition system and Section 5 presents the experimental results. Finally, we draw some conclusions in Section 6. 1a and1b.

The images illustrate considerable differences in the image quality even within a single manuscript.

Saint Gall Database

The Saint Gall database is based on a medieval Latin manuscript from the 9th century that contains the hagiography Vita sancti Galli by Walafrid Strabo. The Abbey Library of Saint Gall, Switzerland, holds a manual copy of the work within the Cod. Sang. 562, which was written by a (probably) single experienced hand in Carolingian script with ink on parchment.

The database currently includes 60 manuscript pages as indicated in Table 1. For further details on this database, we refer to [START_REF] Fischer | Transcription alignment of latin manuscripts using hidden Markov models[END_REF].

Parzival Database

The Parzival database is based on a medieval German manuscript from the 13th century that contains the epic poem Parzival by Wolfram von Eschenbach, one of the most significant epics of the European Middle Ages. There exist several manual copies of the poem that differ in writing style and dialect of the language. For the Parzival database, the Cod. 857 is considered, which is held by the Abbey Library of Saint Gall, Switzerland. It was written in Middle High German by several writers with ink on parchment using Gothic minuscules.

The Parzival database currently consists of 47 manuscript pages as indicated in Table 1. For further details on this database, we refer to [START_REF] Fischer | Lexicon-free handwritten word spotting using character HMMs[END_REF].

IMAGE DEGRADATION MODELS

In this section, the three image degradation models are presented that are used to generate synthetic learning samples. They are applied to binary text line images after a series of image preprocessing steps.

The resulting degradation images have been made publicly available.

Preprocessing

The first step in preprocessing scanned manuscript images for handwriting recognition consists of layout analysis and text line extraction. In order to focus on handwriting recognition, we work directly with pre-segmented text line images in this paper. They were extracted with a semi-automatic procedure detailed in [START_REF] Fischer | Ground truth creation for handwriting recognition in historical documents[END_REF].

Next, the text foreground is determined by means of binarization. For the Saint Gall database, Sauvola's method is employed [START_REF] Fischer | Transcription alignment of latin manuscripts using hidden Markov models[END_REF] and for the Parzival database, local edge enhancement is performed with Difference of Gaussians before applying a global threshold [START_REF] Fischer | Ground truth creation for handwriting recognition in historical documents[END_REF].

Finally, the binary images undergo a series of normalization operations in order to standardize their appearance before recognition. Normalization includes skew correction as well as a vertical division into an upper, middle, and lower region.

An exemplary normalization result is shown in Figures 2b and3b. For further details, we refer to [START_REF] Fischer | Ground truth creation for handwriting recognition in historical documents[END_REF].

Kanungo Noise

The first degradation model was proposed by Kanungo et al.

in [START_REF] Kanungo | Global and Local Document Degradation Models[END_REF] and was validated in [START_REF] Kanungo | A statistical, Nonparametric Methodology for Document Degradation Model Validation[END_REF]. The model distorts character edges by using a non-linear local selection process and a morphological closing operation. First, the selection process flips the value of some foreground and background pixels (foreground pixel changes to background and vice versa) according to a probability function. The probability function depends on a distance transform as follows:

• Foreground pixels: p = α0e -αd + η

• Background pixels: p = β0e -βd + η where d is the distance of a pixel to the nearest character edge. The two input parameters α and β control the number of flipped pixels. The three parameters α0, β0, and η are considered as the initial inputs. Secondly, a morphological closing operation is applied to smooth the character edges. This model works well with binary images and is widely used to simulate the presence of noise for assessing the performances of different document analysis methods.

We consider three degradation levels with increasing distortion effect:

• Level 1: α = β = 7, α0 = β0 = 1, η = 0 • Level 2: α = β = 5.5, α0 = β0 = 1, η = 0 • Level 3: α = β = 4.5, α0 = β0 = 1, η = 0
The effects of the different degradation levels are illustrated in Figures 2c and3c. The number of flipped pixels ranges from 100-200 for level 1, 400-800 for level 2, and 800-1600 for level 3.

Character Degradation

The second degradation model is the character degradation model proposed in [START_REF] Kieu | A Character Degradation Model for Grayscale Ancient Document Images[END_REF], which is specifically designed to As detailed in [START_REF] Kieu | A Character Degradation Model for Grayscale Ancient Document Images[END_REF], this noise can be simulated in three steps. First, the seed-points (centers of noise regions) are selected as the flipped pixels resulting from the non-linear local selection process proposed by Kanungo et al. [START_REF] Kanungo | Global and Local Document Degradation Models[END_REF]. Secondly, at each seed-point, a noise region is defined with three properties: shape, size, and direction. For example, an elliptic noise region has two axes (minor and major) which depend on the gradient vector -→ v at its center and on two input parameters a0 and g. The parameter a0 controls the size of the elliptic region and g its "flatness". The size parameter a0 determines whether noise regions break, connect, or modify character strokes. Finally, in order to obtain degraded regions as realistic as possible, the gray value of the pixels inside each noise region is smoothed with a Gaussian function. Since we are working with binary images rather than grayscale images, this last step is omitted.

Again, we define three degradation levels with increasing distortion effect:

• The effect is shown in Figures 2d and3d. The number of pixels in the elliptic noise regions ranges from 50-100 for level 1, 100-200 for level 2, and 200-400 for level 3.

Geometric Distortion

The third degradation model simulates a geometric surface deformation. It was proposed by J. Liang et al. in [START_REF] Liang | Geometric Rectification of Camera-Captured Document Images[END_REF] to evaluate their restoration algorithm. First, this model maps every point in the plane of the original image to a point on a curved surface. The curved surface is then divided into into planar strips (e.g. quadrilaterals) that are mapped back to the image plane.

We employ this degradation model with different developable surfaces such as sinusoidal and parabolic. The distortion levels are controlled by the amplitude a and wavelength λ of the surface:

• Level 1: sinusoidal surface, a = 10, λ = 0.5

• Level 2: sinusoidal surface, a = 15, λ = 0.25

• Level 3: parabolic surface, a = 20, λ = 0.2
The effect is illustrated in Figures 2e and3e.

RECOGNITION SYSTEM

In this section, we briefly describe the hidden Markov model based recognizer used in the experiments. For a full description, we refer to [START_REF] Marti | Using a statistical language model to improve the performance of an HMM-based cursive handwriting recognition system[END_REF][START_REF] Fischer | Lexicon-free handwritten word spotting using character HMMs[END_REF].

Features

Using a sliding analysis window with a width of 1 pixel moving from left to right, the binary text line images are represented by a sequence of feature vectors. Each window captures nine geometric features including the center of gravity, the second order moment, the fraction of black pixels, the contour positions, the deviation at the contours, and the number of black-white transitions. For more details, we refer to [START_REF] Marti | Using a statistical language model to improve the performance of an HMM-based cursive handwriting recognition system[END_REF].

Training

Morphological character models consist of a series of states arranged in a linear topology.

Recognition

Two goal functions are optimized during recognition with the Viterbi algorithm [START_REF] Rabiner | A tutorial on hidden Markov models and selected applications in speech recognition[END_REF]. First, the best match of the image with the trained morphological character models is established. Secondly, the most likely sequence of words is found with respect to a bigram language model that describes the probability of a word given the preceding word.

The two goal functions are balanced with a grammar scale factor and a word insertion penalty [START_REF] Zimmermann | N-gram language models for offline handwritten text recognition[END_REF].

We estimate the word bigrams on all training and validation text lines and smooth them with the modified Kneser-Ney method [START_REF] Kneser | Improved backing-off for m-gram language modeling[END_REF] to cope with unseen word bigrams. Furthermore, we use a closed vocabulary that contains all words of the database, even those of the test set. It is a common practice to establish comparable results without out-ofvocabulary words. For the Saint Gall database, the recognition vocabulary contains 5, 762 words and for the Parzival database, the vocabulary contains 4, 934 words.

EXPERIMENTAL EVALUATION

We have conducted an experimental evaluation of the proposed image degradation models on the Saint Gall and Parzival data sets (see Section 2).

All three models, namely Kanungo noise, character degradation, and geometric distortion, are applied to the training set individually to create one distorted version of each learning sample. In effect, the size of the training set is doubled.

It is used to train a new recognition system, which is then In addition to an individual model evaluation, we also combine up to two noise models by adding their degradation images to the same training set, which is then three times larger than the original training set.

Setup

First, the text lines of both data sets are split into disjoint sets for training the HMM, validation of system parameters, and testing the final performance. The number of text lines are indicated in Table 2 for each set.

For the image degradation models, the best degradation level L ∈ {1, 2, 3} is determined on the validation set. Also, the best pair of degradation models for the model combination is found with respect to the validation accuracy.

System parameters of the HMM recognizer, which are optimized on the validation set, include the number of training iterations I ∈ {2, 3, 4, 5} per training epoch and the final number of Gaussian mixtures G ∈ {5, 10, . . . , 30}. The number of HMM states and the language model parameters have been adopted from previous experiments [START_REF] Wüthrich | Language model integration for the recognition of handwritten medieval documents[END_REF].

Results

The word accuracy results are indicated in Table 3 for each individual image degradation model as well as for the best combination of two degradation models. Next to the word accuracy, the improvement over the reference system is indicated. With the exception of geometric distortions on the Saint Gall database, 7 out of 8 improvements over the reference system are statistically significant (t-test, α = 0.05). Table 4 lists the relative word error reduction for all models. The reductions are also shown in Figure 4. The synthetic learning samples realize an error reduction of up to 16.53% on the Saint Gall database, achieving an accuracy of 90.81%. On the Parzival database, an error reduction of up to 20.05% is achieved and an accuracy of 87.12% is reported. These re- sults clearly demonstrate the benefit of generating synthetic learning samples using image degradation.

When comparing the individual models, the character degradation model and the Kanungo noise tend to outperform the geometric distortion. Table 5 lists the optimal degradation levels of the different models. In 4 out of 6 cases, the lowest degradation level was pefered.

For both data sets, the combination of noise models achieves the best result. For the Saint Gall database, the optimal pair of models was (Character, Geometric) and for the Parzival database, the optimal pair was (Kanungo, Character).

The combination results demonstrate that multiple sources of degradations are desirable for the synthetic generation of learning samples.

CONCLUSIONS

In this paper, we have investigated image degradation models to generate synthetic learning samples for historical handwriting recognition. The three proposed degradation models include Kanungo noise, character degradation, and geometric distortion. They have been applied to distort binary text line images for training a hidden Markov model based recognizer.

On the medieval Saint Gall and Parzival data sets, significant error reductions could be achieved with the synthetic learning samples. For closed-vocabulary recognition, the best word accuracy results were 90.81% for the Saint Gall database and 87.12% for the Parzival database. This corresponds to an error reduction over the baseline system of 16.53% and 20.05%, respectively. For both data sets, the best results were achieved by a combination of different degradation models.

In the context of historical manuscripts, the importance of learning samples cannot be overemphasized. Especially when bootstrapping a new handwriting recognition system for a new kind of historic script, every learning sample counts. Image degradations have shown a promising potential to incorporate more variability in the training set, which leads to a better generalization capability of the recognizer.

Future work includes improvements in the degradation models, in particular the investigation of realistic 3-dimensional distortions of the writing support. The goal is to create realistic manuscript images that incorporate several types of degradations.

2

 2

1

  Figure 1: Data Sets

Figure 2 :

 2 Figure 2: Image Degradation Saint Gall

Figure 3 :

 3 Figure 3: Image Degradation Parzival

5 •

 5 Level 1: α = β = 8.5, α0 = β0 = 1, η = 0 g = 0.6, a0 = Level 2: α = β = 7, α0 = β0 = 1, η = 0 g = 0.6, a0 ∈ {3, 4, . . . , 7} (random selection) • Level 3: α = β = 7, α0 = β0 = 1, η = 0 g = 0.6, a0 ∈ {3, 4, . . . , 10} (random selection)

Figure 4 :

 4 Figure 4: Error Reduction.

Table 1 :

 1 Data Sets Two medieval data sets of the publicly available IAM-HistDB 1 are used in this paper, namely the Saint Gall database and the Parzival database. Data set statistics are listed in Table 1 and sample images are shown in Figures

		Saint Gall Parzival
	Century Language	9th Latin	13th German
	Writers Pages Text Lines	1 60 1,410	3 47 4,477
	2. DATA SETS

  Each state captures the distribution of the sliding window features with a mixture of Gaussians with diagonal covariance matrices. Training is performed with the Baum-Welch algorithm [19] based on known transcriptions of text line images. We employ a training procedure that increases the number of Gaussians G incrementally. Starting with G = 1, the training set is processed I times in the first training epoch. Then, the number of Gaussians is incremented to G = 2 by splitting and the training set is again processed I times in this second training epoch. This procedure is continued until a maximum number of Gaussians is reached.

Table 2 :

 2 Database Setup

		Saint Gall Parzival
	Training Validation	468 235	2237 912
	Test	707	1328
	Table 3: Recognition Accuracy
		Saint Gall	Parzival
	Reference	88.99	83.89
	Kanungo Character 90.42 +1.43 85.37 +1.48 90.15 +1.16 86.95 +3.06
	Geometric 89.25 +0.26 85.66 +1.77
	Combined 90.81 +1.82 87.12 +3.23
	compared with the reference system trained on the original
	training set.

Table 4 :

 4 Error Reduction

		Saint Gall Parzival
	Kanungo Character	10.54 12.99	18.99 09.19
	Geometric	02.36	10.99
	Combined	16.53	20.05
	Table 5: Degradation Levels
		Saint Gall Parzival
	Kanungo	1	1
	Character Geometric	2 1	1 3

ACKNOWLEDGMENTS

This work has been supported by the Swiss National Science Foundation fellowship project PBBEP2 141453.