作者znmkhxrw (QQ)
看板Math
标题[其他] 电脑浮点数实作两函数问题
时间Wed Jul 19 00:28:59 2023
想请问一下如图这两个函数f跟g如何化简让电脑的误差最小
https://www.desmos.com/calculator/sjxxpwsfwt
其中f(x) = (2^x-(x+1))/(x(x-1)), 0<=x<=1
g(x) = (x^2+1-2^x)/(x(x-1)), 0<=x<=1
可以看到当x越接近0或是1时, f跟g都会遇到0/0的问题, 很吃浮点数的精度
但是由数学理论值可以知道这两边的极限值都存在
因此感觉可以化简成一个让电脑可以不用面临极端值运算的形式
尝试过把2^x做泰勒展开, 但是後续无穷项还是面临项数问题...
再请板友帮忙, 谢谢!
--
※ 发信站: 批踢踢实业坊(ptt.cc), 来自: 59.102.225.191 (台湾)
※ 文章网址: https://webptt.com/cn.aspx?n=bbs/Math/M.1689697741.A.853.html
1F：推 jchin : 在x=0或1, 分子分母取微分? 07/19 06:02
2F：→ yhliu : 当遇到分母为 0 或绝对值低於某值，检查分子是否亦 07/19 07:57
3F：→ yhliu : 接近 0, 如是，则 x 代以另一近似值计算之；否则ERR 07/19 07:59
4F：→ yhliu : 分子分母微分法的好处是此例只要 x 靠近 0 或 1 都 07/19 08:14
5F：→ yhliu : 能得到正确结果。而上面我提出的方法好处是适用一般 07/19 08:15
6F：→ yhliu : 0/0 不定式，而不需改函数式，缺点一当然在此例就是 07/19 08:17
7F：→ yhliu : x 靠近 0 或 1 会遇到除以 0 或溢值问题，另外就是 07/19 08:18
8F：→ yhliu : 结果也许不精确，例如第一个函数在 x<=0.1e9 时开始 07/19 08:21
9F：→ yhliu : 不正确。 07/19 08:21
谢谢j大跟y大的idea, 感觉不管微分还是近似值计算, 都还是要决定一个threshold?
即大於某个threshold采取原式, 小於某个threshold采取近似式?
原本是不想要有这个"tunable parameter", 但是目前看起来只能这样了QQ
好像跟伪逆矩阵一样, 多小的eigenvalue视作0, 否则视作倒数
10F：推 j0958322080 : 试试 bug number algo. 07/19 13:41
11F：→ j0958322080 : bug-->big 07/19 13:41
谢谢j大关键字
12F：→ simonjen : 你可以自己写一个运算用字串去存 07/20 08:45
13F：推 LPH66 : 楼上那就是上面提的 big number (大数) 演算法 07/21 17:11
我查big number algorithm没看到特指哪一种算法, 以及s大说的"用字串去存"指的是什麽
有Reference可以提供吗? 谢谢!
14F：→ PPguest : 我估狗搜寻一下,看起来是指用其他方法如字串,阵列来 07/21 22:17
15F：→ PPguest : 存例如100位有效数字的数 07/21 22:19
16F：推 LPH66 : 大数的基本概念就是这样没错, 用大阵列去模拟小学时 07/22 17:30
17F：→ LPH66 : 做基本加减乘除的演算法 07/22 17:30
18F：→ LPH66 : 例如 GMP (GNU 多重精度运算库) 就是一个用 C 语言 07/22 17:32
19F：→ LPH66 : 写的这类型函式库 07/22 17:32
喔喔 意思就是多开buffer去存更多资料就是了?
20F：推 j0958322080 : 不一定要字串，你也可以开个 int64_t 再注意进位 07/22 20:01
21F：推 wohtp : f(x)在 x = 0 附近改写一下 07/22 22:13
22F：→ wohtp : f(x) = -1/(x-1) + (exp(x ln 2) -1)/x * 1/(x-1) 07/22 22:16
23F：→ wohtp : 主要问题是 (exp(x ln 2) -1)/x 这项 07/22 22:17
24F：→ wohtp : 我不知道你用什麽，但是大部分函式库应该都有expm1 07/22 22:18
25F：→ wohtp : 改用这个就不会有误差问题 07/22 22:20
嗨w大, 第一次看到expm1, 然後查一下还有logp1来增加精度的case, 我再参考, 谢谢!
26F：→ wohtp : f(x)拆成两截，x<0.5的时候用这个，另一头改用y = 1 07/22 22:22
27F：→ wohtp : -x 当变数，然後expm1比照办理 07/22 22:22
了解~
28F：推 wohtp : 伤浮点精度的是大数相减剩小数，不是0/0 07/22 22:26
f(x)/g(x)时, 假设数学上是lim_{x→0}f(x) = lim_{x→0}g(x) = 0
且lim_{x→0} f(x)/g(x) = 1
我就是怕x→0时, f跟g这两个函数在library上的精度不同, 或是函数本身递增速率不同
导致某个很小的x代入g(x)时, 如果电脑已经因为精度问题return true 0
但是f(x)还没, 那f(x)/g(x)就是+inf
反之如果是f(x)先return true 0, 那f(x)/g(x)就是+0
我就是怕以上这种情况才会问这篇, 还是说这个并不会发生?
29F：→ PPguest : 原来有expm1(笔记) 07/22 22:36
30F：→ PPguest : 想请问一下,若想知道实际误差的状况,准确的计算值要 07/22 22:37
31F：→ PPguest : 如何得到？ 07/22 22:37
32F：推 wohtp : 这个问题不就是「怎麽让数值运算更精确」 07/22 22:45
33F：→ wohtp : 开100位浮点数去算啊，不然还能怎样？ 07/22 22:46
34F：推 wohtp : 喔对了，回到原po的问题，x = 0 or 1的时候要直接让 07/22 22:52
35F：→ wohtp : 函数给极限值，所以其实需要四个case 07/22 22:52
我真实要应用的case是(1,0), (s,log_2(s)), (p,log_2(p)), (2,1)
然後1, s, p这三个可能都相等, 可能都相异, 可能两个相等
然後我去观察图形时哪个case, 都还是三次多项式
因此要拆的段数, 要算的极限值好多...
36F：推 wohtp : 所以你的函数最复杂就log？那你怕什麽？ 07/23 01:37
37F：→ wohtp : 这些常见的函数实作精度大概都逼近machine precisio 07/23 01:37
38F：→ wohtp : n。 07/23 01:37
39F：→ wohtp : 你给的例子里面，精度直落的唯一原因是大数相减，完 07/23 01:37
40F：→ wohtp : 全是程式员的锅。 07/23 01:37
41F：推 wohtp : 不然就像上面建议的，开bigfloat硬干吧。你的时间应 07/23 01:47
42F：→ wohtp : 该比cpu时间贵……吧？ 07/23 01:47
我应用的例子是不会有大数相减, 单纯就是一堆0+/0+
如果照w大你所说的这些函数的精度几乎没问题, 那我就放心了
另外问个浮点数"可能出问题"的问题, 今天假设我避开大数相减...等等会大幅丧失精度
的事情, 设计出一个C函数"float f(float x)"
然後我确定其
数学理论式中, x!=0的话f(x)就没问题
但是我还是怕浮点运算中即便x!=0, f(x)也可能会出问题(如这篇讨论的问题)
我除了
暴力跑IEEE754的所有denormal number然後print出来检查, 是不是没其他方式了
换言之, 有没有怎样的
数学式的coding style可以完全避开遇到这些精度问题?
问工程的朋友目前没有general的答案, 都要看函数长相
43F：推 LPH66 : 大减大消光这有个名字叫 catastrophic cancellation 07/23 02:42
44F：→ LPH66 : 一般都是只能建议改动式子不要有让相近数相减的地方 07/23 02:42
45F：→ LPH66 : 但变形方式如你所问到的都得看函数长相才知道 07/23 02:43
46F：→ LPH66 : 你的 0+/0+ 应该是踩到浮点数指数下限所以下溢了 07/23 02:45
47F：→ LPH66 : (denormal number 出现的原因就是缓和接近下限时的 07/23 02:45
48F：→ LPH66 : 下溢状况, 但它只对部份运算有些许缓和作用而已) 07/23 02:46
还有catastrophic cancellation这专有名词, 谢谢提供
然後02:42说的"不要有让相近数相减的地方", 是不论大数或是小数相减都应该避免?
大数相减误差大(浮点在大数精度很低), 小数相减又有下溢问题...
另外重复的问题是, 如果我自认为某个函数实作我都
尽量避免掉这些问题
但是我想要把信心度提高到100%, 就是扫整个浮点数的范围?
49F：推 wohtp : 我说大数小数是相对的啊。要说得更准确些是「相近数 07/23 09:44
50F：→ wohtp : 相减剩下比原数小得多的差」。 07/23 09:44
51F：→ wohtp : 拿你的f(x)为例，2^x的精度不会有问题，但是2^x - 1 07/23 09:44
52F：→ wohtp : 在零附近就烂掉了。出现这个组合是你自己的锅，不是 07/23 09:44
53F：→ wohtp : library精度不行。 07/23 09:44
这个我在吸收一下, 感觉我没搞懂
"单一数的精度问题" VS "相减结果的精度问题"
numpy文件举例exp(10^-10)-1跟expm1(10^-10)的差异确实差很多
但是我自己举例: float x = (1.[x_1~x_23])_2 * 2^0, where x_1~x_22=0, x_23=1
float y = (1.[y_1~y_23])_2 * 2^0, where y_1~y_23=0
其中()_2代表用二进位表示, []里面装着0 or 1的数
可以知道y就是数学上的1.0, x是float中比1.0大的下一个数, 也就是
最接近且大於y的
而观察z=x-y= (1.[z_1~z_23])_2 * 2^(-23), where z_1~z_23=0
可以知道相减後的结果仍是float, 并且
毫无误差
目前有点混乱, 不知道是我举例错误还是刚好举到特例, 还是根本跟expm1的案例不同类
而且如果自己刻比double更大的精度, exp(10^-10)-1的精度也会上升直逼expm1(10^-10)
所以感觉
绿色的VS不太能完全独立? 彼此有关系
54F：推 wohtp : 然後，要保证实作数值方法的品质，你当然对原函数的 07/23 09:58
55F：→ wohtp : 性质要有足够理解。解析方法手算再每个点都挑掉是唯 07/23 09:58
56F：→ wohtp : 一做法，就算像你说扫整个浮点数范围，意义其实也只 07/23 09:58
57F：→ wohtp : 是辅助判别有麻烦的值，最後你还是要自己把函数改写 07/23 09:58
58F：→ wohtp : 成没有精度问题的形式。 07/23 09:58
那可以归纳成 "尽量所有的浮点coding都避开相近数相减" 就好了?
比如我在实作浮点二次多项式P(x)=x^2-x, 我都直接打x*x-x
今天如果我没有要拿P(0+)去除以一个会接近0的分母, 我维持x*x-x就好
但是如果有要除法, 我就要改写x*x-x了?
59F：推 PPguest : 原po担心的0+,就像LPH66讲的,应是怕函数算一算结果 07/23 10:50
60F：→ PPguest : 超出浮点数能表示的范围 07/23 10:50
61F：→ PPguest : 双精度大约比5e-324小一半就会变0 07/23 10:54
62F：→ PPguest : 感觉就是注意实际运算中不会发生这种情况,如果会发 07/23 11:04
63F：→ PPguest : 生,就代表双精度不够用 07/23 11:04
64F：→ PPguest : 例如f(x)=x^2, x有需要用到例如1e-170那麽小吗？ 07/23 11:29
65F：→ PPguest : 若否,那没事;若有需要,那双精度就不够用 07/23 11:32
嗨P大, 如果单纯是f(x)=x^2并且没有f(0+)/0+
那麽f(0+)的精度我确实也觉得不重要(在我的应用中), f(0+)=1e-170 or 0.0都没差
只是目前我对於w大的例子" expm1 " & "2^x-1"
有感觉
w大又提"2^x-1"的问题并不是2^x的精度问题, 但是精度提升不是会更精准吗
而自己举的float x, y, z的相减例子又觉得没问题
目前有点混乱, 是不是要先固定某个精度後(比如float)才能谈相减的误差
有种...某个例子会这样看起来是某个问题, 但是另外一个例子看起来又不是那个问题XD
66F：→ PPguest : 主要是 07/22 22:26 後面你说的那种情况 07/23 16:21
67F：→ PPguest : f(x)=x^2, g(x)=x, x=1e-170 07/23 16:22
68F：→ PPguest : f(x)/g(x)=0, g(x)/f(x)=Inf 07/23 16:22
了解~
69F：推 LPH66 : Catastrophic cancellation 的问题只要你有一个中间 07/23 18:00
70F：→ LPH66 : 使用的精度, 那就有可能发生 07/23 18:00
71F：推 LPH66 : 例如 x=1e-10, 你直接求 2^x 再减 1 的话 07/23 18:05
72F：→ LPH66 : 2^x 得 1.00000000006931 (拿十五位十进位有效举例) 07/23 18:06
73F：→ LPH66 : 减去 1 之後得到 6.931e-11 这个数只有四位有效 07/23 18:06
74F：→ LPH66 : 而如果用其他方式去直接求得 2^x-1 的话 07/23 18:07
75F：→ LPH66 : 或许你能得到 6.93147180583968e-11 这样十五位有效 07/23 18:08
76F：→ LPH66 : 这个"其他方式"就是运用数学关系尝试找出没有这种减 07/23 18:09
77F：→ LPH66 : 的式子 (以 2^x-1 来说可能是泰勒展开後消掉 1) 07/23 18:09
78F：→ LPH66 : 上面这个例子你就算中间用了三十位有效 07/23 18:21
79F：→ LPH66 : 减 1 对消後的有效位数就是不足三十位 07/23 18:21
80F：→ LPH66 : 这就是它被叫做 catastrophic (灾难性) 对消的原因 07/23 18:22
L大这个举例我了解了, 我如此解释:
(1) 2^(1e-10)=1.00000000006931, 这并不是精度不足
这个值就是在固定精度(假设系统是十五位十进位的系统)下所求出来的最准确值了
理论值就是1.00000000006931abcdefg...
这也是为什麽w大说不是library 2^x的精度问题
(2) 2^(1e-10) - 1的理论值是0.00000000006931abcdefg
而这个理论值在系统是十五位十进位是会写成6.931abcdefg...*10^(-11)
也就是说, 这个6.931abcdefg...*10^(-11)这个值才是2^(1e-10) - 1
在这个系统的最准确值
但是今天如果照(1)的算法, 我们却得出6.931*10^(-11)这个不准确的值
这样理解应该是正确的吧?
---------------------------------------
另外请问各位一个实作问题, 我一直以来都在嵌入式系统写定点(int32)程式码
最近第一次有浮点可以用的IC, 写code的速度飙快, 因为之前算出的数学式就直接刻上去
不用像定点还要控制数值范围卡overflow之类的
结果最近这几个需要在意精度的浮点公式让我有种
浮点比定点更难计算与控制误差的感觉
所以想请板上浮点高手给些建议怎样的浮点coding style能少踩一点坑
比如之前说的少用相近数相减
81F：→ PPguest : 原po 07/23 09:44 之後红色那段的例子 07/23 22:45
82F：→ PPguest : 基本上都是10进位可以准确换算成2进位的情况 07/23 22:47
83F：→ PPguest : 例如0.5, 0.25, 0.125。像0.1用二进位的角度来看,在 07/23 22:49
84F：→ PPguest : 浮点数的系统只能用近似值来表达 07/23 22:49
这个我理解, 所以我09:44那段特地取刚好可以表达的浮点数x,y来减少讨论的变因
不过这样看起来会有catastrophic cancellation就是不会发生在这种刚好可以表达的
85F：→ PPguest : 很多数用二进位的角度来看,都是无穷小数,然後搭配 07/23 22:53
86F：→ PPguest : 原po对L大举例的理解,应该就没问题了 07/23 22:54
了解~
87F：→ PPguest : 用比double更大的精度, exp(10^-10)-1 大概会看到更 07/23 22:58
88F：→ PPguest : 多小数,但大概有一定比例的小数其实是不准的 07/23 22:59
如果发生这种事可以先归类在函数实作不精准吼?
前提应该是基於在某个精度系统下, 函数的实作准度是准到尾数的bits数
89F：→ PPguest : 不过小数准的位数应该会变多 07/23 23:01
90F：→ PPguest : "是不是要先固定某个精度......才能谈相减的误差" Y 07/23 23:05
了解~
91F：推 PPguest : 多项式用"nested"的方式求值除了降低计算量,好像也 07/23 23:10
92F：→ PPguest : 算是减少误差 07/23 23:11
93F：推 PPguest : 07/23 10:50 那段，原po实际应用的情况大概不用担心 07/23 23:17
94F：→ PPguest : 会发生.s和p在1,2之间,如果不同,再怎麽近大概是 07/23 23:20
95F：→ PPguest : 2e-16等级的差距,要超出浮点数可表示的范围应该不太 07/23 23:22
96F：→ PPguest : 可能 07/23 23:22
谢谢idea, 放心不少~
97F：→ PPguest : 07/23 22:59"有一定比例的小数...不准" 我用词不好 07/24 11:15
98F：→ PPguest : 单纯是相近的数相减造成不精准 07/24 11:16
99F：→ PPguest : 如果原本有15位有效数字,相减後剩5位有效,後面10位 07/24 11:19
100F：→ PPguest : 不准.用比double更大的精度,计算前例如有30位有效数 07/24 11:20
101F：→ PPguest : 字,相减後还剩20位有效,後面10位不准 07/24 11:21
102F：→ PPguest : 用expm1,可能输入是15位有效数字,输出是14位有效, 07/24 11:23
103F：→ PPguest : 1位不准;用比double更大的精度,输入30位有效数字, 07/24 11:24
104F：→ PPguest : 输出29位有效,1位不准.大概是这种感觉 07/24 11:24
这个举例理解~谢谢P大~
105F：推 LPH66 : 初次接触浮点数运算的话冼镜光老师这篇可以看看 07/24 14:33
106F：→ LPH66 : blog.dcview.com/article.php?a=Az0HYgNrBDU%3D 07/24 14:33
107F：→ LPH66 : 包含了我们这里提过的跟一些其他也很重要的运算概念 07/24 14:34
谢谢L大的reference~
※ 编辑: znmkhxrw (59.102.225.191 台湾), 07/25/2023 01:53:08