Hot answers tagged

65

Except for code which does a significant number of floating-point operations on data that are held in cache, most floating-point intensive code is performance limited by memory bandwidth and cache capacity rather than by flops. $v$ and the products $Av$ and $Bv$ are all vectors of length 2000 (16K bytes in double precision), which will easily fit into a ...


21

Your code is limited by memory bandwidth. For trivial math, it's often better to count memory accesses rather than flops. You'll get the following table: operation memory reads/writes matrix + matrix 3n² matrix * vector 2n²+n (if vector is not cached) matrix * vector n²+2n (if vector is only read once) vector + vector ...


20

Going from MATLAB to Python does introduce quite a bit of syntax overhead. One way to quantify it is the nice QuantEcon cheatsheet which showcases how there's a lot of extra "stuff" going on when trying to write simple linear algebra commands in Python. The verbose NumPy syntax is really just a symptom of how it was not developed as a technical ...


18

The most important aspect of interpolation and curve fitting is to understand why high order polynomial fits can be an issue and what the other options are and then you can understand when they are/are not a good choice. A few issues with high order polynomials: Polynomials are naturally oscillatory functions. As the order of the polynomial increases, the ...


18

Short answer, you want to have the leftmost index on the innermost loop. In your example, the loop indices would go k, j, i and the array indices would be i, j, k. This has to do with how MATLAB stores the different dimensions in memory. For more, see #13 of this reddit post.


16

Note that $\pi/2$ is represented in double precision format in a way that is not exactly equal to $\pi/2$. It's only accurate to about 15 digits. Thus you're starting every so slightly away from the equilibrium position. Since the equilibrium is unstable, it will eventually start moving.


15

The function error('error message'); will exit your program and print the error message to the console.


15

The question has two very different subquestions. I will address the first one only. Matlab's version runs on average 24 times faster than my python equivalent! The second one is subjective. I would say that letting know the user that there is some problem with the integral is a good thing and this SciPy behavior outperforms the Matlab`s one to keep it ...


15

Computing the term $e^x$ is definitely significantly more expensive than computing a lower-order polynomial -- say $x^4$. But it may be ten to 100 times more expensive at most, not "crazy" expensive. So I suspect that if Matlab takes forever to compute something, then that's because the character of your ODE changes significantly. For example, the presence ...


15

I think the two main points have already been made by Brian and Ertxiem: your initial value is an unstable equilibrium and the fact that your numerical computations are never really exact provides the small perturbation that will make the instability kick in. To give a bit more detail how this plays out, consider your problem in the form of a general ...


14

Here is R1, as computed in MATLAB: 1.0e+07 * -7.382605957465515 -9.599867106092937 -2.830412177259742 -0.000000000002830 -0.000000000002830 -1.230434326244253 -1.599977851015490 -0.471735362876624 -0.000000000000472 -0.000000000000472 3.691302978732758 4.799933553046468 1.415206088629871 0.000000000001415 0.000000000001415 -5....


14

First, see Mark L. Stone's answers, which is completely correct. Second, realize that this is the reason why people told you to use relative errors in your numerical analysis class. :) Third, the real question here is why the results do not coincide exactly, since both languages call some BLAS library functions for their computations. There are several very ...


13

There are libraries that you can use in Python that will give you all (or at least nearly all) of the functionality of MATLAB. For example, scipy.integrate.solve_ivp() supports a number of methods for ODE integration that are comparable to what you can get with the various odexxx() functions in MATLAB. So no, you wouldn't have to write your own ODE ...


11

MATLAB has a couple of "exact" functions for this, cond and rcond, with the latter returning a reciprocal of the condition number. Matlab approximate function condest is more fully described below. Often estimates of the condition number are generated as by-products of the solution of a linear system for the matrix, so you might be able to piggyback the ...


11

The cost of standard operations like find-and-replace, maxtrix-matrix multiplications and elementwise matrix-matrix multiplications is well-known: Check out the corresponding wiki site. If you'd really like to know what your MATLAB installation does, you can always go ahead and measure the time, e.g., t = zeros(1,100); for n = 1:100 A = rand(n,n); ...


11

For eigenvalues, simply take $k$ largest or smallest eigenvalues of $T$. They are good approximations of $A$, provided that the number of Lanczos iterations is large compared to $k$. Things are a little trickier if we want eigenvectors as well. The simplest way is to multiply each eigenvector $\mathbf{u}_i$ of $T$ by $V$ to the left, where $V$ is, as you ...


11

MATLAB's \ (aka mldivide) command does not blindly compute the inverse of the matrix. Instead, it uses one of several algorithms based on the type of matrix (see the "Algorithms" section of http://www.mathworks.com/help/matlab/ref/mldivide.html). In the case of a triangular matrix, MATLAB will use a triangular solver which is at least as good as yours in ...


11

And, it is my understanding that the 4 and the 5 are for the order of the global and local error, respectively. Your understanding is wrong. The local error of a Runge–Kutta method of order $n$ is proportional to $h^n$. What ode45 does is to estimate the solution (of one step) with two Runge–Kutta methods with local orders of 4 and 5, respectively (hence ...


11

You have an ill-conditioned eigenvalue problem. Consider a perturbation $\delta A$ in your matrix $A$—with double-precision floats this is $O(10^{-16})$. It turns the eigendecomposition from $$ X^{-1}A X = \Lambda $$ into (approximately keeping $X$ fixed) $$ X^{-1}(A+\delta A)X = \Lambda + \delta \Lambda. $$ This means that the change in eigenvalues is ...


11

A somewhat longer answer that explains why it's more efficient to have the left most index varying most rapidly. There are two key things that you need to understand. First, MATLAB (and Fortran, but not C and most other programming languages) stores arrays in memory in "column major order." e.g. if A is a 2 by 3 by 10 matrix, then the entries will be ...


10

Here I have an example: x = linspace(-5,5,100); y = linspace(-5,5,100); z = linspace(-5,5,100); [X, Y, Z] = meshgrid(x, y, z); Ex = sin(2*pi/5*Z); Ey = 0*X; Ez = 0*X; [Bx, By, Bz, V] = curl(X, Y, Z, Ex, Ey, Ez); Eplot = 0*x; Bplot = 0*x; for i=1:100 %% Integration-like procedure Eplot(i) = mean(mean(Ex(:,:,i),1),2); Bplot(i) = mean(mean(By(:,:,...


10

Unfortunately, convergence of GMRES does not have a clear dependence on the distribution of eigenvalues. It was proved by Greenbaum, Ptak and Strakos in 1996 that you can construct examples with an arbitrary spectrum and an arbitrary convergence history: that is, give me any $n$ nonzero complex numbers, and any decreasing sequence $\|r_k\|$, and I can ...


9

I think you should have a look at PyTables. Especially the tutorial given at PyData 2012. PyTables combines hierarchical datasets with a computational engine. It uses the Blosc compresser to avoid I/O bottlenecks and an optimized evaluator for expressions tables.Expr (based on Numexpr).


9

Before trying to find all of the roots of this function in MATLAB I think it's worth understanding that it has infinitely many roots due to the inclusion of the $\cos()$ term. Additionally, it is easy to find the roots of the function analytically in this case: The roots are defined by $$ \cos(7x)\cdot \exp(-2x^2)\cdot (1-2x^2) = 0. $$ So we have $$ \cos(7x)...


9

This is a linear integral equation (because the right hand side is linear in $h$) and presumably the problem is ill-posed (for one, because of the multiplication by $\sin(\theta)$ but also because the "integral kernel" is smooth). Methods to solve this (approximately) are general Galerkin methods or collocation methods. For Galerkin methods you take a ...


9

Given a nonlinear eigenvalue problem of the form $A(\lambda)x = 0$, reducing it to a real equation $\det(A(\lambda))=0$ is known to be a poor method for just the reason you've discovered yourself. The determinant is too ill-behaved for this to work on non-trivial problems. Even for an ordinary eigenvalue problem $Ax=\lambda x$, if you have a good ...


8

Matlab 6 used to have a function flop to count them, but it was removed in later versions. The main reason was technical (they switched to LAPACK as the linear algebra core, and it did not return a flop count). Today this option isn't available anymore; there is no easy equivalent, so if you wish to have one you'll have to count them yourself using the ...


8

The Matlab command help eps says the following: D = EPS(X), is the positive distance from ABS(X) to the next larger in magnitude floating point number of the same precision as X. X may be either double precision or single precision. In other words, if $\varepsilon_\mathsf{mach}$ is the relative error due to floating point, as defined in the Wikipedia ...


8

The benchmarks on the Julia website 1 2 include R and Matlab as competitors. Note that these are benchmarks focusing on testing the pure speed of the language, not the quality of the underlying linear algebra or FFT libraries. The speed for operations that are outsourced to these libraries (such as a large matrix multiplication) can vary a lot depending on ...


8

When you want to say that you want $u_2$ to be the best approximation of $u_1$ on mesh ${\cal T}_2$, then you have to define what you mean by "best". Let's assume you define it as that function that has the least $L_2$ error, i.e., $$ u_2 = \arg\min_{\phi_h \in V_2} \|\phi_h-u_1\|_{L_2(\Omega)}. $$ Then indeed $u_2$ is the $L_2$ projection of $u_1$ onto ...


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