Generalized implementation of invariant coordinate selection with positive semi-definite scatter matrices
Aurore Archimbaud







1]Aurore Archimbaud

[1]TBS Business School, 1 Place Alphonse Jourdain, 31000 Toulouse, France
[mycorrespondingauthor]Corresponding author. Email address: <a.archimbaud@tbs-education.fr>

§ ABSTRACT 
Invariant coordinate selection (ICS) is an unsupervised multivariate data transformation useful in many contexts such as outlier detection or clustering. It is based on the simultaneous diagonalization of two affine equivariant and positive definite scatter matrices. Its classical implementation relies on a non-symmetric eigenvalue problem (EVP) by diagonalizing one scatter relatively to the other. In case of collinearity, at least one of the scatter matrices is singular and the problem cannot be solved. To address this limitation, three approaches are proposed based on: a Moore-Penrose pseudo inverse (GINV), a dimension reduction (DR), and a generalized singular value decomposition (GSVD). Their properties are investigated theoretically and in different empirical applications. Overall, the extension based on GSVD seems the most promising even if it restricts the choice of scatter matrices that can be expressed as cross-products. In practice, some of the approaches also look suitable in the context of data in high dimension low sample size (HDLSS). 


 
Dimension reduction Generalized Eigenvalue Problem High-dimension Pseudo-inverse Singular scatters Singular value decomposition 
[2020] Primary 62H99 
Secondary 62-08 
Tertiary 65F99 








§ INTRODUCTION

Invariant Coordinate Selection (ICS) is a powerful unsupervised multivariate method designed to identify the structure of multivariate datasets on a subspace. It relies on the joint diagonalization of two affine equivariant and positive definite scatter matrices _1 and _2 and is particularly relevant as a dimension reduction tool prior to clustering <cit.> or outlier detection <cit.>. It goes beyond the well-known Principal Components Analysis (PCA) method by not maximizing the inertia but optimizing a generalized kurtosis. More precisely, some theoretical results <cit.> proved that under some elliptical mixture models, the subspace spanned by the first and/or last components carries the information regarding the multivariate structure and recovers the Fisher discriminant subspace, whatever the choice of scatter matrices is. 


The goal of ICS is to find the p× p matrix B=( b_1,…, b_p)^⊤ of eigenvectors, which simultaneously diagonalizes two scatter matrices _1  ∈𝒫_p and _2 ∈𝒫_p, with 𝒫_p be the set of all symmetric positive definite matrices of order p and where ^⊤ denotes the transpose operation. Typically, this simultaneous diagonalization corresponds to a generalized eigenvalue problem (GEP) or solving a linear matrix pencil, which can be simplified to a non-symmetric eigenvalue problem (EVP) by diagonalizing one scatter relatively to the other. To fix the order of the components and their normalization, we follow the standard definition of <cit.> by  diagonalizing _2 relatively to _1:

    B _1  B^⊤ =  I_p      B _2  B^⊤ =  D,

where D is a diagonal matrix with decreasing diagonal elements ρ_1≥…≥ρ_p>0, which correspond to the eigenvalues of _1^-1_2 and B contains the corresponding eigenvectors as its rows. 

This problem can be  re-written as: 

    _2  b_i = ρ_i _1  b_i 
    ⇔_1^-1_2  b_i = ρ_i   b_i,  for  i∈{1,…,p}.

with the following normalization:

	
  * b_i^⊤_1 b_j=0 for i≠ j and b_j^⊤_1 b_j=1 for i = j, with i, j ∈{1,…,p},
	
  * b_i^⊤_2 b_j=0 for i≠ j and b_j^⊤_2 b_j=ρ_j for i =  j, with i, j ∈{1,…,p}.

Equivalently, as stated in <cit.>, the eigenvalues ρ_i, for i ∈{1,…,p} and the eigenvectors b_1, …, b_p can also be sequentially defined by solving the successive maximization or minimization problems of the ratio: 

    𝒦( b) =  b^⊤_2  b/ b^⊤_1  b,

where ρ_1 is the maximal possible value of 𝒦( b) over b ∈ℝ^p which is achieved in the direction of the eigenvector b_1. 
The so-called invariant coordinates or components are then obtained as: _n =(_n -  1_n T(_n)^⊤) B(_n)^⊤, where _n =(_1,…,_n)^⊤∈ℝ^n × p is a p-variate sample with n observations, 1_n denotes an n-variate vector full of ones, and T(_n) denotes a location estimator, usually the one that goes along with _1.

In practice, the requirement of positive definiteness for the two scatter matrices is limiting 
and we focus on the case where at least one of these scatter matrices is singular. This is a common case since the variance-covariance matrix can be semi-definite positive as soon as some variables are collinear or in a high dimension low sample size context (HDLSS), i.e. if the number of variables exceeds the number of observations. Nowadays, more and more data are easily collected and collinearity issues arrive more frequently even when the number of observations is still higher than the number of observations. In this case, performing ICS is very challenging as we might not be able to compute: (i) one or the two scatter matrices and (ii)  the inverse of _1 as required to solve the GEP <ref>. A simple idea is to perform a variable selection, if many variables are known to be non-relevant. However, this procedure can induce the deletion of a substantial number of variables to obtain a convenient number of dimensions on which both scatter estimators can be defined. And so, it can also lead to a potential loss of information. 

The collinearity is a long-standing issue in multivariate analysis since quite a lot of methods rely on the simultaneous diagonalization of two or more scatter matrices<cit.>, for which the non-singularity of a scatter matrix is required. One of the most well-known methods that solves a GEP of two scatter matrices is 
the classical Fisher Linear Discriminant Analysis (LDA) which maximizes the separation ratio of the between and the within-group covariance matrices Σ_B and Σ_W. 
In the case of collinearity, the rank of the between-class covariance matrix Σ_B  is equal to the number of groups, leading to a singular matrix. So, the maximization problem cannot be performed by solving the  Σ_W^-1Σ_B eigenvalue problem anymore. To overcome such an issue, <cit.> and <cit.> review some solutions even in case the within-group covariance matrix  Σ_W is also singular. Among others, the proposed approaches exploit the Moore-Penrose pseudo-inverse or the Generalized Singular Value Decomposition (GSVD). Another well-used approach is the dimension reduction through a singular value decomposition like in <cit.> or <cit.>. This rank reduction is also used for the LDA method in the HDLSS context, as explained by <cit.> and as a pre-processing step of ICS in <cit.>. 
The objective of this paper is to adapt some of those approaches to generalizing ICS to the singularity issue, to investigate theoretically and practically their properties and to provide R <cit.>  implementations as well.

The structure of this paper is as follows. In Section <ref>, we introduce a definition of ICS with semi-definite positive scatter matrices and we present the challenges associated with such a context. In Section <ref> we propose three approaches to adapting ICS to the case of semi-definite positive scatter estimates based on the Moore-Penrose pseudo inverse (GINV), the dimension reduction (DR) and the generalized singular value decomposition (GSVD). We also investigate theoretically their properties in terms of: (i)  the criterion to optimize, (ii) the affine invariance of the scores and (iii) the symmetry of the roles of the two scatter estimates _1 and _2. Section <ref> focuses on different empirical applications to infirm or confirm the theoretical aspects. Finally, Section <ref> concludes the paper and discusses further perspectives.



§ ICS WITH SEMI-DEFINITE POSITIVE SCATTER MATRICES

In Subsection <ref>, some scatter matrices are detailed and a more general definition is given in case it is semi-definite positive. Subsection <ref> generalizes ICS to positive semi-definite scatter matrices,  while Subsection <ref> presents the main challenges associated with it.




 §.§ Scatter matrices

Let 𝒫_p be the set of all symmetric positive definite matrices of order p, 𝒮𝒫_p be the set of all symmetric positive semi-definite matrices of order p. Generally, a scatter matrix is defined as a p × p scatter matrix (_n) ∈𝒫_p which is affine equivariant in the sense that: 

    (_n  A +  1_n γ^⊤) =  A^⊤(_n)  A,


where A is a full rank p × p matrix, γ a p-vector and 1_n an n-vector full of ones. Among the most common ones there are the regular covariance matrix:

    (_n)= 1/n-1∑_i=1^n (_i-)(_i-)^⊤,


where  denotes the empirical mean, and the so-called scatter matrix of fourth moments:

    _4(_n) = 1/(p+2)n∑_i=1^n r_i^2 (_i-)(_i-)^⊤,

where r_i^2 = ( x_i - x̅)^⊤(_n)^-1( x_i - x̅) is the classical squared Mahalanobis distance. It is well known that those two scatter matrices are not robust against non-normality and the presence of outliers. One of the most widely-used robust alternatives is the minimum covariance determinant estimator () <cit.>. For a tuning parameter α∈ [0.5, 1], the  selects out of the n observations those n_α = ⌈α n ⌉ observations _i_1, …, _i_n_α for which the sample covariance matrix has the smallest determinant: 

    _α(_n) = c_α1/n_α∑_j=1^n_α (_i_j-_α,n)(_i_j-_α,n)^⊤,
 
where _α,n is the sample mean of the selected set of observations and c_α is a consistency factor. To increase efficiency, it is often combined with a reweighting step, see for example <cit.> for more details.

 

Theoretically, we can consider scatter matrices (_n) which are only symmetric positive semi-definite. If _n is not of full rank, then (_n) ∈𝒮𝒫_p. For robust scatter matrices based on a subset of observations lying on a subspace of lower dimension than the entire space, (_n) belongs to 𝒮𝒫_p. This arises for example in the presence of outliers in the orthogonal complement subspace (OC outliers). In this context, we can extend the definition of a scatter matrix (_n) ∈𝒮𝒫_p which is said to be affine equivariant in the sense that:

    (_n  A +  1_n γ^⊤) =  A^⊤(_n)  A,

for any A and γ as previously defined. <cit.> proves that if _n lies in some r ≤ p-dimensional hyperplane, then  
affine equivariant location and scatter statistics are essentially affine equivariant statistics defined on this hyperplane. For  a p × (p-r) matrix M of rank p-r and m ∈ℝ^p-r, let us defined the hyperplane ℋ(M,m)={∈ℝ^p |  M^⊤= m}, then for any affine equivariant scatter matrix (_n): M^⊤(_n)  M =  0_(p-r) × (p-r), where 0_j × k denote the j × k matrix of all zeroes. In addition the lemma states that if L is any p × r matrix such that A =  [ L M ] is non singular then:

    (_n  A ) = [     ^(r)(_n  L)     0_r × (p-r);     0_(p-r) × r 0_(p-r) × (p-r) ],

where  ^(r)( Y) is a scatter matrix for n × r data matrices Y.  <cit.> also note that if the data is in general position[Data is in general position if there is no subset of k observations lying on a subspace of dimension k-2, with k≤ p+1 and p denotes the number of variables.], then all  affine equivariant scatter statistics, symmetric in the observations[i.e. V( Q _n )= V(_n) for any permutation matrix Q of order n and _n ∈ℝ^n × p, the initial data containing n observations, characterized by p variables.], are proportional to the variance-covariance matrix. In practice, if the data is perfectly multicollinear but with n > p then the data is not in general position. If n ≤ p, the situation depends on the data themselves but there are examples, mainly in the automotive field <cit.>, where the data is not in general position.

Another challenge is that many of the well-known robust affine equivariant scatter statistics such as the M-estimators <cit.> or the  are not well-defined when the data contains collinear variables. Usually, we can only define the variance-covariance matrix or the projection-based estimators <cit.> such as the Stahel-Donoho estimator <cit.>. To overcome this issue, regularized estimators of scatter matrices have been proposed: the regularized M-estimators <cit.> or the minimum regularized covariance determinant estimator <cit.> () among others. 


The scatter matrix based on the fourth moments 4 can also be defined in a more general context if the inverse of  is replaced by its pseudo-inverse ^+:

    _4(_n) = 1/(p+2)n∑_i=1^n r_i^2 (_i-)(_i-)^⊤,

where r_i^2 = ( x_i - x̅)^⊤^+(_n)( x_i - x̅) is the classical squared Mahalanobis distance. All those p× p scatter matrices belong to 𝒮𝒫_p if _n is not full rank and are no longer necessary affine equivariant. This is a well-known consequence and it is common to relax the affine invariance of the multivariate methods in case of singularity. For example, <cit.> propose to focus on some “weak invariance" based on the relative ranking of the outlierness measures of the observations instead of requiring the same score value after an affine transformation.

For convenience, for the rest of the paper, the dependence on _n is dropped from the different scatter matrices (_n) when the context is obvious.



 §.§ ICS as a generalized eigenvalue problem

With the common definition of the ICS method, the two scatter matrices _1 and _2 should be definite positive to find a finite and nonzero eigenvalue ρ to the eigenproblem (<ref>): _2  b = ρ_1  b ⇔_1^-1_2  b = ρ b. The positive definiteness of _1 is required to compute its inverse whereas the positive definiteness of _2 ensures nonzero eigenvalues. 
If  _1 is singular and not proportional to _2, then the equivalence (<ref>) is not true anymore since _1 is not invertible. The problem can not be simplified to a non-symmetric eigenvalue problem (EVP) anymore and  so, we have to solve the initial Generalized Eigenvalue Problem  (GEP): 

    _2  b_i = ρ_i _1  b_i      for  i ∈{1,…,p},
 
with _1 ∈𝒮𝒫_p and _2 ∈𝒮𝒫_p which are not necessarily of full ranks. If _1 and/or _2 is singular, that means that the null spaces of the two scatter matrices are not empty and they do not necessarily span the same subspace. Concretely, in this context, solving  the GEP of _1 and  _2, leads to consider the following cases:

	
  * if b ∈(_1) ∩(_2) then ρ∈ℝ^+*,
	
  * if b ∈(_2) -   (_1) then  ρ=0,
	
  * if b ∈(_1) -   (_2) then  ρ=∞,
	
  * if b ∈(_1) ∩(_2) then any ρ∈ℝ is a solution of the GEP. The corresponding eigenvectors are not well-defined and might cause stability issues for some algorithms. However, the structure of the data is not associated with those directions and so, we do not need to consider them further.

So, contrary to the classical ICS,  the directions b associated with infinite or zero eigenvalues should also be analyzed as they might highlight some of the structure of the data. Let us illustrate the challenges of considering semi-definite scatter matrices for ICS on one artificial data.



 §.§ Challenges with singular scatter matrices: an illustrative example


Let X=(X_1,…,X_p)^⊤ be a p-multivariate real random vector and assume the distribution of  X is a mixture of two Gaussian distributions with different covariance matrices:


    X∼ (1-ϵ)   N( 0_p, [ W_1   0;   0   0 ]) +  ϵ  N( 0_p, [ W_1   0;   0 W_2 ]),

with ϵ <1/2, W_1 ∈𝒫_r_1 and W_2 ∈𝒮𝒫_p-r_1 with ( W_2) = r_2 ≤ p-r_1.

Such a distribution illustrates a model containing two clusters: the majority of the data and a group that can be identified as outlying observations. The first cluster follows a Gaussian distribution such that the majority of the data is contained in a r_1-dimensional subspace spanned by the range of W_1. The observations from the second cluster behave the same as previously on the r_1-dimensional subspace but they are also present in r_2 directions not spanned by the majority of the data. The goal of the ICS method is to find this r_2-dimensional subspace where the observations of the second cluster are outlying. Here this subspace is spanned by the range of W_2 which is the orthogonal complement of the range of W_1, i.e. the null space of W_1.

Let us try to recover this subspace using the ICS method with a theoretical “perfectly robust” scatter functional  _1 = [ W_1   0;   0   0 ] and a theoretical “non-robust” scatter functional, the covariance of , _2= [   W_1     0;     0 ϵ W_2 ]. We have _1 ∈𝒮𝒫_p and _2 ∈𝒮𝒫_p, with (_1)=r_1<(_2)≤ p. In addition, (_1) = (  W_1) and (_2) = ( W_1) ⊕(  W_2).

Several of the aforementioned cases arise on this example. First, the intersection of the spaces spanned by the two scatter functionals _1 and _2 corresponds to the r_1-dimensional subspace spanned by W_1, so r_1 nonzero eigenvalues should be found. Then, since (_1) -   (_2) ≠{0},  a new direction associated to an ∞ eigenvalue should also be analyzed. In fact, this is the one that reveals the outliers. Finally, if (_2) < p, then the two scatter functionals share a part of their null subspaces. This subspace is not important since it contains no structure. However, we consider this phenomenon in our analysis because it is common in practice that the data is not of full rank. In addition, this feature could also make some algorithms unstable. 

Practically, to illustrate the model (<ref>), we generate 1000 observations with exactly 20 outliers, W_1= I_2 and W_2 = (2,0). On the left scatterplot matrix of the Figure <ref>, we can see that the outliers represented by some blue triangles behave differently than the majority of the data only on the third variable. The subspace spanned by this third variable is the only one of interest to identify these observations as outliers.
This example can be seen as tricky since the outliers are well-identified on the third variable and no observations lie on the fourth one. So, we apply an affine transformation based on a non-singular p × p particular Toeplitz matrix A:

    A = 
    [     1 p-1/p p-2/p   1/p; p-1/p     1 p-1/p p-2/p; p-2/p p-1/p     1 p-1/p;   1/p p-2/p p-1/p     1 ],

 to transform the initial data  to ^* =  A. We can notice on the right scatterplot matrix of Figure <ref>, that the outliers are no longer as well separated on the third transformed variable as they were initially. In addition, we are no longer able to see that the observations lie in a three-dimensional subspace. However, the structure of outlierness of the data is still contained in one dimension only. The challenge is to be able to recover the direction spanned by the outliers with ICS.

 	
	


§ IMPLEMENTATION OF ICS FOR SEMI-DEFINITE POSITIVE SCATTER MATRICES

Several methods exist to solve a Generalized Eigenvalue Problem. Among others, it exists the well-known QZ-algorithm introduced by <cit.> or the procedure described by <cit.>. However, in practice, solving a GEP of two scatter matrices with a common null space from a numerical point of view is particularly challenging. Indeed, the presence of this null space makes procedures like the well-known QZ-algorithm very unstable. So far, the available algorithms are not satisfactory as they might lead to complex and negative eigenvalues. For that reason, GEP is not directly solved in general and surrogate approaches are used in <cit.>. In this section, we investigate three theoretical and practical implementations of ICS based on  Moore-Penrose pseudo-inverse (GSVD) in Subsection <ref>, dimension reduction (DR) in Subsection <ref> and 
generalized singular value decomposition (GSVD) in Subsection <ref>. We focus on collinearity issue with n>p.



 §.§ ICS with a Moore-Penrose pseudo-inverse

If we assume that _1 and _2 can be defined and computed on _n but that _1 ∈ with (_1)=r_1<p, then it is not possible to solve _1^-1_2 since _1 is not invertible. Instead, we can replace the inverse of _1 with its Moore-Penrose pseudo-inverse _1^+ and solve:

    _1^+ _2  b = ρ b,


where _1^+ =  P Λ^+  P^⊤ = P_1 Λ_r_1^-1 P_1^⊤, with  Λ_r_1^-1 containing only the inverse of the  r_1 nonzero eigenvalues of _1, 
P is an orthogonal matrix containing the eigenvectors of _1. P can be partitioned as P = [ P_1  P_2], with P_1, the p× r_1 matrix containing the first r_1 eigenvectors associated to the  r_1 nonzero eigenvalues of _1, which is an orthonormal basis for the range space of _1. Similarly, the p× p-r_1 matrix P_2 spans the null space of _1. P_1 and P_2  are semi-orthogonal matrices such that:  P_1^⊤ P_1 =  I_r_1 and P_2^⊤ P_2 =   I_p-r_1. 

 Solving _1^+ _2 restricts the direction b associated to the largest eigenvalue ρ_1 to b =  v_1 + v_0 with v_1 ∈(_1), v_0 ∈(_1) only onto the subspace spanned by _1 and expressed by v_1:

    v_1=     b ∈ℝ^p,  b ≠ 0 b^⊤__1_2__1 b/ b^⊤ P_1 Λ_r_1 P_1^⊤ b,

where __1=  P_1  P_1^⊤. The roles of _1 and  _2 are not exchangeable anymore.


 
Instead of solving the non-symmetric EVP _1^+ _2, we transform it to _1^+1/2_2  _1^+1/2, which is symmetric:

    _1^+ _2  b = ρ b ⇔_1^+1/2_2  _1^+1/2 b^* = ρ b^*
      ⇔
    (Λ_r_1^-1/2 P_1^⊤_2  P_1 Λ_r_1^-1/2-ρ I_r_1) b^*=0,



with b= P_1 Λ_r_1^-1/2 b^*. By multiplying by P_1Λ_r_1^1/2, and because  P_1 is only semi-orthogonal, the equation (<ref>) can be rewritten as:

    (  P_1  P_1^⊤_2  P_1   P_1^⊤ -ρ P_1 Λ_r_1 P_1^⊤) b=0,

which leads to the following modified ICS criterion for the eigenvector associated with the largest eigenvalue:


    b ∈ℝ^p,  b ≠ 0max    b^⊤ P_1  P_1^⊤_2  P_1   P_1^⊤ b/ b^⊤ P_1 Λ_r_1 P_1^⊤ b,

with P_1  P_1^⊤= __1, an orthogonal projection matrix onto the (_1), as P_1 is an orthonormal basis for (_1). So ( P_1  P_1^⊤_2  P_1  P_1^⊤) ⊆(_1) and ( P_1 Λ_r_1 P_1^⊤) =  (_1). In addition,  ℝ^p can be decomposed such that: ℝ^p = (_1) ⊕(_1), thus the solution b of the criterion (<ref>) can be expressed as:

    b =  v_1 + v_0   with   v_1 ∈(_1)   and   v_0 ∈(_1),
		

    with   v_1=     b ∈ℝ^p,  b ≠ 0 b^⊤__1_2__1 b/ b^⊤ P_1 Λ_r_1 P_1^⊤ b.

As (_1) ⊆( __1_2  __1), optimizing the new criterion (<ref>) restricts the solutions to directions b only onto the subspace spanned by _1 and expressed by _1. 





If the structure of the data is only visible onto the subspace spanned by _2, in the null space of _1, then it is not possible to highlight it and recover the outlying observations: if b ∈(_1) - (_2), then  ρ=∞.
The roles of  _1 and  _2 are not exchangeable anymore. Indeed, the directions found only span the range of the inverted scatter matrix. So, the ranks of the null spaces of _1 and  _2 are now important. The results remain if _2 is singular or not.


Equivalence with the classical ICS.
If _1 ∈  then solving the _1^-1_2 eigenvalue problem or using the Moore-Penroe pseudo-inverse of _1 is equivalent because _1^+=_1^-1.


If _2(_n) is not defined, the standard ICS algorithm is not applicable but computing _2 on the whitened data might be possible: _2(_n _1^+1/2).




    Going back to our artificial example <ref>, using the Moore-Penrose pseudo-inverse of _1 leads to optimize the following criterion:

    b ∈ℝ^p,  b ≠ 0max    b^⊤_ W_1_2 _ W_1 b/ b^⊤ W_1  b    =     b ∈ℝ^p,  b ≠ 0max    b^⊤ W_1    b/ b^⊤ W_1  b=1.

Clearly, in this case, any b ∈ℝ ^p is a solution of the maximization which implies that the structure of outlierness contained in W_2 cannot be highlighted. We obtain two eigenvalues equal to one since  _1 is two-dimensional and two others equal to zero. The projection of the data onto the eigenvectors space is illustrated in Figure <ref>. Definitely, the outliers cannot be identified because the eigenspace is restricted to the subspace spanned by _1 which does not contain the structure of outlierness defined by W_2. So, the pseudo-inverse of _1 does not always give the correct solution to the singularity issue of the scatter matrices.






If the roots ρ_1,…,ρ_p are all distinct, then for the orthogonal transformation _n^* =  _n  A + 1_n γ^⊤, with A being non-singular and γ∈ℝ^p, the coordinates
_n^* =(_n^* -  1_n T(_n^*)^⊤) B(_n^*)^⊤ and 
_n =(_n -  1_n T(_n)^⊤) B(_n)^⊤, 
 then

    _n^* = _n  J,

where J is a p × p diagonal matrix with diagonal elements ± 1, which means the coordinates Z_n^* and Z_n  are invariant up to their signs through an orthogonal transformation.


 
Let _n^* =  _n  A +  1_n γ^⊤, with A being non-singular and γ∈ℝ^p and _1(_n) ∈ with (_1)<p. By definition of a scatter matrix: _1(_n^*) =  A^⊤_1(_n)  A and if  A is an orthogonal matrix, then: 
_1^+(_n^*) =  A^-1_1^+(X_n) ( A^⊤)^-1 =   A^⊤_1^+(X_n)  A. 
Following the computations detailed in the proof of Property <ref>:   
    _1^+(_n^*) _2(_n^*) b̃ = ρb̃⇔ (  A  P_1  P_1^⊤_2(_n)  P_1   P_1^⊤ A^⊤ -ρ A  P_1 Λ_r_1 P_1^⊤ A^⊤ ) b̃  =0,

with b= A  P_1 Λ_r_1^-1/2b̃ and by multiplying by A  P_1Λ_r_1^1/2.
This leads to the following modified ICS criterion:

    b ∈ℝ^p,  b ≠ 0max    b^⊤ A  __1(_n)_2(_n)  __1(_n) A^⊤ b/ b^⊤ A  P_1 Λ_r_1 P_1^⊤ A^⊤ b.

Compared to the criterion (<ref>), the eigenvectors are rotated by A^⊤ and so projecting the transformed data _n^* onto B A^⊤ or projecting _n onto B leads to the same coordinates _n^* and _n up to their signs. 



If _1 ∈𝒮𝒫_p then the ICS coordinates are not necessarily invariant by an affine transformation since the assumption of orthogonality is required in the proof (see the next counter-example <ref>).



We consider the simulated data transformed by the non-singular matrix A (<ref>). In this case, the structure of outlierness of the data is still contained only in one dimension. So, if the two scatter matrices _1 and _2 are of full ranks then, doing ICS on the initial data _n or on the transformed _n^*=_n  A should lead to the same coordinates. However, if (_1)<p and if we use the pseudo-inverse _1^+ then we lose this affine invariance property of the Invariant Components (IC). Indeed, in the simulated example, we obtain two different eigenvalues, ρ_1=1.1237 and ρ_2=1 instead of the two equal to one, as illustrated on the right panel on Figure <ref>. Obviously,  projecting the data onto the eigenvectors' space leads to new scores, and the affine invariance of the coordinates is lost. 



To conclude, using the generalized inverse shows some differences. First, the initial ICS criterion (<ref>) may be modified to the criterion (<ref>), which leads to finding directions only on the subspace spanned by _1 and so the structure contained in the space spanned by (_1) -(_2) cannot be highlighted. Second, if we use a generalized inverse, the coordinates are invariant up to an orthogonal transformation as for PCA but no longer to an affine transformation. This is unfortunate since an additional choice is required: standardize the data or not.
Finally, the two scatter functionals _1 and  _2 are not exchangeable anymore. Indeed, the directions found only span the range of the inverted scatter matrix. The results remain if _2 is singular or not.



 §.§ ICS with a dimension reduction as pre-processing

Another well-known approach consists of getting rid of the singularity issues by doing a reduction of dimension (DR) first, hoping that no information about the data structure will be lost.  The idea is to perform a Singular Value Decomposition (SVD) of the initial data and to project it onto the right-singular vectors associated with the non-zero singular values. Among others, <cit.> or <cit.> use this pre-processing step before applying their outlier detection algorithms based on PCA or  Mahalanobis distances. This rank reduction is also used for the LDA method in the HDLSS context, as explained by <cit.>.  However, the performance of the preprocessing for the LDA method relies on the rank of the covariance matrix which has to fall into a specific range to ensure that the new within-covariance becomes non-singular. Another pitfall is noted by <cit.> who advised against using an SVD before a robust PCA in the presence of OC outliers.

Considering the reduced data  _n^* =_n P_1 ∈ℝ^n × r__n of rank r__n, we have to solve the GEP of _1(_n^*) and _2(_n^*):

    _2(_n^*)  b = ρ_1(_n^*)  b,

where _n =  U  D  P^⊤ with U and P two n  ×  n and p  ×  p orthogonal matrices, D = [     Δ^1/2_r__n× r__n     0_r__n× (p-r__n);     0_(n-r__n)× r__n 0_(n-r__n)× (p-r__n) ]. The elements of the diagonal matrix Δ^1/2 are the square roots of the positive eigenvalues of _n _n^⊤ and _n^⊤_n. The columns of P are also the eigenvectors of _n^⊤_n and the columns of U are the eigenvectors of _n _n^⊤. U and  P can be partitioned as U=[  U_1   U_2] where  U_1 is n  ×  r__n, U_2 is n  ×  (n-r__n) and P=[  P_1  P_2] where P_1 is p × r__n and  P_2 is p  ×  (p-r__n).  U _1 and P_1 are both semi-orthogonal matrices and U_1 (resp. P_1) are orthonormal basis for the column space (resp. the row space) of _n. 

If _n is of full rank.

If _n is of full rank then performing an SVD as a preprocessing step before  ICS leads to the same components as if we directly compute the invariant coordinates (IC) from the initial data _n. Indeed, in this case,  (_n)=r= p, the data is transformed by a P non-singular orthogonal p × p matrix and it is known that the invariant coordinates are invariant by an orthogonal transformation. However, from a computational point of view, some numerical discrepancies can arise due to the additional step.



If (_1(_n^*))= (_n^*) but  (_2(_n^*)) -   (_n^*) ≠{0}.

If (_1(_n^*))= (_n), then  the GEP (<ref>) of _1(_n^*) and _2(_n^*) can be simplified to the classical EVP: _1(_n^*)^-1_2(_n^*)  b = ρ b.
However,  doing the first step of dimension reduction and transforming the data onto _n^* does not ensure the non-singularity of _2(_n^*). So, it is possible to choose a _2(_n^*) of lower rank than _n. In this case, the solution can lead to some directions associated to a zero eigenvalue of multiplicity potentially greater than one.


If (_1(_n^*)) -   (_n^*) ≠{0}.

If (_1(_n^*)) -   (_n^*) ≠{0}, it means that _1(_n^*) is still singular and doing the first step of dimension reduction does not solve the problem.


 

If _1 = and  _2 is any scatter matrix as defined in (<ref>), then performing ICS with the Moore-Penrose pseudo-inverse of _1 or running ICS on the reduced data leads to the same coordinates up to their signs.



 
Let us start with the initial EVP in <ref>:  _2(_n)  b = ρ_1(_n)  b. By multiplying by P^⊤ and with b =  P b̃ the equation can be rewritten as:

    P^⊤_2(_n)  P b̃ = ρ P^⊤_1(_n)  P b̃⇔_2(_n  P)b̃ = ρ_1(_n P)  b̃,


with 
(_n  P ) = [      ^(r__n)(_n  P_1)      0_r__n× (p-r__n);     0_(p-r__n) × r__n 0_(p-r__n) × (p-r__n) ] and ^(r__n)(_n  P_1) = (_n^*). Using the notations introduced in  Subsection <ref> for the Moore-Penrose pseudo-inverse: _1(_n^*)=Λ_r_1. This is because the right eigenvectors from the singular value decomposition of _n are the ones of _n^⊤_n=_1(_n) =  P_1 Λ_r_1 P_1^⊤ and _1(_n^*)= P_1^⊤_1(_n)  P_1 with P_1 being a p × r__n semi-orthognal matrix and so r__nik<p. 

Finally, we solve:

    [                _2(_n^*)       0_r__n× (p-r__n );          0_(p-r) × r__n 0_(p-r__n ) × (p-r__n ) ]b̃ =  ρ[                  Λ_r__n         0_r × (p-r__n );          0_(p-r) × r__n 0_(p-r__n ) × (p-r__n ) ]b̃,

with b̃ =  v_1 + v_0   with   v_1 ∈(_n)   and   v_0 ∈(X_n),
and so if we restrict the solutions only onto the subspace spanned by (_n^*):

    _2(_n^*)  v_1  = ρ_1(_n^*)  v_1 
    ⇔_2(_n^*)  v_1  = ρΛ_r__n v_1 
    ⇔  v_1 ∈ℝ^p,  v_1 ≠ 0 v_1^⊤_2(_n^*) v_1/ v_1^⊤Λ_r_1 v_1.

If we compare to the modified criterion <ref> obtained using the generalized inverse:


    a ∈ℝ^p,  a ≠ 0max    a^⊤ P_1  P_1^⊤_2(_n)  P_1   P_1^⊤ a/ a^⊤ P_1 Λ_r_1 P_1^⊤ a⇔ a ∈ℝ^p,  a ≠ 0max    a^⊤ P_1 _2(_n^*)  P_1^⊤ a/ a^⊤ P_1 _1(_n^*)  P_1^⊤ a.

with a=  a_1 + a_0 with a_1 ∈(_1(_n)) =  (_n) and a_0 ∈(_1(_n)). So, after projecting, the new coordinates are the same up to their signs. 



So in this case, doing the pre-processing leads to an additional step to the method which is not needed and which implies the same drawbacks as doing ICS with a generalized inverse.



 From a practical point of view, estimating the rank of _n might be very challenging as illustrated in Subsection <ref> and lead to a loss of information regarding the structure of the data.


To conclude, the preprocessing step of dimension reduction does not fulfill all its promises. First, it cannot guarantee that it solves the singularity issues of the scatter matrices. Then, even if it does, if we choose _1 as the variance-covariance matrix, we recover exactly the same modified criterion to solve as when we use the generalized inverse and so the same drawbacks. Finally, if we choose _2 as the variance-covariance matrix, we might be unable to recover the structure of the data if it is only contained on the subspace spanned by _1.



 §.§ ICS with a generalized singular value decomposition

In this section, we focus on an implementation based on a Generalized Singular Value Decomposition (GSVD) as proposed by 
<cit.> and <cit.> in the LDA context. More specifically,  they use a GSVD for computing eigenvectors to define the Fisher's discriminant subspace, when the between Σ_B  and the within-group Σ_W covariance matrices, are susceptible to be singular. The only requirement with this method is to express Σ_B and Σ_W as cross-product matrices, which is easily obtained by their definition. 

This procedure, which uses a GSVD to solve a GEP, can be applied to other scatter matrices which can be expressed as crossproducts.





However, defining a stable algorithm for the GSVD is very challenging and a lot of research was done regarding this topic, such as <cit.> among others. In this section, we present the GSVD procedure as it is given in <cit.>, restricted to the case of real matrices. We retain this definition since it is already implemented into LAPACK and can be used directly in R through the geigen <cit.> package.


Let us define __1∈ℝ^n × p s.t. __1^⊤__1 =_1=_1( X_n) and X_V_2∈ℝ^n × p s.t. X_V_2^⊤ X_V_2 =V_2=_2( X_n). _1, _2   ∈𝒮𝒫_p  with (__1) = (_1) = r_1 ≤ p and (__2) = (_2) = r_2 ≤ p.
The GSVD of __1 and __2 allows us to define the generalized eigenvalues and eigenvectors of the pencil __2^⊤__2-ρ__1^⊤__1: 
		
    B __1^⊤__1 B^⊤ =  B _1(_n)  B^⊤ = [         0         0;         0 D_1^⊤ D_1 ]   and    B __2^⊤__2 B^⊤ =  B _2(_n)  B^⊤ = [         0         0;         0 D_2^⊤ D_2 ],

where __1 =  U  D_1 [  0   R]  Q^⊤, __2 =  V  D_2 [  0   R]  Q^⊤, U and V are n × n, Q is p × p,  U, V and Q are orthogonal. R is r × r, upper triangular and nonsingular with r= ([__1^⊤, __2^⊤]^⊤). [  0   R] is r × p (in other words, the 0 is an r× (p - r ) zero matrix). D_1 and D_2 are n × r. Both are real, nonnegative, and diagonal, satisfying D_1^⊤ D_1 + D_2^⊤ D_2= I_r. Write D_1^⊤ D_1 = (α_1^2,…, α_r^2) and D_2^⊤ D_2 = (β_1^2,…, β_r^2), the ratios α_j/β_j for j = 1, …, r are called the generalized singular values.
B^⊤ =  Q [ I_n-r     0;     0  R^-1 ]. Note that the normalization is not the same as the one presented in the standard definition <ref> but it can easily be adapted. 

Vectorially, it is equivalent to solving a modified version of the GEP (<ref>): _2  b_i = ρ_i _1  b_i, for i=1,…,p:

    β_i^2  _2  b_i = α_i^2 _1  b_i 
    ⇔_2  b_i = ρ_i _1  b_i,  for  i=1,…,p.

where ρ_i  =α_i^2/β_i^2 is real, nonnegative and possibly infinite. 
The rows of B are the eigenvectors of __2^⊤__2-ρ__1^⊤__1 or equivalently of _2 -ρ_1, and the “nontrivial” eigenvalues are the squares of the generalized singular values: ρ_i = α_i^2/β_i^2, for i=p-r+1,…,p. The “trivial” eigenvalues are those corresponding to the leading p -r rows of B, which span the common null space of __1^⊤__1 and __2^⊤__2. These eigenvalues are not well defined and are not of interest. All the cases of interest are summarized in Table <ref>. 



Re-writting the GEP (<ref>) as in (<ref>)  presents some advantages compared to the other two methods. First, it allows to find all the directions which can reveal some structure of the data in the general case where  _1 ∈𝒮𝒫_p and  _2 ∈𝒮𝒫_p, as summarized in the Table <ref>. Second, it is clear that _1 and _2 play a symmetric role. This is important since the other methods can miss the structure of the data if it is contained into the subspace spanned by _2 and in the null space of
_1 in particular. Third, this formulation is still equivalent to the classical EVP (<ref>) if the two scatter matrices are of full ranks with ρ_i=α_i^2/β_i^2. In addition to these nice characteristics,  the invariant coordinates remain invariant by affine transformation.

Affine invariance property.

For two affine equivariant scatter matrices _1 and _2, and using the eigenvectors defined by (<ref>), the invariant coordinates are invariant by an affine transformation.


 
Adaptation of the proofs from <cit.>, appendix A.1, for distinct and multiple roots detailed in Appendix <ref>.





 
Let us take the same example as previously from Subsection <ref>. Using the GSVD of __1 and __2 to solve the GEP (<ref>) leads to investigate four different cases for the direction b as illustrated in the left panel of Figure <ref>:

	
  * if b ∈(_1) ∩(_2) = ( W_1), then the direction b is restricted to the subspace spans by W_1 as when we use  the Moore-Penrose pseudo-inverse and we obtain two eigenvalues equal to one,
	
  * if b ∈(_2) -   (_1) = {0}, then no direction b exists,
	
  * if b ∈(_1) -   (_2) = ( W_2), then  ρ=∞ because β^2=0 and so the direction b can highlight the structure of outlierness contained into the ( W_2), which is not the case when we use the  Moore-Penrose pseudo-inverse,
	
  * if b ∈(_1) ∩(_2) = (_2), then ρ is a “trivial” eigenvalue and any direction b ∈ℝ^p is a solution.

However, only the “non-trivial” eigenvalues, corresponding to the first three cases, are interesting for highlighting the structure of the data. More precisely, in this example, only the eigenvector b ∈(_1) -   (_2) = ( W_2) associated with the infinite eigenvalue, contains the structure of outlierness of the data. This is clearly visible in Figure <ref> which illustrates the projection of our simulated data onto the eigenvectors space. So, using the GSVD outperforms the use of a Moore-Penrose pseudo-inverse because it recovers the structure of outlierness of the data. In addition, in the right panel of Figure <ref>, we can see we obtain the same coordinates and eigenvalues if we transform the initial data by the non-singular matrix A.

	




To conclude, solving the GEP of _1 and _2 through the GSVD of __1 and __2 presents three major advantages. First, it solves the possible singularity issues of _1 and/or _2 by searching in all directions, and remains equivalent to the EVP of _1^-1_2 if the scatter matrices are of full ranks. In addition, it leads that _1 and _2 play a symmetric role since the symmetry of the problem is kept. The affine invariance property of the scores continues to be valid in the general case of semi-definite positive scatter matrices. Finally, it is interesting to note that this GSVD procedure is already implemented into the geigen R package. However, in practice, it is difficult to define another affine equivariant scatter estimator than the variance-covariance matrix. 



Overall, all three approaches present some advantages and limits as summarized in Table <ref>. The next section investigates if in practice those properties are kept.





§ EMPIRICAL APPLICATIONS

In this section, we illustrate the characteristics of the different approaches on different empirical applications. We restrict our analysis to ICS with a generalized inverse (GINV), pre-processed by a dimension reduction (DR) through a singular value decomposition and using the generalized value decomposition (GSVD). We exclude the direct GEP approach since it might result in complex and negative eigenvalues. First, Subsection <ref> analyses the consequences of exchanging the role of _1 and _2 on a simulated correlated mixture of Gaussian distributions. Subsection <ref> focuses on two examples for which estimating the rank of the data is challenging. Subsection <ref> evaluates the impact of transforming the data through an affine transformation on a collinear industrial data set.  Finally, Subsection <ref> investigates if those approaches are applicable also in case of high dimension low sample size (HDLSS) with n<p.




 §.§ Comparison of the three approaches: collinear clustering application

Let X = (_1, …, _d)^⊤ be a d-variate real random vector distributed according to a mixture of two Gaussian distributions such that:

    X∼ϵ_1   N(μ_1,  I_d) +  ϵ_2   N(μ_2,  I_d) ,

with ϵ_1 + ϵ_2 =1, μ_1 =  0_d, μ_2 = (δ, 0, …, 0)^⊤  where δ=10. 
We generate n=1000 observations on d=3 variables for two balanced groups with ϵ_1 = ϵ_2 =0.5. Two collinear variables are added: _4 = _2-3_3 and _5 = _3 +5_4.




  §.§.§ -_4

First of all, we compare the different ICS methods for the scatter pair -_4. For GSVD, we consider _4 since it is not possible to compute _4(_n) as _n is not full rank and so (_n) is singular and it cannot be inverted. Figure <ref> shows the scatterplots matrix of the IC resulting of ICS with the generalized inverse of  (GINV on 1^st column), after a dimension reduction (DR on  2^nd column) and with a generalized singular value decomposition (GSVD on 3^rd column). The second row illustrates the same results when we exchange _1 and _2.

It is interesting to note that depending on the method we do not obtain the same number of components: 5 with GINV and 3 for the other two. Indeed, with a DR the rank of the data is estimated to be 3 and so only three dimensions are kept. For GSVD, three non-trivial eigenvalues are also detected. For GINV, 5 components are illustrated but the last two are associated with almost zero eigenvalues: 1.5e^-15 and -2.5e^-18. Clearly, this indicates numerical issues and the last two components should be disregarded. Now, if we focus on the first three, we can notice that all the methods allow us to easily identify the two clusters on IC_3. In addition, the components are the same between GINV and DR as mentioned in Property <ref>. Finally, if we exchange _1 and _2 then the clustering structure is shown on the first component for the three methods. So here, on an example of simple collinearity it appears that the three methods lead to similar results. The only point of attention is with GINV, where some trivial eigenvalues are estimated and should be put away.






  §.§.§ _0.5-

Focusing on a different scatter pair based on a more robust scatter matrix such as the _0.5 raises several issues. Indeed, it is not possible to compute the _0.5 on our data because “More than half of the observations lie on a hyperplane”, so GINV and GSVD are not applicable. Instead, we can perform ICS with the  on the reduced data or use its regularized version _0.5 as shown in Figure <ref>. As previously, with the DR approach, only three dimensions are kept and the clusters are identifiable on the third IC. This is also true with the  but five eigenvalues are estimated for _0.5- or -_0.5. In addition, two of those eigenvalues are really small or high: 3.5e^-15, -5.3e^-15 and 6.7e^+14, 1.9e^+14. So in each case with _0.5 some eigenvalues need to be disregarded and it is an additional step to take into account.  





 §.§ Challenging estimation of the rank

This section investigates the difficulty of correctly estimating the rank of two empirical applications: a nearly singular industrial data in Subsection <ref> and some simulated data with OC outliers in Subsection <ref>.



  §.§.§ HTP3: nearly singular industrial data

We consider the HTP3 data set, analyzed by <cit.> and available in the R package ICSOutlier<cit.>. It describes n=371 high-tech parts designed for consumer products and characterized by p=33 tests. The part 32 showed defects in use and is considered as an outlier. Here, the data set contains tests in different units which lead to nearly singularity and so the classical ICS algorithm returns an error. 

Using the ICSQR implementation presented in <cit.> solves the issue for a combination of scatter based on a one-step M-scatter matrix and . In Figure <ref>, we compute the so-called squared ICS distances <cit.>, denoted ICSD^2, of the k selected components for different approaches to perform ICS. On the first plot, we use ICSQR with only the first component and the defective part (in orange) is clearly identified as having a high distance and so being an outlier compared to the other observations. We obtain similar results if we use the GSVD approach with -_4 or _4- as illustrated in Appendix <ref>. However, GINV is not working in this context as we can never compute _4. For the DR approach, a new challenge arises regarding the estimation of the rank of the SVD. As mentioned in <cit.>, it is common practice to use a relative rule to estimate the rank based on the first eigenvalue λ_i with i=1,…,p such as: (i)  λ_i/λ_1 < √(()ν) with the epsilon machine ν = 2.2^-16, (ii) λ_i/λ_1 < max(n,p)ν or (iii) (∑_i=1^lλ_i^2)/(∑_i=1^pλ_i^2) < 0.99, to explain at least 99% of the inertia as with PCA for example. Here, we obtain respectively a rank of 23, 33 or 3 based on the different criteria, meaning that we do not reduce the dimension in the second case. For the others, as we can see in the second column of Figure <ref>, the defective part is identified with =23 only if we take two components, and it is not detectable in case of =3. Considering different scatter pairs as the _0.5 is tricky because the scatter matrix cannot be computed on the reduced data. In this case, it is necessary to consider -_0.5 and not _0.5-, but the results are not improved as visible in the Appendix <ref>. So on real datasets, the estimation of the rank for the DR approach can be very challenging and not the best approach.






  §.§.§ OC outliers

This issue for estimating the rank is even more problematic in case of the presence of OC outliers. We generate n=100 observations following the projected mean-shift outlier model presented by <cit.> without noise: _n =  U  D  V^⊤ +(1μ^⊤+S) V^⊤ with random orthogonal U, V, p=5, r=3, D = {1000, 400, 200}, μ=0, and the row outlier matrix S has first O rows as L*[1,…,1] and 0 otherwise, O=4 and L=3.5.

In this context, the true rank of the data set is equal to 4 and so the classical ICS returns an error. With a DR step first, if we estimate the rank to 4 then ICS with -4 identifies the outliers on IC_1 as illustrated in Figure <ref>. However, it is not possible to compute a more robust scatter matrix like the _0.5. In this context, we can estimate the rank based on 95% of explained variance, leading to two dimensions but the information about the OC outliers is lost no matter the scatter pair. However, using a GSVD or GINV approach directly works fine as visible in the second and third plots of Figure <ref> in Appendix <ref>. With - the situation is a bit tricky because the outliers are found on IC_2 as illustrated in Figure <ref>.







 §.§ Impact of affine transformation: HTP2 - collinear industrial data

We consider another industrial data set also analyzed by <cit.> and available in the R package ICSOutlier <cit.> called HTP2 to evaluate the impact of an affine transformation on the data such as the classical standardization. It contains n=149 tests for p=457 high-tech parts with a defective part at number 28. This data set is ill-conditioned and so the classical ICS algorithm returns an error as well as the ICSQR implementation and the GINV approach.


The rank estimation for a DR is quite challenging and unstable. It is estimated to 138 with the first criterion mentioned in Subsection <ref>, 141 with the second and even to 1 with the third one based on the inertia. In addition, if the data is standardized then the two criteria lead to a rank of 141 and 51 for the third one. If we focus on the case of a rank of 138, then we cannot compute a robust scatter matrix like the _0.5. In Figure <ref>, we compute the ICSD^2 based on IC_1 and -_4 and the defective part is weirdly detectable as the observation having the smallest distance instead of the highest. However, doing it on the standardized data with 141 dimensions allows identifying it (see Figure <ref> in Appendix <ref>). This behavior shows that the DR approach is sensible to the estimation of the rank and the standardization of the data. On the contrary, if we perform ICS of -_4 with GSVD on the initial data (2^nd plot) or on the standardized one (3^rd plot), then the outlier is revealed and its ICSD^2 is stable between the two cases. It is noteworthy that the GSVD estimates 141 non-trivial eigenvalues. Finally, the regularized approach based on the  does not identify the outlier (see Figure <ref> in Appendix <ref>).






 §.§ High Dimension Low Sample Size (HDLSS) case

To go further, we generate some data in an HDLSS context with more variables p than observations n but not lying in general position. Following <cit.>, we generate n=50 observations on p=100 variables such as in Subsection <ref>. In this context, we retrieve similar results as when n>p: the classical ICS and GINV are returning an error, the results depend on the estimated rank and GSVD is working fine as we can see in Figure <ref> for -_4 and -_4. For the robust scatter matrix _0.5 it is not possible to compute it on reduced data of 4 components and the outliers are not visible if we kept only two dimensions (see Figure <ref> in Appendix <ref>). The last plot shows that - identifies the outliers, but on IC_2 which is not what is expected. 






§ CONCLUSION

Following previous ideas mainly used in  LDA context, we proposed three new ways of generalizing ICS to the case of positive semi-definite scatter matrices based on generalized inverse, dimension reduction or GSVD. We also investigated their theoretical properties (summarized in Table <ref>) and provided implementations in R. 
Theoretically, the approach based on GSVD looks the most appealing by keeping all the nice properties of the classical ICS. In practice, this method can only deal with a specific type of scatter matrices that has to be expressed as crossproducts and being affine equivariant. Relaxing this affine invariance property of ICS, empirical results showed interesting and stable results for different situations for GSVD with -_4, even on standardized data. For DR, estimating the rank of the data appears to be the main challenge with quite some sensibility on the results. In addition, the method is only orthogonal invariant and does not ensure that some robust scatter matrices can be computed on the reduced data. Finally, GINV seems the least suitable as it requires that V_2 can be computed on the initial data.

Overall, those methods allow to generalize ICS to the context of positive semi-definite scatter matrices and might be also helpful in the HDLSS case as long as the data are not in general position. Depending on the approach, it is important to think about which scatter should be the first one and if some OC outliers are present. In practice, it might be useful to perform ICS multiple times, exchanging _1 and _2 and trying the different implementations to compare the results or performing localized projection pursuit after ICS as suggested by <cit.>. In the future, another idea could be to penalize or regularize ICS, as <cit.> or <cit.> do for LDA for example.













§ COMPUTATIONAL DETAILS

All computations are performed with R version 4.3.3 <cit.> and uses the  packages ICS<cit.> for ICS, ICSClust<cit.>, ICSOutlier<cit.>, rrcov<cit.> for the MCD scatter matrix and geigen<cit.> for computing GSVD.
Replication files are available upon request.



§ ACKNOWLEDGEMENTS


This work is a generalization of the research conducted during my PhD under the supervision of Professor Anne Ruiz-Gazen whom I deeply thank for her guidance and insightful remarks on this topic.







§ APPENDIX A.   CALCULATION DETAILS OF THE MOORE-PENROSE PSEUDO-INVERSE FOR SUBSECTION <REF>

 (_n^*)^+  has to satisfy the four conditions to be a Moore-Penrose pseudo-inverse:
	
		
  * Condition 1: (_n^*) (_n^*) ^+ (_n^*)  = (_n^*).		
		
  * Condition 2: (_n^*)^+(_n^*) (_n^*)^+ = (_n^*)^+.
		
  * Condition 3: ((_n^*) (_n^*)^+)^⊤ = (_n^*) (_n^*)^+.
		
  * Condition 4: ((_n^*)^+ (_n^*))^⊤ = (_n^*)^+ (_n^*). 
	
The proof of conditions 1 and 2 can be generalized to any matrix A but conditions 3 and 4 rely on the assumption of orthogonality of A. 
Indeed, for condition 3, we have ((_n^*) (_n^*)^+)^⊤ = ( A^⊤)^-1_1(_n) _1(_n)^+  A^⊤. Since A is orthogonal ( A^⊤)^-1 =  A and _1(_n)^+ =  A_1(_n)^+  A^⊤,  so we obtained the desired equality:

    ( A^⊤)^-1_1(_n) _1(_n)^+  A^⊤=  A _1(_n)  A^⊤ A _1(_n)^+  A^⊤= (_n^*) (_n^*)^+.

 The proof is similar for the condition 4.




§ APPENDIX B.   AFFINE INVARIANCE FOR ICS WITH A GENERALIZED SINGULAR VALUE DECOMPOSITION FOR SUBSECTION <REF>



 
(i) Adaptation of the proof from <cit.>, appendix A.1, for distinct roots.

Let _n^* = _n  A + 1_n γ^⊤, with γ∈ℝ^p. Then _1(_n^*) =  A^⊤_1(_n)  A and _2(_n^*) =  A^⊤_2(_n)  A.

By definition of ICS we have, for i = 1,…,p:

 
    α_i^2(_n^*) _2(_n^*)  b_i(_n^*)     = β_i^2(_n^*) _1(_n^*)  b_i(_n^*),   
    α_i^2(_n^*)   A^⊤_2(_n)   A  b_i(_n^*)     = β_i^2(_n^*)   A^⊤_1(_n)   A  b_i(_n^*).
    Multiplying by ( A^⊤)^-1:  α_i^2(_n^*)  _2(_n)   A  b_i(_n^*)     = β_i^2(_n^*)  _1(_n)   A  b_i(_n^*).

If α_i^2(_n)/β_i^2(_n) is a distinct root, then α_i^2(_n)/β_i^2(_n) = α_i^2(_n^*)/β_i^2(_n^*) and b_i(_n) ∝ A  b_i(_n^*), so b_i(_n^*) ∝ A^-1 b_i(_n) :
					
    α_i^2(_n)  _2(_n)   b_i(_n) = β_i^2(_n)  _1(_n)   b_i(_n).

Projection onto b_i(_n^*):
					z_i^* =  b_i(_n^*)^⊤ x^* = ( A^-1 b_i(_n))^⊤ A^⊤ x =  b_i(_n)^⊤ x =  z_i.

(ii) Proof from <cit.>, appendix A.1, for multiple roots.

In case of a multiple root of multiplicity p_l, the eigenvectors are not uniquely defined and can be chosen as any linearly independent vectors spanning the corresponding p_l-dimensional eigenspace. However the roots are still the same and so the subspace spanned by the corresponding p_l-dimensional eigenspace is still the same. 




The case of multiple roots may appear  when _1 ∈𝒮𝒫_p and/or _2 ∈𝒮𝒫_p as it means than (_1) ≠{0} and/or (_2) ≠{0}. For example, if we only focus on the cases where B ∈(_2) -   (_1) or  B ∈(_1) -   (_2), if ((_2) -   (_1)) >1 and/or ((_1) -   (_2)) >1 then 0 and/or ∞ are multiple roots.



§ APPENDIX C.   ADDITIONNAL RESULTS FOR THE EMPIRICAL APPLICATIONS IN SECTION <REF>

We display additional results regarding Subsection <ref>
in Figure <ref>, Subsection <ref> in Figure <ref>, 
Subsection <ref> in Figure <ref> and Subsection <ref> in Figure <ref>.












myjmva



