How to find set of directions in Stahel-Donoho outlyingness measure? Currently I’m trying to understand and use the Stahel-Donoho outlyngness measure.
But unfortunately I’ve got a problem in the part where one is taking the maximum over the set of directions. I found papers and tutorials about the Stahel-Donoho estimator, but their explanation is not good enough for me. In one paper,  it has been suggested that one way to find this set of directions is by taking $p$ directions through 2 randomly chosen observations. In another paper [0], the authors use many directions through 2 randomly chosen observations. The original definition itself uses as direction $d$ the direction orthogonal to the hyperplane through $p$ randomly chosen data points. Then, the full set of directions is obtained by repeating  this procedure for many random $p$ subsets (see page 3, 2nd paragraph from the bottom of this paper).
Unfortunately I cannot get any idea how to deal with this problem. Should I use random projection? The problem I'm having is that I don’t know how to obtain the directions $d$ and projections.  I really stuck at the moment. Could someone kindly help me to understand or provide me a good reference for better understanding?


*

*[0] Hubert, M., Rousseeuw P. J. and  Vanden Branden, K. (2005). ROBPCA: A New Approach to Robust Principal Component Analysis. Technometrics
Volume 47, Issue 1.

 A: Summary: the only 'correct' definition is the original one [0][1], the other ones were designed to solve a problem with it that happens in a specific context and would probably have been best called 'pseudo-Stahel-Donoho' distances because they don't, in general, have  the same interpretation.
Now, I will show you how to obtain the standard SD distances in R, explain why they are computed the way they are, then explain what these pseudo-SD 'distances' are and why the need for them arose.
Given an $n$ by $p$ data matrix $X=\{x_1,\ldots,x_n\}$ whose entries lie in general (linear) position in $\mathbb{R}^p$, the original definition is to obtain a $p$-vector $d$ as the vector of coefficients of the hyperplane through $p$ data points chosen randomly out of your sample of $n$ observations. In R this is done so:
moar_directions<-function(x){
    n<-nrow(x)
    p<-ncol(x)
    P<-x[sample(n,p),,drop=FALSE]
    G<-rep(NA,p)
    E<-matrix(1,p,1)
    if((qrP<-qr(P))$rank==p) G<-solve(qrP,E)
    return(G)
}
n<-100
p<-5
x<-matrix(rnorm(n*p),nc=p)
moar_directions(x)

(see for example line 13 of robustbase::adjOutlyingness here, which you can check by tipping:
library(robustbase)
body(robustbase::adjOutlyingness) 

and [2]).
Then, projecting the data unto the direction normal to $d$ is done so:
x%*%moar_directions(x)

The SD outlyingness of $x_i$ with respect to $X=\{x_1,\ldots,x_n\}$ on a single projection unto $d$ is ([0][1]):
$$O(x_i,d,X)=\frac{|x_i'd-\mbox{med}_i(x_i'd)|}{\mbox{mad}_i(x_i'd)}\;\;\;(0)$$
which in R is obtained as:
w<-x%*%moar_directions(x)
abs(w-median(w))/mad(w)

And computing the SD outlyingness as the maximum outlyingness over all the members of $B_p^n$, where 
$B_p^n$ is the set of all ${n\choose p}$ such directions $d$ in a $n$ by $p$ data matrix whose entries lie in general (linear) position in $\mathbb{R}^p$. Often, $|B_p^n|={n\choose p}$ will be too large, and it will not be possible to consider all its members in which case one can sample randomly $K$ directions  from it to compute an approximation to the SD outlyigness of $x_i$ w.r.t. $X$:
$$O(x_i,X)=\underset{\{d_k\}_{k=1}^K\in B_p^n}{\max}\frac{|x_i'd_k-\mbox{med}_i(x_i'd_k)|}{\mbox{mad}_i(x_i'd_k)}$$
which is R is obtained as:
library(matrixStats)
K<-100
AllThoseProj<-matrix(NA,n,K)
for(i in 1:K){
    w<-moar_directions(x)
    AllThoseProj[,i]<-abs(w-median(w))/mad(w)
}
rowMaxs(AllThoseProj)

(I will use for a couple of lines K=100 as my definition of 'many' then latter explain what is meant by many in this context)
To motivate the original definition, notice[3] that:
$$d(x_i,X)=\underset{d\in B_p}{\max}\frac{|x_i'd_k-\mbox{mean}_i(x_i'd_k)|}{\mbox{sd}_i(x_i'd_k)}\;\;(1)$$
where
$$d^2(x_i,X)=(x_i-\mbox{mean}_i(x_i))'\mbox{Cov}_i(x_i)^{-1}(x_i-\mbox{mean}_i(x_i))$$
is the vector of squared Mahalanobis distances. So the original SD outlyingness index was designed as a way to compute 
a consistent estimator of $d^2(x_i)$ with 50% breakdown point (see [1]) and 
equality (1) only holds when the directions $d$ are defined as above.
Now, I will delve a bit more on what is meant by 'many' in ' Often, $|B_p^n|$ will be too large, and ... one can sample randomly many directions'
The notion of many that should be used here depends on the objective being pursued. If you are using the SD to approximate the vector of $d(x_i,X)$'s, then $K$ in the high hundreds should already give a good approximation. If you are using the SD as a robust alternative to   the vector of $d(x_i,X)$'s, then a much higher value of $K$ will be necessary (see page 13 of [0] for a discussion of this). 
Now, I will delve a bit more on these pseudo-SD distances.
The problem with the way $d$ is defined above is that it is only uniquely defined in settings where $n>p$. In the last decade a lot of research started to be done on high dimensional robustness. The idea built up incrementally to search for outliers through projection pursuit (as was done with the SD distances) but using a type of projection that would also 'work' when $p>n$. There are many such notions, but, gradually, a consensus emerged around the the idea of using directions through two points (see [4] for a late implementation of this idea). Directions through two points are defined as (see [4], and line 52 of body(rrcov:::extradir) in package rrcov for example):
moar_directions_2points<-function(x){
        n<-nrow(x)
        p<-ncol(x)
        P<-x[sample(n,2),,drop=FALSE]
        G<-rep(NA,p)
        E<-G
        E<-P[1,]-P[2,]
        N<-sqrt(crossprod(E))
        if(N>1e-8)  G<-E/N        
        return(G)
}

Now, if we substitute in equation (0) the original definition of the directions $d$ by the ones above, the new outlyingness index (let's call it pseudo-SD) is no longer (except when $p=2$) a consistent estimator of  the vector of $d(x_i,X)$'s, so it is a bit hard to interpret what it is. On the other hand, it can still be computed when $p>n$ in which case  the vector of $d(x_i)$'s is not even defined anyway. The fact that we lost consistency by using directions through 2 data points also means we no longer have a target (a quantity we are trying to estimate, like we did with the  the vector of $d(x_i,X)$'s in the case of the SD) so the question of what is meant by 'many' here is not really addressable. rrcov for example uses 'many' (in the case of directions through two points) to mean 250.
References:


*

*[0] Stahel W. (1981).
Breakdown of Covariance Estimators.
Research Report 31, Fachgrupp fur Statistik, E.T.H. Zurich.

*[1] Donoho. D.L.  (1982). 
Breakdown properties of multivariate location estimators
Ph.D. Qualifying Paper Harvard University.

*[2] Hubert, M. and Van der Veeken, S. (2007). Outlier detection for skewed data.     Journal of chemometrics vol:22 issue:3-4 pages:235-246.

*[3] Rousseeuw, P.J. and Leroy, A.M. (1987). 
Robust Regression and Outlier Detection. Wiley, New York.

*[4] Hubert, M., Rousseeuw P. J. and  Vanden Branden, K. (2005). ROBPCA: A New Approach to Robust Principal Component Analysis. Technometrics
Volume 47, Issue 1. 

