Skip to main content
51 votes
Accepted

The easiest way to find intersection of two intervals

We can define a solution to this problem in the following way. Assume the input intervals can be defined as $I_{a} = [a_s, a_e]$ and $I_{b} = [b_s, b_e]$, while the output interval is defined as $I_{o}...
spektr's user avatar
  • 4,278
36 votes

stupid + stupid = brilliant in scientific computing

In W. Kahan, "Interval arithmetic options in the proposed IEEE floating point arithmetic standard". In: Karl L. E. Nickel (ed.), Interval Mathematics 1980, New York: Academic Press 1980, pp. ...
njuffa's user avatar
  • 1,895
28 votes
Accepted

How do compression algorithms compress data so fast?

The short answer to your question is this: If your goal is speed (as it is in typical applications of data compression), then (i) you need to choose a programming language that allows you to write ...
Wolfgang Bangerth's user avatar
28 votes

stupid + stupid = brilliant in scientific computing

Consider the constrained optimization problem $$\min_x f(x) \quad\text{s.t. } g(x) = 0$$ where, to make things nice, we'll assume $f$ is convex and $g$ has convex level sets. There are two bad ways to ...
25 votes
Accepted

Numerically stable way of computing angles between vectors

(I have tested this approach before, and I remember it worked correctly, but I haven't tested it specifically for this question.) As far as I can tell, both $\|\mathbf{v}_1\times \mathbf{v}_2\|$ and $...
Kirill's user avatar
  • 11.4k
25 votes

More stable algorithm to calculate `sqrt(a^2 + b^2) - abs(a)` in MatLab

A manipulation that may help is the following. Assume for simplicity $a>0$. We have the identity $$b^2 = (a^2+b^2)-a^2 = (\sqrt{a^2+b^2}-a)(\sqrt{a^2+b^2}+a),$$ hence $$ \sqrt{a^2+b^2}-a = \frac{b^...
Federico Poloni's user avatar
23 votes
Accepted

How do I find the minimum-area ellipse that encloses a set of points?

Theory The 1997 paper "Smallest Enclosing Ellipses -- Fast and Exact" by Gärtner and Schönherr addresses this question. The same authors provide a C++ implementation in their 1998 paper &...
Richard's user avatar
  • 3,981
19 votes

stupid + stupid = brilliant in scientific computing

Disclaimer: stupid is in the eyes of the beholder. For decades, it seemed like what would now be known as deep learning was combining stupid with stupid. Stupid 1: over parameterization. A rule of ...
Cliff AB's user avatar
  • 291
18 votes
Accepted

More stable algorithm to calculate `sqrt(a^2 + b^2) - abs(a)` in MatLab

You can break down the domain of your function into three distinct cases: $|a|\gg |b|$: In this case, $\sqrt{a^2+b^2} \approx |a|$ and a naive application of the formula will likely result in poor ...
Wolfgang Bangerth's user avatar
13 votes

Numerically stable way of computing angles between vectors

The efficient answer to this question is, not too surprisingly, in another note by Velvel Kahan: $$\alpha=2\arctan\left(\left\|\frac{\mathbf v_1}{\|\mathbf v_1\|}+\frac{\mathbf v_2}{\|\mathbf v_2\|}\...
J. M.'s user avatar
  • 3,135
13 votes
Accepted

Is there an algorithm or graph theory that allows me to not need to store an intermediate matrix when calculating AT*Y1*A + BT*Y2*B?

BLAS may not have a function to compute what you are asking for, but the product $$ Y_N = A^TY_AA + B^T Y_B B $$ means that the entries $(Y_N)_{ij}$ are defined by $$ (Y_N)_{ij} = \sum_{k,l} (A^T)...
Wolfgang Bangerth's user avatar
12 votes
Accepted

Accurate and efficient computation of the inverse Langevin function

The inverse Langevin function $\mathcal{L}^{-1}(x)$ is an odd function. Therefore one needs to consider only approximation on the interval $[0, 1]$; the negative half-plane is treated via symmetry ...
njuffa's user avatar
  • 1,895
11 votes
Accepted

Markov (Chain) image generators?

I've implemented this recently, basically it counts how many times each specific colour borders another colour to make up a frequency table. To generate an image, a random colour and position are ...
Jonno_FTW's user avatar
  • 226
11 votes

Cheap recalculation of eigenvalues and eigenvectors for a low-rank update of the matrix

Unfortunately, I don't think there is a good algorithm to do this efficiently. Given the eigendecomposition $\mathbf A = \mathbf X \mathbf D \mathbf X^T$, one is tempted to project $\mathbf v$ onto ...
rchilton1980's user avatar
  • 4,946
11 votes

How do compression algorithms compress data so fast?

I have written compression software in Rust. The answer is not simple, but compression algorithms are usually designed to be reasonably quick. RFC 1951 has two steps, the first stage is to find ...
George Barwood's user avatar
10 votes

How do I find the minimum-area ellipse that encloses a set of points?

With your formulation of the constraints, you can only produce ellipses where the semi-major and semi-minor axes are aligned with the coordinate axes, but it's clear from the figure that you attached ...
Daniel Shapero's user avatar
9 votes

Should benchmarkings be done at all? What is the point?

Yes, benchmarking should be done. I make this claim as an Editor, Author, and Reviewer. Below, I represent these roles' stances slightly hyperbolically. But let me strongman your argument first. In ...
Richard's user avatar
  • 3,981
9 votes

Is it really necessary to solve a system of linear equations in the Finite Element Method?

I think your question is actually pretty fundamental and deserves a thoughtful answer. Paraphrasing a bit, your question is perhaps motivated by the observation that engineering design is often ...
rchilton1980's user avatar
  • 4,946
9 votes
Accepted

Non-uniform Gaussian spaced vector

What you need is the Gaussian quantile function. The quantile function outputs the value of a random variable such that its probability is less than or equal to an input probability value. The ...
Scientist's user avatar
  • 206
8 votes
Accepted

Methods for solving $Ax=b$, small and sparse A

By direct substitution, trivially. $$ \begin{bmatrix} 0\\f \end{bmatrix} = \begin{bmatrix} M & -I\\ I & 0 \end{bmatrix} \begin{bmatrix} y\\z \end{bmatrix} = \begin{bmatrix} My-z\\ y \end{...
Federico Poloni's user avatar
8 votes

Could we train an AI to find (only) Mersenne primes and beat the current record?

The fundamental difficulty with prime numbers is that for all practice purposes, they are randomly distributed. Most theorems (and open problems) about prime numbers -- say, that there are infinitely ...
Wolfgang Bangerth's user avatar
7 votes
Accepted

Applying the result of Cuthill-McKee in SciPy

The reverse Cuthill-McKee algorithm produces a reordering that applies to both the rows and columns. This is because it works by considering matrices as graphs of (undirected) connected nodes. ...
Nick C.'s user avatar
  • 188
7 votes
Accepted

Eigenvector with maximum overlap

The following paper suggests that the Jacobi-Davidson method can be used to target eigenvectors based on "any property that can be computed from the eigenvector", which would seem to include overlap ...
deemaregee's user avatar
7 votes
Accepted

MD Simulation: Reference for the Neighbor's List Method

I'd recommend "The Art of Molecular Dynamics Simulation" by D. C. Rapaport. The code samples are written in C. I'm not a huge fan of the programming style of the book, but at least it's not ...
lr1985's user avatar
  • 687
7 votes
Accepted

How to justify using available code (in different language) for comparing algorithms

So, you are comparing a generally slower Matlab implementation of algorithm A to a generally faster C++ implementation of algorithm B, and still getting the advantage for A. I would say, ...
Anton Menshov's user avatar
  • 8,702
7 votes
Accepted

Fast algorithm for computing cofactor matrix

So, a cofactor matrix is a transpose of an adjugate matrix. I know of the following paper: G. W. Stewart, "On the adjugate matrix," Lin. Alg. Appl., vol. 283, no. 1–3, pp. 151–164, Nov. ...
Anton Menshov's user avatar
  • 8,702
7 votes

Should benchmarkings be done at all? What is the point?

Benchmarks are useful, but no benchmark tells the whole story. There are many useful benchmarks. For example, the Julia microbenchmarks are an interesting case of an isolating benchmark: it tries to ...
Chris Rackauckas's user avatar
7 votes
Accepted

Solving the time dependent Schrödinger equation with leapfrog integration in 1D

I read through the paper you linked and they give the stability condition for this method to be (eq. A6) $$ \frac{-2}{\Delta t} \le V \le \frac{2}{\Delta t} - \frac{2}{m \Delta r^2} $$ This has to be ...
helloworld922's user avatar
7 votes

How do compression algorithms compress data so fast?

There are actually two questions that can be answered: Why is decompression so much faster than compression? The key point to understand here is, that compression typically involves searching for ...
cmaster - reinstate monica's user avatar

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