Geometry & Recognition

Median은 주어진 값들의 중간에 해당하는 값을 의미한다. 1차원 데이터가 주어졌을 때 median은 주어진 데이터를 크기 순서로 정렬을 한 후 중간에 해당하는 값을 취하면 된다.  데이터를 1차원 직선 위에 한 점으로 대응시킬 때 median을 찾는 것은 각 점으로부터 거리의 합이 최소인 데이터 위치를 찾으면 된다. 즉 1차원 median 찾기는 

$$ C(y)= \sum _i  | x_i - y|$$

로 정의되는 비용함수 $C(y)$를 최소로 만드는 $y$를 주어진 데이터 속에서 찾는 것으로 바뀐다. 

증명: \begin{align}  C (y) = \sum _{x_i < y} (y- x_i)+ \sum _{x_i > y} (x_i -y ) \\  \to ~~~\frac{dC}{dy} = \sum_{x_i < y} 1  - \sum _{x_i > y} 1=0  \end{align}

이어야 하므로 $y=\text{median}\{ x_i \}$임을 알 수 있다.

2차원 이상의 다차원 공간에서도 이 비용함수를 확장하면 다차원 공간에 분포하는 점들의 기하학적인 median을 정의할 수 있다. 

$$ \text{median}( \{ \vec{x}_i\}) = \underset{y}{\text{argmin}} \sum _{i} || \vec{x}_i  - \vec{y}||$$

1차원의 경우는 점들을 크기 순서대로 정렬을 하여 구할 수 있지만, 2차원 이상일 때 닫힌 해는 존재하지 않고, 반복적인 반복을 써서 구할 수 있음이 알려져 있다( Weiszfeld Algorithm). 이는 비용함수에 경사강하법(gradient descent method)을 적용한 알고리즘으로  볼 수 있다. 주어진 단계에서 근사적인 median이 $\vec{y}_{(t)}$일 때 다음 단계에서는 median update 식은 

$$ \vec{y}_{(t+1)} = \frac{\sum_i w_i \vec{x}_i }{ \sum_i w_i}, \quad w_i = || \vec{x}_i - \vec{y}_{(t)}|| $$

로 계산한다. 이는 비용함수의 극소를 찾기 위해서 경사강하법을 적용할 때 step size (또는 learning rate)을

$$ \vec{y}_{(t+1)}=  \vec{y}_{(t)} - \lambda \nabla C \\ \lambda = \frac{1}{\sum_i w_i} $$

로 선택한 결과이다. Weiszfeld algorithm에서는 강하 방향과 step size가 닫힌 형태로 계산되므로 일반적인 경사강하 알고리즘보다도 더 빠르게 근사해를 찾을 수 있다.

참고:  ndevilla.free.fr/median/median/src/optmed_bench.c

/*
 * Optimized median search on 9 values
 */
#define PIX_SORT(a,b) { if ((a) > (b)) PIX_SWAP((a), (b)); }
#define PIX_SWAP(a,b) { BYTE temp = (a); (a) = (b); (b) = temp; }
BYTE opt_med9(BYTE p[9]) {
    PIX_SORT(p[1], p[2]) ; PIX_SORT(p[4], p[5]) ; PIX_SORT(p[7], p[8]) ; 
    PIX_SORT(p[0], p[1]) ; PIX_SORT(p[3], p[4]) ; PIX_SORT(p[6], p[7]) ; 
    PIX_SORT(p[1], p[2]) ; PIX_SORT(p[4], p[5]) ; PIX_SORT(p[7], p[8]) ; 
    PIX_SORT(p[0], p[3]) ; PIX_SORT(p[5], p[8]) ; PIX_SORT(p[4], p[7]) ; 
    PIX_SORT(p[3], p[6]) ; PIX_SORT(p[1], p[4]) ; PIX_SORT(p[2], p[5]) ; 
    PIX_SORT(p[4], p[7]) ; PIX_SORT(p[4], p[2]) ; PIX_SORT(p[6], p[4]) ; 
    PIX_SORT(p[4], p[2]) ; return(p[4]) ;
}

BYTE opt_med25(BYTE p[]) { 
    PIX_SORT(p[0], p[1]) ;   PIX_SORT(p[3], p[4]) ;   PIX_SORT(p[2], p[4]) ;
    PIX_SORT(p[2], p[3]) ;   PIX_SORT(p[6], p[7]) ;   PIX_SORT(p[5], p[7]) ;
    PIX_SORT(p[5], p[6]) ;   PIX_SORT(p[9], p[10]) ;  PIX_SORT(p[8], p[10]) ;
    PIX_SORT(p[8], p[9]) ;   PIX_SORT(p[12], p[13]) ; PIX_SORT(p[11], p[13]) ;
    PIX_SORT(p[11], p[12]) ; PIX_SORT(p[15], p[16]) ; PIX_SORT(p[14], p[16]) ;
    PIX_SORT(p[14], p[15]) ; PIX_SORT(p[18], p[19]) ; PIX_SORT(p[17], p[19]) ;
    PIX_SORT(p[17], p[18]) ; PIX_SORT(p[21], p[22]) ; PIX_SORT(p[20], p[22]) ;
    PIX_SORT(p[20], p[21]) ; PIX_SORT(p[23], p[24]) ; PIX_SORT(p[2], p[5]) ;
    PIX_SORT(p[3], p[6]) ;   PIX_SORT(p[0], p[6]) ;   PIX_SORT(p[0], p[3]) ;
    PIX_SORT(p[4], p[7]) ;   PIX_SORT(p[1], p[7]) ;   PIX_SORT(p[1], p[4]) ;
    PIX_SORT(p[11], p[14]) ; PIX_SORT(p[8], p[14]) ;  PIX_SORT(p[8], p[11]) ;
    PIX_SORT(p[12], p[15]) ; PIX_SORT(p[9], p[15]) ;  PIX_SORT(p[9], p[12]) ;
    PIX_SORT(p[13], p[16]) ; PIX_SORT(p[10], p[16]) ; PIX_SORT(p[10], p[13]) ;
    PIX_SORT(p[20], p[23]) ; PIX_SORT(p[17], p[23]) ; PIX_SORT(p[17], p[20]) ;
    PIX_SORT(p[21], p[24]) ; PIX_SORT(p[18], p[24]) ; PIX_SORT(p[18], p[21]) ;
    PIX_SORT(p[19], p[22]) ; PIX_SORT(p[8], p[17]) ;  PIX_SORT(p[9], p[18]) ;
    PIX_SORT(p[0], p[18]) ;  PIX_SORT(p[0], p[9]) ;   PIX_SORT(p[10], p[19]) ;
    PIX_SORT(p[1], p[19]) ;  PIX_SORT(p[1], p[10]) ;  PIX_SORT(p[11], p[20]) ;
    PIX_SORT(p[2], p[20]) ;  PIX_SORT(p[2], p[11]) ;  PIX_SORT(p[12], p[21]) ;
    PIX_SORT(p[3], p[21]) ;  PIX_SORT(p[3], p[12]) ;  PIX_SORT(p[13], p[22]) ;
    PIX_SORT(p[4], p[22]) ;  PIX_SORT(p[4], p[13]) ;  PIX_SORT(p[14], p[23]) ;
    PIX_SORT(p[5], p[23]) ;  PIX_SORT(p[5], p[14]) ;  PIX_SORT(p[15], p[24]) ;
    PIX_SORT(p[6], p[24]) ;  PIX_SORT(p[6], p[15]) ;  PIX_SORT(p[7], p[16]) ;
    PIX_SORT(p[7], p[19]) ;  PIX_SORT(p[13], p[21]) ; PIX_SORT(p[15], p[23]) ;
    PIX_SORT(p[7], p[13]) ;  PIX_SORT(p[7], p[15]) ;  PIX_SORT(p[1], p[9]) ;
    PIX_SORT(p[3], p[11]) ;  PIX_SORT(p[5], p[17]) ;  PIX_SORT(p[11], p[17]) ;
    PIX_SORT(p[9], p[17]) ;  PIX_SORT(p[4], p[10]) ;  PIX_SORT(p[6], p[12]) ;
    PIX_SORT(p[7], p[14]) ;  PIX_SORT(p[4], p[6]) ;   PIX_SORT(p[4], p[7]) ;
    PIX_SORT(p[12], p[14]) ; PIX_SORT(p[10], p[14]) ; PIX_SORT(p[6], p[7]) ;
    PIX_SORT(p[10], p[12]) ; PIX_SORT(p[6], p[10]) ;  PIX_SORT(p[6], p[17]) ;
    PIX_SORT(p[12], p[17]) ; PIX_SORT(p[7], p[17]) ;  PIX_SORT(p[7], p[10]) ;
    PIX_SORT(p[12], p[18]) ; PIX_SORT(p[7], p[12]) ;  PIX_SORT(p[10], p[18]) ;
    PIX_SORT(p[12], p[20]) ; PIX_SORT(p[10], p[20]) ; PIX_SORT(p[10], p[12]) ;
    return (p[12]);
}

/*---------------------------------------------------------------------------
   Function :   kth_smallest()
   In       :   array of elements, # of elements in the array, rank k
   Out      :   one element
   Job      :   find the kth smallest element in the array
   Notice   :   use the median() macro defined below to get the median. 
                Reference:
                  Author: Wirth, Niklaus 
                   Title: Algorithms + data structures = programs 
               Publisher: Englewood Cliffs: Prentice-Hall, 1976 
    Physical description: 366 p. 
                  Series: Prentice-Hall Series in Automatic Computation 
 ---------------------------------------------------------------------------*/
#define SWAP(a, b) {BYTE t = (a); (a) = (b); (b) = t; }
BYTE kth_smallest(BYTE a[], int n, int k) {
    int l = 0, m = n - 1 ;
    while (l < m) {
        BYTE x = a[k] ;
        int i = l, j = m ;
        do {
            while (a[i] < x) i++ ;
            while (x < a[j]) j-- ;
            if (i <= j) {
                SWAP(a[i], a[j]) ;
                i++ ; j-- ;
            }
        } while (i <= j) ;
        if (j < k) l = i ;
        if (k < i) m = j ;
    }
    return a[k] ;
}
#define median(a, n) kth_smallest(a, n, (((n) & 1) ? ((n) / 2) : (((n) / 2) - 1)))