IGNG — инкрементальный алгоритм растущего нейронного газа +24



При написании статьи о разработке детектора аномалий я реализовывал один из алгоритмов, который называется "Инкрементальный растущий нейронный газ".


В советской литературе российском сегменте Интернета эта тема освещена достаточно слабо, и нашлась только одна статья, да и то с прикладным применением данного алгоритма.


Итак, что же такое — алгоритм инкрементального растущего нейронного газа?


Введение


IGNG, как и GNG, является адаптивным алгоритмом кластеризации.
Сам алгоритм описан в статье Прудента и Еннаджи за 2005-й год.


Также как и в GNG, существует множество векторов данных $X$, либо генерирующая функция $f(t)$, которая предоставляет вектора из произвольно распределённых данных (параметр $t$ — время, либо номер сэмпла в выборке).


Дополнительных ограничений на данные алгоритм не накладывает.
Зато, внутри сильно отличается от GNG.


Данный алгоритм также любопытен ещё и тем, что он чуть точнее, чем GNG моделирует нейрогенез.


Описание алгоритма


Алгоритм разбивает множество данных на кластеры.
По сравнению с GNG, его преимуществом является более высокая скорость сходимости.


Идеи, на которых основан алгоритм:


  • Теория адаптивного резонанса: ?сначала ищется ближайший нейрон, и если разница не превышает порога ("параметра бдительности"), производится корректировка весов или, иначе, изменение координаты нейрона в пространстве данных. Если порог не был преодолен, создаются новые нейроны, лучше приближающие значение сэмпла данных.
  • Как связи, так и нейроны имеют параметр возраст (у GNG — только связи), который сначала равен нулю, но увеличивается по мере обучения.
  • Нейрон появляется не сразу: сначала появляется эмбрион (или зародышевый нейрон), возраст которого увеличивается с каждой итерацией, пока он не созреет. После обучения, в классификации принимают участие только зрелые нейроны.

Основной цикл


Работа начинается с пустого графа. Параметр $\sigma$ инициализируется среднеквадратическим отклонением по обучающей выборке:


$\sigma = \sqrt {\frac{1}{N}\sum\limits_{i = 1}^N {\left( {x_i - \bar x} \right)^2 } }$


Где: $\bar x$ — среднее между координатами по выборке.


Основной цикл на каждом шаге уменьшает значение $\sigma$, которое является порогом близости, и высчитывает разницу между предыдущим уровнем качества кластеризации и уровнем, который был получен после кластеризации процедурой IGNG.



Код диаграммы.
@startuml
start
:TrainIGNG(S);
:<latex>\sigma = \sigma_S,x,y \in S</latex>;
:<latex>IGNG(1, \sigma, age_{mature}, S)</latex>;
:<latex>old = 0</latex>;
:<latex>calin = CHI()</latex>;
while (<latex>old - calin \leq 0</latex>)
  :<latex>\sigma=\sigma - \sigma / 10</latex>;
  :<latex>IGNG(1, \sigma, age_{mature}, S)</latex>;
  :<latex>old = calin</latex>;
  :<latex>calin = CHI()</latex>;
endwhile
stop
@enduml

CHI — это индекс Калинского-Харабаза, показывающий качество кластеризации:


$CHI = \frac{B/(c - 1)}{W/(n - c)}$


Где:


  • $n$ — количество сэмплов данных.
  • $c$ — количество кластеров (в данном случае, количество нейронов).
  • $B$ — матрица внутренней дисперсии (сумма квадратов расстояний между координатами нейронов и средним по всем данным).
  • $W$ — матрица внешней дисперсии (сумма квадратов расстояний между всеми данными и ближайшим нейроном).

Чем больше значение индекса, тем лучше, потому если разность между индексами после кластеризации и до неё будет отрицательной, возможно считать, что индекс стал положительным и превысил предыдущий, т.е. кластеризация успешно выполнена.


Процедура IGNG


Это основная процедура алгоритма.


Она делится на три взаимоисключающих этапа:


  • Нейронов не найдено.
  • Найден один, удовлетворяющий условиям нейрон.
  • Найдено два, удовлетворяющих условиям нейрона.

Если одно из условий проходит, другие этапы не выполняются.


Сначала ищется нейрон лучше всего приближайющий сэмпл данных:


$c_1 = min(dist(\xi, \omega_c))$


Здесь $dist(x_\omega, x_\xi)$ — функция расчёта расстояния, которая обычно является Евклидовой метрикой.


Если нейрон не найден, либо он слишком далеко от данных, т.е. не удовлетворяет критерию приближённости $dist(\xi, \omega_c) \leq \sigma$, создаётся новый эмбриональный нейрон с координатами равными координатам сэмпла в пространстве данных.



Если же проверка на приближённость прошла, производится поиск второго нейрона аналогичным способом и проверка его на приближённость к сэмплу данных.
Если второй нейрон не найден, он создаётся.



Если были найдены два нейрона, удовлетворяющие условию близости к сэмплу данных, их координаты корректируются по следующей формуле:


$ \epsilon(t)h_{c,c_i} = \begin{cases} \epsilon_b,\,если\,c = c_i\\ \epsilon_n,\,если\,есть\,связь\,между\,c = c_i\\ 0,\, в\,ином\,случае \end{cases} $


$ \Delta\omega_c = \epsilon(t)h_c,_{c1}\parallel\xi - \omega_c\parallel\\ \omega_c = \omega_c + \Delta\omega_c $


Где:


  • $\epsilon(t)$ — шаг адаптации.
  • $c_i$ — номер нейрона.
  • $h_c,_{c1}$ — функция соседства нейрона $c$ с нейроном победителем (в данном случае, она вернёт 1 для прямых соседей, 0 в ином случае, потому шаг адаптации, для расчёта $\omega$ будет ненулевым только для прямых соседей).

Иными словами, координата (вес) нейрона-победителя изменяются на $\epsilon_b * \Delta\omega_{i}$, а всех его прямых соседей (тех, которые с ним связаны одним ребром графа) на $\epsilon_n * \Delta\omega_{i}$, где $\omega_i$ — координата соответствующего нейрона до изменения.


Затем создаётся связь между двумя нейронами-победителями, а если она уже была создана, её возраст обнуляется.
Возраст всех других связей увеличивается.


Все связи, возраст которых превысил константу $age_{max}$, удаляются.
После чего удаляются все изолированные (те, у которых нет связи с другими) зрелые нейроны.


Возраст непосредственных нейронов-соседей нейрона-победителя увеличивается.
Если возраст какого-либо из зародышевых нейронов превышает $age_{mature}$, он становится зрелым нейроном.


Окончательный граф содержит только зрелые нейроны.


Условием завершения процедуры IGNG ниже возможно считать фиксированное количество циклов.


Ниже показан алгоритм (рисунок кликабелен):



Код диаграммы.
@startuml
skinparam nodesep 10
skinparam ranksep 20
start
:IGNG(age, sigma, <latex>a_{mature}</latex>, S);
while (Критерий останова достигнут) is (Нет)
   -[#blue]->
  :Выбрать входной вектор e из S;
  :Найти ближайший нейрон c<sub>1</sub>;
  if (Граф пустой или\n<latex>dist(\xi, \omega_{c_1}) \leq \sigma</latex>) then (Да)
    :Вставить новый эмбрион с весом <latex>\omega_{new} = \xi</latex>;
  else (Нет)
    -[#blue]->
    :Найти второй ближайший нейрон;
    if (Если есть только один нейрон или\n <latex>dist(\xi, \omega_{c_2}) \leq \sigma</latex>) then (Да)
      :Вставить новый эмбрион с весом <latex>\omega_{new} = \xi</latex>;
      :Создать связь между <latex>c_1</latex> и <latex>c_2</latex>;
      note
      В оригинальной статье
      тут опечатка, соединяться
      должны именно два нейрона
      end note
    else (Нет)
      -[#blue]->
      :Увеличить возраст всех связей,\nисходящих из <latex>c_1</latex>;
      :<latex>\omega_{c_1} = \omega_c + \epsilon_b(\xi - \omega_{c_1})</latex>;
      :<latex>\omega_n = \omega_n + \epsilon_n(\xi - \omega_n)</latex>;
      note
      n - это все непосредственные
      соседи <latex>c_1</latex>
      (т.е. между ними одно ребро графа)
      end note
      if (c<sub>1</sub> и c<sub>2</sub> связаны) then (Да)
        :Обнулить возраст связи: <latex>age_{c_1 -> c_2} = 0</latex>;
      else (Нет)
        -[#blue]->
        :Создать связь между c<sub>1</sub> и c<sub>2</sub>;
      endif
      :Увеличить возраст всех\nнепосредственных соседей c<sub>1</sub>;
      note
      Заметьте, что увеличивается возраст нейронов,
      а не связей.
      end note
    endif

    repeat
      if (<latex>age(c) \geq a_{mature}</latex>) then (Да)
        :Сделать эмбрион $<!-- math>c</math -->$ зрелым нейроном;
      else (Нет)
        -[#blue]->
      endif

    repeat while (Ещё остались нейроны?)
  endif
  :Удалить связи, возраст которых превышает максимальный;
  :Удалить не связанные нейроны;
  note
  Этот пункт и предыдущий отсутсвуют
  в оригинальном описании процедуры IGNG,
  но они нужны, также как и для GNG.
  Об этом говорится в статье.
  endnote
endwhile (Да)

stop
@enduml

Реализация


Реализация сети выполнена на Python, с использованием библиотеки графов NetworkX. Вырезка кода из прототипа в предыдущей статье приведена ниже. Там же есть краткие пояснения к коду.


Если кого-то интересует полный код, вот ссылка на репозиторий.


Пример работы алгоритма:



Основная часть кода
class NeuralGas():
    __metaclass__ = ABCMeta

    def __init__(self, data, surface_graph=None, output_images_dir='images'):
        self._graph = nx.Graph()
        self._data = data
        self._surface_graph = surface_graph
        # Deviation parameters.
        self._dev_params = None
        self._output_images_dir = output_images_dir
        # Nodes count.
        self._count = 0

        if os.path.isdir(output_images_dir):
            shutil.rmtree('{}'.format(output_images_dir))

        print("Ouput images will be saved in: {0}".format(output_images_dir))
        os.makedirs(output_images_dir)

        self._start_time = time.time()

    @abstractmethod
    def train(self, max_iterations=100, save_step=0):
        raise NotImplementedError()

    def number_of_clusters(self):
        return nx.number_connected_components(self._graph)

    def detect_anomalies(self, data, threshold=5, train=False, save_step=100):
        anomalies_counter, anomaly_records_counter, normal_records_counter = 0, 0, 0
        anomaly_level = 0

        start_time = self._start_time = time.time()

        for i, d in enumerate(data):
            risk_level = self.test_node(d, train)
            if risk_level != 0:
                anomaly_records_counter += 1
                anomaly_level += risk_level
                if anomaly_level > threshold:
                    anomalies_counter += 1
                    #print('Anomaly was detected [count = {}]!'.format(anomalies_counter))
                    anomaly_level = 0
            else:
                normal_records_counter += 1

            if i % save_step == 0:
                tm = time.time() - start_time
                print('Abnormal records = {}, Normal records = {}, Detection time = {} s, Time per record = {} s'.
                      format(anomaly_records_counter, normal_records_counter, round(tm, 2), tm / i if i else 0))

        tm = time.time() - start_time

        print('{} [abnormal records = {}, normal records = {}, detection time = {} s, time per record = {} s]'.
              format('Anomalies were detected (count = {})'.format(anomalies_counter) if anomalies_counter else 'Anomalies weren\'t detected',
                     anomaly_records_counter, normal_records_counter, round(tm, 2), tm / len(data)))

        return anomalies_counter > 0

    def test_node(self, node, train=False):
        n, dist = self._determine_closest_vertice(node)
        dev = self._calculate_deviation_params()
        dev = dev.get(frozenset(nx.node_connected_component(self._graph, n)), dist + 1)
        dist_sub_dev = dist - dev
        if dist_sub_dev > 0:
            return dist_sub_dev

        if train:
            self._dev_params = None
            self._train_on_data_item(node)

        return 0

    @abstractmethod
    def _train_on_data_item(self, data_item):
        raise NotImplementedError()

    @abstractmethod
    def _save_img(self, fignum, training_step):
        """."""
        raise NotImplementedError()

    def _calculate_deviation_params(self, distance_function_params={}):
        if self._dev_params is not None:
            return self._dev_params

        clusters = {}
        dcvd = self._determine_closest_vertice
        dlen = len(self._data)
        #dmean = np.mean(self._data, axis=1)
        #deviation = 0

        for node in self._data:
            n = dcvd(node, **distance_function_params)
            cluster = clusters.setdefault(frozenset(nx.node_connected_component(self._graph, n[0])), [0, 0])
            cluster[0] += n[1]
            cluster[1] += 1

        clusters = {k: sqrt(v[0]/v[1]) for k, v in clusters.items()}

        self._dev_params = clusters

        return clusters

    def _determine_closest_vertice(self, curnode):
        """."""

        pos = nx.get_node_attributes(self._graph, 'pos')
        kv = zip(*pos.items())
        distances = np.linalg.norm(kv[1] - curnode, ord=2, axis=1)

        i0 = np.argsort(distances)[0]

        return kv[0][i0], distances[i0]

    def _determine_2closest_vertices(self, curnode):
        """Where this curnode is actually the x,y index of the data we want to analyze."""

        pos = nx.get_node_attributes(self._graph, 'pos')
        l_pos = len(pos)
        if l_pos == 0:
            return None, None
        elif l_pos == 1:
            return pos[0], None

        kv = zip(*pos.items())
        # Calculate Euclidean distance (2-norm of difference vectors) and get first two indexes of the sorted array.
        # Or a Euclidean-closest nodes index.
        distances = np.linalg.norm(kv[1] - curnode, ord=2, axis=1)
        i0, i1 = np.argsort(distances)[0:2]

        winner1 = tuple((kv[0][i0], distances[i0]))
        winner2 = tuple((kv[0][i1], distances[i1]))

        return winner1, winner2

class IGNG(NeuralGas):
    """Incremental Growing Neural Gas multidimensional implementation"""

    def __init__(self, data, surface_graph=None, eps_b=0.05, eps_n=0.0005, max_age=5,
                 a_mature=1, output_images_dir='images'):
        """."""

        NeuralGas.__init__(self, data, surface_graph, output_images_dir)

        self._eps_b = eps_b
        self._eps_n = eps_n
        self._max_age = max_age
        self._a_mature = a_mature
        self._num_of_input_signals = 0
        self._fignum = 0
        self._max_train_iters = 0

        # Initial value is a standard deviation of the data.
        self._d = np.std(data)

    def train(self, max_iterations=100, save_step=0):
        """IGNG training method"""

        self._dev_params = None
        self._max_train_iters = max_iterations

        fignum = self._fignum
        self._save_img(fignum, 0)
        CHS = self.__calinski_harabaz_score
        igng = self.__igng
        data = self._data

        if save_step < 1:
            save_step = max_iterations

        old = 0
        calin = CHS()
        i_count = 0
        start_time = self._start_time = time.time()

        while old - calin <= 0:
            print('Iteration {0:d}...'.format(i_count))
            i_count += 1
            steps = 1
            while steps <= max_iterations:
                for i, x in enumerate(data):
                    igng(x)
                    if i % save_step == 0:
                        tm = time.time() - start_time
                        print('Training time = {} s, Time per record = {} s, Training step = {}, Clusters count = {}, Neurons = {}, CHI = {}'.
                              format(round(tm, 2),
                                     tm / (i if i and i_count == 0 else len(data)),
                                     i_count,
                                     self.number_of_clusters(),
                                     len(self._graph),
                                     old - calin)
                              )
                        self._save_img(fignum, i_count)
                        fignum += 1
                steps += 1

            self._d -= 0.1 * self._d
            old = calin
            calin = CHS()
        print('Training complete, clusters count = {}, training time = {} s'.format(self.number_of_clusters(), round(time.time() - start_time, 2)))
        self._fignum = fignum

    def _train_on_data_item(self, data_item):
        steps = 0
        igng = self.__igng

        # while steps < self._max_train_iters:
        while steps < 5:
            igng(data_item)
            steps += 1

    def __long_train_on_data_item(self, data_item):
        """."""

        np.append(self._data, data_item)

        self._dev_params = None
        CHS = self.__calinski_harabaz_score
        igng = self.__igng
        data = self._data

        max_iterations = self._max_train_iters

        old = 0
        calin = CHS()
        i_count = 0

        # Strictly less.
        while old - calin < 0:
            print('Training with new normal node, step {0:d}...'.format(i_count))
            i_count += 1
            steps = 0

            if i_count > 100:
                print('BUG', old, calin)
                break

            while steps < max_iterations:
                igng(data_item)
                steps += 1
            self._d -= 0.1 * self._d
            old = calin
            calin = CHS()

    def _calculate_deviation_params(self, skip_embryo=True):
        return super(IGNG, self)._calculate_deviation_params(distance_function_params={'skip_embryo': skip_embryo})

    def __calinski_harabaz_score(self, skip_embryo=True):
        graph = self._graph
        nodes = graph.nodes
        extra_disp, intra_disp = 0., 0.

        # CHI = [B / (c - 1)]/[W / (n - c)]
        # Total numb er of neurons.
        #ns = nx.get_node_attributes(self._graph, 'n_type')
        c = len([v for v in nodes.values() if v['n_type'] == 1]) if skip_embryo else len(nodes)
        # Total number of data.
        n = len(self._data)

        # Mean of the all data.
        mean = np.mean(self._data, axis=1)

        pos = nx.get_node_attributes(self._graph, 'pos')

        for node, k in pos.items():
            if skip_embryo and nodes[node]['n_type'] == 0:
                # Skip embryo neurons.
                continue

            mean_k = np.mean(k)
            extra_disp += len(k) * np.sum((mean_k - mean) ** 2)
            intra_disp += np.sum((k - mean_k) ** 2)

        return (1. if intra_disp == 0. else
                extra_disp * (n - c) /
                (intra_disp * (c - 1.)))

    def _determine_closest_vertice(self, curnode, skip_embryo=True):
        """Where this curnode is actually the x,y index of the data we want to analyze."""

        pos = nx.get_node_attributes(self._graph, 'pos')
        nodes = self._graph.nodes

        distance = sys.maxint
        for node, position in pos.items():
            if skip_embryo and nodes[node]['n_type'] == 0:
                # Skip embryo neurons.
                continue
            dist = euclidean(curnode, position)
            if dist < distance:
                distance = dist

        return node, distance

    def __get_specific_nodes(self, n_type):
        return [n for n, p in nx.get_node_attributes(self._graph, 'n_type').items() if p == n_type]

    def __igng(self, cur_node):
        """Main IGNG training subroutine"""

        # find nearest unit and second nearest unit
        winner1, winner2 = self._determine_2closest_vertices(cur_node)
        graph = self._graph
        nodes = graph.nodes
        d = self._d

        # Second list element is a distance.
        if winner1 is None or winner1[1] >= d:
            # 0 - is an embryo type.
            graph.add_node(self._count, pos=copy(cur_node), n_type=0, age=0)
            winner_node1 = self._count
            self._count += 1
            return
        else:
            winner_node1 = winner1[0]

        # Second list element is a distance.
        if winner2 is None or winner2[1] >= d:
            # 0 - is an embryo type.
            graph.add_node(self._count, pos=copy(cur_node), n_type=0, age=0)
            winner_node2 = self._count
            self._count += 1
            graph.add_edge(winner_node1, winner_node2, age=0)
            return
        else:
            winner_node2 = winner2[0]

        # Increment the age of all edges, emanating from the winner.
        for e in graph.edges(winner_node1, data=True):
            e[2]['age'] += 1

        w_node = nodes[winner_node1]
        # Move the winner node towards current node.
        w_node['pos'] += self._eps_b * (cur_node - w_node['pos'])

        neighbors = nx.all_neighbors(graph, winner_node1)
        a_mature = self._a_mature

        for n in neighbors:
            c_node = nodes[n]
            # Move all direct neighbors of the winner.
            c_node['pos'] += self._eps_n * (cur_node - c_node['pos'])
            # Increment the age of all direct neighbors of the winner.
            c_node['age'] += 1
            if c_node['n_type'] == 0 and c_node['age'] >= a_mature:
                # Now, it's a mature neuron.
                c_node['n_type'] = 1

        # Create connection with age == 0 between two winners.
        graph.add_edge(winner_node1, winner_node2, age=0)

        max_age = self._max_age

        # If there are ages more than maximum allowed age, remove them.
        age_of_edges = nx.get_edge_attributes(graph, 'age')
        for edge, age in iteritems(age_of_edges):
            if age >= max_age:
                graph.remove_edge(edge[0], edge[1])

        # If it causes isolated vertix, remove that vertex as well.
        #graph.remove_nodes_from(nx.isolates(graph))
        for node, v in nodes.items():
            if v['n_type'] == 0:
                # Skip embryo neurons.
                continue
            if not graph.neighbors(node):
                graph.remove_node(node)

    def _save_img(self, fignum, training_step):
        """."""

        title='Incremental Growing Neural Gas for the network anomalies detection'

        if self._surface_graph is not None:
            text = OrderedDict([
                ('Image', fignum),
                ('Training step', training_step),
                ('Time', '{} s'.format(round(time.time() - self._start_time, 2))),
                ('Clusters count', self.number_of_clusters()),
                ('Neurons', len(self._graph)),
                ('    Mature', len(self.__get_specific_nodes(1))),
                ('    Embryo', len(self.__get_specific_nodes(0))),
                ('Connections', len(self._graph.edges)),
                ('Data records', len(self._data))
            ])

            draw_graph3d(self._surface_graph, fignum, title=title)

        graph = self._graph

        if len(graph) > 0:
            draw_graph3d(graph, fignum, clear=False, node_color=(1, 0, 0), title=title,
                         text=text)

        mlab.savefig("{0}/{1}.png".format(self._output_images_dir, str(fignum)))
        #mlab.close(fignum)




К сожалению, не доступен сервер mySQL