작성 중 ...
PCA L1 tracker의 메인 함수인 demo.m에서 호출하는 두 주요 함수는
estwarp_condens_PCAL1.m
sklm.m
이다.
function param = estwarp_condens_PCAL1(frm, tmpl, param, opt)
%% function param = estwarp_condens_PCAL1(frm, tmpl, param, opt)
%% frm: Frame image;
%% tmpl: PCA model;
%% -tmpl.mean: PCA mean vector
%% -tmpl.basis: PCA basis vectors
%% -tmpl.eigval: The eigenvalues corresponding to basis vectors
%% -tmpl.numsample: The number of samples
%% param:
%% -param.est: The estimation of the affine state of the tracked target
%% -param.wimg: The collected sample for update
%% opt:
%% -opt.numsample: The number of sampled candidates
%% -opt.condenssig: The variance of the Guassian likelihood function
%% -opt.ff: Forgotten factor;
%% -opt.bacthsize: The number of collected samples for update
%% -opt.affsig: The variance of affine parameters
%% -opt.tmplsize: The size of warpped image patch
%% -opt.maxbasis: The maximum number of basis vectors
%% -opt.srParam: The parameters of L1 minimization
%% -opt.threshold: The threshold for model update
%%DUT-IIAU-DongWang-2012-05-10
%%Dong Wang, Huchuan Lu, Minghsuan Yang, Online Object Tracking with Sparse
%%Prototypes, IEEE Transaction On Image Processing
%%http://ice.dlut.edu.cn/lu/index.html
%%wangdong.ice@gmail.com
%%
%%************************1.Candidate Sampling************************%%
%%Sampling Number
n = opt.numsample; % 600
%%Data Dimension
sz = size(tmpl.mean); % 32x32
N = sz(1)*sz(2); % 1024
%%Affine Parameter Sampling
% 600개의 affine 파라메터를 랜덤하게 생성
param.param = repmat(affparam2geom(param.est(:)), [1,n]);
% param.est를 600개 복사, size(param.param)=6x600
% 여기서 주의할 점은 param.est를 그대로 사용안하고
% affparam2geom으로 변환 후에 600개를 발생한다.
% affparam2geom의 출력은 [tx,ty,ramda1(수직길이),theta, ramda2/ramda1(길이비),phi]이다.
% 따라서 param.est는 [tx,ty,a11,a12,a21,a22]이다.
randMatrix = randn(6,n); % 6x600, mean(randMatrix,2)=(0.023 0.042 ...-0.036)
param.param = param.param + randMatrix.*repmat(opt.affsig(:),[1,n]);
%%Extract or Warp Samples which are related to above affine parameters
wimgs = warpimg(frm, affparam2mat(param.param), sz);
% size(wimgs)=32x32x600, size(wimgs(:,:,600))=32x32->600개 중 하나의 영상
%%************************1.Candidate Sampling************************%%
%%*******************2.Calucate Likelihood Probablity*******************%%
%%Remove the average vector or Centralizing
diff = repmat(tmpl.mean(:),[1,n]) - reshape(wimgs,[N,n]);
% size(tmpl.mean(:))=1024x1, repmat(..,[1,600])은 이것을 600개의 열로 복사
% 따라서 repmat(...)의 크기는 1024x600, 같은 이미지를 600개의 열에 복사
% size(reshape(wimgs,[N,n]))=1024x600
% size(diff)=1024x600, 변형에 따라 약간씩 평균과 다른 이미지 600개
%
if (size(tmpl.basis,2) > 0)
%%(1)PCA_L1
if (size(tmpl.basis,2) == opt.maxbasis) % maxbasis=16
%(1.1)Calucate representation coefficients of all candidates
alpha = zeros(N+opt.maxbasis,n); % N(1024), opt.maxbasis(16), n(600)
% alpha=1040x600
for num = 1:n
alpha(:,num) = pca_L1(diff(:,num), tmpl.basis, opt.srParam);
% size(diff(:,1))=1024x1, tmpl.basis=1024x16,
% opt.srParam(lambda(0.05), L0(0.05), maxLoopNum(20), tol(10e-3))
end
coeff = alpha(1:size(tmpl.basis,2),:); %%The coefficients of PCA basis vectors
% size(coeff)=16x600
err = alpha(size(tmpl.basis,2)+1:end,:); %%The coefficients of trivial templates
% size(err) = 1024x600
%%(1.2)Calucate observation likelihood via Eq.12
diff = diff-tmpl.basis*coeff-err; %%Reconstruction, diff=1024x600
diff = diff.*(abs(err)<opt.srParam.L0); % abs(err)<opt.srParam.L0 = 1024x600
param.conf = exp(-(sum(diff.^2)
+opt.srParam.L0*sum(abs(err)>=opt.srParam.L0))./opt.condenssig)';
% sum(diff.^2)=1x600, sum은 각 열의 합.
% 두 항이 다 작아야 되는데, 앞의 항은 err가 작은 항들의 합이므로 가림이 없는 부분
% 뒤 항은 err가 큰 항들의 합이니 가림 부분
% 즉, 비 가림부분은 잔차가 작아야 하고, 가림 부분은 err의 합이 작아야(sparse) 해야 함.
%%(2)Traditional PCA
% diff는 이미 평균과 임의로 발생된 600개 이미지와의 차이 이미지이다.
% 그런데 basis를 이용해서 diff를 변환한다.
%
% basis가 16개가 안되면 현재의 basis로 각 600개의 차 이미지를 표현하고
% 600개 각각에서 pca로 표현된 이미지를 뺀다.
else % basis가 16개가 안될 때 수행
coef = tmpl.basis'*diff;
% coff=""5x600"", basis=1024x5, diff=1024x600
% coff(:,1)->5차원의 축에서 1번째 영상의 좌표 (값 5개로 1번째 영상 표현)
% coff(:,2)->5차원의 축에서 2번째 " (값 5개로 2번째 영상 표현)
% size(tmpl.basis)=1024x5 (basis축에서 좌표 값으로 diff를 projection)
% basis는 하나하나가 좌표축의 방향이다.
diff = diff - tmpl.basis*coef;
% tmpl.basis*coef(:,1): 5개의 계수로 1번째 영상을 복원
% tmpl.basis*coef(:,2): 5개의 계수로 2번째 영상을 복원
param.conf = exp(-sum(diff.^2)./opt.condenssig)';
end
else % wimgs가 5개 이상이 되면 basis추출 연산이 수행되므로 wimgs가 5개가 되기 전 단계에서 수행
%%(3)Square Error
param.conf = exp(-sum(diff.^2)./opt.condenssig)'; %condenssig=0.25
% size(diff.^2)=1024x600, size(sum(diff.^2))=1x600
% param.conf=600개, 비슷한 것이 큰 값임.
end
%%*******************2.Calucate Likelihood Probablity*******************%%
%%*****3.Obtain the optimal candidate by MAP (maximum a posteriori)*****%%
param.conf = param.conf ./ sum(param.conf); % pdf를 위한 정규화
[maxprob,maxidx] = max(param.conf); % maxprob=0.2233(130번째)
param.est = affparam2mat(param.param(:,maxidx));
% param.est에 저장은 [tx,ty,a11,a12,a21,a22]형태로 한다.
%%*****3.Obtain the optimal candidate by MAP (maximum a posteriori)*****%%
%%************4.Collect samples for model update(Section III.C)***********%%
wimg = wimgs(:,:,maxidx);
if (size(tmpl.basis,2) == opt.maxbasis)
err = abs(alpha(size(tmpl.basis,2)+1:end,maxidx)); % err=1024x1
%Compute the ratio
errRatio = sum(err>opt.srParam.L0)/length(err); % 7/1024
%Full update: occlusion이나 background clutter부분이 아주 작아 전체 영역을 갱신
if (errRatio < opt.threshold.low) % 0.1
param.wimg = wimg; % 600개중의 최적의 이미지를 현 단계의 추적 이미지로 할당
return;
end
%No update:
%Occlusion이나 background clutter부분이 너무 많아 갱신을 하지 않음
if (errRatio > opt.threshold.high) % 0.6
param.wimg = [];
return;
end
%Partial update: 패치 내의 일부분의 영역만 occlusion이나 clutter임
%이 경우는 가림 영역과 가려지지 않은 영역으로 분리하여
%가려진 영역은 기존 mean영상의 부분에서 가져오고,
%가려지지 않은 부분은 현재 프레임의 최적 데이터에서 가져와 합함.
if ((errRatio>opt.threshold.low) && (errRatio<opt.threshold.high))
param.wimg = (1-(err>opt.srParam.L0)).*wimg(:) +
(err>opt.srParam.L0).*tmpl.mean(:);
% 식의 앞 부분 (1-(err>opt.srParam.L0)): occlusion이 있는 부분은 모두 0으로
% 만들고, 0이 아닌 부분만 wimg(:)의 데이터를 취함
% 뒷 부분 (err>opt.srParam.L0): 가림이 있는 부분만 취함. 이 부분은 모두 1이므로
% tmpl.mean에서 데이터를 가져옴
% 가림이 없는 부분은 현 프레임의 데이터, 있는 부분은 누적 평균 데이터에서 취함
%
% 이 함수의 호출후 메인(demo.m)으로 리턴되면 param.wimg는 누적 됨.
% 5 프레임까지는 누적 프레임의 평균을 구해 tmpl.mean에 저장(sklm.m에서 mu).
% 누적 프레임이 5개가 되면 리셋됨.
param.wimg = reshape(param.wimg,size(wimg));
return;
end
else
param.wimg = wimgs(:,:,maxidx);
end
%%************4.Collect samples for model update(Section III.C)***********%%
2014년 1월 29일 수요일
2014년 1월 22일 수요일
EKF SLAM Tutorial
SLAM
SLAM은 Simultaneous Localization And Mapping의 약자이다. 로봇이나 자율 이동체가 미지의 주변 환경을 탐색하면서 움직인다고 가정하자.
이 때 로봇은 스스로 이동하기 위한 자신의 제어 입력과 주변 환경을 계측하기 위한 센서를 가지고 있다. 제어 입력과 센서 계측 정보를 안다면 이를 이용하여 공간(환경) 내에서 자기의 위치를 인식하고 주변 환경의 지도를 동시에 작성하는 것이 가능하다. 이것이 SLAM이다.
자기 위치를 알면 환경 지도를 만들 수 있고, 환경 지도가 있으면 그 속에서 자신의 위치를 알 수 있다. 그러나, 2개의 문제가 닭과 달걀 문제처럼 서로 연관되어 있으며, 동시에 풀어야 하는 문제이다. 이것이 SLAM을 어렵게 만드는 원인이다.
SLAM에 관련된 이론은 Kalman filter, Monte Carlo localization, Particle filter등이 있으며 대표적 SLAM 알고리즘 중에서 EKF-SLAM과 FAST-SLAM에 대해 차례로 살펴 보기로 한다. 먼저 본 문서에서 EKF(Extended Kalman Filter)-slam에 대해 살펴본다.
SLAM process
SLAM을 수행하는 절차를 몇 개의 그림을 통해서 개념적으로 설명한다. 그림에서 삼각형은 로봇이고 별은 공간에 고정된 랜드 마크(표식)를 나타낸다.
로봇이 미지 환경에 놓여져 있다 하자. 이 때 장착된 센서를 이용하여 처음으로 3개의 랜드 마크의 위치를 측정하였다.
로봇이 제어 입력을 사용하여 이동한다. 제어 입력을 가한 후, 주행 거리계(odometry) 등을 통해 어느 방향으로 얼마나 로봇이 움직였는지 측정이 가능하다. 따라서 로봇은 이동 후 자신의 위치를 짐작할 수 있다.
로봇이 다시 센서를 이용하여 표식의 위치를 측정한다. 그런데, 로봇이 생각하기에 표식이 있어야 하는 위치에 있지 않다는 것을 발견한다. 이것은 다시 말하면 로봇이 자신이 생각한 위치에 있지 않다는 것이다.
로봇이 odometry 보다 센서 데이터를 더 신뢰한다고 하자(실제 로봇 바퀴는 슬립 등으로 인해 제어 입력만큼 이동이 안될 수도 있다).
따라서 자신의 위치를 더 정확하게 파악하기 위해 표식의 위치 정보를 사용하여 로봇 위치를 보정할 수 있다. 그림에서 점선 표시는 로봇이 원래 자신이 있어야 한다고 생각한 위치이다.
진짜 로봇의 위치는 실선으로 표시된 곳이다. 보정된 위치와 약간의 오차가 있다.
센서들은 완전하지 않기 때문에 로봇이 정밀하게 자신의 위치를 알 수는 없다. 그러나 상기한 방법의 위치 추정은 odometry 정보에만 의존하는 것보다 휠씬 정확하게 로봇 위치 추정을 가능하게 한다.
EKF SLAM
EKF process 요약
랜드 마크(LM: Landmark)가 추출되고, 관찰된 LM가 기존 것인지 새로운 것인지를 분석하는 데이터 연관성 해석(DA: data association)이 되었다면, SLAM과정은 3개 단계의 반복을 통해 수행된다:
Step 1: 로봇을 움직이게 하는 제어입력을 사용해서 상태 추정치를 갱신
.이전의 로봇 상태에 제어 입력을 주면 로봇의 이동(상태 변화)이 발생한다.
.만일 이전의 robot 상태가 (x, y, theta)라면, 가장 간단한 제어 입력의 한 예는 (dx, dy, dt)이다. 제어 입력 추가는 로봇이 새로운 상태 값 (x+dx, y+dy, t+dt)을 갖게 한다.
.본 문서에서는 로봇 양 바퀴의 이동 (delta sr, delta sl)이 제어 입력이다.
Step 2: 재 관찰된 LM으로부터 상태 추정치를 갱신
.재 관찰(즉, 이전에 한번 보였던 LM을 다시 관찰) LM이 상태와 상태의 공분산의 갱신에 사용된다.
.Step 1의 갱신 상태 추정치를 사용한다면 재 관찰된 LM이 어디에 있는지에 대한 계산이 가능하다. 실제 센서에서 측정된 값과 계산된 값으로부터 두 값 사이에 차이가 있다면 이 값은 innovation이 된다. 따라서 innovation이란 추정된 로봇 위치와 실제 로봇 위치 사이의 차이가 된다.
.또 Step 2에서는 상태 변수의 공분산도 최근 변화를 반영하여 갱신된다.
Step 3: 현 상태변수에 새로운 LM을 추가
.관찰된 LM이 기존에 관찰된 것이 아니라면 이것은 새로운 LM이고 상태변수에 추가된다.
시스템 상태 변수:
2차원 평면상에서 로봇이 움직인다 가정하면, 어떤 global 기준 축에 대해, 로봇의 위치와 자세는 (xr, yr, thetar)을 가진다.
또한 상태는 각 LM의 위치 (xi,yi)를 포함한다. 상태변수의 크기는 (3+2n)크기의 벡터이다. LM의 수는 n이다.
분산행렬: P
두 변수의 공분산(covariance)이란 두 변수가 상호 얼마나 강하게 상관성이 있느냐에 대한 척도를 제공한다. 이것은 변수들 사이의 선형의존성(linear dependency)을 계산하는 correlation으로 측정한다.
행렬 P는 여러 부분 행렬로 구성된다. 먼저 로봇 위치 변수 사이의 공분산, LM 사이의 공분산, LM과 로봇위치 사이의 공분산 등이 결합되어 구성된다.
위 그림은 공분산 행렬 P를 보여준다. P의 크기는 (3+2n)x(3+2n)이다.
A: 3x3, 로봇 변수(x,y,theta)에 대한 공분산
B: 2x2, 첫번째 LM에 대한 공분산. LM은 각도 정보 없이 위치만 가진다. 각 LM에 대해 대각 방향으로 마지막 LM의 공분산인 C에 도달할 때까지 반복된다.
D: 2x3, 로봇 상태와 첫번째 LM사이의 공분산을 포함한다.
E: 3x2, 첫번째 LM과 로봇 상태에 대한 공분산으로 D와 전치(transpose) 관계이다.
F: 2x2, 마지막 LM과 첫번째 LM 사이의 공분산이다. G는 F와 전치이다.
초기에는 로봇이 어떤 LM도 관찰하지 않았으므로 P는 행렬 A만 포함한다. 초기 행렬은 대각방향으로 어떤 기본 값을 가지도록 초기화 될 수 있다. 이 값은 초기 위치에 대한 불확실성을 반영한다 (불확실성이 없더라도 초기화가 필요하다. 따라오는 여러 행렬 연산의 singularity를 방지한다)
칼만게인: K
칼만게인은 관찰된 LM을 얼마나 신뢰할 수 있는지, 즉 관찰된 데이터를 얼마나 반영해야 하는지의 신뢰도를 나타낸다.
센서는 정확하지 않고 불확실성을 가지며, 로봇의 odometry도 정확하지 않다. 따라서, 상태 갱신을 위해서 두 데이터를 절충할 수 있다.
만일 센서(레인지 측정장치)가 로봇의 odometry에 비해 부정확하다면 센서 데이터는 신뢰하기 어렵고 칼만게인은 낮아진다. 측정 장치가 우수하고 정확하다면 칼만게인은 높아질 것이다. 칼만게인 행렬은 (3+2n)x2의 크기로 다음과 같다:
행렬의 첫번째 행은 상태변수의 첫번째 행에 기여한다. 즉, 로봇 위치변수 x에 innovation을 얼마나 반영해야 하는지를 나타낸다. 레이저 스캐너 센서의 경우 innovation은 2x1 크기인 (LM 거리, LM 각도)'이므로 게인 행렬의 첫 행의 첫번째 열(r: range)은 LM거리 값을 얼마나 반영해야 할지, 두번째 열(b: bearing)은 LM각도를 얼마나 반영해야 할지를 나타낸다. 둘 다는 상태변수의 첫번째 행에 대한 것이므로 로봇 위치 x에 대한 것이다.
센서모델 및 자코비안: H
먼저 센서 모델은 다음과 같다:
여기서 lambda x는 LM의 x방향의 좌표 위치, x는 로봇 위치의 추정치이다. 이 식은 로봇 위치에서 측정했을 때, LM에 대한 거리(range)와 각도(bearing)에 대한 예측치를 준다. 거리는 로봇과 LM의 직선거리이고, 각도는 기준축(수평축)에 대한 LM의 각도에서 로봇의 자세각을 뺀 값이다. v는 노이즈이다.
이 센서 모델의 x, y, theta에 대한 자코비안은 다음과 같다:
H행렬은 (x,y,theta)가 변할 때, 얼마나 많이 range와 bearing이 변할지를 보여준다.
1행의 첫번째 요소는 x축의 변화에 대한 range의 변화이다. 두번째 요소는 y축에서 변화에 대한 range의 변화에 대한 것이다.
1행의 마지막 요소는 theta의 변화에 대한 것이다. 물론 로봇이 회전할 때 range변화는 없으므로 이것은 0이다.
실제 SLAM에서는 상태 변수에 LM가 포함되어 있으므로, H행렬은 LM에 대한 항들도 포함해야 한다.
한 예로 위 행렬은 우리가 관찰한 센서 정보가 LM 2번일 때의 h및 자코비안 H행렬이다.
HLM2행렬의 첫 3열은 바로 앞에서 나온 H행렬이다.
그리고 LM 하나당 두개의 열이 추가된다. LM 2에 대한 H행렬은 해당 LM 위치인 6과 7열에서만 값을 가지고, 나머지는 모두 0이다. 즉, hLM2가 LM 2에 대한 정보만을 포함하므로 LM 2를 제외한 나머지 LM 상태변수에 대한 h의 도함수는 모두 0이다.
로봇의 이동 모델
제어 입력에 대한 로봇의 이동 모델은 다음과 같다.
여기서 delta sl, sr은 로봇의 왼쪽 바퀴와 오른쪽 바퀴가 회전하여 이동한 양이다. b는 좌우 바퀴 사이의 거리이다. 로봇의 이동 후 위치는 xr(k)=xr(k-1)+f(xr,u)이다.
로봇 상태에 대한 이동 모델의 자코비안은 다음과 같다:
프로세스 노이즈: Q
제어 입력의 노이즈는 좌우 바퀴의 입력에 대한 불확실성이다.
kl과 kr는 좌, 우 바퀴의 에러 상수이다.
제어 입력에 대한 이동 모델의 자코비안: W
이동 모델의 제어 입력에 대한 1차 도함수 행렬로 구성된다:
센서 측정 오차에 대한 공분산 행렬: R
range와 bearing을 측정하는 레이저스캐너는 오차를 가진다. 상수 c, d는 센서의 정확도를 나타내는 값으로 range값이 1cm의 분산을 가진다면 c=0.01이다. 만일 각도가 1도의 분산을 가진다면 d=1이다.
EKF SLAM의 세 단계
Step 1: odometry 데이터를 사용한 상태의 갱신
실제 SLAM에서 상태 변수는 LM을 포함하므로 상태 추측 모델(이동 모델)은 로봇 변수와 LM을 모두 포함하며 다음과 같다:
LM의 위치는 stationary model이다.
공분산 행렬 P의 갱신에 필요한 상태 추측 모델의 상태변수에 대한 자코비안을 구하면 다음과 같다:
행렬 Ak의 크기는 (3+2n)x(3+2n)이다.
제어 입력 (delta sr, delta sl)'에 대한 이동 모델의 자코비안은 다음과 같다:
공분산 행렬 P의 갱신 식은 다음과 같다:
Step 2: 재 관찰된 LM로부터 상태 갱신
로봇이 제어 입력에 의해 이동 후 LM가 하나 관찰되었다면 이것을 이용하여 상태 변수와 상태의 공분산 행렬을 갱신할 수 있다. 즉, 로봇에 설치된 센서에서 LM와의 거리(r: range)와 각도(b: bearing) 값의 벡터 z를 계측할 수 있다.
그런데 이 LM이 이미 한번 관측한 것이라면 로봇의 이동 후 (r, b)의 예측도 가능하다. 예측 값을 벡터 h라 하자. 그러면, z와 h의 차이를 이용해서 상태값과 공분산 행렬을 갱신하는 것이 가능하다. 단, 이 때 칼만 게인을 통해 로봇의 이동(예측) 모델과 센서의 측정 모델의 불확실성을 반영해 준다.
먼저 칼만 게인 식은 다음과 같다:
여기서, 게인 행렬 계산에 사용된 P, H, R 등의 행렬 크기가 표시되었고, S 및 K행렬도 크기를 나타내었다. 게인 계산에 사용되는 행렬들은 모두 앞 부분에 유도 되어 있다. R이 커지면 게인은 작아진다.
일단 칼만 게인이 계산되면 상태와 공분산 행렬의 갱신이 가능하다:
Step 3: 현 상태에 새로운 LM의 추가
만일 관찰한 LM가 기존에 관찰된 LM가 아니라면 이 LM는 새로운 것으로 상태변수에 추가된다. 또한 상태의 크기가 커지면 공분산 행렬도 커진다.
새로운 LM 위치를 y라 하자. y는 LM가 로봇에 의해 관찰된 시점의 값이다. 새로운 상태와 공분산 행렬은 다음과 같다:
여기서, r은 측정된 LM의 range값이고, b는 bearing 값이다.
Jxr는 y의 로봇 상태에 대한 변화율(자코비안)이고, Jz는 y의 센서 측정값 z에 대한 변화율이다.
센서와 로봇은 각각 불확실성을 나타내는 분산 값 R과 Pxx를 가지므로 이 값들이 LM 위치의 불확실성으로 전파된 값이 공분산 행렬에 추가된다.
.Jz: 센서 측정치 (r, b)의 변화에 대한 LM y의 자코비안
.Jxr: 로봇 상태 (xr,yr,thetar)의 변화에 대한 LM y의 자코비안
.(dy/dxr)Pxx: 불확실성 Pxx는 xr의 것이고 이것이 y에 영향을 줌
.(dy/dz)R: 불확실성 R은 z의 것이고 이것이 y에 영향을 줌
(c.f.)
(df/du)Q: 불확실성 Q는 제어입력 u에 대한 것이고 로봇 이동모델 f에 영향을 줌
JxrPxxJxr'는 로봇 공분산이 LM 공분산으로 전파된 것이고, JzRJz'는 센서의 불확실성이 LM 공분산으로 확산된 것을 보여준다. LM 값에 로봇 위치와 센서 계측이 둘 다 영항을 주므로 두 전파된 공분산의 합이 대각방향에 추가되는 공분산 값이다.
JxrPx(x~yn)은 추가 설명...
참고 자료
[1] SLAM for dummies, MIT Open Course ware.
[2] Open Robotics blog, SLAM part.
[3] Cyrill Stachniss, EKF SLAM 자료
SLAM은 Simultaneous Localization And Mapping의 약자이다. 로봇이나 자율 이동체가 미지의 주변 환경을 탐색하면서 움직인다고 가정하자.
이 때 로봇은 스스로 이동하기 위한 자신의 제어 입력과 주변 환경을 계측하기 위한 센서를 가지고 있다. 제어 입력과 센서 계측 정보를 안다면 이를 이용하여 공간(환경) 내에서 자기의 위치를 인식하고 주변 환경의 지도를 동시에 작성하는 것이 가능하다. 이것이 SLAM이다.
자기 위치를 알면 환경 지도를 만들 수 있고, 환경 지도가 있으면 그 속에서 자신의 위치를 알 수 있다. 그러나, 2개의 문제가 닭과 달걀 문제처럼 서로 연관되어 있으며, 동시에 풀어야 하는 문제이다. 이것이 SLAM을 어렵게 만드는 원인이다.
SLAM에 관련된 이론은 Kalman filter, Monte Carlo localization, Particle filter등이 있으며 대표적 SLAM 알고리즘 중에서 EKF-SLAM과 FAST-SLAM에 대해 차례로 살펴 보기로 한다. 먼저 본 문서에서 EKF(Extended Kalman Filter)-slam에 대해 살펴본다.
(from wikipedia)
SLAM process
SLAM을 수행하는 절차를 몇 개의 그림을 통해서 개념적으로 설명한다. 그림에서 삼각형은 로봇이고 별은 공간에 고정된 랜드 마크(표식)를 나타낸다.
로봇이 미지 환경에 놓여져 있다 하자. 이 때 장착된 센서를 이용하여 처음으로 3개의 랜드 마크의 위치를 측정하였다.
로봇이 제어 입력을 사용하여 이동한다. 제어 입력을 가한 후, 주행 거리계(odometry) 등을 통해 어느 방향으로 얼마나 로봇이 움직였는지 측정이 가능하다. 따라서 로봇은 이동 후 자신의 위치를 짐작할 수 있다.
로봇이 다시 센서를 이용하여 표식의 위치를 측정한다. 그런데, 로봇이 생각하기에 표식이 있어야 하는 위치에 있지 않다는 것을 발견한다. 이것은 다시 말하면 로봇이 자신이 생각한 위치에 있지 않다는 것이다.
로봇이 odometry 보다 센서 데이터를 더 신뢰한다고 하자(실제 로봇 바퀴는 슬립 등으로 인해 제어 입력만큼 이동이 안될 수도 있다).
따라서 자신의 위치를 더 정확하게 파악하기 위해 표식의 위치 정보를 사용하여 로봇 위치를 보정할 수 있다. 그림에서 점선 표시는 로봇이 원래 자신이 있어야 한다고 생각한 위치이다.
진짜 로봇의 위치는 실선으로 표시된 곳이다. 보정된 위치와 약간의 오차가 있다.
센서들은 완전하지 않기 때문에 로봇이 정밀하게 자신의 위치를 알 수는 없다. 그러나 상기한 방법의 위치 추정은 odometry 정보에만 의존하는 것보다 휠씬 정확하게 로봇 위치 추정을 가능하게 한다.
EKF SLAM
EKF process 요약
랜드 마크(LM: Landmark)가 추출되고, 관찰된 LM가 기존 것인지 새로운 것인지를 분석하는 데이터 연관성 해석(DA: data association)이 되었다면, SLAM과정은 3개 단계의 반복을 통해 수행된다:
Step 1: 로봇을 움직이게 하는 제어입력을 사용해서 상태 추정치를 갱신
.이전의 로봇 상태에 제어 입력을 주면 로봇의 이동(상태 변화)이 발생한다.
.만일 이전의 robot 상태가 (x, y, theta)라면, 가장 간단한 제어 입력의 한 예는 (dx, dy, dt)이다. 제어 입력 추가는 로봇이 새로운 상태 값 (x+dx, y+dy, t+dt)을 갖게 한다.
.본 문서에서는 로봇 양 바퀴의 이동 (delta sr, delta sl)이 제어 입력이다.
Step 2: 재 관찰된 LM으로부터 상태 추정치를 갱신
.재 관찰(즉, 이전에 한번 보였던 LM을 다시 관찰) LM이 상태와 상태의 공분산의 갱신에 사용된다.
.Step 1의 갱신 상태 추정치를 사용한다면 재 관찰된 LM이 어디에 있는지에 대한 계산이 가능하다. 실제 센서에서 측정된 값과 계산된 값으로부터 두 값 사이에 차이가 있다면 이 값은 innovation이 된다. 따라서 innovation이란 추정된 로봇 위치와 실제 로봇 위치 사이의 차이가 된다.
.또 Step 2에서는 상태 변수의 공분산도 최근 변화를 반영하여 갱신된다.
Step 3: 현 상태변수에 새로운 LM을 추가
.관찰된 LM이 기존에 관찰된 것이 아니라면 이것은 새로운 LM이고 상태변수에 추가된다.
시스템 상태 변수:
2차원 평면상에서 로봇이 움직인다 가정하면, 어떤 global 기준 축에 대해, 로봇의 위치와 자세는 (xr, yr, thetar)을 가진다.
또한 상태는 각 LM의 위치 (xi,yi)를 포함한다. 상태변수의 크기는 (3+2n)크기의 벡터이다. LM의 수는 n이다.
분산행렬: P
두 변수의 공분산(covariance)이란 두 변수가 상호 얼마나 강하게 상관성이 있느냐에 대한 척도를 제공한다. 이것은 변수들 사이의 선형의존성(linear dependency)을 계산하는 correlation으로 측정한다.
행렬 P는 여러 부분 행렬로 구성된다. 먼저 로봇 위치 변수 사이의 공분산, LM 사이의 공분산, LM과 로봇위치 사이의 공분산 등이 결합되어 구성된다.
위 그림은 공분산 행렬 P를 보여준다. P의 크기는 (3+2n)x(3+2n)이다.
A: 3x3, 로봇 변수(x,y,theta)에 대한 공분산
B: 2x2, 첫번째 LM에 대한 공분산. LM은 각도 정보 없이 위치만 가진다. 각 LM에 대해 대각 방향으로 마지막 LM의 공분산인 C에 도달할 때까지 반복된다.
D: 2x3, 로봇 상태와 첫번째 LM사이의 공분산을 포함한다.
E: 3x2, 첫번째 LM과 로봇 상태에 대한 공분산으로 D와 전치(transpose) 관계이다.
F: 2x2, 마지막 LM과 첫번째 LM 사이의 공분산이다. G는 F와 전치이다.
초기에는 로봇이 어떤 LM도 관찰하지 않았으므로 P는 행렬 A만 포함한다. 초기 행렬은 대각방향으로 어떤 기본 값을 가지도록 초기화 될 수 있다. 이 값은 초기 위치에 대한 불확실성을 반영한다 (불확실성이 없더라도 초기화가 필요하다. 따라오는 여러 행렬 연산의 singularity를 방지한다)
칼만게인: K
칼만게인은 관찰된 LM을 얼마나 신뢰할 수 있는지, 즉 관찰된 데이터를 얼마나 반영해야 하는지의 신뢰도를 나타낸다.
센서는 정확하지 않고 불확실성을 가지며, 로봇의 odometry도 정확하지 않다. 따라서, 상태 갱신을 위해서 두 데이터를 절충할 수 있다.
만일 센서(레인지 측정장치)가 로봇의 odometry에 비해 부정확하다면 센서 데이터는 신뢰하기 어렵고 칼만게인은 낮아진다. 측정 장치가 우수하고 정확하다면 칼만게인은 높아질 것이다. 칼만게인 행렬은 (3+2n)x2의 크기로 다음과 같다:
행렬의 첫번째 행은 상태변수의 첫번째 행에 기여한다. 즉, 로봇 위치변수 x에 innovation을 얼마나 반영해야 하는지를 나타낸다. 레이저 스캐너 센서의 경우 innovation은 2x1 크기인 (LM 거리, LM 각도)'이므로 게인 행렬의 첫 행의 첫번째 열(r: range)은 LM거리 값을 얼마나 반영해야 할지, 두번째 열(b: bearing)은 LM각도를 얼마나 반영해야 할지를 나타낸다. 둘 다는 상태변수의 첫번째 행에 대한 것이므로 로봇 위치 x에 대한 것이다.
센서모델 및 자코비안: H
먼저 센서 모델은 다음과 같다:
여기서 lambda x는 LM의 x방향의 좌표 위치, x는 로봇 위치의 추정치이다. 이 식은 로봇 위치에서 측정했을 때, LM에 대한 거리(range)와 각도(bearing)에 대한 예측치를 준다. 거리는 로봇과 LM의 직선거리이고, 각도는 기준축(수평축)에 대한 LM의 각도에서 로봇의 자세각을 뺀 값이다. v는 노이즈이다.
이 센서 모델의 x, y, theta에 대한 자코비안은 다음과 같다:
H행렬은 (x,y,theta)가 변할 때, 얼마나 많이 range와 bearing이 변할지를 보여준다.
1행의 첫번째 요소는 x축의 변화에 대한 range의 변화이다. 두번째 요소는 y축에서 변화에 대한 range의 변화에 대한 것이다.
1행의 마지막 요소는 theta의 변화에 대한 것이다. 물론 로봇이 회전할 때 range변화는 없으므로 이것은 0이다.
실제 SLAM에서는 상태 변수에 LM가 포함되어 있으므로, H행렬은 LM에 대한 항들도 포함해야 한다.
한 예로 위 행렬은 우리가 관찰한 센서 정보가 LM 2번일 때의 h및 자코비안 H행렬이다.
HLM2행렬의 첫 3열은 바로 앞에서 나온 H행렬이다.
그리고 LM 하나당 두개의 열이 추가된다. LM 2에 대한 H행렬은 해당 LM 위치인 6과 7열에서만 값을 가지고, 나머지는 모두 0이다. 즉, hLM2가 LM 2에 대한 정보만을 포함하므로 LM 2를 제외한 나머지 LM 상태변수에 대한 h의 도함수는 모두 0이다.
로봇의 이동 모델
제어 입력에 대한 로봇의 이동 모델은 다음과 같다.
여기서 delta sl, sr은 로봇의 왼쪽 바퀴와 오른쪽 바퀴가 회전하여 이동한 양이다. b는 좌우 바퀴 사이의 거리이다. 로봇의 이동 후 위치는 xr(k)=xr(k-1)+f(xr,u)이다.
로봇 상태에 대한 이동 모델의 자코비안은 다음과 같다:
프로세스 노이즈: Q
제어 입력의 노이즈는 좌우 바퀴의 입력에 대한 불확실성이다.
제어 입력에 대한 이동 모델의 자코비안: W
이동 모델의 제어 입력에 대한 1차 도함수 행렬로 구성된다:
센서 측정 오차에 대한 공분산 행렬: R
range와 bearing을 측정하는 레이저스캐너는 오차를 가진다. 상수 c, d는 센서의 정확도를 나타내는 값으로 range값이 1cm의 분산을 가진다면 c=0.01이다. 만일 각도가 1도의 분산을 가진다면 d=1이다.
EKF SLAM의 세 단계
Step 1: odometry 데이터를 사용한 상태의 갱신
실제 SLAM에서 상태 변수는 LM을 포함하므로 상태 추측 모델(이동 모델)은 로봇 변수와 LM을 모두 포함하며 다음과 같다:
공분산 행렬 P의 갱신에 필요한 상태 추측 모델의 상태변수에 대한 자코비안을 구하면 다음과 같다:
행렬 Ak의 크기는 (3+2n)x(3+2n)이다.
제어 입력 (delta sr, delta sl)'에 대한 이동 모델의 자코비안은 다음과 같다:
공분산 행렬 P의 갱신 식은 다음과 같다:
Step 2: 재 관찰된 LM로부터 상태 갱신
로봇이 제어 입력에 의해 이동 후 LM가 하나 관찰되었다면 이것을 이용하여 상태 변수와 상태의 공분산 행렬을 갱신할 수 있다. 즉, 로봇에 설치된 센서에서 LM와의 거리(r: range)와 각도(b: bearing) 값의 벡터 z를 계측할 수 있다.
그런데 이 LM이 이미 한번 관측한 것이라면 로봇의 이동 후 (r, b)의 예측도 가능하다. 예측 값을 벡터 h라 하자. 그러면, z와 h의 차이를 이용해서 상태값과 공분산 행렬을 갱신하는 것이 가능하다. 단, 이 때 칼만 게인을 통해 로봇의 이동(예측) 모델과 센서의 측정 모델의 불확실성을 반영해 준다.
먼저 칼만 게인 식은 다음과 같다:
여기서, 게인 행렬 계산에 사용된 P, H, R 등의 행렬 크기가 표시되었고, S 및 K행렬도 크기를 나타내었다. 게인 계산에 사용되는 행렬들은 모두 앞 부분에 유도 되어 있다. R이 커지면 게인은 작아진다.
일단 칼만 게인이 계산되면 상태와 공분산 행렬의 갱신이 가능하다:
Step 3: 현 상태에 새로운 LM의 추가
만일 관찰한 LM가 기존에 관찰된 LM가 아니라면 이 LM는 새로운 것으로 상태변수에 추가된다. 또한 상태의 크기가 커지면 공분산 행렬도 커진다.
새로운 LM 위치를 y라 하자. y는 LM가 로봇에 의해 관찰된 시점의 값이다. 새로운 상태와 공분산 행렬은 다음과 같다:
Jxr는 y의 로봇 상태에 대한 변화율(자코비안)이고, Jz는 y의 센서 측정값 z에 대한 변화율이다.
센서와 로봇은 각각 불확실성을 나타내는 분산 값 R과 Pxx를 가지므로 이 값들이 LM 위치의 불확실성으로 전파된 값이 공분산 행렬에 추가된다.
.Jz: 센서 측정치 (r, b)의 변화에 대한 LM y의 자코비안
.Jxr: 로봇 상태 (xr,yr,thetar)의 변화에 대한 LM y의 자코비안
.(dy/dxr)Pxx: 불확실성 Pxx는 xr의 것이고 이것이 y에 영향을 줌
.(dy/dz)R: 불확실성 R은 z의 것이고 이것이 y에 영향을 줌
(c.f.)
(df/du)Q: 불확실성 Q는 제어입력 u에 대한 것이고 로봇 이동모델 f에 영향을 줌
JxrPxxJxr'는 로봇 공분산이 LM 공분산으로 전파된 것이고, JzRJz'는 센서의 불확실성이 LM 공분산으로 확산된 것을 보여준다. LM 값에 로봇 위치와 센서 계측이 둘 다 영항을 주므로 두 전파된 공분산의 합이 대각방향에 추가되는 공분산 값이다.
JxrPx(x~yn)은 추가 설명...
참고 자료
[1] SLAM for dummies, MIT Open Course ware.
[2] Open Robotics blog, SLAM part.
[3] Cyrill Stachniss, EKF SLAM 자료
2014년 1월 17일 금요일
기계학습에서 Data cleaning 기술
앞에서는 기계학습에서 발생하는 imbalanced data 문제를 해결하기 위해 여러가지의 샘플링 기반의 데이터 합성 기술을 살펴 보았다.
데이터를 합성하면 overlapping 문제가 생기는데, 이것은 정,부의 데이터들이 서로 뒤 섞여서 학습 성능을 저하시키는 것을 말한다.
따라서 Tomek link와 같은 데이터 삭제 기술을 적용하여 샘플링에 의한 합성 데이터로 생긴 overlapping을 제거할 수 있다.
Tomek link는 정의를 위해 샘플사이의 거리 값을 이용한다. 단, 서로 다른 부류가 연결된 것이며, 상호간의 최소 거리를 가진 두 샘플의 쌍(pair)을 말한다.
예를 들면, (xi, xj)의 쌍에 대해 xi는 Smin(minority of S)에 속하며, xj는 Smaj(majority of S)에 속한다.
또한 d(xi, xj)는 xi xj사이의 거리이며, 만일 d(xi,xk)<d(xi,xj)이거나 d(xj,xk)<d(xi,xj)가 되는 샘플 xk가 없다면 이것을 Temek link라 부른다.
만일 어떤 두 샘플이 Tomek link를 형성한다면, 두 샘플의 어떤 하나는 노이즈이거나 둘 다가 부류의 경계근처에 있는 것으로 볼 수 있다.
따라서, 데이터 합성 후, 원치 않는 부류들 사이의 overlapping을 제거하기 위해 모든 Temek link를 찾아 제거한다. 제거 과정은 모든 최소 거리 연결 쌍이 같은 부류에 있을 때까지 계속된다.
overlapping 샘플들을 제거함에 의해 학습 샘플들은 잘 정의된 부류의 cluster를 얻을 수 있고, 이것은 결국 개선된 성능으로 나타난다.
그림 (a)는 overlapping이 있는 imbalanced된 원 데이터를 보여준다. 그림 (b)는 SMOTE에 의한 데이터 합성 후의 분포를 보여준다. SMOTE에 의해 생긴 overlapping 증가가 나타난다. 그림 (c)에서 Tomek link는 점선 박스로 표시되고, 그림 (d)는 Tomek link를 제거한 결과를 보여준다. 잘 정의된 클래스의 cluster가 나타난다.
참고 논문
[1] Haibo He, et. al., Learning from Imbalanced data, IEEE Trans. Knowledge and data Eng., 21(9), 2009
[2] I, Tomek, Two modifications of CNN, IEEE Trans. System, Man, Cybernetics, 6(11), 1976.
데이터를 합성하면 overlapping 문제가 생기는데, 이것은 정,부의 데이터들이 서로 뒤 섞여서 학습 성능을 저하시키는 것을 말한다.
따라서 Tomek link와 같은 데이터 삭제 기술을 적용하여 샘플링에 의한 합성 데이터로 생긴 overlapping을 제거할 수 있다.
Tomek link는 정의를 위해 샘플사이의 거리 값을 이용한다. 단, 서로 다른 부류가 연결된 것이며, 상호간의 최소 거리를 가진 두 샘플의 쌍(pair)을 말한다.
예를 들면, (xi, xj)의 쌍에 대해 xi는 Smin(minority of S)에 속하며, xj는 Smaj(majority of S)에 속한다.
또한 d(xi, xj)는 xi xj사이의 거리이며, 만일 d(xi,xk)<d(xi,xj)이거나 d(xj,xk)<d(xi,xj)가 되는 샘플 xk가 없다면 이것을 Temek link라 부른다.
만일 어떤 두 샘플이 Tomek link를 형성한다면, 두 샘플의 어떤 하나는 노이즈이거나 둘 다가 부류의 경계근처에 있는 것으로 볼 수 있다.
따라서, 데이터 합성 후, 원치 않는 부류들 사이의 overlapping을 제거하기 위해 모든 Temek link를 찾아 제거한다. 제거 과정은 모든 최소 거리 연결 쌍이 같은 부류에 있을 때까지 계속된다.
overlapping 샘플들을 제거함에 의해 학습 샘플들은 잘 정의된 부류의 cluster를 얻을 수 있고, 이것은 결국 개선된 성능으로 나타난다.
그림 (a)는 overlapping이 있는 imbalanced된 원 데이터를 보여준다. 그림 (b)는 SMOTE에 의한 데이터 합성 후의 분포를 보여준다. SMOTE에 의해 생긴 overlapping 증가가 나타난다. 그림 (c)에서 Tomek link는 점선 박스로 표시되고, 그림 (d)는 Tomek link를 제거한 결과를 보여준다. 잘 정의된 클래스의 cluster가 나타난다.
참고 논문
[1] Haibo He, et. al., Learning from Imbalanced data, IEEE Trans. Knowledge and data Eng., 21(9), 2009
[2] I, Tomek, Two modifications of CNN, IEEE Trans. System, Man, Cybernetics, 6(11), 1976.
2014년 1월 16일 목요일
기계학습에서 Imbalanced learning
교사학습(supervised learning)이란 학습 샘플과 이 샘플이 속하는 부류(class) 정보가 함께 주어질 때 수행하는 기계학습을 말한다.
교사학습에서 imbalanced learning은 정(positive) 학습 샘플과 부(negative) 학습 샘플 수가 비슷하지 않거나 한 쪽의 수가 다른 쪽의 수보다 아주 작을 때 발생한다.
보통, negative 샘플 수는 아주 많고, positive 샘플 수는 상대적으로 작은 경우가 많다.
m개의 샘플을 가진 학습 데이터 집합을 S라 하자.
|S| = m.
이 때, S={(xi,yi)}, i=1,...,m이고 샘플 벡터 xi는 n차원의 특징 공간 X에 속한다.
X={f1,f2,...,fn}
즉, xi는 n개의 특징값을 가진 벡터이다. yi는 Y에 속하는데 Y={1,...,C}이다. 여기서, C는 xi에 대한 class label이다. C=2인 경우는 binary 분류 문제이다.
Imbalanced data에서 S는 두개의 부분 집합으로 나눌 수 있는데 Smin은 샘플 수가 작은(minority) 클래스의 샘플 집합이고, Smaj는 샘플 수가 큰(majority) 클래스의 샘플 집합이다. Smin과 Smaj의 교집합은 공집합이고, 합집합은 전체집합 S이다.
정, 부 샘플 수의 균형을 마추기 위해 샘플 수가 부족한 클래스에 합성 샘플(synthetic example)을 생성하여 샘플 수를 증가시키는 전략을 살펴보자.
SMOTE
SMOTE(synthetic minority oversampling technique)는 여러 응용분야에서 성공적인 결과를 보여준 방법이다.
Smin에 속하는 어떤 샘플 xi에 대해 K-NN(Nearest Neighbors)인 Smin 내의 샘플을 생각하자. K-NN은 기준 샘플 xi에 대해 xi에 가장 가까운 거리를 가진 K개의 샘플을 말한다. K는 어떤 정수 값이고, 메트릭은 특징 공간에서 euclidian거리이다.
합성 샘플을 만들기 위해 K-NN 중에서 한 샘플을 랜덤하게 선택한다. 합성 샘플은 다음과 같이 계산한다.
xnew = xi + (xi_h - xi) * delta
여기서 xi는 Smin에 속하는 기준 샘플이고 xi_h는 xi에 대한 K-NN의 하나이다. xi_h또한 Smin에 속한다. delta는 0과 1사이의 랜덤 수이다.
그림에서 원은 negative, 별은 positive샘플이다. 그림 (a)를 보면 기준 샘플 xi에 대해 가장 가까운 샘플 6개(K-NN)가 선택되었다. 이 중에서 임의로 선택된 샘플이 xi_h(hat of xi)이고, 위 수식에 따라 xi와 xi_h을 연결하는 직선 위에 존재하는 합성 샘플이 하나 추가 된다.
Borderline-SMOTE
SMOTE에서는 과도하게 발생하는 일반화(over generalization) 문제가 있다. 즉, SMOTE는 주변 샘플에 대한 고려 없이 Smin에 속하는 모든 개별 샘플에 대해 동일한 방법으로 합성 데이터를 생성시킨다. 이는 부류들 사이에서 중첩 발생(overlapping occurrence)를 증가시킨다.
Borderline-SMOTE는 이러한 문제점을 해결하기 위해 제안된 적응 샘플링 기법(adaptive sampling method)이다.
먼저, 각각의 Smin에 속하는 xi에 대해 m-NN을 구한다. 단, SMOTE와의 차이는 m-NN을 구할 때, Smin에 속하는 주변 샘플이 아니라 Smin, Smaj를 가리지 않고 가까운 m-NN을 구한다.
다음, m-NN 중에서 Smaj에 속하는 샘플의 수를 확인한다. 이 수를 Nh라고 하자.
최종으로, Nh가 m/2보다 크거나 같고 m보다 작은 조건을 만족하는 xi를 선택한다.
이러한 방법은 위 그림에서 보듯이, Smin보다 Smaj에 속하는 이웃을 더 많이 가진 xi가 선택되게 한다. 위 그림에서 DANGER set에 해당한다.
DANGET set에 속하는 샘플은 SMOTE 알고리즘에 의해 합성 샘플을 생성한다.
만일 m-NN의 모두가 Smaj라면 (즉, 위 그림에서 C에 해당한다), xi는 노이즈로 취급되고, 합성데이터는 발생되지 않는다.
Borderline-SMOTE는 SMOTE에 비해 경계(borderline)에 더 가까운 Smin의 샘플에서 합성 데이터를 생성시키는 효과가 있다.
ADASYN
ADASYN(adaptive synthetic sampling)도 SMOTE의 문제를 해결하기 위해 제안된 적응형 방법의 하나이다.
ADASYN은 주위 데이터의 분포에 따라 발생시킬 합성 데이터의 수를 좀 더 체계적으로 조절하는 방법이다.
먼저, Smin에 대해 발생할 합성 데이터의 수를 계산한다:
G = (|Smaj| - |Smin|) * beta
여기서 beta는 0과 1사이의 값으로 합성 데이터를 생성한 후에, 두 그룹간의 데이터 수의 균형을 정하기 위해 사용된다.
다음, Smin에 속하는 각 샘플 xi에 대해 euclidean거리에 따라 K-NN을 찾고, Gama(i)를 계산한다:
Gama(i) = (delta(i)/K) / Z, i=1,...,|Smin|
여기서 delta(i)는 xi의 K-NN 중에서 Smaj에 속하는 샘플의 수이다. Z는 Gama(i)가 p.d.f(probability density function)가 되도록 하는 정규화 상수(Sigma(Gama(i)) = 1)이다.
다음, Smin에 속하는 각각의 xi에 대해 발생될 필요가 있는 합성 데이터 샘플의 수를 결정한다.
g(i) = Gama(i) * G.
최종적으로, 각각의 xi에 대해 SMOTE에 따라 g(i)개의 합성 데이터를 생성한다.
위 식에서 delta(i)가 크면 많은 데이터를 발생 시킨다. 즉 Smaj에 속하는 샘플의 수가 많은 xi가 많은 데이터를 발생시킨다.
ADYSYN의 주요 아이디어는 합성 데이터의 수를 자동적으로 결정하기 위한 표준으로 p.d.f인 Gama를 사용하는 것이다.
참고 논문
[1] Haibo He, et. al., Learning from Imbalanced data, IEEE Trans. Knowledge and data Eng., 21(9), 2009
[2] N.V. Chawla, et. al., SMOTE: Synthetic Minority Over-Sampling Technique, J. Artificial Intelligence Research, 16, 2002.
[3] H. Han, et.al., Borderline-SMOTE: A New Over-Sampling Method in Imbalanced Data Sets Learning, Proc. ICIC, 2005.
[4] H. He, et. al., ADASYN: Adaptive Synthetic Sampling Approach for Imbalanced Learning, Proc. Int. J. Comf. Neural Networks, 2008.
교사학습에서 imbalanced learning은 정(positive) 학습 샘플과 부(negative) 학습 샘플 수가 비슷하지 않거나 한 쪽의 수가 다른 쪽의 수보다 아주 작을 때 발생한다.
보통, negative 샘플 수는 아주 많고, positive 샘플 수는 상대적으로 작은 경우가 많다.
m개의 샘플을 가진 학습 데이터 집합을 S라 하자.
|S| = m.
이 때, S={(xi,yi)}, i=1,...,m이고 샘플 벡터 xi는 n차원의 특징 공간 X에 속한다.
X={f1,f2,...,fn}
즉, xi는 n개의 특징값을 가진 벡터이다. yi는 Y에 속하는데 Y={1,...,C}이다. 여기서, C는 xi에 대한 class label이다. C=2인 경우는 binary 분류 문제이다.
Imbalanced data에서 S는 두개의 부분 집합으로 나눌 수 있는데 Smin은 샘플 수가 작은(minority) 클래스의 샘플 집합이고, Smaj는 샘플 수가 큰(majority) 클래스의 샘플 집합이다. Smin과 Smaj의 교집합은 공집합이고, 합집합은 전체집합 S이다.
정, 부 샘플 수의 균형을 마추기 위해 샘플 수가 부족한 클래스에 합성 샘플(synthetic example)을 생성하여 샘플 수를 증가시키는 전략을 살펴보자.
SMOTE
SMOTE(synthetic minority oversampling technique)는 여러 응용분야에서 성공적인 결과를 보여준 방법이다.
Smin에 속하는 어떤 샘플 xi에 대해 K-NN(Nearest Neighbors)인 Smin 내의 샘플을 생각하자. K-NN은 기준 샘플 xi에 대해 xi에 가장 가까운 거리를 가진 K개의 샘플을 말한다. K는 어떤 정수 값이고, 메트릭은 특징 공간에서 euclidian거리이다.
합성 샘플을 만들기 위해 K-NN 중에서 한 샘플을 랜덤하게 선택한다. 합성 샘플은 다음과 같이 계산한다.
xnew = xi + (xi_h - xi) * delta
여기서 xi는 Smin에 속하는 기준 샘플이고 xi_h는 xi에 대한 K-NN의 하나이다. xi_h또한 Smin에 속한다. delta는 0과 1사이의 랜덤 수이다.
그림에서 원은 negative, 별은 positive샘플이다. 그림 (a)를 보면 기준 샘플 xi에 대해 가장 가까운 샘플 6개(K-NN)가 선택되었다. 이 중에서 임의로 선택된 샘플이 xi_h(hat of xi)이고, 위 수식에 따라 xi와 xi_h을 연결하는 직선 위에 존재하는 합성 샘플이 하나 추가 된다.
Borderline-SMOTE
SMOTE에서는 과도하게 발생하는 일반화(over generalization) 문제가 있다. 즉, SMOTE는 주변 샘플에 대한 고려 없이 Smin에 속하는 모든 개별 샘플에 대해 동일한 방법으로 합성 데이터를 생성시킨다. 이는 부류들 사이에서 중첩 발생(overlapping occurrence)를 증가시킨다.
Borderline-SMOTE는 이러한 문제점을 해결하기 위해 제안된 적응 샘플링 기법(adaptive sampling method)이다.
먼저, 각각의 Smin에 속하는 xi에 대해 m-NN을 구한다. 단, SMOTE와의 차이는 m-NN을 구할 때, Smin에 속하는 주변 샘플이 아니라 Smin, Smaj를 가리지 않고 가까운 m-NN을 구한다.
다음, m-NN 중에서 Smaj에 속하는 샘플의 수를 확인한다. 이 수를 Nh라고 하자.
최종으로, Nh가 m/2보다 크거나 같고 m보다 작은 조건을 만족하는 xi를 선택한다.
이러한 방법은 위 그림에서 보듯이, Smin보다 Smaj에 속하는 이웃을 더 많이 가진 xi가 선택되게 한다. 위 그림에서 DANGER set에 해당한다.
DANGET set에 속하는 샘플은 SMOTE 알고리즘에 의해 합성 샘플을 생성한다.
만일 m-NN의 모두가 Smaj라면 (즉, 위 그림에서 C에 해당한다), xi는 노이즈로 취급되고, 합성데이터는 발생되지 않는다.
Borderline-SMOTE는 SMOTE에 비해 경계(borderline)에 더 가까운 Smin의 샘플에서 합성 데이터를 생성시키는 효과가 있다.
ADASYN
ADASYN(adaptive synthetic sampling)도 SMOTE의 문제를 해결하기 위해 제안된 적응형 방법의 하나이다.
ADASYN은 주위 데이터의 분포에 따라 발생시킬 합성 데이터의 수를 좀 더 체계적으로 조절하는 방법이다.
먼저, Smin에 대해 발생할 합성 데이터의 수를 계산한다:
G = (|Smaj| - |Smin|) * beta
여기서 beta는 0과 1사이의 값으로 합성 데이터를 생성한 후에, 두 그룹간의 데이터 수의 균형을 정하기 위해 사용된다.
다음, Smin에 속하는 각 샘플 xi에 대해 euclidean거리에 따라 K-NN을 찾고, Gama(i)를 계산한다:
Gama(i) = (delta(i)/K) / Z, i=1,...,|Smin|
여기서 delta(i)는 xi의 K-NN 중에서 Smaj에 속하는 샘플의 수이다. Z는 Gama(i)가 p.d.f(probability density function)가 되도록 하는 정규화 상수(Sigma(Gama(i)) = 1)이다.
다음, Smin에 속하는 각각의 xi에 대해 발생될 필요가 있는 합성 데이터 샘플의 수를 결정한다.
g(i) = Gama(i) * G.
최종적으로, 각각의 xi에 대해 SMOTE에 따라 g(i)개의 합성 데이터를 생성한다.
위 식에서 delta(i)가 크면 많은 데이터를 발생 시킨다. 즉 Smaj에 속하는 샘플의 수가 많은 xi가 많은 데이터를 발생시킨다.
ADYSYN의 주요 아이디어는 합성 데이터의 수를 자동적으로 결정하기 위한 표준으로 p.d.f인 Gama를 사용하는 것이다.
참고 논문
[1] Haibo He, et. al., Learning from Imbalanced data, IEEE Trans. Knowledge and data Eng., 21(9), 2009
[2] N.V. Chawla, et. al., SMOTE: Synthetic Minority Over-Sampling Technique, J. Artificial Intelligence Research, 16, 2002.
[3] H. Han, et.al., Borderline-SMOTE: A New Over-Sampling Method in Imbalanced Data Sets Learning, Proc. ICIC, 2005.
[4] H. He, et. al., ADASYN: Adaptive Synthetic Sampling Approach for Imbalanced Learning, Proc. Int. J. Comf. Neural Networks, 2008.
2014년 1월 7일 화요일
L1 minimization
L1-min 문제는 부족 제한(under-constrained)된 선형 시스템 b=Ax의 해를 구하는 방법 중의 하나이다[1]. 관찰 벡터 b의 요소 수가 미지수 벡터 x의 요수 수보다 작기 때문에 A, b가 주어질 때, x를 구하는 것은 non-trivial linear solver 문제이다.
만일 x가 충분히 sparse하다면(즉, 표준 축(canonical coordinate)에서 요소의 대부분이 0이라면) x값은 minimum l1-norm을 계산함에 의해 얻어진다(주1):
min ||x||1, subject to b=Ax
x
실제로는 observation b는 noise를 포함하므로 equality condition은 다소 완화되고 constrained BPDN (basis pursuit denoising)문제가 된다:
min ||x||1, subject to ||b-Ax||2 <= epsilon
x
여기서, epsilon>0 은 미리 정해진 noise의 레벨을 나타낸다. 이 식은 가중치 ramda를 도입하면 unconstrained BPDN문제(주2)로 바뀐다:
min (1/2)||b-Ax||^2 + ramda.||x||1.
x
위 식들에서 ||.||1은 l1-norm, ||.||2는 l2-norm이다. l2 norm은 Euclidean거리이고, l1 norm은 벡터 x의 요소들의 절대치의 합이다. 목표 x는 가능하면 sparse하면서 b=Ax를 만족하는 것이 정답이된다.
An example of L1-min solver
unconstrained BPDN을 풀어내는 Shrinkage 알고리즘을 살펴 보자. 가림이 존재하는 물체를 L1 tracker를 사용하여 추적하는 예에서 Objective function은 다음과 같다:
min L(z,e), where L(z,e)=(1/2)||y-Uz-e||^2 + ramda||e||1
z,e
여기서, y는 관찰(observation) 벡터, U는 PCA 기저(basis), z는 기저에 대한 계수, e는 trivial template에 대한 계수 벡터이다. 이 식에 대한 closed-form 풀이법은 없다. 다음의 2단계로 이것을 푼다[2].
step 1: Given eopt, zopt는 zopt = U'(y-eopt)로 얻어진다.
step 2: Given zopt, eopt = Sr (y-Uzopt).
여기서 Sr(.)는 shrinkage operator이다. 식의 유도와 증명은 [2]를 참조한다. 알고리즘은 다음과 같다.
아래의 Matlab code를 호출하기 위해서는 W에 PCA기저, srParam에는 미리 정해진 상수 값들을 넘겨 주어야 한다. Wang[2]의 논문에서는 32x32크기의 window patch에 대해 16개의 PCA 기저를 사용하므로 W는 1024x16크기의 행렬이다.
srParam값은 다음과 같다:
srParam.(lambda(0.05), L0(0.05), maxLoopNum(20), tol(10e-3))
data에는 물체를 포함하는 32x32크기의 window patch에서 평균값을 제거한 1024x1 크기의 벡터 값이 넘어가야 한다. Wang의 코드에서는 이전 프레임까지 추적된 affine계수값의 근방에서 정규분포로 샘플을 600개 발생시킨다. 다음 프레임의 영상에 이렇게 발생된 600개의 affine 계수를 차례로 적용시켜 만들어낸 600개의 이미지에서 평균을 제거하고 남은 값을 차례로 data값으로 넘겨준다.
600번의 pca_L1함수 호출의 결과로 얻어낸 600개의 z,e값을 가지고 각 particle의 confidence를 평가한다.
function alpha = pca_L1(data, W, srParam)
%% function alpha = pca_L1(data, W, srParam)
%% Object representation via sparse prototypes
%% (PCA basis vectors +trivial templates)
%% Input:
%% data: A normalized data vector
%% W: PCA basis vectors
%% srParam: Parameters for L1 minimization
%% -srParam.lambda
%% -srParam.maxLoopNum
%% -srParam.tol
%% Output:
%% alpha: The representation coefficient
%%
%%1.Initialization:
coeff = zeros(size(W,2),1); % The coefficients of PCA basis vectors
% size(coeff) = 16x1, size(W)=1024x16 (16 basis for 32x32 window patch)
err = zeros(size(data)); % The coefficients of trivial templates
% size(err) = 1024x1 for 32x32 window patch
objValue = zeros(1,srParam.maxLoopNum); % The values of objective functions
% size(objValue) = 1x20 for max iteration # 20
%%2.Iterative solution:
for num = 1:srParam.maxLoopNum % opt.srParam=20
%(2.1)Fix 'err', slove 'coeff' --(Lemma 1 in Section III; Step 3 in Table)
coeff = W'*(data-err); % U'(y-e)
%(2.1)Fix 'coeff', slove 'err'--(Lemma 2 in Section III; Step 4 in Table)
y = data-W*coeff; % Sr(y-Uz)
err = max(abs(y)-srParam.lambda, 0).*sign(y); % size(abs(y))=1024x1
% abs(y) 값 중에서 lambda(0.05)보다 작은 잔차를 가진 경우(즉, 가림이나
% background clutter가 아니라고 판단)는 err 값을 0으로 만듬.
% 따라서, err가 살아 있는 부분은 가림이나 clutter에 의한 오차 부분이다.
if num == 1
alpha = [coeff; err]; % Save the result
continue;
end
objValue(num) = (y-err)'*(y-err) + srParam.lambda*sum(abs(err));
% objValue는 가림은 최소화하면서 data를 잘 근사하는 coeff와 err
% 값을 구하는 것이 목표이다.
% 식의 앞 부분에서 y=data-W*coeff는 후보 윈도 패치의 잔차이고, 이 중에서
% 가림이 존재하는 부분은 err를 빼줌에 의해 제거되어 보상된다.
% 뒤의 항은 err벡터 내부에 값을 가진 항을 최소화해서 err가 패치 내에서
% 작용하는 픽셀의 수를 줄이도록 해준다.
if abs(objValue(num)-objValue(num-1)) < srParam.tol
break;
else
alpha = [coeff; err]; % Save the result
end
end
% 첫번째 호출에서 objVaue=1.4406 1.4247 1.4181 1.4159 1.4151이고
% 5번 반복 후에 수렴하였음. objective func. 값은 점차 작아져서 수렴함.
% z, e <-- alpha for input 'data'
위 함수를 600번 호출하여 만든 600개의 z,e 집합을 confidence값으로 바꾸어 주어야 한다. 먼저 target template를 표현하는 16개의 target coefficient coeff와 trivial template를 표현하는 1024개의 trivial coefficient err를 얻어 낸다. 또, 원 데이터인 data에서 z,e로 재 생성한 데이터를 빼서 잔차 diff를 만든다.
600개 중에서 하나만 평가해 보자. err 내의 600개 벡터 중의 하나를 err_i라 하자. 먼저, err_i값을 어떤 임계값(threshold)을 중심으로 분리한다.
confidence가 높아지기 위해서는 err_i 요소 중에서 th보다 큰(즉, occluded되었다고 가정된) 요소의 합이 작아져야 한다. 즉, sparse해야 한다는 것이다. 또한, err_i 요소 중 th보다 작은 요소의 값은 occlusion이나 background clutter가 없는 픽셀이므로 이 경우는 잔차 diff가 가능하면 작아야 한다.
아래의 코드는 이러한 아이디어를 반영하고 있다.
coeff = alpha(1:size(tmpl.basis,2),:); %The coefficients of PCA basis vectors
% size(coeff)=16x600, tmpl.basis = 16
err = alpha(size(tmpl.basis,2)+1:end,:); %The coefficients of trivial templates
% size(err) = 1024x600
%% Calucate observation likelihood
% diff는 600개의 각각의 후보 패치에서 평균 템플릿을 제거한 영상 패치.
% 하나의 후보 패치는 32x32=1024x1이므로 전체는 1024x600이다.
diff = diff-tmpl.basis*coeff-err; %Reconstruction, diff=1024x600
diff = diff.*(abs(err)<opt.srParam.L0); % abs(err)<opt.srParam.L0 = 1024x600
param.conf = exp(-(sum(diff.^2)+opt.srParam.L0
*sum(abs(err)>=opt.srParam.L0))./opt.condenssig)'; % 0.05
% sum(diff.^2)=1x600, sum은 각 열의 합.
% 두 항이 다 작아야 되는데, 앞의 항은 err가 작은 항들의 합이므로 가림이 없는 부분
% 뒤 항은 err가 큰 항들의 합이니 가림 부분
% 즉, 비 가림부분은 잔차가 작아야 하고, 가림 부분은 err의 합이 작아야(sparse) 해야 함.
(주1) 원래 Ax=b의 선형시스템 방정식에서 관찰 값 b내의 요소 수가 미지수 x 내의 요소 수보다 작다면 이 식은 부족 제한(under-constrained)이며, 해의 수가 무한개가 되므로 풀 수 없게 된다.
단, x가 sparse하다면 이 식은 BP(basis pursuit)문제가 되고 제한적으로 풀 수 있다.
Ax = b
에서, 0이 아닌 x의 요소는 A(보통 측정행렬이라 부름)내 어떤 열(atom)에 대응한다. x를 결정하는 것이란 신호(measurement) b를 재생성하기 위해 필요한 A행렬 내의 열의 요소를 선택하는 것이다.
이러한 x 내의 0이 아닌 요소의 집합을 basis라고 부른다. 따라서 x 내의 0이 아닌 요소를 찾아 내는 것을 BP(basis pursuit) 라고 부른다.
L1-놈의 사용은 A행렬의 각 atom에 어떤 비용을 할당하는 것이다. 이 비용의 나열이 x벡터이다. 단, 이 비용의 합은 최소화되도록 한다.
(주2) BPDN의 해를 찾는 문제는 convex quadratic problem이다. quadratic은 2차 방정식을 의미하고, convex의 의미는 어떤 제한된 구간 내에서 임의의 직선을 그었을 때, 2차 식과 만나는 두 점 아래의 그래프의 어떤 점도 이 직선 위에 있지 않을 때 convex라 부른다[7]. convex한 형태의 방정식은 gradient-descendant법으로 최적해를 찾았을 때, 항상 전역 최적값을 보장하기 때문에 최적화에 유리하다.
BPDN 문제는 2차 미분 값이 항상 양수이므로 convex이고 전역 최적값을 보장한다[3, 5].
참고 논문
[1] A. Y. Yang, et. al., Fast l1-minimization algorithms for robust face recognition, IEEE Trans. IP, 2013
[2] D.Wang, et.al., Online object tracking with sparse prototypes, IEEE Trans. IP, 2011.
[3] Basis pursuit, Wikipedia
[4] BPDN, Wikipedia
[5] Convexity of quadratic functions in two variables, Webpage
[6] Kristen Michelle Cheman, OPTIMIZATION TECHNIQUES FOR SOLVING BASIS PURSUIT PROBLEMS, Ph.D. thesis
[7] Convex function, Wikipedia
만일 x가 충분히 sparse하다면(즉, 표준 축(canonical coordinate)에서 요소의 대부분이 0이라면) x값은 minimum l1-norm을 계산함에 의해 얻어진다(주1):
min ||x||1, subject to b=Ax
x
실제로는 observation b는 noise를 포함하므로 equality condition은 다소 완화되고 constrained BPDN (basis pursuit denoising)문제가 된다:
min ||x||1, subject to ||b-Ax||2 <= epsilon
x
여기서, epsilon>0 은 미리 정해진 noise의 레벨을 나타낸다. 이 식은 가중치 ramda를 도입하면 unconstrained BPDN문제(주2)로 바뀐다:
min (1/2)||b-Ax||^2 + ramda.||x||1.
x
위 식들에서 ||.||1은 l1-norm, ||.||2는 l2-norm이다. l2 norm은 Euclidean거리이고, l1 norm은 벡터 x의 요소들의 절대치의 합이다. 목표 x는 가능하면 sparse하면서 b=Ax를 만족하는 것이 정답이된다.
An example of L1-min solver
unconstrained BPDN을 풀어내는 Shrinkage 알고리즘을 살펴 보자. 가림이 존재하는 물체를 L1 tracker를 사용하여 추적하는 예에서 Objective function은 다음과 같다:
min L(z,e), where L(z,e)=(1/2)||y-Uz-e||^2 + ramda||e||1
z,e
여기서, y는 관찰(observation) 벡터, U는 PCA 기저(basis), z는 기저에 대한 계수, e는 trivial template에 대한 계수 벡터이다. 이 식에 대한 closed-form 풀이법은 없다. 다음의 2단계로 이것을 푼다[2].
step 1: Given eopt, zopt는 zopt = U'(y-eopt)로 얻어진다.
step 2: Given zopt, eopt = Sr (y-Uzopt).
여기서 Sr(.)는 shrinkage operator이다. 식의 유도와 증명은 [2]를 참조한다. 알고리즘은 다음과 같다.
아래의 Matlab code를 호출하기 위해서는 W에 PCA기저, srParam에는 미리 정해진 상수 값들을 넘겨 주어야 한다. Wang[2]의 논문에서는 32x32크기의 window patch에 대해 16개의 PCA 기저를 사용하므로 W는 1024x16크기의 행렬이다.
srParam값은 다음과 같다:
srParam.(lambda(0.05), L0(0.05), maxLoopNum(20), tol(10e-3))
data에는 물체를 포함하는 32x32크기의 window patch에서 평균값을 제거한 1024x1 크기의 벡터 값이 넘어가야 한다. Wang의 코드에서는 이전 프레임까지 추적된 affine계수값의 근방에서 정규분포로 샘플을 600개 발생시킨다. 다음 프레임의 영상에 이렇게 발생된 600개의 affine 계수를 차례로 적용시켜 만들어낸 600개의 이미지에서 평균을 제거하고 남은 값을 차례로 data값으로 넘겨준다.
600번의 pca_L1함수 호출의 결과로 얻어낸 600개의 z,e값을 가지고 각 particle의 confidence를 평가한다.
function alpha = pca_L1(data, W, srParam)
%% function alpha = pca_L1(data, W, srParam)
%% Object representation via sparse prototypes
%% (PCA basis vectors +trivial templates)
%% Input:
%% data: A normalized data vector
%% W: PCA basis vectors
%% srParam: Parameters for L1 minimization
%% -srParam.lambda
%% -srParam.maxLoopNum
%% -srParam.tol
%% Output:
%% alpha: The representation coefficient
%%
%%1.Initialization:
coeff = zeros(size(W,2),1); % The coefficients of PCA basis vectors
% size(coeff) = 16x1, size(W)=1024x16 (16 basis for 32x32 window patch)
err = zeros(size(data)); % The coefficients of trivial templates
% size(err) = 1024x1 for 32x32 window patch
objValue = zeros(1,srParam.maxLoopNum); % The values of objective functions
% size(objValue) = 1x20 for max iteration # 20
%%2.Iterative solution:
for num = 1:srParam.maxLoopNum % opt.srParam=20
%(2.1)Fix 'err', slove 'coeff' --(Lemma 1 in Section III; Step 3 in Table)
coeff = W'*(data-err); % U'(y-e)
%(2.1)Fix 'coeff', slove 'err'--(Lemma 2 in Section III; Step 4 in Table)
y = data-W*coeff; % Sr(y-Uz)
err = max(abs(y)-srParam.lambda, 0).*sign(y); % size(abs(y))=1024x1
% abs(y) 값 중에서 lambda(0.05)보다 작은 잔차를 가진 경우(즉, 가림이나
% background clutter가 아니라고 판단)는 err 값을 0으로 만듬.
% 따라서, err가 살아 있는 부분은 가림이나 clutter에 의한 오차 부분이다.
if num == 1
alpha = [coeff; err]; % Save the result
continue;
end
objValue(num) = (y-err)'*(y-err) + srParam.lambda*sum(abs(err));
% objValue는 가림은 최소화하면서 data를 잘 근사하는 coeff와 err
% 값을 구하는 것이 목표이다.
% 식의 앞 부분에서 y=data-W*coeff는 후보 윈도 패치의 잔차이고, 이 중에서
% 가림이 존재하는 부분은 err를 빼줌에 의해 제거되어 보상된다.
% 뒤의 항은 err벡터 내부에 값을 가진 항을 최소화해서 err가 패치 내에서
% 작용하는 픽셀의 수를 줄이도록 해준다.
if abs(objValue(num)-objValue(num-1)) < srParam.tol
break;
else
alpha = [coeff; err]; % Save the result
end
end
% 첫번째 호출에서 objVaue=1.4406 1.4247 1.4181 1.4159 1.4151이고
% 5번 반복 후에 수렴하였음. objective func. 값은 점차 작아져서 수렴함.
% z, e <-- alpha for input 'data'
위 함수를 600번 호출하여 만든 600개의 z,e 집합을 confidence값으로 바꾸어 주어야 한다. 먼저 target template를 표현하는 16개의 target coefficient coeff와 trivial template를 표현하는 1024개의 trivial coefficient err를 얻어 낸다. 또, 원 데이터인 data에서 z,e로 재 생성한 데이터를 빼서 잔차 diff를 만든다.
600개 중에서 하나만 평가해 보자. err 내의 600개 벡터 중의 하나를 err_i라 하자. 먼저, err_i값을 어떤 임계값(threshold)을 중심으로 분리한다.
confidence가 높아지기 위해서는 err_i 요소 중에서 th보다 큰(즉, occluded되었다고 가정된) 요소의 합이 작아져야 한다. 즉, sparse해야 한다는 것이다. 또한, err_i 요소 중 th보다 작은 요소의 값은 occlusion이나 background clutter가 없는 픽셀이므로 이 경우는 잔차 diff가 가능하면 작아야 한다.
아래의 코드는 이러한 아이디어를 반영하고 있다.
coeff = alpha(1:size(tmpl.basis,2),:); %The coefficients of PCA basis vectors
% size(coeff)=16x600, tmpl.basis = 16
err = alpha(size(tmpl.basis,2)+1:end,:); %The coefficients of trivial templates
% size(err) = 1024x600
%% Calucate observation likelihood
% diff는 600개의 각각의 후보 패치에서 평균 템플릿을 제거한 영상 패치.
% 하나의 후보 패치는 32x32=1024x1이므로 전체는 1024x600이다.
diff = diff-tmpl.basis*coeff-err; %Reconstruction, diff=1024x600
diff = diff.*(abs(err)<opt.srParam.L0); % abs(err)<opt.srParam.L0 = 1024x600
param.conf = exp(-(sum(diff.^2)+opt.srParam.L0
*sum(abs(err)>=opt.srParam.L0))./opt.condenssig)'; % 0.05
% sum(diff.^2)=1x600, sum은 각 열의 합.
% 두 항이 다 작아야 되는데, 앞의 항은 err가 작은 항들의 합이므로 가림이 없는 부분
% 뒤 항은 err가 큰 항들의 합이니 가림 부분
% 즉, 비 가림부분은 잔차가 작아야 하고, 가림 부분은 err의 합이 작아야(sparse) 해야 함.
단, x가 sparse하다면 이 식은 BP(basis pursuit)문제가 되고 제한적으로 풀 수 있다.
Ax = b
에서, 0이 아닌 x의 요소는 A(보통 측정행렬이라 부름)내 어떤 열(atom)에 대응한다. x를 결정하는 것이란 신호(measurement) b를 재생성하기 위해 필요한 A행렬 내의 열의 요소를 선택하는 것이다.
이러한 x 내의 0이 아닌 요소의 집합을 basis라고 부른다. 따라서 x 내의 0이 아닌 요소를 찾아 내는 것을 BP(basis pursuit) 라고 부른다.
L1-놈의 사용은 A행렬의 각 atom에 어떤 비용을 할당하는 것이다. 이 비용의 나열이 x벡터이다. 단, 이 비용의 합은 최소화되도록 한다.
(주2) BPDN의 해를 찾는 문제는 convex quadratic problem이다. quadratic은 2차 방정식을 의미하고, convex의 의미는 어떤 제한된 구간 내에서 임의의 직선을 그었을 때, 2차 식과 만나는 두 점 아래의 그래프의 어떤 점도 이 직선 위에 있지 않을 때 convex라 부른다[7]. convex한 형태의 방정식은 gradient-descendant법으로 최적해를 찾았을 때, 항상 전역 최적값을 보장하기 때문에 최적화에 유리하다.
BPDN 문제는 2차 미분 값이 항상 양수이므로 convex이고 전역 최적값을 보장한다[3, 5].
참고 논문
[1] A. Y. Yang, et. al., Fast l1-minimization algorithms for robust face recognition, IEEE Trans. IP, 2013
[2] D.Wang, et.al., Online object tracking with sparse prototypes, IEEE Trans. IP, 2011.
[3] Basis pursuit, Wikipedia
[4] BPDN, Wikipedia
[5] Convexity of quadratic functions in two variables, Webpage
[6] Kristen Michelle Cheman, OPTIMIZATION TECHNIQUES FOR SOLVING BASIS PURSUIT PROBLEMS, Ph.D. thesis
[7] Convex function, Wikipedia
2014년 1월 5일 일요일
Symbolic에 의한 Jacobian 연산
물체 추적이나 필터 설계에서 Jacobian을 계산해야 할 경우가 많다. 자세 행렬의 Jacobian을 구할 경우 회전 행렬을 각도 변수에 대해 편미분해야 하며 손으로 계산하기에는 복잡하다. Matlab의 symbolic연산을 이용하여 Qurternion을 이용하는 회전 행렬의 편미분을 계산해 보자.
% Symbolic 연산 테스트
% Qurternion의 Jacobian 계산
% 01/2014, funMV
clear; clc;
syms alpa; % 변수 alpa를 symbol로 선언
syms beta;
syms gama;
cos_alpa = cos(alpa/2); sin_alpa=sin(alpa/2);
cos_beta = cos(beta/2); sin_beta=sin(beta/2);
cos_gama = cos(gama/2); sin_gama=sin(gama/2);
q1=cos_alpa*cos_beta*sin_gama-sin_alpa*sin_beta*cos_gama;
q2=cos_alpa*sin_beta*cos_gama+sin_alpa*cos_beta*sin_gama;
q3=sin_alpa*cos_beta*cos_gama-cos_alpa*sin_beta*sin_gama;
q4=cos_alpa*cos_beta*cos_gama+sin_alpa*sin_beta*sin_gama;
%{
q1 =
cos(alpa/2)*cos(beta/2)*sin(gama/2) - cos(gama/2)*sin(alpa/2)*sin(beta/2)
%}
jacobian(q1,alpa)
% -(cos(alpa/2)*cos(gama/2)*sin(beta/2))/2-(cos(beta/2)*sin(alpa/2)*sin(gama/2))/2
jacobian(q1,beta)
% -(cos(beta/2)*cos(gama/2)*sin(alpa/2))/2-(cos(alpa/2)*sin(beta/2)*sin(gama/2))/2
jacobian(q1,gama)
% (sin(alpa/2)*sin(beta/2)*sin(gama/2))/2+(cos(alpa/2)*cos(beta/2)*cos(gama/2))/2
jacobian(q2,alpa)
% (cos(alpa/2)*cos(beta/2)*sin(gama/2))/2 - (cos(gama/2)*sin(alpa/2)*sin(beta/2))/2
jacobian(q2,beta)
% (cos(alpa/2)*cos(beta/2)*cos(gama/2))/2 - (sin(alpa/2)*sin(beta/2)*sin(gama/2))/2
jacobian(q2,gama)
% (cos(beta/2)*cos(gama/2)*sin(alpa/2))/2 - (cos(alpa/2)*sin(beta/2)*sin(gama/2))/2
% product of two quaternion b and q
syms b1;
syms b2;
syms b3;
syms b4;
qq1 = b1*q1-[b2 b3 b4]*[q2; q3; q4];
qq2 = (b1*[q2 q3 q4]+q1*[b2 b3 b4])+(cross([b2 b3 b4]',[q2 q3 q4]'))';
%{
qq1 =
b1*(cos(alpa/2)*cos(beta/2)*sin(gama/2) - cos(gama/2)*sin(alpa/2)*sin(beta/2))
- b2*(cos(alpa/2)*cos(gama/2)*sin(beta/2) + cos(beta/2)*sin(alpa/2)*sin(gama/2))
- b3*(cos(beta/2)*cos(gama/2)*sin(alpa/2) - cos(alpa/2)*sin(beta/2)*sin(gama/2))
- b4*(sin(alpa/2)*sin(beta/2)*sin(gama/2) + cos(alpa/2)*cos(beta/2)*cos(gama/2))
qq2(1)= b1*(cos(alpa/2)*cos(gama/2)*sin(beta/2) + cos(beta/2)*sin(alpa/2)*sin(gama/2))
+ b2*(cos(alpa/2)*cos(beta/2)*sin(gama/2) - cos(gama/2)*sin(alpa/2)*sin(beta/2))
+ b3*(sin(alpa/2)*sin(beta/2)*sin(gama/2) + cos(alpa/2)*cos(beta/2)*cos(gama/2))
- b4*(cos(beta/2)*cos(gama/2)*sin(alpa/2) - cos(alpa/2)*sin(beta/2)*sin(gama/2)),
qq2(2)= b1*(cos(beta/2)*cos(gama/2)*sin(alpa/2) - cos(alpa/2)*sin(beta/2)*sin(gama/2))
- b2*(sin(alpa/2)*sin(beta/2)*sin(gama/2) + cos(alpa/2)*cos(beta/2)*cos(gama/2))
+ b3*(cos(alpa/2)*cos(beta/2)*sin(gama/2) - cos(gama/2)*sin(alpa/2)*sin(beta/2))
+ b4*(cos(alpa/2)*cos(gama/2)*sin(beta/2) + cos(beta/2)*sin(alpa/2)*sin(gama/2)),
qq2(3)= b1*(sin(alpa/2)*sin(beta/2)*sin(gama/2) + cos(alpa/2)*cos(beta/2)*cos(gama/2))
+ b2*(cos(beta/2)*cos(gama/2)*sin(alpa/2) - cos(alpa/2)*sin(beta/2)*sin(gama/2))
- b3*(cos(alpa/2)*cos(gama/2)*sin(beta/2) + cos(beta/2)*sin(alpa/2)*sin(gama/2))
+ b4*(cos(alpa/2)*cos(beta/2)*sin(gama/2) - cos(gama/2)*sin(alpa/2)*sin(beta/2))
%}
jacobian(qq1, b1)
%cos(alpa/2)*cos(beta/2)*sin(gama/2) - cos(gama/2)*sin(alpa/2)*sin(beta/2)
jacobian(qq1, b2)
%- cos(alpa/2)*cos(gama/2)*sin(beta/2) - cos(beta/2)*sin(alpa/2)*sin(gama/2)
jacobian(qq1, b3)
%cos(alpa/2)*sin(beta/2)*sin(gama/2) - cos(beta/2)*cos(gama/2)*sin(alpa/2)
jacobian(qq1, b4)
%- sin(alpa/2)*sin(beta/2)*sin(gama/2) -%cos(alpa/2)*cos(beta/2)*cos(gama/2)
jacobian(qq1, alpa)
%{
b4*((cos(beta/2)*cos(gama/2)*sin(alpa/2))/2 - (cos(alpa/2)*sin(beta/2)*sin(gama/2))/2)
- b2*((cos(alpa/2)*cos(beta/2)*sin(gama/2))/2 - (cos(gama/2)*sin(alpa/2)*sin(beta/2))/2)
- b3*((sin(alpa/2)*sin(beta/2)*sin(gama/2))/2 + (cos(alpa/2)*cos(beta/2)*cos(gama/2))/2)
- b1*((cos(alpa/2)*cos(gama/2)*sin(beta/2))/2 + (cos(beta/2)*sin(alpa/2)*sin(gama/2))/2)
%}
jacobian(qq2, b1)
% cos(alpa/2)*cos(gama/2)*sin(beta/2) + cos(beta/2)*sin(alpa/2)*sin(gama/2)
% cos(beta/2)*cos(gama/2)*sin(alpa/2) - cos(alpa/2)*sin(beta/2)*sin(gama/2)
% sin(alpa/2)*sin(beta/2)*sin(gama/2) + cos(alpa/2)*cos(beta/2)*cos(gama/2)
jacobian(qq2, alpa)
%{
b1*((cos(alpa/2)*cos(beta/2)*sin(gama/2))/2 - (cos(gama/2)*sin(alpa/2)*sin(beta/2))/2)
- b2*((cos(alpa/2)*cos(gama/2)*sin(beta/2))/2 + (cos(beta/2)*sin(alpa/2)*sin(gama/2))/2)
- b3*((cos(beta/2)*cos(gama/2)*sin(alpa/2))/2 - (cos(alpa/2)*sin(beta/2)*sin(gama/2))/2)
- b4*((sin(alpa/2)*sin(beta/2)*sin(gama/2))/2 + (cos(alpa/2)*cos(beta/2)*cos(gama/2))/2),
b1*((sin(alpa/2)*sin(beta/2)*sin(gama/2))/2 + (cos(alpa/2)*cos(beta/2)*cos(gama/2))/2)
+ b2*((cos(beta/2)*cos(gama/2)*sin(alpa/2))/2 - (cos(alpa/2)*sin(beta/2)*sin(gama/2))/2)
- b3*((cos(alpa/2)*cos(gama/2)*sin(beta/2))/2 + (cos(beta/2)*sin(alpa/2)*sin(gama/2))/2)
+ b4*((cos(alpa/2)*cos(beta/2)*sin(gama/2))/2 - (cos(gama/2)*sin(alpa/2)*sin(beta/2))/2),
b2*((sin(alpa/2)*sin(beta/2)*sin(gama/2))/2 + (cos(alpa/2)*cos(beta/2)*cos(gama/2))/2)
- b1*((cos(beta/2)*cos(gama/2)*sin(alpa/2))/2 - (cos(alpa/2)*sin(beta/2)*sin(gama/2))/2)
- b3*((cos(alpa/2)*cos(beta/2)*sin(gama/2))/2 - (cos(gama/2)*sin(alpa/2)*sin(beta/2))/2)
- b4*((cos(alpa/2)*cos(gama/2)*sin(beta/2))/2 + (cos(beta/2)*sin(alpa/2)*sin(gama/2))/2)
%}
Object tracking with sparse representation (L1 tracker)
Sparse 표현은 최근 cv 분야에서 폭 넓게 인용되고 연구되는 주제이다. 물체 추적 문제에서 가림(occlusion) 문제의 해결은 여전히 어려운데, sparse표현이 이를 다루고 있다.
Mei[1]는 L1 tracker를 제안하여 물체 추적 시 부분 가림을 trivial template를 이용하여 해결 하였다. 목표 물체를 표현하는 방법을 몇 가지 살펴 보자.
(1) PCA에 의한 관찰 물체 표현
그림은 목표 물체의 영상을 표현하기 위해 pca를 적용한 경우를 보여준다. 주축(pca basis) 몇 개와 주축에서의 좌표 값(target coefficient)이 곱해져서 목표 영상을 생성한다:
y = Uz +e
여기서, y는 관찰 벡터, U는 기저이고, z는 좌표 값, e는 오차이다. y를 안다면, z는 z=U'y로 계산 가능하다. 잔차(residual)은 ||y-UU'y||^2이다.
(2) L1 tracker
Mei는 그림에서 보인 바와 같이 sparse표현을 제안 했다. 물체의 가림을 다루기 위해 trivial template를 도입하였다. 물체에서 가림 없이 보이는 부분은 target template가 나타내고, 가려지는 부분은 trivial template에 곱해지는 trivial 계수가 표현한다. 예를 들면, 얼굴의 크기가 32x32라면 trivial 계수는 1024의 요소를 가진 벡터이다. target 계수는 여러 템플레이트 중에서 목표 물체와 가장 유사한 템플릿의 인덱스를 나타낸다.
y=Az+e = [A a][z e]' = Bc
여기서, y는 관찰 벡터, A는 템플릿들의 행렬, z는 템플릿을 선택하는 계수, e는 trivial template에 대한 계수 벡터이다. 행렬 a는 32x32 크기의 trivial 템플릿(즉, 1024크기의 벡터)가 1024개 있는 행렬로(즉, 1024x1024이다) 하나의 벡터에는 해당 픽셀에 대응되는 위치에 1이 저장되어 있고 나머지가 모두 0인 벡터이다. trivial 템플릿의 계수벡터 e는 1024 벡터이다. 계수 벡터의 값이 occlusion을 표현하게 된다.
수식의 풀이는 L1 minimization을 통한다.
여기서, ||.||2는 l2놈이고, ||.||1은 l1놈이다. l1놈은 모든 계수의 절대치의 합의 값으로 생각하면 된다. 식의 의미를 살펴보면, Bc는 가능하면 관찰 벡터 y와 일치해야 하고, c벡터는 sparse(l1놈이 작아지려면 벡터 내에 값이 듬성듬성 들어 있어야 한다)해야 한다는 의미이다. 즉, occlusion은 최소화하면서 y와 잘 일치하도록 계수치들을 결정하면 된다. 단점은 연산량이 크다.
(3) PCA 기저를 이용한 L1 tracker
Wang[2]은 SKL과 L1 tracker를 결합하였다. 그림 (c)는 위 그림 (b)의 target 템플릿을 pca 기저(basis)로 바꾼 것이다. 여러 템플릿 영상을 pca의 기저로 표현하면 템플릿의 누적, 갱신 등이 효과적이다. 즉, SKL(sequential KL)을 적용할 수 있다.
y = Uz+e = [U a][z e]'
여기서, y는 관찰 벡터, U는 pca기저, z는 기저에 대한 계수, e는 trivial 템플릿에 대한 계수벡터이다. 이 식을 다시 쓰면
이다. 위 그림에서 (b)와 (c)를 비교해 보면, (b)의 경우는 물체나 trivial 템플릿 계수 벡터나 둘 다 sparse하다. 즉, target template들은 모두 하나의 물체를 나타내는 영상이므로, 이들 사이에는 상관성이 존재한다. (c)의 경우는 pca기저 벡터가 서로 독립(orthonormal)하기 때문에 기저 계수 벡터는 sparse하지 않고, trivial 계수만 sparse하다. 일반적으로 pca기저는 몇개 안되며(wang논문에서는 16개 사용), trivial template수는 아주 많다(wang에서는 1024개).
식에서 y-Uz는 y가 Uz(pca 기저에 의해 표현되는 템플릿 패치)에 가까와야 한다는 것을 표현하고, y-Uz의 값이 큰 부분(즉, 가림이나 background clutter)은 e에 의해 제거(보상)된다. 단, e가 패치 내의 모든 픽셀에 대해 작용하지 않도록 y-Uz값을 분석해서 잔차가 큰 부분에만 선택적으로 e가 작용하게 한다.
식에서 y-Uz는 y가 Uz(pca 기저에 의해 표현되는 템플릿 패치)에 가까와야 한다는 것을 표현하고, y-Uz의 값이 큰 부분(즉, 가림이나 background clutter)은 e에 의해 제거(보상)된다. 단, e가 패치 내의 모든 픽셀에 대해 작용하지 않도록 y-Uz값을 분석해서 잔차가 큰 부분에만 선택적으로 e가 작용하게 한다.
참고 논문
[1] X.Mei, H. Ling, Robust visual tracking using L1 minimization, CVPR 2009.
[2] D.Wang, et.al., Online object tracking with sparse prototypes, IEEE Trans. IP, 2011.
피드 구독하기:
글 (Atom)