31

I think the (first order) right thing to do is look at the ratio of flops to bytes needed in the algorithm, which I call $\beta$. Let $F_{\mathrm{max}}$ be the maximum flop rate of the processor, and $B_{\mathrm{max}}$ the maximum bandwidth. If $\frac{F_{\mathrm{max}}}{\beta} > B_{\mathrm{max}}$, then the algorithm will be bandwidth limited. If $B_{\...


27

There is a wide variety of algorithms; Barnes Hut is a popular $\mathcal{O}(N \log N)$ method, and the Fast Multipole Method is a much more sophisticated $\mathcal{O}(N)$ alternative. Both methods make use of a tree data structure where nodes essentially only interact with their nearest neighbors at each level of the tree; you can think of splitting the ...


26

Are these integers or floating point numbers? Assuming it's floating point, I would go with the first option. It's better to add the smaller numbers to each other, then add the bigger numbers later. With the second option, you'll end up adding a small number to a big number as i increases, which can lead to problems. Here's a good resource on floating point ...


24

animal_magic's answer is correct that you should add the numbers from smallest to largest, however I want to give an example to show why. Assume we are working in a floating point format that gives us a staggering 3 digits of accuracy. Now we want to add ten numbers: [1000, 1, 1, 1, 1, 1, 1, 1, 1, 1] Of course the exact answer is 1009, but we can't get ...


22

I don't see why one has to be the "winner"; this isn't a zero-sum game, where flop counts and memory access have to drown the other out. You can teach both of them, and I think they both have their uses. After all, it's hard to say that your $O(N^4)$ algorithm with $O(N)$ memory accesses is necessarily going to be faster than your $O(N \log N)$ algorithm ...


20

$n \log\log n$ is between $n$ and $n \log n$, and is a relatively common one to find in the wild.


19

The closest positive answers to your question that I could find is for sparse diagonal perturbations (see below). With that said, I do not know of any algorithms for the general case, though there is a generalization of the technique you mentioned for scalar shifts from SPD matrices to all square matrices: Given any square matrix $A$, there exists a Schur ...


19

I will try to answer your question considering that you are asking for Python specifically. I will describe my own method of tackling a simulation problem. Strategies for faster simulations are given in this description. First, I prototype new simulations in Python. Of course, I try to make use of NumPy and SciPy as much as I can. Whereas NumPy provides a ...


18

Adding arbitrary floating point numbers will usually give some rounding error, and the rounding error will be proportional to the size of the result. If you calculate a single sum and start by adding the largest numbers first, the average result will be larger. So you would start adding with the smallest numbers. But you get better result (and it runs ...


18

(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 $\mathbf{v}_1\cdot \mathbf{v}_2$ can suffer from catastrophic cancellation if they are almost parallel/perpendicular—atan2 can't give you good accuracy if ...


16

Encapsulation and data hiding are extremely important for extensible libraries in scientific computing. Consider matrices and linear solvers as two examples. A user just needs to know that an operator is linear, but it may have internal structure such as sparsity, a kernel, a hierarchical representation, a tensor product, or a Schur complement. In all cases, ...


16

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} = [o_s, o_e]$. We can find the intersection $I_{o} = I_{a} \bigcap I_{b}$ doing the following: if ( $b_s \gt a_e$ or $a_s \gt b_e$ ) { return $\emptyset$ } ...


15

For the general case, I'd use compensated summation (or Kahan summation). Unless the numbers are already sorted, sorting them will be much more expensive than adding them. Compensated summation is also more accurate than sorted summation or naive summation (see the previous link). As for references, What every programmer should know about floating-point ...


15

Quasi-Newton methods construct an approximate Hessian for an arbitrary smooth objective function $f(x)$ using values of $\nabla f$ evaluated at the current and previous points. At each iteration of the method the quasi-Newton approximate Hessian is updated using the gradient evaluated at the latest iterate, $x^{(k)}$. These approximate Hessians aren't ...


14

It is cumbersome for the user to specify every aspect of an algorithm. If the algorithm allows nested components, then no finite number of options would be sufficient. Therefore, it is critical that options do not necessarily "bubble up" to the top level, as in the case of explicit arguments or template parameters. This is sometimes called the "configuration ...


14

Well, everything is a trade-off of some sort or another. For random number generators, I group them into 3 basic categories: Good enough for homework. Good enough to bet your company on. Good enough to bet your country on. Linear congruential PRNGs (the method generally implemented in most libraries) are solidly in category 1. Both Fortuna and ...


14

It depends a lot on the size of your matrix, in the large-scale case also on whether it is sparse, and on the accuracy you want to achieve. If your matrix is too large to allow a single factorization, and you need high accuracy, the Lanczsos algorithm is probably the fastest way. In the nonsymmetric case, the Arnoldi algorithm is needed, which is ...


13

Look at the fast multipole method. It is highly scalable and $O(n)$. It allows trading off between precision and cost. Here's an example where it is run at 42 Tflops on a GPU cluster.


13

It's a rather rough algorithm, but I'd use the following procedure for a crude estimate: if, as you say, the purported $f(x)$ that represents your $(x_i,y_i)$ is already almost linear as $x$ increases, what I'd do is to take differences $\dfrac{y_{i+1}-y_i}{x_{i+1}-x_i}$, and then use an extrapolation algorithm like the Shanks transformation to estimate the ...


13

You look at how deal.II (http://www.dealii.org/) does it -- there, dimension independence lies at the very heart of the library, and is modeled as a template argument to most data types. See, for example, the dimension-agnostic Laplace solver in the step-4 tutorial program: http://www.dealii.org/developer/doxygen/deal.II/step_4.html See also https://...


13

Having done numerical software for 15 years, I can unambiguously state the following: Encapsulation is important. You do not want to pass around pointers to data (as you suggest) since it exposes the data storage scheme. If you expose the storage scheme, you can never change it again because you will access the data all over the entire program. The only way ...


13

A fairly simple method would be to choose a basis in function space and convert the integral transformation to a matrix. Then you can just invert the matrix. Mathematically, here's how that works: you need some set of orthonormal basis functions $T_i(x)$. (You can get away without them being normalized too, but it's easier to explain this way.) Orthonormal ...


12

If $(D_{i} + A)$ is diagonally dominant for each $i$, then recent work by Koutis, Miller, and Peng (see Koutis' website for work on symmetric diagonally dominant matrices) could be used to solve each system in $\mathcal{O}(n^2 \log(n))$ time (actually $\mathcal{O}(m\log(n))$ time, where $m$ is the maximum number of nonzero entries in $(D_{i} + A)$ over all $...


12

The fastest method will likely depend upon the spectrum and normality of your matrix, but in all cases Krylov algorithms should be strictly better than power iteration. G.W. Stewart has a nice discussion of this issue in Chapter 4, Section 3 of Matrix Algorithms, Volume II: Eigensystems: The power method is based on the observation that if $A$ has a ...


12

The question highlights that most "plain" programming languages (C, Fortran, at least) do not allow you to do this cleanly. An added constraint is that you want notational convenience and good performance. Therefore, instead of writing a dimension-specific code, consider writing a code that generates a dimension-specific code. This generator is dimension-...


12

Knowing how other CAS do this might help you. To my knowledge, Mathematica uses a variation on the following basic algorithm for plotting a one-variable function $f(x)$ or a parametric curve $(x(t), y(t))$ (I'm going to assume $f(x)$ for this description). Start with a regularly spaced grid of points on the plotting domain. (In Mathematica, there's a ...


12

Yes, but the answer is a little more complicated than you were probably hoping for. In short: one color in general corresponds to a distribution of different frequencies, or wavelengths. You have three types of cone cells in your eyes [+ one type of rod cell]. Color-blind people have less than three types (or their response curves are closer together). ...


12

Adams-Moulton method is significantly more stable. The analogy used when I was taught the difference is the same as extrapolation and interpolation. Interpolation is relatively safe numerically. Extrapolation can blow up if you happen to have an asymptote or some other odd feature. For instance, solving the ode $y'(t) = -y(t)$ with $y(0) = 1$ using ...


11

I've faced this problem several times when developing my own simulation codes from scratch: which parameters should go in an input file, which should be taken from the command line, etc. After some experimenting, the following turned out to be efficient. (It is not as advanced as PETSc.) Instead of writing an experimental simulation 'program', I'm more ...


11

Boost Graph Library and LEMON As Daniel mentions in his comprehensive answer, the most full-featured general C++ library is the Boost Graph Library. There is a new distributed-memory extension capable of doing some basic algorithms such as breadth-first and depth-first search, minimum spanning trees, and connected components search, but I am not very ...


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