20 votes
Accepted

Inequality caused by float inaccuracy

In typical floating point implementations, the result of a single operation is produced as if the operation was performed with infinite precision, and then rounded to the nearest floating-point number....
user avatar
  • 25.2k
17 votes

How to prove that matrix inversion is at least as hard as matrix multiplication?

If you want to multiply two matrices $A$ and $B$ then observe that $$\begin{pmatrix}I_n&A&\\&I_n&B\\&&I_n\end{pmatrix}^{-1}= \begin{pmatrix}I_n&-A&AB\\&I_n&-B\\&...
user avatar
  • 1,116
14 votes

Implementation of Naive Bayes

The usual trick to avoid this underflow is to compute with logarithms, using the identity $$ \log \prod_{i=1}^n p_i = \sum_{i=1}^n \log p_i. $$ That is, instead of using probabilities, you use their ...
user avatar
13 votes

Multi-point evaluations of a polynomial mod p

No, $O(n \lg q)$ running time is not achievable. It takes $\Omega(q)$ space even just to write out the answer, so any algorithm will necessarily have running time at least $\Omega(q)$. However, you ...
user avatar
  • 141k
13 votes
Accepted

Fastest way to solve a system of linear equations

A LU decomposition of a $n \times n$ matrix can be computed in $O(M(n))$ time, where $M(n)$ is the time to multiply two $n \times n$ matrices. Therefore, you can find a solution to a system of $n$ ...
user avatar
  • 141k
9 votes
Accepted

Where can I find an original reference for this integer square root algorithm

The wikipedia article on Methods of computing square roots: base 2 presents a strikingly similar snippet of C-code [a], but the link to the source is dead. Let's try to do better. The snippets from ...
8 votes
Accepted

Program transformations for numeric stability

There actually is some research on improving the numerical stability of floating point expressions, the Herbie project. Herbie is a tool to automatically improve the accuracy of floating point ...
user avatar
8 votes

Arbitrary precision integer square root algorithm?

One of the problems with Newton's method is that it requires a division operation in each iteration, which is the slowest basic integer operation. Newton's method for the reciprocal square root, ...
user avatar
  • 18.9k
7 votes

Is everything in CS either a numeric method or a symbolic method?

No. You could make a decent argument that every algorithm is in some sense a symbolic manipulation, though not in the sense normally used for symbolic computation - they are not explicitly algebraic, ...
user avatar
7 votes
Accepted

Project to nearest point in convex polytope

A quadratic program is an optimization problem where the goal is to minimize $y^T Q y + c^T y$ subject to $A y \leq b$. If $Q$ is positive definite, then this is a convex quadratic program and we can ...
user avatar
  • 116
6 votes

Why are transcendental functions of large numbers inaccurate on computers?

There's nothing fundamentally hard about computing $\sin(10^{99})$. You simply compute $x = 10^{99} \bmod 2\pi$, then compute $\sin(x)$. (Why is this valid? It's because $\sin(x)=\sin(y)$ if $x\...
user avatar
  • 141k
6 votes
Accepted

The stability of log(1+x)

Consider the case that x is small. (1 + x) has a rounding error; the result that you get is not (1 + x) but (1 + x') for some x' close to x. If x is very small, the relative difference between x' and ...
user avatar
  • 25.2k
6 votes

Why aren't computables used for numerical calculations?

It sounds extremely inefficient compared to floating point. We have a very good understanding of how to control the errors in floating-point calculations (e.g., adding small numbers before large ones, ...
user avatar
6 votes

Why is adding log probabilities considered "numerically stable"?

A "numerically stable" method calculates a result in a way that will not produce excessive rounding errors. Given a small number x, let's say x = 0.00079, it is possible to calculate log (1 + x) ...
user avatar
  • 25.2k
6 votes

fast and stable x * tanh(log1pexp(x)) computation

With some algebraic manipulation (as pointed out in @orlp's answer), we can deduce the following: $$f(x) = x \tanh(\log(1+e^x)) \tag{1}$$ $$ = x\frac{(1+e^x)^2 - 1}{(1+e^x)^2 + 1} = x\frac{e^{2x} + 2e^...
user avatar
  • 265
5 votes
Accepted

Monetary computations theory (manual/textbook)

Let me illustrate one problem which could happen, and one way to solve it. You want to distribute a given amount of cents $N$ into $k$ piles, in proportions $p_1,\ldots,p_k$, where $p_1,\ldots,p_k \...
user avatar
5 votes

Do recursive algorithms generally perform better than their for-loop counterpart?

The answer will depend on the compiler. As @vonbrand wrote, "Given a good enough compiler, you might even get the very same object code." In particular, good compilers will do tail-call elimination. ...
user avatar
  • 141k
5 votes

Calculating decimal digits of pi, using something similar to a Bailey–Borwein–Plouffe formula

The Bailey–Borwein–Plouffe formula only works in hexadecimal. There might be other formulas for other bases, but I'm not aware of a decimal-based formula. If you want to obtain the $N$th decimal digit,...
user avatar
5 votes
Accepted

Arbitrary precision integer square root algorithm?

You can use Newton's method or any of a number of other methods for finding approximations to roots of the polynomial $p(x) = x^2 -c$. The rate of convergence for Newton's method will be quadratic, ...
user avatar
  • 141k
5 votes
Accepted

Computing Von Neumann Entropy Efficiently

A paper Computing the Entropy of a Large Matrix by Thomas P. Wihler, Bänz Bessire, André Stefanov suggests approximating $x \lg x$ with a polynomial. Then you can use the trace of powers of the matrix ...
user avatar
  • 5,722
5 votes
Accepted

Numerical stability of linear interpolation

"Numerical stability" is a much vaguer term than most people realise. We typically use it when referring to an approximation method, such as some kind of linear analysis, or numeric quadrature, or ...
user avatar
  • 18.9k
5 votes
Accepted

Are there parallel matrix exponentiation algorithms that are more efficient than sequential multiplication?

If you have multiple processors that can work in parallel, then you can calculate any power up to the power (2^k) in k steps. For example: To calculate $M^{15}$, you calculate: Stage 1: Calculate $M^...
user avatar
  • 25.2k
5 votes
Accepted

Check if a given polynomial is primitive

In order to check that a degree $n$ polynomial $P$ over $GF(2)$ is primitive, you first need to know the factorization of $2^n-1$ (you can look it up in tables, or use a CAS). Then, you test that $x^{...
user avatar
5 votes

Inequality caused by float inaccuracy

Java uses IEEE 754 binary floating point representation, which dedicates 23 binary digits to the mantissa, that is normalized to begin with the first significant digit (omitted, to save space). $0....
user avatar
5 votes

Inequality caused by float inaccuracy

The binary floating point format supported by computers is essentially similar to decimal scientific notation used by humans. A floating-point number consists of a sign, mantissa (fixed width), and ...
user avatar
  • 842
5 votes

Fastest way to solve a system of linear equations

There is what you want to achieve, and there is reality, and sometimes they are in conflict. First you check if your problem is a special case that can be solved quicker, for example a sparse matrix. ...
user avatar
  • 25.2k
5 votes
Accepted

numerically stable log1pexp calculation

Let $0 < \varepsilon \lll 1$ be the relative error bound of the floating-point system—$2^{-53}$ in IEEE 754 binary64 arithmetic. First, the naive formula ...
user avatar
4 votes

Why are transcendental functions of large numbers inaccurate on computers?

Taking the sine of large numbers is a numerically unstable operation. Considering an argument like $10^{99}$, you can get a completely different value of the sine by adding, say $1$ to it. Think that ...
user avatar
  • 3,912
4 votes
Accepted

Efficient algorithm to compute the minimum of multiple piecewise linear functions

This is basically an instance of the line segment intersection problem. One standard approach is to use a sweep line algorithm. For instance, the Bentley-Ottman algorithm would be a reasonable ...
user avatar
  • 141k
4 votes
Accepted

Approximate a float using a minimal fraction

The partial convergents of the continued fraction of $x$ consists of all the best rational approximations of $x$; see Wikipedia, for example. A best rational approximation of $x$ is a rational number $...
user avatar

Only top scored, non community-wiki answers of a minimum length are eligible