R-CNN(Rich feature hierarchies for accurate object detection and semantic segmentation)

convolution, CNN이 중요하고 많이 쓰인다는 것은 알겠지만 최근 네트워크가 왜, 어떤 방향으로 발전되어오고 있는지에 대해 알려고 하지 않았던 것 같다. 따라서 YOLO에 다다르기 까지 여러 논문을 리뷰하고자 한다. 리뷰할 논문은 R-CNN이다.

R-CNN 구조

R-CNN
Object detection with R-CNN

R-CNN의 obect tetection은 3개의 모듈로 구성되어있다.

generates category-independent region proposals

large CNN extracts a fixed-length feature vector from region

set of class specific SVMs

Region proposals

Region proposal 이란 이미지에서 object가 있을 만한 영역을 찾는 것이다. 가장 Naive 하게 접근한다고 하면 sliding window 방식으로 이미지를 탐색하는 것이다.

Sliding Window

하지만 딱 보면 알 수 있듯이 굉장히 오래 걸릴것이라고 예상할 수 있다. 따라서 이미지의 모든 영역에 대해 탐색하는 것이 아닌, object가 있을만한 '후보 영역'을 탐색 빠르게 하기위해 제안된 것이 region proposal이다.

Region Proposal

Region proposal의 방법으로는 objectness, selective search, category-independent object proposals 등등 다양한 방법이 있으나 R-CNN은 selective search를 이용한다고 한다. Selective Search 도 알면 좋겠지만 논문에서 말하고자 하는 것은 아니니 일단 넘어간다.

Feature extraction

Feature Extraction

Region Proposal 이후, 후보 영역에 대해 4096-dimensional feature vector를 추출한다. 이 부분은 이해가 잘 되지 않았는데 깔끔하게 정리된 블로그가 있어서 참고했다. 링크 Region proposal에서 2천개 정도의 후보 영역을 모두 227x227(fixed-size CNN)로 warp한 후, 4096 차원의 feature vector를 추출한다.

Classification

추출된 feature vector는 SVM에 통과시켜 각 영역에 점수를 매긴다. 각 영역에 Score가 나오면 greedy non-maximum suppression을 진행한다 Non-Maximum Suppression이란 1)동일한 object에 대해 여러개의 box가 겹쳐있을 때, 2)가장 스코어가 높은 박스를 제외하고 IoU(Intersection-over-Union)가 일정값을 넘어가면 나머지 박스들을 제거하는 작업이다.

IoU

논문에선 IoU가 0.5보다 클 경우 동일 object로 결정한다.

R-CNN Limitations

R-CNN의 한계는 무엇일까?

후보 영역의 크기, 비율에 상관없이 fixed-sized CNN에 맞추기 위한 warping 으로 localization에 취약

많은 후보 영역에 대한 CNN inference로 시간이 오래걸림

Selective Search, SVM 은 GPU에 적합하지 않음

수행한 Computation Share이 없음

이정도가 R-CNN의 큰 틀이라고 한다. 정리하자면 R-CNN은 CNN을 이용한 최초의 Object Detection 방법이고 정확도, 속도를 크게 향상시켰다 이후 많은 후속 연구들이 진행되어 많이 발전되었다고 한다.

c++ smart pointer

지금까지 c를 위주로 코딩했더니 c++을 쓸줄 안다고 하면서도 생각보다 모르는 부분이 많다는 것을 깨달았다. 그래서 마침 코드 읽을 일도 생겼을 겸 따로 공부좀 해야 될 것 같다.

smart pointer

c/c++로 프로그래밍을 하다보면 메모리 이슈가 많이 생긴다. 할당한 메모리를 해제해주지 않아 생기는 메모리 누수문제도 이에 해당한다. 물론 프로그래머가 모든 경우의 수를 생각하고 프로그래밍을 할 수 있다면 큰 문제가 아닐 수 있다. 하지만 사람은 항상 실수를 하기 마련이다.

일반적으로 객체를 생성하면 함수를 빠져나가면서 소멸자의 호출로 인해 소멸된다. 하지만 포인터를 이용할 경우 객체가 이니기 때문에 소멸자가 호출되지 않는다. 그래서 생각한 것이 스마트 포인터라는 개념이다. 포인터'객체'를 만들어 자신이 가리키고 있는 데이터를 같이 delete 된다.

unique_ptr

스마트 포인터 중 1개인 unique_ptr은 메모리 관련 이슈 중 double free 버그를 해결하는 포인터 객체다. 딱 이 문제를 해결한다기 보단 객체의 소유권을 부여한다고 한다.

void function(){
    std::unique_ptr\ p_a(new A());
}

위 예제 코드에서 함수를 빠져나오게 되면 클래스 A에서 메모리할당이 일어나도 소멸자로 인해 메모리가 해제된다.

void function(){
    std::unique_ptr<A> p_a(new A());
    std::unique_ptr<A> pp_a = p_a; // compile error!
    std::unique_ptr<A> pp_a = move(p_a); // correct!
}

따라서 위 예제코드처럼 move라는 소유권 이전을 할순 있지만 pp_a에 p_a를 복사할 순 없다. 이런 기능을 이용해 프로그래머의 실수를 줄일 수 있다.

shared_ptr

shared pointer는 unique pointer와 같이 smart pointer 종류 중 하나이다. 다만 차이점은 unique_ptr 와는 다르게 두개 이상의 포인터가 단일 객체를 가리킬 수 있다

void function(){
    std::unique_ptr<A> p_a(new A());
    std::unique_ptr<A> p_a2(p_a); // compile error

    std::shared_ptr<B> p_b(new B());
    std::shared_ptr<B> p_b2(p_b);
}

unique_ptr을 이용할 경우 유일한 소유권만 인정되기 때문에 컴파일 오류가 발생한다. 한지만 shared_ptr의 경우 단일 객체를 여러 포인터가 가리킬 수 있다. 다만 객체의 소멸자는 reference count 가 0이 되어야만 호출되기 때문에 주의해야한다.

ML(머신러닝)과 inference(추론) 최대한 빨리 공부해보기(5) (수정중)

개인적으로 텐서플로우를 이용해보고 싶었지만 이런저런 이유로 MATLAB의 예제부터 진행을 해보도록 한다.

MNIST Example

% Data load
digitDatasetPath = fullfile(matlabroot,'toolbox','nnet','nndemos', ...
    'nndatasets','DigitDataset');
% Create imageData object
imds = imageDatastore(digitDatasetPath, ...
    'IncludeSubfolders',true,'LabelSource','foldernames');

일단 실험에 쓸 데이터가 필요하다. 저 fullfile이라는 함수는 해당인자들을 순차적으로 타고 들어가는 폴더 경로를 반환한다. fullfile이 반환하는 경로에 가보니 'digitTest'와 'digitTrain'이라는 excel파일과 폴더'0,1, ... 9' 가 존재했다. imageDatastore은 이미지들을 저장할 객체를 생성해주는데 굉장히 많은 일을 한번에 해주는 함수다. digitDatasetPath 하위폴더의 image들을 저장하는데 폴더이름과 같은 이름의 label을 달아준다. 폴더 이름이 '0,1, ..., 9'이기 때문에 각 이미지에 0 부터 9 까지의 label이 달린다.

figure;
perm = randperm(10000,20);
for i = 1:20
    subplot(4,5,i);
    imshow(imds.Files{perm(i)});
end

그 이후에 sample로 몇개의 이미지들을 가시화하는데 사실 큰 의미는 없는듯하다. figure창을 4x5로 분할하고 분할된 영역에 random한 20개의 파일을 띄운다.

labelCount = countEachLabel(imds)

img = readimage(imds,1);
size(img)

그리고 이미지저장소 객체에서 각 label에 해당하는 이미지가 몇개인지 센다. 해당 예제에는 각 1000개씩 존재한다. 그 후로 첫번째 이미지를 저장하고 사이즈(28x28)를 출력하는데 이것도 별 의미 없는 것 같다.

numTrainFiles = 750;
[imdsTrain,imdsValidation] = splitEachLabel(imds,numTrainFiles,'randomize');

그 후 전체 데이터들 중 train 데이터와 validation 데이터를 분리한나. 750개의 기존 imds에서 다시 두개의 이미지저장소 객체를 만드는데, train 저장소엔 각 레이블마다 750개의 이미지를 저장하고 validation 저장소에 나머지 데이터를 저장한다. 기존에 주어진 코드는 75 : 25 로 Train : Validation 셋을 나누는데 보통 train이 더 적지 않나?? 3:7로 줄였더니 정확도가 내려간다.

layers = [
    imageInputLayer([28 28 1])
    
    convolution2dLayer(3,8,'Padding','same')
    batchNormalizationLayer
    reluLayer
    
    maxPooling2dLayer(2,'Stride',2)
    
    convolution2dLayer(3,16,'Padding','same')
    batchNormalizationLayer
    reluLayer
    
    maxPooling2dLayer(2,'Stride',2)
    
    convolution2dLayer(3,32,'Padding','same')
    batchNormalizationLayer
    reluLayer
    
    fullyConnectedLayer(10)
    softmaxLayer
    classificationLayer];

그리고 드디어 네트워크의 구조를 만들어주는데 대체 문법이 이상한건지 콤마도 없고.. 그래도 한개씩 보자.

inputlayer : 28x28 이미지이고 1채널이다. gray scale이란 뜻이겠지.

convolution2dLater : filter size : 3x3, filter num : 8, padding을 이용해 output size를 input과 같게 만듬

batchNormalizationLayer :모르겠다. 나중에 보자.

ReLU Layer : activation function으로 RELU를 쓴다.

maxPooling2dLayer : down-sampling을 위해 max pooling을 이용한다. size는 2x2

fullyConnectedLayer : output 10 짜리 fully connected layer.

options = trainingOptions('sgdm', ...
    'InitialLearnRate',0.01, ...
    'MaxEpochs',4, ...
    'Shuffle','every-epoch', ...
    'ValidationData',imdsValidation, ...
    'ValidationFrequency',30, ...
    'Verbose',false, ...
    'Plots','training-progress');

다음으로는 학습을 어떻게 진행할지 옵션을 준다.

net = trainNetwork(imdsTrain,layers,options);

YPred = classify(net,imdsValidation);
YValidation = imdsValidation.Labels;

accuracy = sum(YPred == YValidation)/numel(YValidation) 

마침내 학습을 하고 학습된 모델을 이용해 classification을 진행한다. 만약 최적화를 한다고 치면 trainNetwork, classify를 직접 까보거나 대체하기 위한 함수를 만들 것이다. 정확도를 올리기 위해서는 네트워크의 구조를 변경해야 할 것이다. 예제를 한번 봤는데 실제로 의미있는 부분은 함수로 전부 가려져있다. 좀더 봐야할 듯 하다.

ML(머신러닝)과 inference(추론) 최대한 빨리 공부해보기(4) (수정중)

드디어 Convolution Neural Network 까지 왔다. 정말 많이들어본 녀석 중 하나지만 대충 이미지 좋다는 정도만 알고있었다. 정확히 어떤 기능을 하는지, 무엇이 현재 이 분야의 이슈인지 빨리 알아봐야 할 것이다. 수많은 논문중에 눈이 가는 제목 2개를 뽑아 조금 훑어봤는데 꽤나 흥미로운 것 같다.

Zhijian Liu et al., "Point-Voxel CNN for Efficient 3D Deep Learning" CVPR '19

Yue Wang et al., "Dynamic Graph CNN for Learning on Point Clouds" CVPR '19

두번째껀 읽지 않았지만 CNN을 통해 3D Deep Learning을 한다는 내용인것 같다. 사실 2D에서 3D로 확장될 수 있다는 것은 어찌보면 당연한 것일지도 모르지만 이쪽에 아예 조예가 없는 나로선 꽤나 흥미로운 것 같다. 잡소리는 여기까지.

Convolution Neural Network

이미지를 컴퓨터로 인식할 수 있을까. 물론 있다. 내가 잘모를 뿐. 아무튼 기본적인 아이디어는 이런 내용인 것 같다. 우리가 이미지를 네트워크로 학습시킬 때 각각의 픽셀을 네트워크의 인풋으로 넣어준다고 가정하자.

픽셀 단위로 학습을 진행한다면, 위와 같이 같은 A글자더라도 왼쪽으로 이동한다거나 회전된다거나 할 경우, 인식하는데 문제가 발생할 것이라는 것을 직관적으로 이해할 수 있다. 더군다나 이미지라는 것은 관측자(카메라)의 위치, 각도, 조도 등 많은 요인에 의해 영향을 받기 때문에 좀더 안정적이고 효율적인 방안이 필요하다. 그렇기 때문에 나온것이 CNN이라는 것이다. 그러면 CNN은 뭐가 어떻게 더 좋은 것일까.

Convolution Layer

일단 Convolution 이란 연산을 지칭하는 단어인데, 수학적으로 정확한 정의는 모르지만 적어도 우리가 다루려고 하는 것은 행렬조각과 Kernel이라 불리는 작은 행렬의 원소간 곱셈의 총합을 구하는 것이다.

위와 같이 이미지에 필터를 슬라이딩 시키면서 연산을 진행한다. 연산을 진행하고 결과로 나온 행렬을 feature map이라 부른다. 이때 궁금한 것은 filter 역시 학습을 통해 구해지는 것일까? 레이어마다 나오는 feature의 개수, 즉 네트워크 구조는 그냥 때려 맞추는 것일까? 아니면 합리적 결정, 혹은 학습을 통해 결정될까. 넘나 어렵다.

Pooling Layer

이렇게 결과로 나온 Feature map은 다시 Pooling Layer를 거치는데 Pooling은 또 뭘까. 찾아보니 네트워크를 학습할때 parameter가 많으면 over-fitting이 발생할 확률이 높다고한다. 이를 방지하는 것이 Pooling Layer라고 한다. Pooling 연산? 작업? 은 2개가 존재하는데 Max Pooling과 Avg Pooling이다.

그림과 같이 연산 자체는 굉장히 간단하다. 연산이 중요한 것이 아니고 데이터의 사이즈를 줄여 over-fitting을 방지하는 것 같은데 해당 내용은 over-fitting이 무엇인지 좀더 알아보면 좀더 명확히 이해할 수 있지 않을까.

CNN 구조

결국 CNN이란 위에서 언급된 Convolution layer와 Pooling layer들이 반복되고 마지막에 Fully connected layer가 붙은 형태를 띈다. 네트워크마다 어떤 이유로 차이가 있을텐데 이것들은 나중에 알아보도록 하자.

대충 이런 내용들을 훑어봤는데 역시 디테일을 이해하기 위해선 실습이 필요하지 않을까. 사실 얼마전에 Anaconda 설치했는데 뭔가 익숙하지 않아서... Matlab 예제로 mnist를 해볼지, python으로 해볼지 고민해봐야겠다.

ML(머신러닝)과 inference(추론) 최대한 빨리 공부해보기(3) (수정중)

자율주행, 혹은 robotics에서 어떤 기술이 필요할까 생각하면서 바로 CNN, 이미지/영상 처리를 볼까 생각하다가 그래도 Neural Network한번 보고 넘어가야겠다고 생각을 고쳐먹었다. 계속 이전 포스팅을 정리하지 못하고 넘어가는 느낌이지만 어쩔수 없다 생각한다.

Neural Network
Non-linear Hypothsis

지금까지 우리는 hypothesis가 linear하다는 가정하에 진행을 해왔다. 하지만 어떤 복잡한 문제는 linear하지 않을것이다. 예를들어 집값을 결정할 때 집값에 영향을 주는 feature가 100개라 가정하자. 이때 hypothesis가 linear하다는 전제가 없으면 feature space 가 너무 넓어진다. 기존 linear hypothesis를 구하는 방식을 2차, 3차,...n차로 확장하는 것은 연산량의 증가만 야기하고 복잡한 함수에 맞추다가 overfitting 문제가 발생할 수 있다.

그래서 Neural Network는 뭐가 다른다는 것일까? 유튜브 강의를 보다가 뭔가 intuition이 온 것 같다. 사실 예전에 학부강의에서 XOR classification에 대해 설명 듣긴 했지만 한귀로 듣고 흘려보냈는데...

위 그림과 같이 데이터가 분포할 때 함수로서 데이터들을 구분하기는 어렵다. 그에 비해 AND 나 OR의 경우 쉽게 구할 수 있는데 예시를 보도록하자.

만약 우리가 위처럼 AND, OR의 weight들을 찾았다 가정해보자.

결국 아이디어는 간단한 분류기들을 잘 합성한다면 복잡한 형태의 분류기를 만들 수 있다는 것이다. 이는 우리가 어떤 문제에 접근할 때 말도안되는 함수를 찾아내기보단, 네트워크의 weight을 구하는 것이 효율적이다라는 관점을 보여주는 것 같다(맞나?).

Machine Learning With Neural Network

NN을 이용해 학습을 한다고 하자. 결굴 근본은 cost function을 최소화하는 weight을 찾는 것인데 이제 이 weight들은 네트워크 layer 사이의 간선, 즉 행렬값을 찾는 것과 같아진다.

만약 위와 같은 네트워크에 대한 학습이 끝났을 땐, 5x4, 5x6, 4x5의 크기의 행렬이 구해져있을 것이다. Cost function과 back propagation을 온전히 이해해야 '학습'에 대해 안다고 말할텐데 cost function 이 굉장히 복잡하고 back propagation 역시 직관적으로 와닿지 않으니 일단 보류하고 CNN으로 넘어가보자.

ML(머신러닝)과 inference(추론) 최대한 빨리 공부해보기(2) (수정중)

Linear Regression

시간이 촉박하기 때문에 바로 Linear Regression으로 넘어가보자. 이전 포스팅에서 언급했듯이 Regression은 데이터셋에 fit한 함수 f를 찾아나가는 것이다.

Linear Regression With One Variable
Housing Prices

최대한 구현이 어떻게 될지 생각할 수 있게 예시를 사용하도록 하자. 집값에 영향을 주는 요소를 생각해보자. 1차적으로 면적이 넓을수록 가격이 높을 것이다. 그렇다면 내게 면적-집값 이라는 데이터셋이 주어질때 어떻게 학습을 시킬까?

내게 주어진 데이터셋으로 학습을 해서 어떤 모델(함수)을 만들고, 그 모델(함수)에 집의 면적을 입력했을 때 추정 집값이 나와야한다. 집값이 면적으로만 결정된다면, 즉 1차원 함수라 가정하면 다음과 같이 표현될 것이다.

그렇다면 함수를 찾는다는 것은 결국 w(weight), b(bias)를 구하는 것과 같다. regression에 대한 개념을 다시 생각해보면 실제 데이터값 y와 내가 구한 함수 h(x)과 비슷해야 한다. 즉 다음 수식이 최소가 될때의 w,b를 구해야한다.

아주 깔끔하게 요약된 영상을 캡처해서 올려보도록 한다.
출처 : 유튜브 강의

위 그림의 cost function이 최소가 되는 theta0,theta1(w,b)를 찾아야한다. 초기엔 값을 모르기 때문에 임의의 값(0,0)으로 초기화를 한다. 그리고 이런 과정을 거치길 바라는 것이다.

최초 (0,0)으로 초기화 했을때 비용이 빨간 점에서 시작했다면, theta값들을 조절하며 min이 되는 지점을 찾아가는것.

저 derivative가 의미하는 바를 생각해보자. 1변수 함수라고 가정하면 특정 theta에서의 기울기를 뜻하고 결과적으로 theta는 함수값이 감소하는 방향으로 갱신된다. 2변수 함수로 확장하면 산에서 가장 가파른 길로 내려오는 상황이 만들어지는 것이다.

ML(머신러닝)과 inference(추론) 최대한 빨리 공부해보기

Machine Learning

SLAM 포스팅을 계속 이어갈 생각이었지만 의도치 않게 머신러닝, 특히 추론과정에 대한 기회가 왔기 때문에 잠시 SLAM 포스팅을 접고 새로운 포스팅을 시작한다. 인공지능, 그리고 머신러닝이라는 용어는 컴쟁이 뿐만 아니라 모두가 들어본 기억이 있을텐데 필자도 딱 그 정도이다. 기본적인 개념과 틀은 있지만 디테일과 이 기술이 어디까지 왔는지, 주요 이슈는 모르는 상태다. 최대한 빠르게 진행을 해봐야할 것 같다.

What is Machine Learning?

A computer program is said to learn from experience E with respect to some class of tasks T and performance measure P, if its performance at tasks in T, as measured byP, improves with experience E. - T. Michell
뭔가 멋있게 적혀있지만 짧게 task T를 달성하기 위해 경험 E를 통해 성능 P를 향상시킨다면 이것은 머신러닝이다. 어떤 문제를 인공지능을 통해 해결하려고 한다면 이 3가지 요소를 잘 파악하고, 설계하란 뜻이다.

TASK

Task는 기본적으로 크게 3가지라고 한다.

classification

regression

clustering

classification은 직역한 대로 분류하는 것이다. 가장 대표적인 것으로는 손글씨로 써있는 숫자를 0~9의 숫자로 분류하는 것인데 필자도 연구실 친구가 MNIST 뭐라고 했던 것을 들어본 기억이 있다.

Regression의 경우는 classification 같이 정답이 discrete하게 주어지는 것이 아닌 연속한 값일 때, 즉 입력값과 결과값이 주어졌을 때 이를 도출해내는 함수를 구하는 과정과 같다. 근본적인 차이가 있는지는 모르겠지만 그렇다고 한다...

마지막으로 clustering은 raw data가 주어졌을 때 형성된 군집을 추출하는 과정이라 한다.

PERFORMANCE

우리가 결국 머신러닝을 통해 원하는 것은 컴퓨터가 어떤일을 인간 대신 수행했을 때 인간과 비슷하게, 혹은 보다 더 잘 하길 바라는 것이다. 그러기 위해선 컴퓨터가 학습을 잘 했는지 판단 할 수 있는 metric이 필요하고, 이 metric을 통해 학습을 진행시켜야 할 것이다.
MNIST를 다시 예로 들어보자. 필기체 숫자 4가 주어졌을 때 컴퓨터가 4라고 인식한다면 정답이지만 만약 9라고 인식했다면 에러가 발생한 것이다. 따라서 입력값 x에 대해 도출한 결과 y에 대해 정답/오답을 판별 할 수 있고, 이 정확도를 쉽게 구할 수 있다.
이와 다르게 regression의 경우 실제 데이터에 fit이 맞는 함수를 추출하는 과정이다. 이때 실제 데이터가 이쁘게 100% 어떤 함수 f와 일치하면 좋겠지만 데이터에 노이즈나 오차가 있을 수 있고, 완벽하게 함수로 표현되기 어려울 수도 있다. 따라서 실제 데이터와 함수값이 얼마나 가까운지를 측정해 이를 metric으로 사용하는 것이다.

EXPERIENCE

Experience는 간단하게 데이터라고 받아들이면 된다. 학습하기 위해선(지도학습의 경우) 입력(x)-결과(y) 세트가 있어야 그에 대해 평가할 수 있는 것이 당연하다.

이정도까지가 머신러닝의 정의라고 볼 수 있다. 다음 포스팅에서는 머신러닝이 어떻게 학습되고 추론되는지 알아 볼 수 있으면 좋겠다.

matlab SLAM(Simultaneous Localization And Mapping) example(2)

SLAM Algorithm

저번 포스팅에 이어서 실제 SLAM이 어떻게 수행되는지 바로 알아보자.

for i=1:length(pClouds)
    % Read point clouds in sequence
    pc = pClouds{i};

반복문의 시작. 매회 point clouds 1개의 cell을 로드한다. 즉 simulation이 아니라 실제 로봇이었다면 1회 scan을 진행하는 것.

% Remove invalid points outside the max range and unnecessary points
    ind = (-maxLidarRange < pc(:,1) & pc(:,1) < maxLidarRange ...
        & -maxLidarRange  < pc(:,2) & pc(:,2) < maxLidarRange ...
        & (abs(pc(:,2))>abs(0.5*pc(:,1)) | pc(:,1)>0));
    
    
    pcl = pointCloud(pc(ind,:));

주석을 보면 Lidar sensor의 range 밖의 point등과, 불필요한 point들을 필터링한다고 되어있다. matlab을 처음 접한 나로선 저 :(콜론)을 포함한 코드가 뭘 뜻하는지 찾아봐야한다. 콜론은 matlab에서 가장 유용한 연산자 중 하나로, 벡터나 첨자 배열을 만드는 데 쓴다고 한다. 이런 문법같은 것은 사실 코드를 여러번 작성하면서 익숙해지는 것이 답인데 일단 matlab 공부하는 것이 아니니 코드 뜻만 알아보도록 하겠다. pc(:,1) 은 행렬 A의 1번째 열을 뜻한다. 그런데 지금보니까 matlab은 C/C++에서와 같이 인덱싱을 0부터 하지 않는듯 하다. 그리고 지금보니 Z축이 땅에서 수직이다. 코드를 해석보면 이렇게 그릴 수 있지 않을까.

이 코드는 시뮬레이션 코드이기 때문에 Lidar sensor 범위를 조절하고 편의를 위해 정사각형 영역을 사용하고 있지만 사실은 오른쪽과같이 원과 같은 영역을 탐지할 것이다. 다만 왜 2번째 행이 1번째 행의 두배보다 커야하는걸까? 나는 단순히 1,2,3 열이 x,y,z 좌표를 의미한다고 생각했는데 좀더 말이 되는쪽으로 생각해보니 Camera Coordinate 상의 좌표인듯 하다. 하지만 그래도 2행이 1행보다 2배이상 커야 하는 이유는 모르겠다. 보류하고 넘어가자. 결국 pcl 에는 1회 filtered scan point cloud가 배열로 저장될 것이다.

% Remove the points on the ground
    [~, ~, outliers] = ...
        pcfitplane(pcl, maxDistance,referenceVector,maxAngularDistance);
    pcl_wogrd = select(pcl,outliers,'OutputSize','full');
    
    % Remove the points on the ceiling
    [~, ~, outliers] = ...
        pcfitplane(pcl_wogrd,0.2,referenceVector,maxAngularDistance);
    pcl_wogrd = select(pcl_wogrd,outliers,'OutputSize','full');

솔직히 다른 언어들을 접하는 것에 대해 크게 생각해본적 없는데 생각보다 다르니 찾아봐야하는 것이 많다. ~ 연산자는 논리적 not 이나 출력 제한을 의미한다는데 일단 처음 pcfitplane은 outlier[변수]에 outlier index를 저장한다. 그리고 논리적으로는 pcl_wogrd(설마 pcl without ground??)에 기존 pcl에서 outlier들을 제거하는게 맞는거같다. select함수는

matlab SLAM(Simultaneous Localization and Mapping) example(1)

한 6개월 전부터 SLAM에 관심이 생겨서 여기저기 유튜브강의와 Probabilistic Robotics 책을 통해 독핵을 해보았다.
뭔가 기본지식을 쌓아나가면서도 내가 제대로 이해한 것이 맞나 의문이 한켠에 남는 느낌을 지울 수가 없었다.
왜 why? 첫째로는 이쪽 분야를 처음 접했으니 여러 수식들과 개념들이 바로바로 이해가 안된다. 기본기 부족이랄까...
둘째로는 부족한 이해를 실습으로 보충하면 좋겠지만 일단 로봇이 없다. 그래서 인터넷이 돌아다니는 dataset을 활용해볼까 했지만 역시 기본기 부족으로 제대로 활용하기 어려웠다.
이렇게 처음으로 시도한 C++로 가장 기본적인 Kalman Filtering 구현이 실패해버렸다

MATLAB

SLAM에 대해 배우고 많은 실험을 해보고 싶기 때문에 뒤적거리다 보니 matlab에서 SLAM 을 simulation 할 수 있는 예제코드가 있었다.
나는 matlab을 써본적도 없고 문법도 모르지만 일단 학교 라이센스를 이용해 설치후 예제를 돌려본다.
이 포스팅 자체가 해당 예제에 대한 포스팅이 될 예정..
matlab과 SLAM 둘다 맨땅에 해딩하는 것이나 마찬가지기 때문에 험난한 여정이 예상되지만 어떻게든 되겠지라는 마인드로 시작한다.
하지만 그 만큼 포스팅 자체의 퀄리티는 일기장 수준일 듯 하다.

Perform SLAM Using 3-D Lidar Point Clouds

matlab slam을 검색하면 처음으로 나오는 예제이다.
일단 여기서 다루는 데이터부터 한번 이해해보자.

load pCloud.mat

맨 처음에 하는것이 바로 pCloud.mat이라는 데이터를 로드하는 것인데 이 데이터가 정확히 어떤 데이터를 담고 있는지 궁금하다.
실제 데이터 파일을 열게되면 다음과 같이 1x240 의 cell을 가진 행벡터로 보여진다.

그리고 그 한개의 cell에는 대량 50000x3 의 데이터가 존재하는데 로봇이 움직일때 lidar 센서로 들어온 좌표정보가 아닐까 추측한다(예제 페이지 설명보니 아마 맞는듯).

즉 240x50000개의 x,y,z 정보가 있는셈이다.

referenceVector = [0, 0, 1];    % normal vector
maxDistance = 0.5;              % for removing ceiling/ground
maxAngularDistance = 15;        %

% sampling (sampling ratio : 0.25)
% point cloud downsampled
randomSampleRatio = 0.25;

% voxel grid size in NDT(Normal Distribution Tranformation)
gridStep = 2.5;
distanceMovedThreshold = 0.3;

데이터 로드 후 이것저것 변수들을 초기화 한다. 일단 randomSampleRatio는 말 그대로 point cloud를 sampling 하는 비율이다. 모든 point cloud들을 전부 보기엔 굉장히 오래걸리니까 샘플링을 하는 듯 하다. point cloud registration의 accuracy도 올려준다는데 샘플링을 안하는게 제일 정확한 것이 아닌가??? 일단 넘어가보자. 또 gridStep이라는 변수가 있는데 내가 알기로는 point 하나하나를 이용하기 보다 여러 point들을 묶어 하나의 voxel로 만들어 이용하는 것이 효율적이라 들었다 이 gridStep이 voxelization에 관여하는 변수인듯 하다. Voxelization에 대해 자세하게 알기 위해선 Normal Distribution Transformation에 대해 알아봐야겠지만 이부분도 일단 스킵하자. 마지막으로 distanceMovedThreshold는 0.3 이상 움직인 후 scan을 수행한다는 뜻.

% Loop Closure Estimation Algorithm?
loopClosureSearchRadius = 3;

nScansPerSubmap = 3;
subMapThresh = 50;

annularRegionLimits = [-0.75,0.75];

rmseThreshold = 0.26;

loopClosureThreshold = 150;
optimizationInterval = 2;

Loop Closure가 뭔가 검색해봤다. Loop Closure(루프결합)은 정확한 지도를 작성함에 있어 중요한 역할을 한다. 루프결합이란 로봇이 지도상에 이미 등록한 위치를 다시 방문하였을 때, 현재 측정된 sensor data와 이미 작성된 map data 사이의 일치를 찾고 누적된 에러를 감소시켜 정확도를 향상하는 과정이다. 시간이 지날 수록 error가 축적되는 SLAM 알고리즘과 같은 문제에서 필수적으로 필요한 프로세스로 보인다. 여기까지 알아보고 다시 진행해보자 loopClosureSearchRadius는 동일 위치에 재방문했다는 decision을 내릴 때 현재 나의 추정위치와 지도상 위치의 거리 Threshold 역할을 하는 듯 하다. 나 자신의 이해를 돕기위해 그림을 그려보도록 하자.

로봇이 파란 곡선을 따라 이동했다고 생각해보자. 사실 저 파란 경로 역시 로봇이 자신이 이동했다고 '추정'하는 경로일 뿐이다. 현재 로봇이 파란지점에 도달했다고 해보자. 이때 저 노란 별을 관측한 데이터와 지도상 빨간점에서 관측한 데이터가 유사할 때 로봇은 '아 내가 사실은 아까 그 빨간 점까지 돌아왔구나' 라고 생각하며 경로와 지도를 수정할 수 있다.

그리고 갑자기 Submap 이란 것이 튀어나오는데 일단 정확히 파악한 것은 아니지만 이해한대로 적어본다. 로봇이 움직이면서 지도를 그리고 submap이라는 local information을 유지한다. loop closure decision을 할때 유지한 submap와 비교를 하는데 바로 직전 subMapThresh 개의 submap은 고려하지 않는다 라는 뜻인듯 하다.

어제에 이어 포스팅을 계속해보자. 일단 어제 코드를 보며 생각했던 것은 일단 SLAM 알고리즘 자체가 크게 3개로 이루어져있다.
1. Kalman filter SLAM
2. Particle filter SLAM
3. Graph SLAM
그런데 일단 이 셋중에 어떤 알고리즘이 이용되어있는지 감을 잡기 힘들다는 것. 코드를 계속 보면 감이올지 안올지 모르겠으나 일단 어제 코드에 이어서 계속 봐보도록 하자.

% 3D Posegraph object for storing estimated relative poses
% poseGraph3D object stores information for 3D pose containing nodes
% connected by edges
pGraph = poseGraph3D;

% Default serialized upper-right triangle of 6-by-6 Information Matrix
infoMat = [1,0,0,0,0,0,1,0,0,0,0,1,0,0,0,1,0,0,1,0,1];

% Number of closure edges added since last pose graph optimization and map
% refinement
numLoopClosuresSinceLastOptimization = 0; 
% True after pose graph optimization until the next scan
mapUpdated = false;
% Equals to 1 if the scan is accepted
scanAccepted = 0;

% 3D Occupancy grid object for creating and visualizing 3D map

mapResolution = 2; % cells per meter
omap = occupancyMap3D(mapResolution);

일단 실질적으로 SLAM에서 사용되는 것으로 보이는 poseGraph3D의 객체를 생성하고 information matrix 의 upper triangular matrix(상삼각행렬)을 초기화한다. 그리고 scanAccepted 라는 flag 변수가 있는 것을 보니 모든 scan이 accept되지 않는다는 것도 알 수 있다.

pcProcessed = cell(1,length(pClouds));
lidarScans2d = cell(1,length(pClouds)); 
submaps = cell(1,length(pClouds)/nScansPerSubmap);

pcsToView = cell(1,length(pClouds)); 

% Set to 1 to visualize created map and posegraph during build process
viewMap = 1; 
% Set to 1 to visualize processed point clouds during build process
viewPC = 0;

rng(0);

그 후 메모리를 미리 할당해놓는데 이때 눈에 띄는 것은 submaps 를 할당할 때 length(pClouds)/nScansPerSubmap 의 사이즈로 할당하는데 그럼 이전에 언급했던 240 x 50000 이라고 생각했던 x,y,z 좌표 정보는 사실 240번의 scan 정보에 지나지 않았던 것 아닐까? 대략 1번의 scan에 50000개의 좌표가 들어오나보다. 실제로 몇번의 scan이 add 되었나 count 하는 변수가 229까지 카운팅 된것을 보면 240개의 scan중 11개가 drop되고 229개가 accept 된거같다. 이따 확인해봐야겠다. 이제 실제로 루프를 돌며 실제로 SLAM 알고리즘이 어떻게 돌아가는지 봐야하니 이후내용은 다음 포스팅에 올린다.

GPU에서의 희소행렬 곱셈(spGEMM)(3) 수정중

희소행렬이 무엇인지, 희소행렬 곱셈이 어떻게 이루어지는지는 사실 생각보다 간단하다. 하지만 실제 데이터에 대해 연산을 한다면 이런저런 문제가 발생한다.

GPU architecture

일단 GPU를 이용해 연산한다는 가정을 깔고 있기 때문에 간단하게 복기하자면 GPU의 가장 큰 특징은 코어가 여러개란 것이다. 따라서 병렬처리에 용이하지만 만약 자원을 효율적으로 사용하지 못한다면 CPU를 이용해 계산하는 것만 못할 수 있다. CUDA terminology를 빌려 쓰자면 Thread들의 집합인 Thread block 이 GPU의 SM(Streaming Multiprocessor)에 스케줄링 되고, 블록내의 Thread들은 SP(Streaming Processor)에 의해 수행된다. 이때 32개의 Thread 단위로 같은 instruction을 수행하는데 이를 Warp이라고 부르며 GPU의 기본 수행단위이다. GPU 메모리는 크게 Global memory, L2, L1 cache, shared memory, register file로 나뉘어지는데 최대한 Global memory access를 줄이는 것이 유리하다.

Load balancing problem

희소행렬 곱셈의 수행시간은 사실 인풋행렬의 형태에 의해 크게 좌우된다. 만약 행렬의 banded matrix를 띄거나 랜덤분포를 이루면 큰 문제가 없겠지만 어떤 행렬들은 power-law degree distribution이라는 특성을 갖는다. 간단하게 말하면 행렬의 특정 row나 column에 0이 아닌 원소들이 몰려있는 분포이다. 이럴 경우 어떤 문제가 발생할까?

Row-wise product-based spGEMM의 경우 각 행마다 원소의 개수가 크게 차이나게 되면 Thread 마다 workload가 차이나게 된다. 물론 이를 해결하기 위해 따로 Thread들의 작업량을 균등 분배하기 위해 전처리를 가할 수 있지만 "전처리 시간 + 추가 memory access cost"를 고려한다면 효과가 미미하거나 오히려 손해를 볼 가능성이 있다.

Outer product-based spGEMM의 경우엔 column의 각 element에 곱해지는 B의 벡터가 같기 때문에 Thread간 load imbalance 문제는 발생하지 않는다. 하지만 power-law 특성이 더 강해지면 Outer product에선 문제가 더 심화된다.
Power-law 특성이 강해지면 Row-wise나 Outer product 방식이나 두 방식 모두 load balancing 문제가 발생하게 되지만 Outer-product의 경우 column의 원소가 모두 동일한 벡터와 곱해진다. 따라서 column, row 모두 0이 아닌 원소가 많은 밀도 높은 벡터일 경우 workload amplification이 일어난다. 즉 Thread 간 workload는 동일하지만 Thread block 간 심각한 load imbalance가 발생하는 것이다. 보통의 경우 Thread block이 충분히 많다면 Round-Robin 스케줄링으로 격차가 메워지지만 Power-law가 심한 경우엔 심각한 성능저하가 발생한다. 이를 방지하기 위해 벡터를 쪼개는 방식을 이용한다.

Too Sparse Matrix

인풋 행렬이 Power-law 특성을 가진다는 것은 몇몇 노드들은 0이 아닌 원소가 굉장히 많다는 것을 의미하면서 동시에 대부분의 노드들은 0이 아닌 원소가 거의 없다는 것을 뜻한다. GPU에서 가장 작은 수행단위가 32개의 Thread들의 집합인 warp이기 때문에 특별한 전처리를 하지 않는다면 각 연산마다 적어도 32개의 0이 아닌 원소가 필요하다. 하지만 데이터를 분석해보면 대부분의 노드의 차수(0이 아닌 원소의 개수)가 32도 안되는 경우가 많다. 이 경우 GPU의 자원이 낭비된다. Thread 가 instruction을 수행해야 하지만 issue 가능한 instruction이 없기 때문에 첫째로는 sp 코어가 놀게된다. 또 GPU는 많은 양의 register file로 각 Thread가 필요한 data들을 load 해놓고 stall이 필요할 때 context-switching을 이용하여 latency-hiding을 하지만 이 경우엔 latency-hiding을 할 수 없게된다. 이 문제를 해결하는 방법은 여러가지겠지만 Thread coarsening 이라는 직관적이고 간단한 방법으로 전처리 cost가 거의 없이 어느정도 완화할 수 있다.

Output management

희소행렬 곱셈을 수행할 때 생기는 문제들은 이것저것 있지만 사실 가장 큰 문제는 너무 큰 아웃풋 사이즈, 그리고 아웃풋을 관리하는 것이다. 첫째로 아웃풋의 원소가 몇개생길지 계산해보기 전까진 알 방법이 없다. 따라서 미리 아웃풋의 메모리를 할당받아야 하는 GPU의 경우 최대 얼만큼의 메모리가 필요한지 계산하는 과정이 필요하다. 이 과정을 거쳐야 2번째 포스팅 코드의 offset을 구할 수 있다. 두번째로 곱셈계산을 끝내면 위치가 중복된 원소들이 존재한다. Row-wise 방식의 경우 Output Matrix C의 각 행마다 겹치는 원소가 존재하고 Outer-product 방식의 경우 판단위로 겹치는 원소가 존재하는데 결국 이 원소들을 전부 합쳐서 1개의 원소로 만들어야 곱셈연산이 완벽하게 끝난 것이다. 그럼 두가지 선택을 할 수 있는데 1) 모든 원소가 계산된 후 Merge Process를 수행한다. 와 2) 계산 중 겹치는 원소들을 Merge한다. 가 있다. 개인적으로 Row-wise product의 경우 2), Outer product의 경우 1)로 노선을 잡는게 맞는 것 같다. Row-wise product의 경우는 행단위로 merge process가 이루어지기 때문에 겹치는 원소들이 이미 같은 SM에 존재하지만 Outer product의 경우 판단위 merge process가 필요하기 때문에 겹칠 원소들을 미리 모아놓지 않는다면 SM간 communication cost가 필요하고 굉장히 알고리즘적으로 복잡할 것이기 때문이다. 아웃풋을 어떻게 관리할지는 아직 연구가 필요한 부분이라고 생각되며 뚜렷한 해결책이 생각나지 않는다. 사실 merge process는 결국 memory transaction이 굉장히 많이 일어나기 때문에 최적화의 여지가 별로 없다. merge 될 원소와 merge가 필요 없는 원소를 따로 분리하고 memory 접근을 줄이는 방향으로 생각해봤지만 사실 어떤 원소가 merge될지 안될지는 찍기나 다름없다.

이외에도 matrix내의 locality를 증가시킨다거나 하는 접근법이 있긴하지만 1) 전처리 시간이 오래걸린다는 점과 random 분포에선 사실상 NP problem이기 때문에 이는 고려하지 않겠다. 희소행렬 곱셈에 대해 연구하다보면 이런 문제들이 존재하고 이를 해결할 수 있는 접근법이 필요해진다.

GPU에서의 희소행렬 곱셈(spGEMM)(2)

저번 포스팅에서는 간단하게 희소행렬을 저장할 때 사용되는 포맷에 대해 올렸다. 간단히 복기하자면 메모리 효율성을 위해 COO(Coordinate list ;Triplet) 포맷이나 CSR, CSC(Compressed Sparse Row/Column) 포맷을 이용한다는 것이다.

희소행렬 곱셈 알고리즘

보통 '행렬곱셈' 하면 내적을 이용한 곱셈을 떠올리지만 희소행렬 곱셈에서는 내적방식이 잘 쓰이지 않는다. 만약 행렬이 2D 그 이유는 내적연산을 위해선 두 벡터의 "인덱스가 같은 원소"를 곱해줘야 하지만 희소행렬은 각 행값, 또는 열값이 idx array에 저장되어있다. 따라서 idx array를 읽고 두 index가 같은지 체크해줘야하는 비효율이 생긴다.

Row-wise product

희소행렬 곱셈을 수행할 때 row-wise(row-row) product이 널리 쓰인다. Nvidia에서 제공하는 cuSPARSE, CUSP 역시 이 방식으로 곱셈을 진행한다. 구체적인 알고리즘은 위와 같이 A의 행의 0이 아닌 원소와 그 원소의 열값에 해당하는 B의 행 전체를 곱하고 결과로 C의 행 1개가 생성된다. A, B 모두 원소를 row-wise로 접근하기 때문에 두 행렬 모두 CSR 포맷을 이용한다. 간단히 Naive 한 GPU 코드를 작성해보면 다음과 같다.

// row-wise product based spGEMM
__global__ void spGEMM_expansion(
    int* a_ptr, int* a_idx, double* a_val,
    int* b_ptr, int* b_idx, double* b_val,
    int* c_ptr, int* c_idx, double* c_val)
{
    int bid = blockIdx.x;
    int tid = threadIdx.x;
    
    for(int a_iterator = a_ptr[bid]; a_iterator < a_ptr[bid + 1]; ++a_iterator)
    {
        int a_col = a_iter[a_iterator];
        
        for(int b_iterator = b_ptr[a_col]+tid; b_iterator < b_ptr[a_col + 1]; b_iterator+=blockDim.x)
        {
            // Merge?
            c_idx[offset] = b_idx[b_iterator];
            c_val[offset] = a_val[a_iterator]*b_val[b_iterator];   
        }
    }
}

간단하게 1개의 Thread block이 C의 행 1개를 생성하는 코드이다. 알고리즘 자체는 굉장히 간단하지만 몇가지 문제들이 발생한다. 이슈는 이후에 언급하도록 하겠다.

Outer product

row-wise(row-row) product 방식 이외에는 outer-product(column-row) 방식이 쓰인다. 구체적인 알고리즘은 위와 같이 A의 열과 B의 행의 곱셈으로 C의 일부가 생성된다. row-wise product 방식과 차이점은 연산의 결과가 전체 행렬 한판이라는 것이다. 이는 꽤나 중요한 차이지만 역시 이후에 언급. A의 access pattern은 column-wise, B는 row-wise이기 때문에 A는 CSC 포맷, B는 CSR 포맷을 이용한다. 간단히 Naive 한 GPU 코드를 작성해보면 다음과 같다.

// outer product based spGEMM
__global__ void spGEMM_expansion(
    int* a_ptr, int* a_idx, double* a_val,
    int* b_ptr, int* b_idx, double* b_val,
    int* c_ptr, int* c_idx, double* c_val)
{
    int bid = blockIdx.x;
    int tid = threadIdx.x;
    
    for(int a_iterator = a_ptr[bid]; a_iterator < a_ptr[bid + 1]; ++a_iterator)
    {
        int a_col = a_iter[a_iterator];
        
        for(int b_iterator = b_ptr[bid]+tid; b_iterator < b_ptr[bid + 1]; b_iterator+=blockDim.x)
        {
            // Merge?
            c_idx[offset] = b_idx[b_iterator];
            c_val[offset] = a_val[a_iterator]*b_val[b_iterator];   
        }
    }
}

간단하게 1개의 Thread block이 C행렬의 일부가 되는 C 한 판을 생성하게 된다.

GPU에서의 희소행렬 곱셈(spGEMM)(1)

컴퓨터좀 만져본 사람중 행렬 곱셈 코드를 안짜본 사람은 없다.
나도 이제 깨우치는 중이지만 상당히 많은 데이터가 행렬의 형태를 띄고, 많은 문제들이 수학적으로 모델링 될 때 행렬곱셈 연산이 포함되기 때문에 행렬 곱셈은 컴퓨터공학도에게는 꽤 중요한 연산중 하나다.

희소행렬

최근 데이터가 커지고 네트워크가 희소해지면서 희소행렬에 대한 연산이 중요해지고 있는데 희소행렬은 일반적인 행렬과 달리 2D array의 형태로 저장되지 않는다.
이유는 메모리 효율성 때문이다. 행렬이 N x N 이라면 2D array의 형태로 저장하기 위해선 N²의 메모리공간이 필요하다.

Space complexity for 2D array : \(O(N^{2})\)

따라서 스케일이 크고 희소한 행렬은 크게 두가지 포맷을 이용한다.

COO format
각각의 원소는 (row, col, value)로 표현된다. 위 행렬을 coo 포맷으로 저장하면 다음과 같다.

typedef struct _coo {
    int row;
    int col;
    double val;
}coo;

typedef struct _spmat{
    coo* elem;
}spmat;

int main(){
    spmat A;
    A.elem = (double*)malloc(sizeof(double)*4);
    A.elem[0].row = 1; A.elem[0].col = 0; A.elem[0].val = 5;
    A.elem[1].row = 1; A.elem[1].col = 1; A.elem[1].val = 8;
    A.elem[2].row = 2; A.elem[2].col = 2; A.elem[2].val = 3;
    A.elem[3].row = 3; A.elem[3].col = 1; A.elem[3].val = 6;
    return 0;
}

희소행렬을 COO format으로 저장한다면 행렬원소의 개수가 E일 경우 필요한 메모리 공간은 대략 3E이다.

Space complexity for COO : \(O(3E)\)

COO 포맷은 굉장히 직관적이고 어느정도의 메모리 효율성을 보장하지만 실제 연산에서는 자주 이용되지 않는다. 이유는 각 행이나 열에 direct access 할 수가 없기 때문이다.
실제 연산에서는 CSR(Compressed Sparse Row), CSC(Compressed Sparse Column) 포맷을 이용한다.
Compressed sparse format 는 3개의 array(ptr, idx , val)로 이루어져 있다.
CSR 포맷을 예로 들어보면

ptr array는 각 행의 첫 원소의 위치를 저장하고
idx array는 각 원소의 col을 저장하고
val array는 각 원소의 value를 저장한다.

COO 포맷과 다른 것은 ptr array 뿐이다.

#define ROW_MAJOR 0
#define COL_MAJOR 1
typedef struct _cs {
    int number_of_rows;
    int number_of_cols;
    int number_of_elems;
    int order;
    int* ptr;
    int* idx;
    double* val;
}cs;

typedef struct _spmat{
    coo* elem;
    cs csr;
    cs csc;
}spmat;

int main(){
    spmat A;
    // CSR representation
    A.csr.order = ROW_MAJOR;
    A.csr.number_of_rows = 4;   
    A.csr.number_of_cols = 4;
    A.csr.number_of_elems = 4;
    A.csr.ptr = (int*)malloc(sizeof(int)*(number_of_rows+1));
    A.csr.idx = (int*)malloc(sizeof(int)*(number_of_elems));
    A.csr.val = (double*)malloc(sizeof(double)*(number_of_elems));
    
    A.csr.ptr[0] = 0; // ptr[0] is always zero
    A.csr.ptr[1] = 0; // because there's no nonzero element in row 0
    A.csr.ptr[2] = 2; // two nonzero elements in row 1
    A.csr.ptr[3] = 3; // one nonzero elements in row 2
    A.csr.ptr[4] = 4; // one nonzero elements in row 3
    
    A.csr.idx[0] = 0;
    A.csr.idx[1] = 1;
    A.csr.idx[2] = 2;
    A.csr.idx[3] = 1;

    A.csr.val[0] = 5;
    A.csr.val[1] = 8;
    A.csr.val[2] = 3;
    A.csr.val[3] = 6
    return 0;
}

만약 위 행렬에서 row 1 에 대한 탐색을 하고싶다면 A.csr.ptr[1]을 참조하고 A.csr.ptr[2]-A.csr.ptr[1] 개의 원소를 탐색하면 된다 CSC 포맷은 CSR 포맷과 구조는 같지만 column-major order로 저장되고 ptr는 각 열의 첫번째 원소의 위치, idx array는 행 값을 저장한다.
희소행렬을 CSR/CSC format으로 저장한다면 행렬원소의 개수가 E, 행렬의 크기가 NxN일 경우 필요한 메모리 공간은 N+2E이다.

Space complexity for CSC,CSR : \(O(N+2E)\)

희소행렬 데이터 소스 : https://sparse.tamu.edu/

Matrix Multiplication Acceleration

행렬 곱셈은 단순하지만 상당히 비용이 큰 연산이다.

C = A * B 일때 C의 원소 하나를 구하기 위해선 A의 행벡터와 B의 열벡터를 곱해야 한다. 이를 모든 원소에 대하여 수행해야하기 때문에 행렬의 크기가 클 경우 컴퓨터로 연산을 하더라도 상당한 시간이 걸릴 수 있다.


다만 다행인 것은 행렬의 연산은 모든 원소에 대해 같은 방식, 다른 데이터로 구하는 것이라는 거다. 즉 SIMD task 라는 것이다. 따라서 GPU로 가속한다면 성능향상을 기대할 수 있다.

1. Naive implementation
가장 단순하게 구현하는 방법을 생각해보자. CPU에서의 연산을 그대로 사용하되 이를 GPU를 이용해 병렬적으로 계산하는 것이다.

__kernel
void matrix_multiplication1(const int M, const int N, const int K,
 const __global float* A,const __global float* B, __global float* C ) {
    const int globalRow = get_global_id(1);
    const int globalCol = get_global_id(0);

    float localsum = 0.0f;
    for (int k = 0; k < K; k++) {
        localsum += A[globalRow*K+k]*B[N*k + globalCol];
    }
    C[K*globalRow + globalCol] = localsum;
}

위는 내가 수치 컴퓨팅을 수강하자 마자 집가서 작성해본 코드이다. 간단하지만 그만큼 가속 성능이 좋지 않다. 그 이유론 잦은 Global Memory access를 생각해 볼 수 있다.
GPU에서 Global memory를 access 하는것은 상당히 오래 걸리는 작업이다. 물론 I/O 동안 context switching, caching 이 일어나지만 이로썬 부족하다.

register usage : 6(parameters) + 5(local variables)
memory access : (2K+1)*size^2[GM] for WG

그렇다면 이를 최적화 해서 성능을 끌어올릴 방법은 어떤 것이 있을까?
첫째로 Local Memory 이용을 생각해 볼 수 있다.

2.Implementation using Local Memory
위와 같은 방식이지만 Global memory를 줄이고 Local Memory 를 이용하는 것이다.

#define TS 16

__kernel 
void matrix_multiplication2(const int M, const int N, const int K,
    const __global float* A, const __global float* B, __global float* C) {
    const int row = get_local_id(1);
    const int col = get_local_id(0);
    const int globalRow = TS * get_group_id(1) + row;
    const int globalCol = TS * get_group_id(0) + col; 
 
    __local float Asub[TS][TS];
    __local float Bsub[TS][TS];

    float localsum = 0.0f;
    const int numTiles = K / TS;
    for (int i = 0; i < numTiles; i++) {
        const int tiledRow = TS * i + row;
        const int tiledCol = TS * i + col;
        Asub[row][col] = A[globalRow*K + tiledCol];
        Bsub[row][col] = B[tiledRow*N + globalCol];
  
        barrier(CLK_LOCAL_MEM_FENCE);
        for (int k = 0; k < TS; k++) {
            localsum += Asub[row][k] * Bsub[k][col];
        }
        barrier(CLK_LOCAL_MEM_FENCE);
    }

    C[globalRow*N + globalCol] = localsum;
}

Local Memory를 이용하기 위해 Global Memory에서 데이터를 복사하는 과정이 추가되었다. 또한 Local Memory를 할당받는 것은 하드웨어적 제약이 따르므로 Tiling 을 한다.
그 외 나머지 작업은 동일하다.

register usage : 6(parameters) + 9(local variables)
memory access : (2K+size)*size[GM] + 2K*size^2[LM] for WG

이로써 어느정도 최적화를 했다고 생각할 수 있다.
하지만 우린 파워 공대생이기 때문에 최적화를 멈출 수 없다.

위는 Local Memory를 이용한 행렬곱셈 코드를 PTX assembly code 로 dump 시킨 것이다. C의 원소 계산을 위한 fma instruction당 두번의 Local Memory Access가 필요하단 것이다.
그렇다면 하나의 CU에 더 효율적으로, 더 많은 일을 시킬 수 있을까?

3.More Work per Thread
나의 실력이 부족해 이 방식을 이해하고 직접 구현하는데 정말 많은 고민을 했다.
설명할 자신은 없고 다만 커널에서 원소를 한개가 아니고 여러개 계산한다는 것이 포인트다. 친절한 설명따윈 없으니 고민해 보길 바란다.

#define TS 64
#define WPT 8
#define RTS 8

__kernel 
void matrix_multiplication3(const int M, const int N, const int K,
    const __global float* A, const __global float* B, __global float* C) {
    const int row = get_local_id(1);
    const int col = get_local_id(0);
    const int globalRow = TS * get_group_id(1) + row;
    const int globalCol = TS * get_group_id(0) + col; 

    __local float Asub[TS][TS];
    __local float Bsub[TS][TS];

    float acc[8];
    for (int i = 0; i < WPT; i++)
        acc[i] = 0.0f;

    const int numTiles = K / TS;
    for (int t = 0; t < numTiles; t++) {           
        for (int w = 0; w < WPT; w++) {            
            const int tiledRow = TS * t + row;
            const int tiledCol = TS * t + col;
            Asub[row][col + w * RTS] = A[globalRow * K + (tiledCol + w * RTS)];
            Bsub[row][col + w * RTS] = B[tiledRow * N + (globalCol + w * RTS)];
        }
        barrier(CLK_LOCAL_MEM_FENCE);
  
        for (int w = 0; w < WPT; w++) {
            for (int k = 0; k < TS; k++) {
                acc[w] += Asub[row][k] * Bsub[k][col + w * RTS];
            }
        }
        barrier(CLK_LOCAL_MEM_FENCE);
    }
    for (int w = 0; w < WPT; w++) {
        C[globalRow * N + globalCol + w * RTS] = acc[w];
    }
}

이로써 우리는 clBlas와 비슷한 성능의 코드를 작성할 수 있었다. 하지만 파워 공대인은 멈추지 않는다. 최적화 방법은 상당히 많이 남아있고 아직 그 부분은 구현을 하지 않았다.
내가 생각하고 있는 것은 Transpose를 이용한 memory access pattern을 조금이나마 효율적으로 바꾸는 것과 벡터를 이용한 계산등을 시도해 볼 예정이다.

관심이 있는 사람은 https://cnugteren.github.io/tutorial/pages/page1.html 에서 많은 도움을 얻을 수 있을 것이다.