Vanishing Component Analysisがダメだと思うわけ
少し前に、少し話題になっていたVanishing Component Analysis(VCA)を実装してみたのだけど、よく考えると、これは全然ダメじゃないかという気がしてきた。
Vanishing Component Analysis
http://jmlr.org/proceedings/papers/v28/livni13.pdf
(2014/04/12追記)結論としては、わたしの実装が間違ってたorz下のコードは修正済み(行列Aを転置してなかった。。)。円や楕円程度なら、ほぼ確実にそれっぽい答えが返ってくるっぽい(但し、epsパラメータ/論文のεを小さくし過ぎると怪しいっぽい?)。修正したコードで半径1の円近傍に点をまいて、VCAで得られた二次多項式(これでは、一次元なので、もっと高次の多項式も返ってくるけど)を、Wolfram alphaで描画してみると、こんな感じ。
(2014/04/12追記)以下に書いたような"Huモーメント不変量の推定"みたいな、より複雑で高次の多項式でも、うまくいくのかは、また今度検証
VCAは、アフィン空間上の点の集合を与えた時、その点集合の(近似的な)定義イデアルを決めてやるアルゴリズムらしい。話題になってた時は、全然興味がなかったのだけど、ふと、これが本当にうまく機能するなら、Huモーメント不変量を推定させることができるのでは、ということを思った。つまり、x,yをそれぞれ、2つの画像のモーメントを低次から順にいくつか並べたものとして、これらのモーメントが、互いに拡大縮小・回転・並進を施してうつりあう画像に由来する場合に、d(x)=d(y)となる関数dを発見するという問題を考える。Huモーメント不変量は、条件を満たすけども、これは、有理関数であるので、VCAの適用範囲外である。それに、VCAは定義イデアルを発見するので、上のdを直接求めることはできない。それで少し見方を変えて、モーメントx,yが互いに拡大縮小・回転・並進を施してうつりあう画像に由来する場合、F(x,y)=0となる多項式Fを探すという形で考える。d(x)=P(x)/Q(x)で、P,Qは多項式とすると、F(x,y)=P(x)Q(y)-Q(x)P(y)は条件を満たす多項式となる。それにまぁ、画像を比較する時、最終的に欲しいのは、不変量そのものより、(モーメント間で定義される)距離関数であることを思うと、Fを探すほうが、reasonableな気もする
そんなことを思って、期待して実装して見たのだけど、例えば、半径1の円周近傍に点をばらまいて、x^2+y^2=1を推定してくれるかというと、全然そんな風にはならない。わたしの実装が間違ってるという可能性もあるけど、ある多項式fが、零点集合の定義多項式として、十分よいものであるかどうか判定するのは容易ではなく、この部分に致命的な問題があると思う。要するに、最小二乗法で、y=ax+bでfittingする代わりに、ax+by+c=0でfittingするようなもので、後者の場合、何を最小化すればよいのか適切な基準は定かでない。それだけでなく、a,b,cに同じ定数をかけて(例えば、0.01とか)小さくすると、サンプル点集合の上で、結構小さな値を取るので、何かよい近似に見えてしまったりする。VCAでも、係数が、どんどん小さくなってしまったような答えが返ってきて、こういうことが起きてるのだと思う(具体的には、論文のFindRangeNullで、行列Aを作るのに、サンプル点集合上での多項式の値を使ってるので、これがダメ)
もう一つは、より微妙なことであるけど、論文のTheorem5.2とか見ると、ε= 0でVCAを走らせると、完全な定義イデアルが得られるぜ、と書いてあるのだけど、そもそも、サンプル点集合は、所詮有限集合で0次元なので、その定義イデアルが欲しいわけでもないし、これが本当にアルゴリズムが満たしてほしい性質と言えるか曖昧なとこである(サンプル点集合以外の部分でも、あちこち0になるようでは、論外なので、この定理は、そんなことが起きないと保証しているということなのだろうけど)。半径1の円周上の点列を与えたら、x^2+y^2=1と、その他なんか意味のない高次の多項式が得られる！という風になってれば、使えると言えるけれども、そううまいこといくという保証もない気がする。これは正規表現の学習とも似ている側面がある。有限個の正例p1,...,pnがあれば、それにマッチする正規表現は、p1|....|pnで済むけど、これは、勿論通常欲しい答えとはかけ離れている(あるいは、普通の多項式回帰でも、次数をあげていけば、いくらでもいいfittingができるけど、やはり、これは欲しい答えとはかけ離れる。この場合は、適当な情報量基準でトレードオフを測ったりする)。正規表現の学習の場合には、この問題に対処するために極限同定という概念が考えられたのだった。今回の問題でも、アルゴリズムが満たす性質として、どういうものが望ましいのか、よく考える必要がある気がする(極限同定可能性に類似した概念は、今回も意味がありそうな気はする。円周x^2+y^2=1上の乱数列を与えた時、最初のN個を取って定義イデアルを考えても仕方ないけど、N=∞の極限に於いては、それは円周上の稠密な可算集合となって、それら全体で消える多項式は、x^2+y^2=1のみとなる。勿論、実際にはデータに必ず誤差があるので、もっと複雑である)
そういうわけで、わたしの理解が間違ってる可能性もあるけど、これが、どこだかのbest paperだとかいうのは、みんなちゃんと読んでないんじゃないか的な
一応、Pythonによる実装。WeylSequenceというクラスは、円周x^2+y^2=1に分布する点列を作るためのもので、VCAに必要なものではない。名前の由来は、Weylの無理数回転。
LDS(1)
http://d.hatena.ne.jp/m-a-o/20091212#p3
で実装したものを流用している
import numpy as np
from math import cos,acos,pi
class Monomial(object):
def __init__(self):
self.degs = {}
def __hash__(self):
return tuple(self.degs.items()).__hash__()
def __eq__(self , other):
if len(self.degs)!=len(self.degs):
return False
else:
for k,v in self.degs.iteritems():
if other.degs.get(k,0)!=v:
return False
return True
def __mul__(self , other):
ret = Monomial()
vars = set(self.degs.keys()) | set(other.degs.keys())
for var in vars:
d = self.degs.get(var , 0) + other.degs.get(var , 0)
if d!=0:
ret.degs[var] = d
return ret
def __str__(self):
if len(self.degs)==0:
return str(1)
else:
terms = []
for k,v in self.degs.iteritems():
if v==1:
terms.append( str(k) )
else:
terms.append( "%s**%d" % (k,v) )
return "*".join( terms )
def __repr__(self):
if len(self.degs)==0:
return str(1)
else:
terms = []
for k,v in self.degs.iteritems():
if v==1:
terms.append( str(k) )
else:
terms.append( "%s^%d" % (k,v) )
return "*".join( terms )
class Polynomial(object):
def __init__(self):
self.coeffs = {}
def __eq__(self , other):
if not isinstance(other , Polynomial):
if len(self.coeffs)>1:
return False
elif len(self.coeffs)==1:
k,v = self.coeffs.items()[0]
return (len(k.degs)==0 and v==other)
else:
return (other==0)
elif len(self.coeffs)!=len(other.coeffs):
return False
else:
for k,v in self.coeffs.iteritems():
if other.coeffs.get(k , 0)!=v:
return False
return True
def __ne__(self,other):
return not(self==other)
def __add__(self , other):
if isinstance(other , Polynomial):
ms = set(self.coeffs.keys()) | set(other.coeffs.keys())
ret = Polynomial()
for m in ms:
c = self.coeffs.get(m , 0) + other.coeffs.get(m , 0)
if c!=0:ret.coeffs[m] = c
return ret
else:
ret = Polynomial()
cm = Monomial()
for m,c in self.coeffs.iteritems():
ret.coeffs[m] = c
if other!=0:ret.coeffs[cm] = ret.coeffs.get(cm , 0) + other
return ret
def __radd__(self , other):
return self.__add__(other)
def __sub__(self , other):
if isinstance(other , Polynomial):
ms = set(self.coeffs.keys()) | set(other.coeffs.keys())
ret = Polynomial()
for m in ms:
c = self.coeffs.get(m , 0) - other.coeffs.get(m , 0)
if c!=0:ret.coeffs[m] = c
return ret
else:
ret = Polynomial()
cm = Monomial()
for m,c in self.coeffs.iteritems():
ret.coeffs[m] = c
if other!=0:ret.coeffs[cm] = ret.coeffs.get(cm , 0) - other
return ret
def __neg__(self):
ret = Polynomial()
for k,v in self.coeffs.iteritems():
ret.coeffs[k] = -v
return ret
def __rsub__(self , other):
tmp = self.__sub__(other)
return tmp.__neg__()
def __mul__(self , other):
if isinstance(other , Polynomial):
ret = Polynomial()
for m1,c1 in self.coeffs.iteritems():
for m2,c2 in other.coeffs.iteritems():
m = m1*m2
c = c1*c2
if c!=0:ret.coeffs[m] = ret.coeffs.get(m , 0) + c
return ret
else:
ret = Polynomial()
for m,c in self.coeffs.iteritems():
c0 = c*other
if c0!=0:ret.coeffs[m] = ret.coeffs.get(m , 0) + c0
return ret
def __rmul__(self , other):
return self.__mul__(other)
def __pow__(self , n):
assert(isinstance(n,int))
assert(n>=0)
ret = 1
for _ in xrange(n):
ret *= self
return ret
def __call__(self,**args): #-- 代入操作
vars = [str(v) for v in args.keys()]
ucvars = [("v_%d" % n) for n in xrange(len(vars))]
terms = []
for (m,c) in self.coeffs.iteritems():
if len(m.degs)==0:
terms.append( str(c) )
else:
mterms = []
for (varname , degree) in m.degs.iteritems():
if varname in vars:
mterms.append('%s**%d' % (ucvars[vars.index(varname)] , degree))
else:
mterms.append('Symbol("%s")**%d' % (str(varname) , degree))
terms.append('%s*%s' % (str(c) , '*'.join(mterms)))
if len(terms)==0:
pyexpr = ("lambda %s:%s" % (",".join(ucvars) , "0"))
else:
pyexpr = ("lambda %s:%s" % (",".join(ucvars) , "+".join(terms)))
return apply(eval(pyexpr) , args.values())
def __str__(self):
terms = []
for (m,c) in self.coeffs.iteritems():
if len(m.degs.keys())==0:
terms.append( str(c) )
elif c==1:
terms.append( str(m) )
elif c==-1:
terms.append( "-%s" % str(m) )
else:
terms.append( "%s*%s" % (str(c),str(m)) )
if len(terms)!=0:
return "+".join(terms)
else:
return "0"
def __repr__(self):
"""
strとreprの(主な)違いは、べきをpython流に"**"で表すか"^"で表すか
"""
terms = []
for n,(m,c) in enumerate(self.coeffs.iteritems()):
if n>0 and c>0:
terms.append('+')
if len(m.degs.keys())==0:
terms.append( str(c) )
elif c==1:
terms.append( repr(m) )
elif c==-1:
terms.append( "-%s" % repr(m) )
else:
terms.append( "%s*%s" % (str(c),repr(m)) )
if len(terms)!=0:
return "".join(terms)
else:
return "0"
class Symbol(Polynomial):
def __init__(self , s):
super(Symbol , self).__init__()
m = Monomial()
m.degs[s] = 1
self.coeffs[m] = 1
class Constant(Polynomial):
def __init__(self , c):
super(Constant , self).__init__()
m = Monomial()
self.coeffs[m] = c
class WeylSequence:
def __init__(self, _x):
self.prev = 1
self.x = _x
self.cur = _x
self.mirror = False
self.reflect = False
self.theta = acos(_x)/pi
def __iter__(self):
return self
def next(self):
val = acos(self.cur)/pi
if self.mirror:
ret = 1-val
else:
ret = val
#ret = 1-val if self.mirror else val
#Chebyshev polynomial
next = 2*(self.x)*(self.cur) - (self.prev)
self.prev = self.cur
self.cur = next
if val < self.theta or 1-val < self.theta:
if self.reflect:
self.reflect = False
else:
self.mirror = not(self.mirror)
self.reflect = True
return ret
def VCA(S , varnames , eps=1.0e-3):
def findRangeNull(F,C,S,vars,eps):
k = len(C)
A,FP=[] , []
for f in C:
tmp = f
for g in F:
tmp = tmp - np.sum([apply(f,[],dict(zip(vars,pt)))*apply(g,[],dict(zip(vars,pt))) for pt in S])*g
A.append( [apply(tmp,[],dict(zip(vars,pt))) for pt in S] )
FP.append(tmp)
A = np.asmatrix(np.array(A))
L,D,U = np.linalg.svd(A.transpose())
F1,V1=[],[]
for i in range(k):
g = np.sum([U[i,j]*FP[j] for j in range(k)])
if i < len(D) and D[i]>eps:
val = np.sum([apply(g,[],dict(zip(vars,pt)))**2 for pt in S])
val = 1.0/np.sqrt(val)
F1.append( val*g )
else:
V1.append(g)
return F1,V1
#-- VCA main --
assert(type(S)==list and len(S) > 0),"invalid argument"
n,m=len(varnames) , len(S)
F = [Constant(np.sqrt(1.0/m))]
C = [Symbol(name) for name in varnames]
(F1,V) = findRangeNull(F , C , S , varnames , eps)
F , NF = F+F1,F1
while True:
C = [g*h for g in NF for h in F1]
if C==[]:break
(NF,NV) = findRangeNull(F , C , S , varnames , eps)
F , V = F+NF,V+NV
yield V
import itertools
if __name__=="__main__":
seq = WeylSequence(0.42)
circle = [((np.random.rand()-0.5)*0.01+np.cos(2*t*np.pi),(np.random.rand()-0.5)*0.01+np.sin(2*t*np.pi)) for t in itertools.islice(seq,0,180)]
for n,vc in enumerate( VCA(circle,['x' , 'y'],eps=1.0e-3) ):
if len(vc)>0:print repr(vc[0]);break
以下は、上記の実行の結果得られた二次の多項式の零点を、Wolfram alphaで描画した結果。まぁだめっぽい(別に収束してないとか、そういうことではないと思う)
