2016-06-21

JuliaのFloatとRoundingの奇妙な挙動[追記]

JuliaのFloatとRoundingの奇妙な挙動

with_rounding(Float64, RoundNearest) do
    exp(-0.999)
end

0.36824750461366296

with_rounding(Float64, RoundUp) do
    exp(-0.999)
end

0.3682475046136629

with_rounding(Float64, RoundDown) do
    exp(-0.999)
end

0.3682475046136629

本来、 Down <= Nearest <= Upとならなければならないところが
見た目上は Up < Nearest となってしまっている。

追記: Up <= Nearest を確認してみると成立しているため、計算上の問題はない。表記がおかしいだけで、内部的には影響は小さいようだ

もちろん、Roundingは四則演算と平方根のみなのだが、結果が逆になったのは興味深い。
JuliaのRoundingのドキュメントを参照すると

Waring

Rounding is generally only correct for basic arithmetic functions (+(), -(), *(), /() and sqrt()) and type conversion operations. Many other functions assume the default RoundNearest mode is set, and can give erroneous results when operating under other rounding modes.

とあるので、expはおかしくなってしまったんだろう。issueをさらに調べると
exp(1) != e #3370　　

Since exp(1) just calls openlibm's exp function, this appears to be an openlibm issue.

とあるので、openlibm側の問題のようだ。
当たり前だがexpにはRoundingは効かず、場合によっては妙な挙動を起こすことがわかる。
なお、これは利用者側が気をつけなければならない仕様であるので、Julia側に落ち度はない。
with_roundingやset_rounding内では、四則演算と平方根以外は利用しないほうが良いだろう。
(できれば、四則演算と平方根以外は最近点丸めで計算してくれるとありがたいんだけど…… )

なお、mpfrを利用するBigFloatではこのような問題は起こらない

with_rounding(BigFloat, RoundNearest) do
    exp(BigFloat(-0.999))
end

3.682475046136629215368118018604961977815208270355409127685022306741116356975955e-01

with_rounding(BigFloat, RoundUp) do
    exp(BigFloat(-0.999))
end

3.682475046136629215368118018604961977815208270355409127685022306741116356975998e-01

with_rounding(BigFloat, RoundDown) do
    exp(BigFloat(-0.999))
end

3.682475046136629215368118018604961977815208270355409127685022306741116356975955e-01

見づらいんで、大小比較すると

with_rounding(BigFloat, RoundNearest) do
    global near
    near = exp(BigFloat(-0.999))
end

with_rounding(BigFloat, RoundUp) do
    global up
    up = exp(BigFloat(-0.999))
end

with_rounding(BigFloat, RoundDown) do
    global down
    down = exp(BigFloat(-0.999))
end

3.682475046136629215368118018604961977815208270355409127685022306741116356975955e-01

down <= near <= up

true

となるのでmpfr上では大丈夫
なお、この現象はJulia 0.4.5で確認し、windows 10(64bit)、OS X 10.11.5で確認した

2016-06-14

numpyを使って、行列に対するtwosum, twoproductを高速化する

　柏木先生の考案された改良型twosum, twoproductを行列に適用し、それをnumpy流にベクトル化し、高速化した。
コードと結果は以下
github.com
github.com 環境は

Mac Book Pro 2012 mid (ただし、メモリを16GBに増やしている)

anaconda 4.0.8 (以下anaconda内のパッケージ)

 * python 3.5.1  
 * scipy 0.17.0 np110py35_4  
 * numpy 1.10.4 py35_2  
 * mkl 11.3.3

結果は

%timeit twosum(a, b)

10 loops, best of 3: 80.7 ms per loop

%timeit twosum_old(a, b)

1 loop, best of 3: 2 s per loop

%timeit twoproduct(a, b)

10 loops, best of 3: 76.8 ms per loop

%timeit twoproduct_old(a, b)

1 loop, best of 3: 4.22 s per loop

　特に特殊な手法は用いていないが、ベクトル化するだけで大幅に高速化した。

2016-01-31

分数区間を用いた、浮動小数点数の区間包含についてのアイデア

分数区間を用いた、浮動小数点数の区間包含について

python 3.5.1で動作を確認
区間演算の実装について(2)をpython上解決する場合のアイデア。
いろいろと検証が足りていないが公開はする。
jupyter notebook用のコードもアップロードした。

import pint as pn
from pint import roundfloat as rf
from pint import roundmode as rdm
import fractions

def frac_interval_proto(a):
    aH, aL = rf.split(a)
    if aL ==0:
        answer = pn.interval(a)
    else:
        aS = rf.succ(a)
        aP = rf.pred(a)
        aS_c, aS_p = aS.as_integer_ratio()
        aP_c, aP_p = aP.as_integer_ratio()
        a_c, a_p = a.as_integer_ratio()
        bS_c = a_c * aS_p + aS_c * a_p
        bS_p = a_p * aS_p * 2
        bP_c = a_c  * aP_p + aP_c * a_p
        bP_p = a_p * aP_p * 2
        answer = pn.interval(0.)
        answer.inf = fractions.Fraction(bP_c, bP_p)
        answer.sup = fractions.Fraction(bS_c, bS_p)
    return answer

今回は、代入演算子を用いて区間型を生成した際に、初期値を包含しない区間が発生するのを防ぐ分数区間を生成する。
すなわち

a = 0.1

itv_a = pn.interval(a)

format(itv_a.inf, '.17g')

'0.10000000000000001'

format(itv_a.sup, '.17g')

'0.10000000000000001'

のような区間が発生することを防ぐ。
なお、今回以外の解決方法としては

itv_a.inf = rf.pred(a)
itv_a.sup = rf.succ(a)

format(itv_a.inf, '.17g')

'0.099999999999999992'

format(itv_a.sup,'.17g')

'0.10000000000000002'

のようにsucc predを用いることで解決可能である。
しかし、今回のアルゴリズムのほうが、より区間幅を縮小できる。

まず、上位ビットと下位ビットを分離する

aH, aL = rf.split(a)

もし、下位ビットが存在する場合、else文以降の通りにする。
(下位ビットが存在しない浮動小数点数で点区間を生成した場合、元の初期値を包含する区間になる可能性が高い
元の初期値の下位ビットが全て1で、近似の際にくりあがりで下位ビットが消去された場合、おかしくなるかも)

if aL ==0:
    print(True)
    answer = pn.interval(a)
else:
    print(False)

False

以下はelseの続き

次にaをsucc predし、as_intger_ratio()で近似分数にする

aS = rf.succ(a)
aP = rf.pred(a)

aS_c, aS_p = aS.as_integer_ratio()
aP_c, aP_p = aP.as_integer_ratio()

a_c, a_p = a.as_integer_ratio()

次に、aとsuccしたaの中点と、aとpredしたaの中点を求める。

bS_c = a_c * aS_p + aS_c * a_p
bS_p = a_p * aS_p * 2
bP_c = a_c  * aP_p + aP_c * a_p
bP_p = a_p * aP_p * 2

こうして得た、分数bSとbPとaは以下の性質を満たす

浮動小数点数として評価した場合、同値である。

2016-02-01 例外がたくさん発生した
正しくは bP_c / bP_p <= a <= bS_c / bS_p

bS_c / bS_p == a == bP_c / bP_p

True

分数として評価するとbP <= a <= bSが成立する
(以下は、分母と通分するためbS_pやら、a_pがかけられている。今回は bP < a < bSとなっていることをわかりやすくするため、等号の場合は考慮していない)

a <= bS

a_c * bS_p < bS_c * a_p

True

bP <= a

bP_c * a_p < a_c * bP_p

True

よって、分数で評価すると大小関係が存在するが、浮動小数点数で評価すると値が一致する分数区間が得られた

なお、succ predで得られたitv_aと比較すると

rf.pred(a) < bP_c / bP_p

True

bS_c / bS_p < rf.succ(a)

True

となり、succ predよりも狭い区間幅が得られた。

answer = pn.interval(0.)

answer.inf = fractions.Fraction(bP_c, bP_p)
answer.sup = fractions.Fraction(bS_c, bS_p)

print(answer)

[14411518807585587/144115188075855872,14411518807585589/144115188075855872]

今回は

引数が正の場合しか考慮していない
float.as_intger_ratio()で得られる分数が、最近点丸めで近似されたaと厳密に同じかどうか考慮していない
分数演算の有効桁数に由来する計算誤差を考慮していない(分母、分子ともに巨大な整数になりがちである)
その後の分数の演算方法を考慮していない(作ったは良いが、使い方を考えていない)

などの問題がある
最後の分数の演算に関しては以下のようにするアイデアがあるが、検証はしていない

分数 * 分数: 分数演算をしたあと、丸め制御で適切な浮動小数点数にする
分数 * 浮動小数点数(下位ビットがない): 分子 * 小数を計算後、丸め制御で適切な浮動小数点数にする
分数 * 浮動小数点数(下位ビットがある): 小数を分数に変換したあと、分数 * 分数に帰着させる

2016-01-05

pythonの区間演算ライブラリ進捗2016-01-05

実は大体できた(唐突)。
今回はvmathクラスの紹介

github.com

import pint as pn

from pint import vmath

a = 2.

itv_a = pn.interval(a)

print(itv_a)

[2.0,2.0]

vmath.sqrt(a)

1.4142135623730951

floatを代入するとfloatが返ってくる

vmath.sqrt(itv_a)

[1.414213562373095,1.4142135623730951]

intervalを代入するとintervalが返ってくる
じゃあ、少し変なことをしてみる

itvmath = pn.vmath()

itvmath.arg_type = pn.interval

このitvmathのmath関数を使う時は代入する前にinterval型に変換するように設定

itvmath.sqrt(a)

[1.414213562373095,1.4142135623730951]

返ってくるのはinterval型になる。
itvmath.sqrt(pn.interval(a))と同じ

itvmath.sqrt(itv_a)

[1.414213562373095,1.4142135623730951]

itvmath2 = pn.vmath()

itvmath2.return_type = pn.interval

itvmath2.sqrt(a)

[1.4142135623730951,1.4142135623730951]

return_typeはreturnする直前に変換を行う
return return_type(answer)となるそのため、変換するタイミングが違うitvmath.sqrt(2.)とitvmath2.sqrt(2.)は結果が少し違う
このようにmath関数の振る舞いを自由に変更することができるのがvmathクラスの特徴です。
今はinterval型とビルトインされている型にしか対応してませんが、そのうち増やす予定

2015-12-21

メモ:as_interger_ratio()では正確に元の分数には戻せない

as_interger_ratio()では正確に元の分数には戻せない

a = 0.1 を区間型に代入する際にうまく包含できないので、色々実験してみた。

何が問題なのか?

import pint as pn

x=pn.interval(0.1)

と、2進数小数で表現不可能な数字なので下限と上限は

format(x.inf,'.17g')

'0.10000000000000001'

format(x.sup,'.17g')

'0.10000000000000001'

となり、0.1はうまく包含できない

解決策になりそうだったもの

pythonのfloat型にはas_integer_ratio()という関数があり、

a = 0.1

x,y = a.as_integer_ratio()

print(x)

3602879701896397

print(y)

36028797018963968

x/y

0.1

となり、0.1の分母と分数を取得できる
これを用いれば、さきほどのinterval型への代入問題を解決できそうだが

sup = pn.roundfloat.rddiv(x, y, pn.roundmode.up)

inf = pn.roundfloat.rddiv(x, y, pn.roundmode.down)

format(sup,'.17g')

'0.10000000000000001'

format(inf,'.17g')

'0.10000000000000001'

となり、0.1を包含するように計算できない。
なぜならば a=0.1 とした時点で 0.1 から、0.1よりずれた小数に近似されるため、as_interger_ratio で取得できる分母と分子は 0.1よりずれた小数に合うように取得されるためである。

誤解回避

0.1...となる分数を取得することは可能だが
厳密に0.1になるような分数はas_integer_ratioでは取得できない
(このように厳密に計算するとき以外は大した問題にならないと思われる)

2015-12-21

メモ:twosumアルゴリズムの研究

twosumアルゴリズムの研究

from pint import roundfloat as rf

from pint import roundmode as rdm

a=0.1

b=0.2

誤差が発生するパターン

rf.rdadd(a, b, rdm.nearest)

0.30000000000000004

pythonの区間演算ライブラリ進捗2015-11-16

成果

github.com
ソースコードが長くなりそうだったので分離した。
本当は一つのファイルで行きたかったけど、管理や読み込みなどを考えたら分離せざるを得なくなった。