- 업데이트: 내부적으로 GLPK 를 이용할 수 있는 옵션이 있어서 써보니까 훨씬 빠르다. 'ㅅ'
- 업데이트: GLPK 에 2*n*(n+m) 크기의 행렬을 넘겨줘야 하는데 이게 너무 억울해서.. MathProg 를 쓰면 n*m 크기의 행렬만 넘겨줘도 되니까 빨라지겠지? 헤헷 하고 열심히 MathProg 를 배워서 모델 짜고 데이터 제너레이터 짜서 넣었더니... 안빨라져... orz 패러미터 설정하면서 갖고 놀거나, lp-solve 를 써보면 더 나아질 수도 있겠다 싶지만 귀찮아서 이만.
Convex by Boyd 을 보면, 많은 경우에 L2 optimized estimator 보다 L1 optimized estimator 가 더 robust 하다고 되어 있다. Outlier 에 훨씬 덜 민감하니까 당연한 노릇. 보면 LP 로 바꿔서 풀 수 있다고 되어 있다. 그래서 파이썬용 Convex Optimization 패키지인 cvxopt 로 풀어보았다.
파란 선을 중심으로 랜덤 샘플을 만들고 (가우시안 노이즈 추가) 적절한 개수로 엉뚱한 위치에 아웃라이어를 좀 던져둔다. 초록색 선이 L2 핏, 빨간색 선이 L1 핏이다. 훨씬 원래 분포에 가까운 결과를 보여줌.
cvxopt 샘플에 L1 approximation 이 샘플로 나와 있길래 좋아라 하면서 써먹어 보려 했는데 안돌아가... 행렬 크기가 안맞는대... 좀 검증해 보고 올리지 쩝. ㅠ.ㅠ 미워라 덕분에 무식하게 모델링해서 돌렸다.
속도를 재 보면:
초록색이 처음에 써 본 cvxopt 의 기본 solver 속도. 이건 뭔가 exponential 하게 올라가는 삘인데.. -_-;; 과연 현실적으로 쓸 수 있을까.. 하고 후덜덜했는데, GLPK 써보니 거의 linear 하게 올라가는 듯? 뭐.. 일단 지금은 무지하게 무식하게 짠 점도 있고,뭔가 problem structure 를 이용하면 좀 더 효율적으로 돌릴 수 있을 듯. 그러려면 저 Boyd 책을 봐야 하나.... -_- ..... 어쨌든 소스코드
lang:py
# A and b are numpy arrays.
def l1simple(A, b):
# n = number of samples
# m = number of features
n, m = A.shape
# minimize 1^T t
# subject to
# -t <= Ax - b <= t
# which translates to
# -Ax - t <= -b
# and
# Ax - t <= b
# We optimize [x,t]: m+n variables... oTL
c = matrix(m*[0.0] + n*[1.0])
G = matrix(0.0, (2*n, n+m))
for i in xrange(n):
for j in xrange(m):
G[i,j], G[i+n,j] = -A[i,j], A[i,j]
G[i,m+i], G[i+n,m+i] = -1.0, -1.0
h = matrix(hstack((-b,b)))
sol = solvers.lp(c, G, h)
print sol["x"][:m]
return sol["x"][:m]
def test(samples = 100, outliers = 10):
global A, b
import numpy.random
x = array([1.436743232, -0.2371234])
A = hstack((rand(samples,1) * 100, ones((samples,1))))
b = dot(A, x) + array([numpy.random.normal() * 10 for i in xrange(samples)])
for i in range(outliers):
devX = numpy.random.normal() * 10
devY = numpy.random.normal() * 10
if numpy.random.random() < 0.5:
A[i,0], b[i] = devX, devY + 150
else:
A[i,0], b[i] = 100+devX, devY
clf()
scatter(A[:,0], b)
st = time.time()
x1, _, __, ___ = linalg.lstsq(A,b)
L2 = time.time() - st
st = time.time()
x2 = l1simple(A, b)
L1 = time.time() - st
px = array([min(A[:,0]),max(A[:,0])])
plot(px, px*x[0]+x[1])
plot(px, px*x1[0]+x1[1])
plot(px, px*x2[0]+x2[1])
legend(("underlying","L2","L1"))
print "Solving Least Squares: %.8lf sec, Solving L1 Norm: %.8lf sec" % (L2, L1)
return L1
패키지 자체는 알기 쉬워서 좋은데 입력을 column major order 로 주는 것 때문에 죽도록 헷갈림;;;;;; 대체 왜 이렇게 해놓은거야아아;;;;;
나중에 짜본 GLPK 용 모델. Mathprog 는 심지어 vim syntax file 도 찾기 힘든 마이너 언어...
# L1-norm linear fit
/* sets */
set SAMPLES;
set FEATURES;
/* parameters */
param FeatureValues { i in SAMPLES, j in FEATURES };
param TargetVariable { i in SAMPLES };
/* decision variable X is given as { x, t }
x = weight of each feature
t = residual of each sample */
var x { j in FEATURES };
var t { i in SAMPLES } >= 0;
/* objective function */
minimize obj: sum{ i in SAMPLES } t[i];
/* constraints */
s.t. AB{i in SAMPLES} : -t[i] - sum{j in FEATURES} FeatureValues[i,j] * x[j]
<= -TargetVariable[i];
s.t. BC{i in SAMPLES} : -t[i] + sum{j in FEATURES} FeatureValues[i,j] * x[j]
<= TargetVariable[i];
end;
샘플 데이터 파일.
data;
set SAMPLES := 1 2 3 4 5 6 7 8;
set FEATURES := 1 2;
param FeatureValues: 1 2 :=
1 1 1
2 2 1
3 3 1
4 4 1
5 5 1
6 6 1
7 7 1
8 9 1;
param TargetVariable :=
1 2
2 3
3 4
4 5
5 6
6 7
7 8
8 1;
glpsol 깔고 glpsol -m l1fit.mod -d l1fit.data -w output 식으로 돌린다.
세상에서 제일 무식한 파이썬 드라이버.
lang:py
def l1glpk(A, b):
n, m = A.shape
fp = open("l1glpk.data", "w")
fp.write(
"""data;
set SAMPLES := %s;
set FEATURES := %s;
param FeatureValues: %s :=\n""" % (" ".join(map(str, xrange(1,n+1))),
" ".join(map(str, xrange(1,m+1))),
" ".join(map(str, xrange(1,m+1)))));
for i in xrange(n):
fp.write("%d\t%s\n" % (i+1, " ".join(["%g" % v for v in A[i]])));
fp.write(";\n");
fp.write("param TargetVariable :=\n");
for i in xrange(n):
fp.write("%d\t%g\n" % (i+1, b[i]))
fp.write(";\n");
fp.close()
os.system("glpsol -m l1fit.mod -d l1glpk.data -w out")
return [float(line.split()[1]) for line in open("out").readlines()[2*n+3:2*n+3+m]]
GLPK 도 파이썬 드라이버 있던데 담에 쓸일 있으면 그거나 써야겠다....
