The medcouple kernel

We make the following observations about the kernel function ${\displaystyle h(x_{i}^{+},x_{j}^{-})}$:

1. The kernel function is location-invariant.^[1]: 999 If we add or subtract any value to each element of the sample ${\displaystyle X}$, the corresponding values of the kernel function do not change.
2. The kernel function is scale-invariant.^[1]: 999 Equally scaling all elements of the sample ${\displaystyle X}$ does not alter the values of the kernel function.

These properties are in turn inherited by the medcouple. Thus, the medcouple is independent of the mean and standard deviation of a distribution, a desirable property for measuring skewness. For ease of computation, these properties enable us to define the two sets

${\displaystyle Z^{+}:=\left.\left\{{\frac {x_{i}^{+}-x_{m}}{r}}~\right|~x_{i}^{+}\in X^{+}\right\}}$

${\displaystyle Z^{-}:=\left.\left\{{\frac {x_{j}^{-}-x_{m}}{r}}~\right|~x_{j}^{-}\in X^{-}\right\}}$

where ${\displaystyle r=2\max _{0\leq i\leq n-1}|x_{i}|}$. This makes the set ${\displaystyle Z:=Z^{+}\cup Z^{-}}$ have range of at most 1, median 0, and keep the same medcouple as ${\displaystyle X}$.

For ${\displaystyle Z}$, the medcouple kernel reduces to

${\displaystyle h(z_{i}^{+},z_{j}^{-}):={\begin{cases}\displaystyle {\frac {z_{i}^{+}+z_{j}^{-}}{z_{i}^{+}-z_{j}^{-}}}&{\text{ if }}z_{i}^{+}>z_{j}^{-}\\\operatorname {signum} (p-1-i-j)&{\text{ if }}z_{i}^{+}=0=z_{j}^{-}\end{cases}}}$

Using the recentred and rescaled set ${\displaystyle Z}$ we can observe the following.

3. The kernel function is between -1 and 1,^[1]: 998 that is, ${\displaystyle |h(z_{i}^{+},z_{j}^{-})|\leq 1}$. This follows from the reverse triangle inequality ${\displaystyle |a|-|b|\leq |a-b|}$ with ${\displaystyle a=z_{i}^{+}}$ and ${\displaystyle b=z_{j}^{-}}$ and the fact that ${\displaystyle z_{i}^{+}\geq 0\geq z_{j}^{-}}$.
4. The medcouple kernel ${\displaystyle h(z_{i}^{+},z_{j}^{-})}$ is non-decreasing in each variable.^[1]: 1005 This can be verified by the partial derivatives ${\displaystyle {\frac {\partial h}{\partial z_{i}^{+}}}}$ and ${\displaystyle {\frac {\partial h}{\partial z_{j}^{-}}}}$, both nonnegative, since ${\displaystyle z_{i}^{+}\geq 0\geq z_{j}^{-}}$.

With properties 1, 2, and 4, we can thus define the following matrix,

${\displaystyle H:=(h_{ij})=(h(z_{i}^{+},z_{j}^{-}))={\begin{pmatrix}h(z_{0}^{+},z_{0}^{-})&\cdots &h(z_{0}^{+},z_{q-1}^{-})\\\vdots &\ddots &\vdots \\h(z_{p-1}^{+},z_{0}^{-})&\cdots &h(z_{p-1}^{+},z_{q-1}^{-})\end{pmatrix}}.}$

If we sort the sets ${\displaystyle Z^{+}}$ and ${\displaystyle Z^{-}}$ in decreasing order, then the matrix ${\displaystyle H}$ has sorted rows and sorted columns,^[1]: 1006

${\displaystyle H={\begin{pmatrix}h(z_{0}^{+},z_{0}^{-})&\geq &\cdots &\geq &h(z_{0}^{+},z_{q-1}^{-})\\\geq &&&&\geq \\\vdots &&\ddots &&\vdots \\\geq &&&&\geq \\h(z_{p-1}^{+},z_{0}^{-})&\geq &\cdots &\geq &h(z_{p-1}^{+},z_{q-1}^{-})\end{pmatrix}}.}$

The medcouple is then the median of this matrix with sorted rows and sorted columns. The fact that the rows and columns are sorted allows the implementation of a fast algorithm for computing the medcouple.

Robustness

The breakdown point is the number of values that a statistic can resist before it becomes meaningless, i.e. the number of arbitrarily large outliers that the data set ${\displaystyle X}$ may have before the value of the statistic is affected. For the medcouple, the breakdown point is 25%, since it is a median taken over the couples ${\displaystyle (x_{i},x_{j})}$ such that ${\displaystyle x_{i}\geq x_{m}\geq x_{j}}$.^[1]: 1002

Values

Like all measures of skewness, the medcouple is positive for distributions that are skewed to the right, negative for distributions skewed to the left, and zero for symmetrical distributions. In addition, the values of the medcouple are bounded by 1 in absolute value.^[1]: 998

Naïve algorithm

The naïve algorithm for computing the medcouple is slow.^[1]: 1005 It proceeds in two steps. First, it constructs the medcouple matrix ${\displaystyle H}$ which contains all of the possible values of the medcouple kernel. In the second step, it finds the median of this matrix. Since there are ${\displaystyle pq\approx {\frac {n^{2}}{4}}}$ entries in the matrix in the case when all elements of the data set ${\displaystyle X}$ are unique, the algorithmic complexity of the naïve algorithm is ${\displaystyle O(n^{2})}$.

More concretely, the naïve algorithm proceeds as follows. Recall that we are using zero-based indexing.

function naïve_medcouple(vector X):
    // X is a vector of size n.
    
    // Sorting in decreasing order can be done in-place in O(n log n) time
    sort_decreasing(X)
    
    xm := median(X)
    xscale := 2 * max(abs(X))
    
    // Define the upper and lower centred and rescaled vectors
    // they inherit X's own decreasing sorting
    Zplus  := [(x - xm)/xscale | x in X such that x >= xm]
    Zminus := [(x - xm)/xscale | x in X such that x <= xm]
    
    p := size(Zplus)
    q := size(Zminus)
    
    // Define the kernel function closing over Zplus and Zminus
    function h(i, j):
        a := Zplus[i]
        b := Zminus[j]
        
        if a == b:
            return signum(p - 1 - i - j)
        else:
            return (a + b) / (a - b)
        endif
    endfunction
    
    // O(n^2) operations necessary to form this vector
    H := [h(i, j) | i in [0, 1, ..., p - 1] and j in [0, 1, ..., q - 1]]
    
    return median(H)
endfunction

The final call to median on a vector of size ${\displaystyle O(n^{2})}$ can be done itself in ${\displaystyle O(n^{2})}$ operations, hence the entire naïve medcouple algorithm is of the same complexity.

Fast algorithm

The fast algorithm outperforms the naïve algorithm by exploiting the sorted nature of the medcouple matrix ${\displaystyle H}$. Instead of computing all entries of the matrix, the fast algorithm uses the K^th pair algorithm of Johnson & Mizoguchi.^[4]

The first stage of the fast algorithm proceeds as the naïve algorithm. We first compute the necessary ingredients for the kernel matrix, ${\displaystyle H=(h_{ij})}$, with sorted rows and sorted columns in decreasing order. Rather than computing all values of ${\displaystyle h_{ij}}$, we instead exploit the monotonicity in rows and columns, via the following observations.

     function greater_h(kernel h, int p, int q, real u):
         // h is the kernel function, h(i,j) gives the ith, jth entry of H
         // p and q are the number of rows and columns of the kernel matrix H
         
         // vector of size p
         P := vector(p)
         
         // indexing from zero
         j := 0
         
         // starting from the bottom, compute the [[supremum|least upper bound]] for each row
         for i := p - 1, p - 2, ..., 1, 0:
                  
             // search this row until we find a value less than u
             while j < q and h(i, j) > u:
                 j := j + 1
             endwhile
             
             // the entry preceding the one we just found is greater than u
             P[i] := j - 1
         endfor
         
         return P
     endfunction

     function less_h(kernel h, int p, int q, real u):
     
         // vector of size p
         Q := vector(p)
         
         // last possible row index
         j := q - 1
         
         // starting from the top, compute the [[infimum|greatest lower bound]] for each row
         for i := 0, 1, ..., p - 2, p - 1:
         
             // search this row until we find a value greater than u
             while j >= 0 and h(i, j) < u:
                 j := j - 1
             endwhile
             
             // the entry following the one we just found is less than u
             Q[i] := j + 1
         endfor
         
         return Q
     endfunction

K^th pair algorithm

A visualisation of the fast medcouple algorithm. It begins with a matrix with sorted rows and sorted columns, where darker squares are smaller than lighter squares. At each iteration, the weighted median of row medians is picked, in yellow. It is then compared to the rest of the matrix to produce candidate red upper and blue lower boundaries. The algorithm then selects the boundary which is known to exclude the global matrix median, by considering the number of entries excluded by this boundary (which is equivalent to considering the rank of the yellow entry). The algorithm then proceeds until the yellow weighted median of row medians is exactly the medcouple, or the number of candidate entries is small enough to perform a selection sort amongst the remaining entries.

Putting together these two observations, the fast medcouple algorithm proceeds broadly as follows.^[4]: 148

1. Compute the necessary ingredients for the medcouple kernel function ${\displaystyle h(i,j)}$ with ${\displaystyle p}$ sorted rows and ${\displaystyle q}$ sorted columns.
2. At each iteration, approximate the medcouple with the weighted median of the row medians.^[4]: 148
3. Compare this tentative guess to the entire matrix obtaining right and left boundary vectors ${\displaystyle P}$ and ${\displaystyle Q}$ respectively. The sum of these vectors also gives us the rank of this tentative medcouple.
   1. If the rank of the tentative medcouple is exactly ${\displaystyle pq/2}$, then stop. We have found the medcouple.
   2. Otherwise, discard the entries greater than or less than the tentative guess by picking either ${\displaystyle P}$ or ${\displaystyle Q}$ as the new right or left boundary, depending on which side the element of rank ${\displaystyle pq/2}$ is in. This step always discards at least 1/4 of all remaining entries.
4. Once the number of candidate medcouples between the right and left boundaries is less than or equal to ${\displaystyle p}$, perform a rank selection amongst the remaining entries, such that the rank within this smaller candidate set corresponds to the ${\displaystyle pq/2}$ rank of the medcouple within the whole matrix.

The initial sorting in order to form the ${\displaystyle h(i,j)}$ function takes ${\displaystyle O(n\log n)}$ time. At each iteration, the weighted median takes ${\displaystyle O(n)}$ time, as well as the computations of the new tentative ${\displaystyle P}$ and ${\displaystyle Q}$ left and right boundaries. Since each iteration discards at least 1/4 of all remaining entries, there will be at most ${\displaystyle O(\log n)}$ iterations.^[4]: 150 Thus, the whole fast algorithm takes ${\displaystyle O(n\log n)}$ time.^[4]: 150

Let us restate the fast algorithm in more detail.

function medcouple(vector X):
    // X is a vector of size n
    
    // Compute initial ingredients as for the naïve medcouple
    sort_decreasing(X)
    
    xm := median(X)
    xscale := 2 * max(abs(X))
    
    Zplus  := [(x - xm)/xscale | x in X such that x >= xm]
    Zminus := [(x - xm)/xscale | x in X such that x <= xm]
    
    p := size(Zplus)
    q := size(Zminus)
    
    function h(i, j):
        a := Zplus[i]
        b := Zminus[j]
        
        if a == b:
            return signum(p - 1 - i - j)
        else:
            return (a + b) / (a - b)
        endif
    endfunction
    
    // Begin Kth pair algorithm (Johnson & Mizoguchi)
    
    // The initial left and right boundaries, two vectors of size p
    L := [0, 0, ..., 0]
    R := [q - 1, q - 1, ..., q - 1]
    
    // number of entries to the left of the left boundary
    Ltotal := 0
    
    // number of entries to the left of the right boundary
    Rtotal := p*q
    
    // Since we are indexing from zero, the medcouple index is one
    // less than its rank.
    medcouple_index := floor(Rtotal / 2)
    
    // Iterate while the number of entries between the boundaries is
    // greater than the number of rows in the matrix.
    while Rtotal - Ltotal > p:
        
        // Compute row medians and their associated weights, but skip
        // any rows that are already empty.
        middle_idx  := [i | i in [0, 1, ..., p - 1] such that L[i] <= R[i]]
        row_medians := [h(i, floor((L[i] + R[i])/2) | i in middle_idx]
        weights := [R[i] - L[i] + 1 | i in middle_idx]
        
        WM := weighted median(row_medians, weights)
        
        // New tentative right and left boundaries
        P := greater_h(h, p, q, WM)
        Q := less_h(h, p, q, WM)
        
        Ptotal := sum(P) + size(P)
        Qtotal := sum(Q)
        
        // Determine which entries to discard, or if we've found the medcouple
        if medcouple_index <= Ptotal - 1:
            R := P
            Rtotal := Ptotal
        else:
            if medcouple_index > Qtotal - 1:
                L := Q
                Ltotal := Qtotal
            else:
                // Found the medcouple, rank of the weighted median equals medcouple index
                return WM
            endif
        endif
   
    endwhile
    
    // Did not find the medcouple, but there are very few tentative entries remaining
    remaining := [h(i, j) | i in [0, 1, ..., p - 1],
                            j in [L[i], L[i] + 1, ..., R[i]]
                            such that L[i] <= R[i] ]
    
    // Select the medcouple by rank amongst the remaining entries
    medcouple := select_nth(remaining, medcouple_index - Ltotal)
   
    return medcouple
endfunction

In real-world use, the algorithm also needs to account for errors arising from finite-precision floating point arithmetic. For example, the comparisons for the medcouple kernel function should be done within machine epsilon, as well as the order comparisons in the greater_h and less_h functions.