In the mathematical subfield of numerical analysis, numerical stability is a generally desirable property of numerical algorithms. The precise definition of stability depends on the context. One is numerical linear algebra and the other is algorithms for solving ordinary and partial differential equations by discrete approximation.
Questions tagged [numerical-stability]
104 questions
58
votes
2 answers
Numerically stable way to compute sqrt((b²*c²) / (1-c²)) for c in [-1, 1]
For some real value b and c in [-1, 1], I need to compute
sqrt( (b²*c²) / (1-c²) ) = (|b|*|c|) / sqrt((1-c)*(1+c))
Catastrophic cancellation appears in the denominator when c approaches 1 or -1. The square root probably also does not help.
I was…
jeff
- 623
- 5
- 6
24
votes
4 answers
In Python small floats tending to zero
I have a Bayesian Classifier programmed in Python, the problem is that when I multiply the features probabilities I get VERY small float values like 2.5e-320 or something like that, and suddenly it turns into 0.0. The 0.0 is obviously of no use to…
Pravel
- 295
- 1
- 2
- 8
24
votes
1 answer
Weird numpy.sum behavior when adding zeros
I understand how mathematically-equivalent arithmentic operations can result in different results due to numerical errors (e.g. summing floats in different orders).
However, it surprises me that adding zeros to sum can change the result. I thought…
shx2
- 61,779
- 13
- 130
- 153
22
votes
1 answer
Why does chisq.test sort data in descending order before summation
Why does chisq.test function in R sorts data before summation in descending order?
The code in question is:
STATISTIC <- sum(sort((x - E)^2/E, decreasing = TRUE))
If I was worried about numerical stability due to usage of float arithmetic and…
user824276
- 617
- 1
- 7
- 20
18
votes
5 answers
Robust atan(y,x) on GLSL for converting XY coordinate to angle
In GLSL (specifically 3.00 that I'm using), there are two versions of
atan(): atan(y_over_x) can only return angles between -PI/2, PI/2, while atan(y/x) can take all 4 quadrants into account so the angle range covers everything from -PI, PI, much…
HuaTham
- 7,486
- 5
- 31
- 50
12
votes
5 answers
Quaternions and numerical stability
I'm learning about unit quaternions and how to use them to represent and compose rotations. Wikipedia says they are more numerically stable than matrix representations, but doesn't give a reference. Can anyone explain to me (preferably with some…
rmp251
- 5,018
- 4
- 34
- 46
12
votes
1 answer
c++: strategies for stability of floating point arithmetic
Can anyone recommend any C++ libraries/routines/packages that contain strategies for maintaining the stability of various floating point operations?
Example: suppose you would like to sum across a vector/array of one million long double in the unit…
cmo
- 3,762
- 4
- 36
- 64
10
votes
4 answers
How do I check and handle numbers very close to zero
I have some math (in C++) which seems to be generating some very small, near zero, numbers (I suspect the trig function calls are my real problem), but I'd like to detect these cases so that I can study them in more detail.
I'm currently trying out…
Petriborg
- 2,940
- 3
- 28
- 49
9
votes
1 answer
Java code optimization leads to numerical inaccuracies and errors
I'm trying to implement a version of the Fuzzy C-Means algorithm in Java and I'm trying to do some optimization by computing just once everything that can be computed just once.
This is an iterative algorithm and regarding the updating of a matrix,…
rano
- 5,616
- 4
- 40
- 66
8
votes
2 answers
Numerically Stable Implementation
I need to compute a normalized exponential of a vector in Matlab.
Simply writing
res = exp(V)/sum(exp(V))
overflows in an element of V is greater than log(realmax) = 709.7827.
(I am not sure about underflow conditions.)
How should I implement it…
user25004
- 1,868
- 1
- 22
- 47
7
votes
5 answers
numerically stable inverse of a 2x2 matrix
In a numerical solver I am working on in C, I need to invert a 2x2 matrix and it then gets multiplied on the right side by another matrix:
C = B . inv(A)
I have been using the following definition of an inverted 2x2 matrix:
a = A[0][0];
b =…
Steve
- 8,153
- 9
- 44
- 91
7
votes
1 answer
What are the specific reasons for CVXPY to throw `SolverError` exception?
I am using CVXPY (version 1.0) to solve a quadratic program (QP) and I often get this exception:
SolverError: Solver 'xxx' failed. Try another solver.
which makes my program really fragile. I have tried different solvers, including CVXOPT, OSQP,…
Tony
- 445
- 6
- 13
7
votes
4 answers
Logsoftmax stability
I know how to make softmax stable by adding to element -max _i x_i. This avoids overflow and underflow.
Now, taking log of this can cause underflow. log softmax(x) can evaluate to zero, leading to -infinity.
I am not sure how to fix it. I know this…
Abhishek Bhatia
- 9,404
- 26
- 87
- 142
7
votes
2 answers
How do I implement a numerically stable weighted logaddexp?
What is the most numerically stable way of calculating:
log[(wx * exp(x) + wy * exp_y)/(wx + wy)]
where the weights wx, wy > 0?
Without the weights, this function is logaddexp and could be implemented in Python with NumPy as:
tmp = x - y
return…
Neil G
- 32,138
- 39
- 156
- 257
7
votes
1 answer
Power of number close to 1
I'm guessing there is some standard trick that I wasn't able to find: Anyway I want to compute a large power of a number very close to 1(think 1-p where p<1e-17) in a numerically stable fashion. 1-p is truncated to 1 on my system.
Using the taylor…
Tobias Madsen
- 857
- 8
- 10