2013-03-09

多項式計算を確率的ロジックで

The Synthesis of Robust Polynomial Arithmetic with Stochastic Logic
Bernstein polynomials
利点
- 算術演算が単純な論理回路で構成可能になる
入出力
- 0/1の値を持つ長さnのベクトルXを入力として、0/1の出力Yをする
- Xのベクトルの要素が1をとる確率 $p_i;i=1,2,...,n$ の多項式としてYが1となる確率があらわされる
入出力を繰り返す
- Serial stream,または Parallel stream
任意の多項式を実装するためにBernstein polynomialを導入する
- - $B_i^n(t) = \begin{pmatrix} n\\i \end{pmatrix} t^i(1-t)^{n-i}$
Bernstein polynomial の基底たち
多項式のBernstein polynomials分解
- $f(t) = \sum_{i=0}^n a_i t^i = \sum_{i=0}^n b_i B_i^n(t)$ のとき
- $b_i = \sum_{j=0}^i \frac{\begin{pmatrix}i\\n\end{pmatrxi}}{\begin{pmatrix}n\\j\end{pmatrix}} a_j$
多項式のBernstein polynomials 分解は一意ではなく、power formでn次の多項式はm>=n次のBernstein polynomials 分解ができて、そのときは、次のような関係で、１次ずつ挙げていくことができる
- $b_0^{m+1} = b_0^m$
- $b_i^{m+1} = (1-\frac{i}{m+1})b_i^m + \frac{i}{m+1}b_{i-1}^m; i = 1,2,...,m$
- $b_{m+1}^{m+1} = b_m^m$
加算器(adder)のこと
- ２つのバイナリ値を受け取って、2ケタでその和を表すのが1-bit adder
- 同様に2-bit adder, 3-bit adder, 4-bit adderなどもあり、
- 結局、n個のバイナリ値を受け取って、n個のうちi個が1で残りのn-i個が0であるようなときに、0,1,2,3,...,nに対応するn+1桁の数値(0/1)を使って、i個が1だ、という情報にすることができる(0,1,..,としたときの0から数えてi+1番目を1にして残りを全部0にする)
- この回路は二分岐木構成で深くなっている、そんな回路
Bernstein polynomial表記になっている多項式の値を計算するような回路作成が、上述の回路を使うとできる(ことが示せる)
- ただし、すべての $b_i^n$ は0-1の間の値
- i/n 個が1だ、という長さn+1のboolean ベクトルと、n+1個の値のベクトルZのブーリアン積のORをとると、Zのi番目の要素の値が取り出せる
- n個の確率変数が、いずれも確率tで値1をとるときi+1番目のビットが立つという確率は(n個ののうちi個が1のとき) $\begin{pmatrix}n\\i\end{pmatrix} t^i(1-t)^{n-i}=B_i^n(t)$ とBernstein 多項式になっている
- 今、n+1個の確率があって、それが確率 $b_i^n$ で値１をとるとすると $\sum_{i=0}^n b_i^n B_i^n(t)$ は、この演算ユニットが1を返す確率になっている
- このことは、 $b_i^n;i=0,1,...,n$ という $[0,1]$ の値をとる長さ $n+1$ の任意の値のベクトルを係数とするBernstein polynomial 形式の多項式の $t$ のときの値が確率として(複数試行をして、そのうち１が立つ割合として)計算できることを意味する
- ここで任意のpower formの１変数多項式はBernstein polynomial形式で書けることはすでに示してあるから、どんなpower form多項式も計算できることになる
- ただし、 $b_i^n$ が $[0,1]$ に収まらないとだめなわけで、power formの多項式の係数から、Bernstein polynomial 形式の係数に直したうえで、その次数を上げていって、全部の $b_i^n$ が $[0,1]$ に収まるようにしたらよい
そんなことができるような多項式は、特定の制約があるのか、あるならそれはどんな制約か、というと：
- $g(t) = 0$ , $g(t)=1$
- $0 < g(t) < 1, \forall t \in (0,1); 0 \le g(0),g(1) \le 1$
さて、制約はあるものの、多項式の値を確率的に求める回路の設計ができることは分かった
いわゆる普通に計算機で計算するのと比較する
- 計算精度は普通の計算の場合には、ビット数・桁数で決まる。確率回路の場合は、何セットの乱数を作るかで決まる。その数を相当値とすることで、似たり寄ったりの精度が確保できる
- そのうえで、計算の遅さとかを比較する
- また、入力データにノイズがあるときのパフォーマンスも比較する
- 特殊関数とかは、どのみち級数展開で近似しているので、そんな関数での性能比較もしている

# Bernstein polynomial

Bernstein <- function(n,m,x){
	if(m<0){
		return(rep(0,length(x)))
	}else if(m > n){
		return(rep(0,length(x)))
	}else{
		#tmp <-lgamma(n+1)-lgamma(m+1)-lgamma(n-m+1) +m * log(1-x) + (n-m) * log(x)
		#return(exp(tmp))
		return(gamma(n+1)/gamma(m+1)/gamma(n-m+1)*x^m*(1-x)^(n-m))
	}
}

x <- seq(from=-3,to=4,length=100)

n <- 3
m <- 1
plot(x,Bernstein(n,m,x))



k <- 4
as <- sample(-5:5,k)
# power formをBernstein formへ
Power2Bernstein <- function(a){
	b <- rep(0,length(a))
	n <- length(a)-1
	for(i in 1:length(b)){
		for(j in 1:i){
			tmp.i <- i-1
			tmp.j <- j-1
			b[i] <- b[i] + factorial(tmp.i)*factorial(n-tmp.j)/factorial(tmp.i-tmp.j)/factorial(n)*a[j]
		}
	}
	b
}

bs <- Power2Bernstein(as)
bs

y.p <- rep(0,length(x))
for(i in 1:length(as)){
	y.p <- y.p + x^(i-1) * as[i]
}

plot(x,y.p,type="l")
y.b <- rep(0,length(x))
n <- length(bs)-1
for(i in 1:length(bs)){
	y.b <- y.b + bs[i] * Bernstein(n,(i-1),x)
}

plot(y.b,y.p)
# degree n の Bernstein分解を m >= nに変換する

Bernstein.1plus <- function(b){
	m <- length(b)-1
	ret <- rep(0,length(b)+1)
	ret[1] <- b[1]
	ret[length(ret)] <- b[length(b)]
	for(i in 2:length(b)){
		ret[i] <- (1-(i-1)/(m+1))*b[i]+(i-1)/(m+1)*b[i-1]
	}
	ret
}
Bernstein.m <- function(b,m){
	ret <- b
	n <- length(b)-1
	for(i in 1:(m-n)){
		ret <- Bernstein.1plus(ret)
	}
	ret
}

n <- length(bs)-1

m <- n + 3

bs.m <- Bernstein.m(bs,m)

y.b.m <- rep(0,length(x))
n <- length(bs.m)-1
for(i in 1:length(bs.m)){
	y.b.m <- y.b.m + bs.m[i] * Bernstein(n,(i-1),x)
}

plot(y.b,y.b.m)

- 多項式の近似法の一つ
- [0,1]の範囲で特に有用
- ベジエ曲線の基礎(n点を通るn-1次曲線)(イラストレータとかでも使われている曲線)