추천시스템 Pytorch 데이터셋과 데이터로더

Pytorch 는 Dataset 과 Dataloader 클래스를 상속해 커스텀 데이터셋과 데이터로더를 만들 수 있다. 보통 Dataset을 커스텀하여 학습할 데이터를 텐서화하며 데이터로더는 collate_fn 등을 활용해 미니배치상 부가 절차를 수행한다.

implicit feedback 데이터를 학습할 땐 포지티브, 네거티브 샘플링이 필요하며 유저와 아이템의 상호작용이 있으면 포지티브 샘플, 상호작용이 없는 모든 아이템들은 네거티브 샘플이 될 수 있다.

그간 머신러닝 딥러닝 추천 모델을 구축하면서 point-wise, pair-wise loss 를 범용적으로 적용할 수 있는 데이터셋을 디자인하는데 머리를 썼는데 이를 정리하고자 한다.

커스텀 데이터셋

implicit feedback은 웹 로그 데이터로 쉽게 수집할 수 있으며 웹 접속 세션(시간대)별로 구분하거나 환경에 따라 학습할 시계열만큼 정제할 수 있다.

ml-latest-small 데이터셋을 확인해보자.

유저아이디, 타임스탬프 컬럼으로 정렬한 모습이다. explicit feedback 이지만 평점을 매겼다면 일단 본 것으로 간주(모두 포지티브 샘플)하고 userId, movieId 컬럼만으로 학습 데이터셋을 구성할 수 있다.

explicit 데이터임을 십분 활용해 평점이 1점 이거나, 유저별로 편향을 고려하여 유저가 비교적 낮은 평점을 부여한 영화도 네거티브 샘플로 간주해볼 수 있다.

일단은 모두 포지티브 샘플이라 가정한다.

이 데이터셋의 Id 정보는 정수로 저장되어 있어 그대로 CF 모델을 만들 수도 있지만, 아이템 Id가 unique 사이즈 대비 매우 크기 때문에 sparcity가 학습 불가할 정도로 높아질 수 있어 라벨인코딩을 진행해야 한다.

인코딩 사전 정보는 따로 저장하여 추후 inference 시 인코딩 엔진으로 활용할 수 있다.

uid2idx = {v: k for k, v in enumerate(ratings_df['userId'].unique())}
iid2idx = {v: k for k, v in enumerate(ratings_df['movieId'].unique())}

# 추후 디코딩 용
idx2uid = {k: v for k, v in enumerate(ratings_df['userId'].unique())}
idx2iid = {k: v for k, v in enumerate(ratings_df['movieId'].unique())}

# 인코딩
ratings_df['user_idx'] = ratings_df['userId'].map(uid2idx)
ratings_df['item_idx'] = ratings_df['movieId'].map(iid2idx)

위와 같이 python dict를 활용하여 간단하게 인코딩 할 수 있다.

이 후 학습을 위해 train, valid split 을 진행하는데 이 포스팅에선 세션, 시계열을 나누지 않고 통째로 학습한다. (최근 시계열만 학습하는 방향도 가능하나 그만큼 학습데이터가 줄어들 수 있다. 영화 개봉연도에 따라 최신 영화에 대해서만 학습하는 방향도 있으나 모든 데이터를 학습 후 추론 시 개봉연도 기준으로 topk 를 정제할 수도 있다.)

train, valid set 은 유저별로 아이템을 그룹바이한 상태에서 랜덤하게 추출한다.

valid_df = ratings_df.groupby('user_idx').sample(n=2, random_state=42)
train_df = ratings_df.drop(valid_df.index)

# number of cold-start items
print(len(set(valid_df['item_idx']) - set(train_df['item_idx'])))

위와 같이 valid set 을 나누는 과정에서 cold-start 아이템이 생길 수 있다. CF 계열 모델의 학습데이터가 곧 라벨인 특성 때문으로 평가지표 확인 시 감안해야 한다.

이후 user_idx, item_idx 컬럼으로 구성된 데이터 프레임으로만 pytorch 데이터셋을 구성한다. train_df 에 대한 Dataset 클래스를 구성하며 유저-포지티브 아이템 페어가 나오도록 한다.

추천시스템 엔진 구축에 있어 대부분의 데이터는 유저와 아이템의 상호작용(click, purchase) 여부의 데이터만 존재할 수 있다. 이를 implicit feedback 이라 칭한다.

반대로 유저가 아이템에 매긴 직접적인 평가 점수(평점, 후기 등)를 explicit feedback 이라 칭하며 현실적으로 implicit feedback 대비 데이터를 구하기 어렵고 유저, 아이템 편향 및 신뢰도 확인 등 부가적인 정제 절차가 필요하다.

따라서 대부분의 추천시스템 엔진, 모델은 많은 경우 implicit feedback 데이터로 학습을 진행한다.

이 포스트는 이러한 implicit feedback 학습 시 필요한 Negative Sampling 기본적인 구현 방법을 다루며 속도를 높이는 것에 집중한다.

기본적인 네거티브 샘플링 기법

보통 추천시스템 딥러닝 학습 코드를 보면 미니배치를 돌 때마다 유저 번호를 확인하고, 그 유저의 history (interaction, session etc) 에 없는 아이템에서 뽑아온다.

# psuedo code for one user
import random

...
# in minibatch
neg_sample = []
while len(neg_sample) != n_negs:
    neg_candidate = random.choice(total_items)
    if neg_candidate not in user_history:
        neg_sample.append(neg_candidate)

이 방법은 네거티브 샘플을 필요할 때에만 뽑음으로써 메모리를 가장 절약하여 사용할 수 있다.

하지만 매번 전체 아이템에서 하나씩 랜덤하게 뽑아 유저 히스토리에 포함되는지 일일이 확인하기 때문에 유저 히스토리 길이에 따라 시간이 들쭉날쭉하며, 파이썬 코드 상에선 전반적으로 매우 느리다.

Jagged Array 형태로 네거티브 샘플링 pool 유지 기법

네거티브 샘플링 풀은 유저별로 정해져 있으며, 이를 통해 유저별 네거티브 샘플 테이블을 구성하여 jagged array 형태로 저장해두면 이 array 에서 random.choices 혹은 random.sample 함수를 활용해 보다 빠르게 샘플링이 가능하다.

네거티브 샘플링 풀을 메모리에 띄우기 때문에 메모리 사용량은 늘어날 수 있다.

import pandas as pd

user_negatives = {}
# create user-negative jagged array with pandas.Series
user_positives = train_data.groupby('user_id')['item_id'].apply(set)
total_items = set(train_data['item_id'])

for u, pos_list in user.positives.items():
    # neg sample pool per user
    neg_set = total_item_set - pos_list
    self.user_negatives[u] = list(neg_set)

user_negatives = pd.Serise(user_negatives)

...

neg_batch = self.user_negatives[user_ids].apply(
    lambda neg_list: random.choices(neg_list, k=n_negs)
).tolist()

위와 같이 pandas 의 apply 함수를 통해서 jagged array 를 구성하고, 이후 리스트 인덱싱을 같이 활용해 다수 유저에 대한 빠른 네거티브 샘플링이 가능하다.

모델 구축 및 속도 비교

ml-latest-small 데이터셋을 SGD-MatrixFactorizaion 으로 학습하면서 두 방법의 속도 차이가 어느정도 나는지 확인해보자. pytorch를 활용하겠다.

해당 데이터셋은 600명의 유저와 9700개의 영화로 이루어져 있다. 총 인터렉션은 100,000개 이다. valid_set 은 유저별 히스토리에서 랜덤하게 두 개씩 뽑아서 준비하였다.

import pandas as pd

ratings_df = pd.read_csv('data/ml-latest-small/ratings.csv')
ratings_df.sort_values(by=['userId', 'timestamp'], inplace=True)

valid_df = ratings_df.groupby('userId').sample(n=2, random_state=42)
train_df = ratings_df.drop(valid_df.index)

# cold-start items
len(set(valid_df['movieId']) - set(train_df['movieId']))

train_data = train_df[['user_idx', 'item_idx']]

네거티브 샘플링은 미니배치의 유저별로 시행하며, pytorch dataloader의 collate_fn 인자를 활용해 미니배치 뽑을 때 옆에 붙여주는 형태를 취했다. RNSCollateFn1 클래스는 필요할 때 일일이 뽑는 샘플링이다.

class MlDataset(Dataset):
    def __init__(self, train_data: pd.DataFrame):
        self.duets = train_data.values

    def __len__(self):
        return self.duets.shape[0]

    def __getitem__(self, index):
        return self.duets[index]

class RNSCollateFn1:
    def __init__(self, train_data: pd.DataFrame, n_negs: int):
        self.train_data = train_data
        self.total_items = train_data['item_idx'].unique()
        self.n_negs = n_negs

    def __call__(self, batch):
        batch = np.array(batch)
        neg_batch = []
        for user_idx in batch[:, 0]:
            neg_sample = []
            user_positive = \
                self.train_data[self.train_data['user_idx']==user_idx]\
                ['item_idx'].unique()
            while len(neg_sample) != self.n_negs:
                neg_candidate = random.choice(self.total_items)
                if neg_candidate not in user_positive:
                    neg_sample.append(neg_candidate)
            neg_batch.append(neg_sample)
        neg_batch = np.array(neg_batch)
        return torch.LongTensor(np.concatenate((batch, neg_batch), axis=1))

train_dataset = TrainDataset(train_data=train_data)

collate_fn = RNSCollateFn1(train_data, n_negs=5)
duet_dataloader = DataLoader(
    dataset=train_dataset,
    shuffle=True,
    batch_size=4096,
    collate_fn=collate_fn,
    num_workers=4
)

collate_fn 을 클래스로 구성한 이유는 미니배치 처리 시 train_data 를 통한 부가적인 데이터들을 collate 내부에서 만들고 관리하기 위함이다.

dataloader 가 iterate 할 때 마다 collate 함수를 포인터로써 호출하는데 _call_ 함수를 통해 클래스 객체가 함수처럼 동작하도록 디자인 할 수 있다.

SGD-MF 모델은 다음과 같이 간략히 구성하였고, BCE Loss 를 활용함으로써 SGD-MF 를 implicit feedback 을 통해 학습하는 형태로 디자인 하였다.

class MatrixFactorization(nn.Module):
    def __init__(self, n_users, n_items, embedding_size):
        super().__init__()
        self.user_embedding = nn.Embedding(n_users, embedding_size)
        self.item_embedding = nn.Embedding(n_items, embedding_size)

    def forward(self, users, items):
        user_embed = self.user_embedding(users)
        item_embed = self.item_embedding(items)
        return torch.mul(user_embed, item_embed).sum(dim=1)

배치사이즈 1024 기준 한 epoch당 12초 걸렸다. 네거티브 샘플링이 학습 프로세스보다 오래 걸린다.

네거티브 샘플링 풀을 구성하는 방식으로 collate_fn 을 교체해보자.

class RNSCollateFn2:
    def __init__(self, train_data: pd.DataFrame, n_negs: int):
        train_data = train_data[['user_idx', 'item_idx']]
        user_positives = \
            train_data.groupby('user_idx')['item_idx'].apply(set)
        total_items_set = set(train_data['item_idx'].unique())

        self.user_negatives = {}
        for u, pos_list in user_positives.items():
            neg_set = total_items_set - pos_list
            self.user_negatives[u] = list(neg_set)

        # for list-indexing
        self.user_negatives = pd.Series(self.user_negatives) 
        self.n_negs = n_negs

    def __call__(self, batch):
        batch = np.array(batch)
        neg_batch = np.array(self.user_negatives[batch[:, 0]].apply(
            lambda neg_list: random.choices(neg_list, k=self.n_negs)
        ).tolist())
        return torch.LongTensor(np.concatenate((batch, neg_batch), axis=1))

네거티브 샘플링 아이템 차이 때문에 매 에폭마다 metric 점수가 약간 낮지만 속도는 채 1초가 걸리지 않았다. epoch 을 늘리면 금방 위의 방법을 따라잡는다.

더 큰 데이터 셋으로 학습하면 시간 차이는 훨씬 더 크다. ml-1m 데이터셋은 기본 네거티브 샘플링 시간이 분 단위를 넘어가지만 개선 시 채 10초가 걸리지 않는다.

이 방법은 pytorch 딥러닝 모델 뿐만 아니라 각종 네거티브 샘플링이 필요한 모든 학습에서 응용할 수 있다.

1. 연쇄법칙으로 분해

우리는 이전포스트에서 행렬곱의 미분이 어떤식으로 나타나는지 계속 확인하였다. 그럼 이제 MSE의 미분은 어떻게 나타나는지 알아보자.

2. MSE 미분

이럴수가 말도 안되게 간결해져 버렸다. 실제로 다중회귀분석이 이루어 지는지 코드로 확인하겠다.

3. 최종 결과

import numpy as np

from sklearn.datasets import load_iris
from sklearn.linear_model import LinearRegression as SkLinear

iris = load_iris()
X, t = iris['data'], iris['target'].reshape(-1, 1)
W = np.full((X.shape[1], 1), 0.1)

# add bias
nX = np.concatenate((X, np.ones((X.shape[0], 1))), axis=1)
nW = np.concatenate((W, np.full((1,1), 0.2)), axis=0)

def forward(X, W, t):
    # dot
    Y = np.dot(X, W)
    # mse
    return np.sum(0.5 * ((Y - t) ** 2)) * 0.5

def backward(X, W, t):
    # dmse
    dY = 0.5 * (np.dot(X, W) - t)
    # ddot
    return np.dot(X.T, dY)

iter_num = 100000

for i in range(iter_num):
    loss = forward(nX, nW, t)
    if i % (iter_num//10) == 0:
        print(f'loss: {loss:.5f}, W: {nW.squeeze()}')

    # W update
    nW -= 0.0001 * backward(nX, nW, t)

sk_lr = SkLinear()
sk_lr.fit(X, t)
print('sklearn result:', sk_lr.coef_, sk_lr.intercept_)

loss: 1.74015, W: [-0.11184572 -0.04004399 0.22863744 0.60922891 0.18608863] sklearn result: [[-0.11190585 -0.04007949 0.22864503 0.60925205]] [0.18649525]

이 코드에서는 dot 부분과 mse 부분을 하나의 forward, backward 함수에서 퉁쳐 이루어지도록 했다. *실행시간이 채 10초가 되지 않는데 sklearn의 회귀분석과 거의 일치하는 결과를 낸다. *sklearn은 내부적으로 정규방정식을 사용하는 것으로 보여 더 정확하고 빠르게 계산하는 것 같다.

참고로 정규방정식은 데이터의 특징이 많을 수록 행렬의 크기가 커지기 때문에(유사역행렬 개념) 빅데이터를 다룰 땐 이와 같은 경사하강법을 사용해야 한다.

오차역전파는 연쇄법칙, 사람이 직관적으로 계산해 알 수 있는 도함수 (ex. y = x^2 -> 2x^1), 간결한 행렬 표현을 직접 대입해 활용할 수 있기 때문에 미세값을 직접 일일이 넣는 수치미분보다 훨씬 빠르게 가중치를 업데이트 할 수 있다.

이 시리즈를 포스팅하면서 행렬곱의 미분원리를 알게 되었고 동시에 오차역전파의 전반적인 구조가 어떻게 빠른 학습을 해내는지 알 수 있었다. 여기에선 다중회귀분석의 결과만 빠르게 확인하고자 퉁쳐서 만들었지만, 레이어 단위로 구현하면서 활성화함수, 소프트맥스 함수 등을 같이 구현하면 쉽게 딥러닝 뉴런을 구성할 수 있다.

이후부턴 많은 책과 강의, 블로그에서 역전파 구현을 다루고 있기 때문에 여기까지 마무리하겠다.

1. 계산그래프

계산그래프는 함수에 입력값이 들어올 때 출력을 연결리스트처럼 표현하는 기법으로, 역으로 거슬러 올라가면 그 함수의 이전에 들어왔던 입력값과 출력값에 따른 미분값을 자동으로 계산한다. 이러한 특징이 곧 오차역전파를 구현하는 핵심이 된다. 아래 그림은 다변수의 합(summation)으로 이루어진 다항식 함수, 곱(product)으로 이루어진 다항식 함수를 계산그래프로 표현한 것이다.

원 모양의 노드는 스칼라 값, 네모 모양의 노드는 연산(or 함수)이다. 간선의 윗부분은 forward 부분으로, 왼쪽에서 오른쪽으로 흐르며 원에서 나올 땐 그대로, 연산 노드를 거치면 해당 연산의 값이 흐른다.

반대로 간선의 아랫부분은 backward 부분으로, 오른쪽에서 왼쪽으로 흐르며 연산노드를 거치면 해당 연산노드의 편미분 값들이 흐른다. 주목할 것은 미분값이 흐를 때 이전 미분값을 곱하면서 흐르는데, 이는 미분의 연쇄법칙(합성함수) 때문이다. 다양한 연산들을 엮어낼 경우 가장 오른쪽 연산부터 각각 편미분 값이 흐르며 연쇄법칙을 통해 각 연산의 미분 값이 곱해지면서 지나게 된다.

1부터 흐르게 한 이유는 연산노드의 값을 그 함수에 그대로 미분했다고 가정했기 때문이다.(해당 연산노드가 계산그래프의 최종이면 1을 넣으면서 시작한다.) df(x)/df(x) = 1 이라고 생각하면 된다.

1-1. Summation, Product Node

행렬의 곱은 원소의 곱하기와 더하기 연산만 필요함으로 위 두 연산노드만 있으면 된다. 예시코드는 아래와 같다.

import numpy as np

class ProductNode:
    def __init__(self):
        self.elems = None
        self.grads = None

    # 따로 리스트로 받지 않고 매개변수 나열로 받을 수 있도록 패킹
    def forward(self, *elems):
        self.elems = elems
        return np.prod(self.elems)

    def backward(self, dout):
        self.grads = \
        [dout * np.prod(self.elems[:i] + self.elems[i + 1:]) \
        for i in range(len(self.elems))]
        return self.grads


class SumationNode:
    def __init__(self):
        self.elems = None
        self.grads = None

    def forward(self, *elems):
        self.elems = elems
        return np.sum(self.elems)

    def backward(self, dout):
        self.grads = [dout * 1 for i in range(len(self.elems))]
        return self.grads

처음 구현할 땐 두 변수만 받는 더하기, 곱하기로 구현했지만 summation, product 로 구현하면 더 많은 변수를 한번에 처리할 수 있어 이렇게 구현하였다. summation 노드는 사실 더하기로 이루어진 다항식이기 때문에 편미분 시 모두 1이다. backward 멤버함수는 dout을 받는데, dout을 곱해주어 이전에 흘러오는 미분값을 곱해 연쇄법칙을 구성한다.

1-2. Dot layer

이제 위의 스칼라 계산그래프를 이용해 행렬곱 연산을 만들어보자. 행렬곱 연산은 말그대로 행렬을 다루기 때문에 행렬 곱의 규칙에 따라 출력을 구상해야 한다. 스칼라 연산 노드들을 뭉쳐 만들었으므로 레이어라 칭하자.

2 by 2 행렬끼리의 곱은 위와같은 계산그래프 형태를 갖는다. 복잡해 보이지만 중복되는 부분이 많을 뿐이다. 행렬곱 규칙에 따라 product 노드는 8개, summation 노드는 4개가 필요하다. 입력은 행렬 X를 받고(행렬 W는 멤버변수로 존재) 행렬 Y를 출력하는 forward 함수를 구상하여 행렬 곱을 계산한다. 얼핏 보면 트리 따위의 자료구조를 만들어야 할까 싶지만, 2차원 리스트의 인덱스 접근을 통해 다중 for문으로 해결할 수 있다. 각 연산 노드들을 가지는 2차원 리스트를 만들어 계산할 원소의 위치마다 적용하도록 하면 된다.

class DotLayer:
    def __init__(self, n, W):
        self.X = None
        self.W = W
        self.n = n

        self.Ms = [[[ProductNode() for _ in range(W.shape[0])] \
        for _ in range(W.shape[1])] for _ in range(n)]
        self.Ss = [[SumationNode() for _ in range(W.shape[1])] \
        for _ in range(n)]

        self.dW = None

    def forward(self, X):
        # 행렬곱에 의해 출력될 틀
        ret = np.zeros((self.n, self.W.shape[1]), dtype=np.float32)

        for i in range(self.n):
            for j in range(self.W.shape[1]):
                outs = []
                for k in range(self.W.shape[0]):
                    outs.append(self.Ms[i][j][k].forward(X[i][k], \
                    self.W[k][j]))
                ret[i][j] = self.Ss[i][j].forward(*outs)

        return ret

    def backward(self, dout):
        # 각 원소별 미분값이 출력될 틀
        dX = np.zeros((self.n, self.W.shape[0]), dtype=np.float32)
        self.dW = np.zeros_like(self.W)

        for i in range(self.n):
            for j in range(self.W.shape[1]):
                dout_2s = self.Ss[i][j].backward(dout[i][j])
                for k, dout_2 in enumerate(dout_2s):
                    dx, dw = self.Ms[i][j][k].backward(dout_2)
                    dX[i][k] += dx
                    self.dW[k][j] += dw
        return dX

backward 함수의 경우 각 스칼라 노드로 향해가는 과정에서 결과 행렬을 행렬 X로 혹은 행렬 W로 미분하게 된다. 이전 포스트에서 보았다시피 각각 영향을 미치는 원소들의 편미분 합으로 나타낼 수 있으므로 처음에 미분값 행렬을 0을 채운 행렬로 만든 다음, 각각의 위치와 같은 원소에 backward 로 도달할 때 마다 나온 값을 0을 채운 행렬에 차례로 누적하면 된다. 이전 포스트에선 dW 만 구했지만, dX 까지 구하는 이유는 딥러닝 은닉층에서 활용할 수 있기 때문이다.

Dot 레이어를 만들었으니 결과를 테스트해보자.

X = np.array([[1,2], [3,4]])
W = np.array([[5,6], [7,8]])

al = DotLayer(2, W)
Y = al.forward(X)
dX = al.backward(np.array([[1,1],[1,1]]))

print(f'Y:\n {Y}')
print(f'X.T:\n {X.T}')
print(f'W.T:\n {W.T}')
print(f'dX:\n {dX}')
print(f'dW:\n {al.dW}')

dout으로 들어오는 행렬을 1로 구성된 행렬(Y와 사이즈가 같아야 한다.)을 넣으면 dX는 dout * W.T 와 같고, dW는 X.T * dout 임을 알 수 있다. 행렬이므로 교환법칙이 성립하지 않는다는 것을 상기하자.

1-3. Exponent node

MSE 함수를 구성하기 전에 연산노드가 하나 더 필요하다. 지수함수를 나타내는 exponent node이다. 아래 그림처럼 구성되며, 미분 대상이 밑이든 지수이든 두 경우 모두의 미분값을 계산한다.

class ExponentNode:
    def __init__(self):
        self.base = None 
        self.exponent = None

        self.dbase = None
        self.dexponent = None


    def forward(self, base, exponent):
        self.base = base
        self.exponent = exponent

        return base ** exponent


    def backward(self, dout):

        self.dbase = dout * (self.exponent * (self.base) ** (self.exponent - 1))
        # 지수로 미분 시 밑 0이거나 너무 작으면 안됨
        if self.base < 0.0:
            self.base = -self.base
        log_base = -self.base + 1e-7 if self.base <= 0 else self.base
        self.dexponent = dout * (self.base  ** self.exponent) * np.log(log_base)
        return [self.dbase, self.dexponent]

지수연산 노드가 필요한 이유는 mean Squared error 때문이다. 사실 Product Node를 적절히 활용하면 되지만 향후 시그모이드, 소프트맥스, 크로스 엔트로피 손실함수 등에서 활용할 수 있도록 만들어 두었다. 덧붙여 지수연산 노드는 밑과 지수 두 가지 변수만 받는다.

1-4. MSE layer

대망의 평균제곱오차 손실함수의 레이어를 구상하겠다.

복잡해 보이지만 중복된 곳이 많고 행렬 Y의 원소별 위치에 따라 2차원 리스트로 연산 노드들을 접근하면 된다. 정답행렬 T는 음수로 만들어 Y에 더하고, 각 원소들을 제곱하고 행 단위로 더한뒤, 다시 1/2를 곱하고 모든 열을 다 더하고(...) 데이터 개수만큼 나누어 주면 loss 값이 나온다.

class MSELayer:
    def __init__(self, T):
        self.T = T
        self.Ss_1 = [[SumationNode() for _ in range(T.shape[1])] for _ in range(T.shape[0])]
        self.Es = [[ExponentNode() for _ in range(T.shape[1])] for _ in range(T.shape[0])]
        self.Ss_2 = [SumationNode()for _ in range(T.shape[0])]
        self.Ms = [ProductNode()for _ in range(T.shape[0])]
        self.lastSum = SumationNode()
        self.lastMul = ProductNode()

    def forward(self, Y):
        last_outs = []
        for i in range(Y.shape[0]):
            outs = []
            for j in range(Y.shape[1]):
                out = self.Ss_1[i][j].forward(Y[i][j], -self.T[i][j])
                outs.append(self.Es[i][j].forward(out, 2))

            out = self.Ss_2[i].forward(*outs)
            last_outs.append(self.Ms[i].forward(out, 0.5))

        out = self.lastSum.forward(*last_outs)
        return self.lastMul.forward(out, 1.0 / len(last_outs))

    def backward(self, dout):
        dY = np.zeros_like(self.T, dtype=np.float32)

        dout = self.lastMul.backward(dout)[0]
        dout_2s = self.lastSum.backward(dout)
        for i, dout_2 in enumerate(dout_2s):
            dout_2 = self.Ms[i].backward(dout_2)[0]
            dout_3s = self.Ss_2[i].backward(dout_2)
            for j, dout_3 in enumerate(dout_3s):
                dout_3 = self.Es[i][j].backward(dout_3)[0]
                dY[i][j] = self.Ss_1[i][j].backward(dout_3)[0]

        return dY

MSE 레이어는 생성할 때 정답행렬을 넣어주고, forward 할 때 Y 행렬을 넣어주면 loss 값을 계산할 수 있다. backward는 이 레이어가 최종단에 있기 때문에 단순히 1을 넣어주면 된다.

Dot 레이어와 MSE 레이어를 완성했으니 다중회귀분석 모델을 만들어 테스트하겠다.

2. 오차역전파(스칼라 단위) 경사하강법을 통한 다중회귀분석

유틸리티가 sklearn의 LinearRegression 과 동일하도록 구성했다.

class LinearRegression:
    def __init__(self):
        self.W = None
        self.Y = None

    def get_W(self):
        return self.W

    weights = property(get_W)

    def fit(self, X, t, lrate = 0.001, bias=True, iter=1000, printIter=False):

        # preprocessing
        X = np.array(X, dtype=np.float32)
        if X.ndim != 2:
            raise Exception("X shape is not (N, 1) matrix")

        t = np.array(t, dtype=np.float32).reshape((len(t),-1))

        if len(X) != len(t):
            raise Exception("len(X) != len(t)")

        self.W = np.zeros((X.shape[1], 1))
        if bias:
            X = np.concatenate((X, np.ones((X.shape[0], 1))), axis=1)
            self.W = np.concatenate((self.W, np.zeros((1,1))), axis=0)

        # preparing layers for backprop and grad descent
        af = AffineLayer(X.shape[0], self.W)
        ms = MSELayer(t)

        out = af.forward(X)
        ms.forward(out)

        for i in range(iter):
            dout = ms.backward(1)
            af.backward(dout)

            af.W -= lrate * af.dW

            out = af.forward(X)
            loss = ms.forward(out)

            if printIter:
                if i % (iter//10) == 0 and i != 0:
                    print(f'iter: {i} loss: {loss}')
                    print(af.W.squeeze())

처음에 LinearRegression 객체를 생성하고, fit 함수를 통해 계수를 찾는다. bias의 경우 하단과 같이 두가지 구성을 통해 만들 수 있는데, 행렬에 concat하는 방법이 훨씬 간단함으로 후자의 방법을 사용한다.

sklearn dataset에 있는 붓꽃 패키지를 회귀분석 해보았다. 결과값은 다음과 같으며, 이는 sklearn의 선형회귀 방법과 상당히 비슷한 결과를 내었다. 하지만 시간은 20분이 넘게 걸려서 사실상 수치미분을 통한 경사하강법 보다 느렸다.

만든 모델 -> W = [-0.10406646 -0.03316504 0.2296437 0.60228353], b = 0.12368906

sklearn -> W = [-0.11190585 -0.04007949 0.2286450 0.60925205], b = 0.186495247206253

3. 결론

이번 포스팅에서는 행렬 곱의 미분값이 각 행렬의 전치행렬이면서 곱셈의 위치도 각각 바뀌는 것을 스칼라 단위로 구현하여 실제 그렇게 되는지 확인했던 것에 의의가 있다.

오차역전파가 딥러닝을 구현하는데 사용되는 이유는 실제 수치미분보다 빠르게, 자동으로 미분하는 것에 그 핵심이 있다. 행렬곱은 단순히 전치행렬 곱으로 나타날 수 있다는 것(추가로 MSE 역시 훨씬 간략하게 행렬 곱셈으로 나타낼 수 있다) 는 것을 여기서는 밝히기만 했지 직접 사용하지는 않았다. 다음 포스팅에선 위의 스칼라 단위가 아닌 행렬 단위로 구현해 속도를 높여보겠다.

이 포스팅은 "밑바닥부터 시작하는 딥러닝" 책을 공부하다가

𝜕L/𝜕W = 𝜕Y/𝜕W * 𝜕L/𝜕Y (연쇄법칙) = X.T * 𝜕L/𝜕Y

위 수식을 보고 이해가 안돼서 쓰기 시작한다. 왜 가중치 행렬로 손실함수를 미분하면 데이터 행렬 X의 전치행렬을 맨 앞에 곱하는 형태가 되는지?

이 글을 포스팅하는 순간의 나는 아직 다변수 미적분을 모른다. 그래서 행렬 수식을 이해할 수 있는데까지 풀어보고 이를 스칼라 단위의 계산그래프로 구현해 결과를 확인하고자 한다.

논의를 간단하게 하고자 활성화함수 등은 생략한 다중회귀분석 모형으로 진행하겠다. 어차피 오차역전파 구현 시 그대로 사용할 수 있다. 최종 손실함수는 MSE이다. (교차 엔트로피든 상관 없다.)

데이터 행렬 X는 (3 x 2) 행렬이다. 입력 데이터가 총 3개, 각 데이터는 2개의 특징을 갖는다. 배치처리를 고려한다.

가중치 행렬 W는 (2 x 3) 행렬이다. 열 값이 3임으로 Y(X * W) 행렬은 각 행마다 3개의 스칼라 출력이 나온다.

보통의 다중회귀분석이라면 W는 행렬이 아닌 열벡터로 주어지고, Y도 열벡터가 된다.

딥러닝 은닉층에서도 활용 가능하도록 W를 3개 열을 갖는 행렬로 정했다. Y 행렬에서 Y의 행들은 같은 위치의 X의 행들(각 데이터)에 대한 예측값인데, 이게 벡터로 나온다고 생각하면 된다.

은닉층을 위한 것 뿐 아니라, 다중회귀분석 이더라도 정답데이터가 one-hot 인코딩 방식 등등의 형태라면 적용될 수 있다.

다중회귀분석의 신경망을 그려보았다. Dot 구간에서는 가중치 행렬 W와 행렬곱을 통해 Y가 되고 MSE 구간에서는 정답행렬 T와 Y의 각 행벡터 별 오차제곱합들을 평균내는 과정이다. 참고로 MSE 방식과 회귀분석에서 SSR 등의 개념으로 접근할 때의 방식은 표기에 약간의 차이가 있다. 1/2 를 곱하는 이유는 one-hot 인코딩에서 최대 오차값을 1로 맞추기 위함이다.

MSE는 입력이 벡터인가, 행렬인가에 따라 내부 식 표기가 약간 달라진다. 개념적으로 원소끼리의 차이의 제곱(오차제곱)을 서로 더하고, 평균을 냄으로 행렬이 들어가도 스칼라 값이 나온다. 앞서 언급한 것 처럼 다중회귀분석은 보통 y가 벡터 형태로 나오기 때문에 식이 다르다. 정답 원소가 벡터일 경우 오차 제곱합 후 평균값을 내기 위해 한번 더 합하는 과정이 필요하다.

오차역전파는 최종 스칼라 손실값 L을 출력하는 함수 MSE를 가중치행렬 W로 편미분하는 과정에서 연쇄법칙을 이용하게 된다. L이 값이자 손실함수라고 가정하자.

다음과 같이 곱셈형태로 분할하여 나타낼 수 있다. 주의할 점은 오른쪽으로 푸는게 아니라 왼쪽으로 풀어야 한다.

경사하강법을 위해 구해야할 것은 𝜕L/𝜕W 이고, 이는 가중치 행렬 W에 편미분 행렬을 element-wise로 뺄셈하여(학습률 곱하고) 구현하게 된다.

여기서 짚고 넘어가야 할 것은 분모중심 표현이다. 스칼라 함수를 벡터나 행렬로 미분할 때 그 결과값이 분모로 들어가는 백터나 행렬의 형태를 따르는 것이다.

분자에 스칼라 함수가 아니라 벡터나 행렬이면 전치 된다. 여튼 분모(분모중심 표현에서)의 shape으로 브로드캐스팅 되는 식 인것 같다. 스칼라/벡터, 벡터/스칼라, 스칼라/행렬, 행렬/스칼라, 벡터/벡터 까지는 브로드캐스팅으로 단순하게 2차원을 떠올릴 수 있는데 벡터/행렬, 행렬/행렬은 지금 포스팅 시점에선 모르겠다. 아마도 텐서곱, 야코비 행렬 등을 공부해야 하는것 같다. 후술하겠지만 사실 텐서단위로 올라갈 필요는 없다.

2. 손실함수와 X*W 의 편미분

하여튼 가중치 행렬 W에 그대로 요소별 뺄셈을 하기 때문에 분모중심 표현으로 시작해서 연쇄미분 식에 접근해보자.

MSE를 통해 행렬 Y가 들어오고, 내부적으로 정답행렬 T와 오차제곱합의 평균을 계산한다. MSE(Y) = L 인 형태로 스칼라를 행렬로 미분하는 식이다. 이는 분모의 원소 배열로 브로드캐스팅 되어 간단하게 계산할 수 있다. 다음 단계로 넘어가자.

3. 행렬곱(X*W) 함수와 W의 편미분

이번에는 행렬의 원소가 행렬인 4차원 텐서? 같은 것이 되었다.

3-1. 원소별 편미분

다행히 위의 summation 수식을 적용하면 역시 각 원소행렬의 모든 원소의 합으로 표현된다. (참고 블로그 - 1) 갑자기 왜 편미분 값들을 합하는지 의아할 수 있는데, L값을 계산하는 전체 다항식을 펼쳐보면 이해할 수 있다.

w11은 y11, y21, y31에만 각각 영향을 주기 때문에 편미분의 합이 되는 것이다. 위 그림의 마지막 줄을 보면 행렬 곱의 결과인데, 다항식을 편미분했을 땐 더하는 과정에서 스칼라가 나왔지만, 행렬의 원소를 모두 더하기 직전의 형태를 보이기 위해 행렬의 곱으로 표현해 보았다. 행렬 곱의 경우 1행을 제외하고 모두 0이고 이들을 다 더한 것이, 하나의 미분값 원소가 되는 것을 알 수 있다.

이 신경망 그림에서 y_11 의 밑에 포개진 것들은 y_21, y_31이다. 배치를 그려내고자 이런 그림이 되었다. 그럼 이제 가중치 행렬 W로 미분한 결과를 행렬로 다시 나타내보자.

드디어 결론에 도달했다.

4. 결론

사실 여기까지 오면서 다차원의 미적분의 지식이 필요한 것이 아니었다. 그저 마지막 손실값의 다항식을 기반으로 가중치 행렬의 원소들을 편미분하는 과정에서, 연쇄법칙으로 분해 시 행렬곱으로 나타나는 부분들은 그대로 나타내고, 텐서 형태가 될 경우 다항식을 참고하여 차원을 낮추는 덧셈을 통해 스칼라 값으로 바꾸어 주었던 것이다.

행렬곱으로 나타내는 이유는 컴퓨터 연산이 행렬곱에 특화될 수 있기 때문인 것으로 보인다. 결과적으로 마지막 손실함수와 평균의 과정에서 덧셈으로 차원이 줄어들기 때문에 가능한 것 같다.

다변수 미적분과 선형대수학을 깊게 공부하고 돌아오면 또 다르게 생각할지도 모르겠지만, 적어도 지금에선 이 결론이 최선이다. 행렬로 표현하다 차원이 늘어나면 최종 다항식에 맞게 줄인다. 표기법으로써의 행렬 응용이 맞는 말이다. (참고 블로그 - 2)

길고 긴 다항식의 결과로 나오는 마지막 손실값은 스칼라이고, 그 중간 중간을 연쇄법칙에 맞게 나누어 이를 행렬의 곱으로 표현되도록 한 것이다.

단순히 shape을 같게 한다던가, 적절히 조정해줘야 해서 전치행렬이 된다든가 하는 설명은 선형대수학도 잘 모르는 나한테 전혀 와닿지 않았는데, 이번에 깊게 생각해보면서 조금은 더 이해할 수 있었다. 연립방정식, 다항식을 벡터, 행렬의 곱으로 표현하고 역으로도 쉽게 떠올릴 수 있는 훈련이 필요한 것 같다.

다음 포스트는 실제로 스칼라 단위로 계산그래프를 구현하여 결과를 확인하겠다.

att48 데이터셋, 미국 48개 주 오래된 데이터셋인데 gif 형태로 만들어 보았다. 최적경로와 1% 이내 차이를 낸다.

TSP 최적 경로를 찾는 프로젝트

전체 코드 https://github.com/p-idx/TSP-project

프로젝트 배경

TSP 경로 찾기 문제는 대표적인 NP-hard 문제다.

완전탐색은 시간복잡도 O(n!), 분기한정법 및 동적계획법은 O(2^n*n^2) 으로 다항시간 안에 최적경로를 찾는 알고리즘은 아직 없다.

도시 간 거리테이블이 주어졌을 때, 동적계획법으로 풀면 20개 도시만 해도 상용 cpu에서 3분 이상 걸린다. 30개 이상이면 연 단위로 오래걸린다.

이보다 빠르게 최적경로를 찾는 알고리즘이 아직 없기 때문에, 근사하여 경로를 구하는 방법들이 있다.

아쉽게도 최적해를 보장하지 않는다. 다항 시간안에 풀긴 해도 지도가 어떻게 생겼냐에 따라 영향을 많이 받는 것 같다.

근사알고리즘 외에도 유전알고리즘을 통해 해를 구하기도 하는데, 경로 자체를 일종의 유전자배열처럼 고려하여 부모 유전자 경로의 일정 부분을 가져와 조합하는 방식으로 최적해를 찾아나간다.

유전알고리즘 내 교차, 선택, 돌연변이 연산 등을 잘 만들고 좋은 컴퓨팅 성능이 뒷받침되면 도시가 많은 지도에서 괜찮은 해를 찾아주기도 한다. 시간을 많이 들일수록 좋은 해가 나온다.

여기까지 TSP 알고리즘들을 공부하면서

*“현실적인 시간 안에 어떻게 생긴 지도든 괜찮은 해를 찾을 수 있는가?” *

라고 고민했고 이 프로젝트를 진행하게 되었다.

프로젝트 목표

*1. 현실적인 시간 안에 * 유전 알고리즘 등 근사해 찾는 알고리즘이 사용하는 시간 대비 빠르게

2. 어떻게 생긴 지도든 + 도시 수가 많아도 일본처럼 길게 늘여지면서 대마도 같이 멀리 떨어진 섬도 있는 지도

3. 보기에 나쁘지 않은 경로를 찾자 보통 경로끼리 가로지르지 않는다 최적해가 구해진 지도가 있다면 그 거리에 최대한 근접하도록

프로젝트 내용

이 프로젝트의 해법은 다섯 단계로 이루어진다.

1. Kmeans 클러스터링을 통해 지도를 작은 지도(클러스터)로 분할

클러스터링 기법은 많지만, kmeans 클러스터링이 유클리디안 거리를 이용해 중점을 구하기 때문에 군집들의 크기가 가장 균일하다.(물론 적정한 하이퍼 파라미터 K가 주어졌을 때)

또 Gaussian mixture 모델이 비슷하게 내긴 하지만 그냥 간단명료한 kmeans를 사용하였다. 나중에 해봐야겠다.

K는 얼마나 오래 돌릴거냐에 따라 정해지는데, 여기서는 사용자가 K를 직접 정하지 않고 다음단계에서 알아서 적정 값을 고르도록 한다. 대신 K의 최대값은 받는다.

2. 클러스터들 간 동적계획법 TSP를 수행

클러스터링이 진행되면 각 클러스터들의 중점을 결과값으로 얻을 수 있는데, 이들을 대상으로 동적계획법 TSP를 수행한다. 따라서 클러스터링 할 때 K 값이 20 이하여야 현실적인 시간안에 해를 낼 수 있다.

K 값은 주어진 지도의 도시 수에 따라 정해진다. 최대값은 1단계에서 받은 K 최대값이다. 그보다 작은 값은 ‘각 클러스터 내 최대 도시의 수’ 로 나눌 수 있는 K를 정한다.

예를 들면 50개 도시가 주어지고 각 클러스터 내 최대 도시 수 가 10일 경우 K는 50 / 10 = 5 가 된다.

500개 도시가 주어질 경우, 500 / 10 = 50 이므로 사용자가 입력한 K 최대값 20을 선택한다.

K = min ( city_n / max_city_bound_in_each_cluster, user imput K )

따라서 사용자는 K 최대 바운드 값과 더불어 한 클러스터 내 최대 도시 수도 파라미터로 입력해야 한다. 추후 설명하겠지만 한 클러스터 내 최대 도시 수는 12개를 넘기기 어렵다.

3. 클러스터들을 분열

이 부분 때문에 머리가 아팠다. 원래 이 프로젝트의 큰 틀은 지도를 소지도로 나누고, 소지도 두 개 끼리 도시간 경로(path-TSP 라 한다)를 구하고, 전부 연결하는 것이었다. 아래의 그림과 같다.

근데 클러스터링 진행시 도시가 무수히 많은 지도에서 얼마나 나눠야 할지가 문제였다.

군집간 TSP에서, 시간문제로 군집 개수를 최대 20~25개 정도로 정해버리니 도시가 수백개 수천개인 지도의 군집내 도시가 너무 많아져 버렸다.

A 군집의 도시를 하나 정하고, 다음 B 군집의 도시를 하나 정한뒤, A 군집 내부의 도시를 전부 방문해서 B군집 도시까지 가는 경로를 구하는게 path-tsp 문제인데 (한 지도에서 시작점과 도착점을 정하고 모두 한번만 방문하는 경로)

이 문제는 그냥 TSP 문제보다 어렵다.

시간복잡도가 O(n!) 이므로 12를 초과하면 하루종일 돌려도 안끝난다. 결론적으로 한 군집 내 최대 도시수가 약 12 이하여야 한다.

이 문제를 해결하면서 있던 시행착오들을 정리했다.

초기에는 군집들을 12개 이하 도시를 가지도록 그냥 여러 개로 나누었다. 군집의 개수가 많아져도 이들을 동적계획법 TSP가 아니라 유전알고리즘 TSP로 해결하고자 하였다.

결국 군집들의 개수는 (도시수 / 12) 였기 때문에 천개 만개 도시의 지도는 너무 오래걸렸고 유전알고리즘 TSP 자체도 시간 대비 좋은 성능을 내지 못했다. 유전알고리즘은 내려놓고 다른 방안을 모색하였다.

또한 군집들의 도시수가 제각각인 것도 문제였다. 도시가 많은 군집이 가까우면서 도시가 적은 군집에게 하나씩 주는 것도 고려했다.

나중에 보니 same-size kmeans 클러스터링 이라고 명명된 기법들이 이러한 방식으로 구현되어 있었다. 하지만 이 방법은 군집의 모양새를 나쁘게 만들기도 하고 레이스 컨디션에 빠지기도 하기 때문에 역시 접었다. (서로 인접한 군집이 최대값으로 11개, 10개 도시를 가진다면 자기들끼리 주고받는다) 물론 안그렇게 구현하면 되지만 다른 방법이 더 나을것 같았음

결과적으로 큰 군집을 분열하는 방법을 떠올렸다. 현재 적용한 방법은 군집을 인접 군집 경로를 바탕으로 반절 잘라내는 것이다. 마치 세포분열하듯이 잘라내는데 아래 그림과 같다.

잘라낸 군집들을 연결할 땐 TSP 경로를 바탕으로 가장 짧은 경로를 찾아 다시 이어준다. 가장 도시가 많은 군집을 고르고 분열하는 과정을 최대 도시수가 이하가 될 때까지 반복한다.

이 방법을 통해 위의 문제들을 부분적으로 해결할 수 있었다.

여담) _ _작성하면서 떠올랐는데 군집을 두개로 분열하기보다, _ _분할정복 아이디어로 군집 안에서 다시 클러스터링을 진행하는 것도 나쁘지 않아보인다.

군집안에 군집안에 군집안에 … 부분적으로 구현해놓았다. _ _일단 주석으로 추가만 해놓았고 나중에 다시 시간이 나면 코드 접합해서 이걸로 해봐야겠다.

프로젝트를 할 땐 중간중간 글로 과정을 정리하는게 필요한거 같다. 정리하면서 했으면 좋은 아이디어들을 빨리 골라낼 수 있었을 듯

4. 정리된 클러스터들과 TSP 경로를 바탕으로 path-TSP 수행

한 군집의 시작점을 선택, 다음 군집의 도착점을 선택 후 군집의 모든도시를 돌고 도착점으로 가는 경로를 만들어서 이것들을 싹 붙이면 된다.

참고로 path-tsp 알고리즘은 dfs 를 통해 모든 경로를 체크하는 것으로 했으며 상태공간트리가 아래 그림처럼 그려진다.

여기에 백트래킹 아이디어로 유망한지 판단하는 과정을 넣어 시간을 줄였다. 위의 트리를 전부 순회하는것이 아니라

한번 끝까지 내려가 경로를 정했을 때, 그 거리 값을 전역에 저장해놓고 다른 경로를 탐색할 때 끝까지 내려가지 않았는데도 현재 거리가 전역 값보다 길면 그 밑의 노드 경로는 확인할 필요가 없다.

최선은 n, 최악은 n! 이다. 근데 최선과 최악의 확률은 정말 낮다. 평균적으로 시간을 반절 아낄 수 있었다.

완전탐색으로 10 ~ 11개 도시까지 현실적으로 가능했다면 백트래킹 부분을 추가해 12개까지 가능하도록 한 정도다.

다음 군집의 도착점을 어떤 로직으로 고를 수 있을까?

처음에는 다음 군집의 도시들 중에 이전 군집의 센트로이드와 가장 가까운 도시를 고르는 것이었다. 이는 얼핏 쉬워보이지만 다음 그림처럼 어이없는 문제가 발생한다.

군집의 모양새에 따라 매우 불편하게 도착점을 골라버린다.

이전 군집의 중점과 가까운 것을 고르면 길쭉한 군집에서 그림과 같은 선택을 하기 때문에 path-TSP를 수행할 때 비효율적인 경로를 만든다.

따라서 도착점을 정하는 아이디어를 다른 방법으로 생각했는데, 우선 특징들을 골라내면서 시작한다.

얻을 수 있는 정보는 다음 군집의 모양새, 다음 경로의 방향, 이전 군집이 타고 있는 경로의 방향 등이 있다.

이들을 잘 조합하면 괜찮은 도착점을 잘 고르지 않을까 고민했고 머리가 터지기 직전에 다음과 같은 아이디어로 정리하였다.

1. 다음 군집과 그다음 군집의 중점을 연결하는 직선A를 긋는다.
2. 시작 군집의 시작점과 그 군집의 중점을 연결하는 직선B를 긋는다.
3. 다음 군집에서 직선A와는 최대한 멀고, 직선B와는 최대한 가까운 점을 순위매겨 정한다.
4. 후보가 여러 개라면 시작 군집과 가장 가까운 도시로 정한다.

위 방법은 아래의 그림과 같다.

경우에 따라 비효율적으로 고르기도 하지만 군집 중점에 가까운 도착점 vs 위의 방법의 도착점 비교는 아래와 같다.

타원의 방정식과 상관계수를 활용하거나, 혹은 군집 내 도시들의 주성분 분석 등을 한다던가 하면 더 깔끔하게 군집의 모양새를 수치화할 수 있을 것 같다. 근데 그 정도로 수학이 아직 익숙치 않기 때문에...

5. 윈도우를 설정해 path-tsp를 수행

이 단계는 나중에 추가했는데, 4번의 방법도 꼬이는 경로가 나오는 한계가 있었다.

이걸 어떻게 풀지 생각했는데 근사알고리즘 중에서 꼬인 경로들을 풀어내는 걸 봤던 것 같다.

다른 일이 많아 다음에 보완하기로 하고 대신 단순한 방법으로 접근했다.

12개 이하의 구간을 한칸씩 이동하며 일일이 path-tsp를 수행하면서 풀어주면 된다. 이 단계에서는 클러스터링 정보는 더 이상 필요가 없다.

결론

이 프로젝트의 주된 로직은 군집화와 path-TSP에 있다. path-TSP 는 팩토리얼 경우의 수를 다 뒤져서 완전한 최적해를 찾는다. 핵심은 15! 보다 10! * 100,000번이 더 적다는 것이다.

시간을 더 투자하면 더 괜찮은 경로를 내지만 클러스터링이 초기에 얼마나 잘 군집화를 해주는지, 클러스터 분열이 효과적인지, 도착점 찾기가 효과적인지 등에 따라 달라진다.

** att532 set으로 테스트하였고 20분 이내 최적해보다 10% 이내의 차이를 내었다.**