# Tag Info

61

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 ...

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 computing ...

19

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 ...

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

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 ...

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

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 ...

7

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 ...

6

In general, I agree with Chris's comment that using a compiled language with the allocation of the matrices on the stack can help significantly. Several possibilities if we are limited to Python and numpy: consider np.array vs np.matrix, it might happen that np.matrix is faster than np.array matrix-matrix product (it is unclear what you are using now, and ...

5

There are several things to consider in this experiment: Why Matlab sparse direct might be "so fast": (for your particular test) In 2D (of course, problem-dependent), your matrix $A$ arising after FEM discretization, after some reorderings might appear to be "close to banded" structure. The smaller is the bandwidth of $A$, the more efficient a sparse ...

5

Estimates for GMRES convergence based on eigenvalue distribution often implicitly assume that the matrix is normal. Sometimes the convergence rate is still provable in an asymptotic sense in the non-normal case, but if the matrix is severely non-normal then the "pre-asymptotic" behavior will make such convergence rates never reachable in practice. Your ...

5

The particular set of constraints you have chosen does not prevent a rigid body rotation about node 1. Thus the stiffness matrix is singular, as you have noted. One way to prevent this rigid body rotation is to set the y-displacement at node 2 to zero. You could also constrain the x-displacement at either node 3 or node 4 to prevent the rotation. One way ...

4

The initial assumption was that the initial position was at a stable equilibrium (i.e., a minimum of the potential energy) with zero kinetic energy and the system started moving away from the equilibrium. Since physically it can't happen (if we consider classical mechanics), two things came to my mind: The first one is that maybe the initial position is: ...

4

We can break down the code rhsu = -a*rx.*(Dr*u) + LIFT*(Fscale.*(du)); into the two parts -a*rx.*(Dr*u) and LIFT*(Fscale.*(du)); The first part is simply taking a derivative in reference space by multiplying with the matrix Dr, then transforming into physical space by multiplying with the Jacobian rx and then the coefficient a. Next we need to do the ...

4

Look at the components of the forces calculated in your functions. You will probably find they are never exactly zero, because as other answers have said, you can't represent the value of $\pi$ exactly in computer arithmetic, and the routines that calculate trig functions are not exact either. Eventually, the tiny forces (probably of order $10^{-16}$ at ...

4

First of all, MATLAB's gmres assumes that the preconditioner you use is linear. This is important! Actually it is the main difference between FGMRES and GMRES. Right preconditioned GMRES and FGMRES are exactly the same if you use a linear preconditioner, however, FGMRES allows the use of non-linear preconditioners. What do I mean by a non-linear ...

4

As mentioned in a comment by @AloneProgrammer, it's unlikely that your system has a solution. But to make any kind of progress, it's also useful to rewrite the system in a way that makes its structure simpler. To this end, notice that in the first equation, you always have terms of the form $A_{ij}=B_{ij}/V_i$. So the first of the two equations might be ...

4

Standard examples of PDE to solve with the typically taught basic discretization methods (Crank-Nicolson et al.) are Transport equations, and other first order equations like Burger's, have often explicit solutions and conservation laws that the numerical methods more-or-less satisfy The heat equation with different boundary conditions and source terms is ...

4

Crank-Nicolson is a very good classical approach for parabolic PDE like the heat transfer PDE to which it was originally applied. It is relatively easy to understand and implement so it is often presented in basic courses on numerical methods for PDE. pdepe is also very well-suited to this class of PDE (the second "p" in pdepe stands for parabolic). It has ...

3

As far as I know, QZ decomposition is given for two matrices, so that for $A,B\in \mathbb F^{n \times n}$: $$A=QSZ^*, B=QTZ^* \tag{1} \label{QZ}$$ where $Q,Z^*$ are unitary, and $S,T$ are upper-triangular, and $\mathbb F$ is the field (real $\mathbb R$ or complex $\mathbb C$). QZ decomposition is usually called generalized Schur decomposition. For a ...

3

The situation seems to be: You have some input function $x$ which more or less follows a model $\dot x=-ax+bu$. There may be noise involved, so the values of $x$ are not exact, and simply computing difference quotients will in general not be close to the right side of the differential equation. To filter out the noise some averaging is required. This you do ...

3

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, congratulations, you certainly have a stronger point now, since the "competing" algorithm is given an advantage. It's worth no note though, that implementation in C++ is ...

3

Not sure whether this provides a full answer, but at least it provides some thoughts and hints. Concept I think the expression $\exp(\mathcal{Lt})\rho$ needs clarification. This is not necessarily a matrix exponential, but rather an abstract notation for the solution operator. For ODEs $\partial_t x = A x$, where $x$ is a vector and $A$ is a matrix (such ...

3

You might look up "complete fraction-free factorization" methods. The paper "Generalized fraction-free LU factorization for singular systems with kernel extraction" contains pseudo-code.

3

It is unlikely that you can benefit a lot from rewriting this particular code in C++. Main reasons: you are already assembling a sparse matrix using Matlab-specific framework for sparse matrices the hotspot is the eigenvalue computation, for which Matlab will use a call to a highly-optimized LAPACK implementation. Which would be similar to what you would ...

3

You will need to write your problem such that the unknowns are a single vector, not a matrix. In your example with $N=2$, you will have an unknown vector $x(t)$ of size $20\times 1$ (not a matrix of $10\times 2$). You will solve a problem of the following shape \dot{x}(t) = A(t)\, x(t), \mathrm{ with }\ x \in \mathbb{R}^{n},\ A(t) \in \mathbb{R}^{n\times n}...

3

It's almost impossible to say whether a direct solver will outperform an iterative solver or vice-versa without knowing more specific information about the sparse matrix. The key problem with direct solvers is fill-in, which happens during factorization phase and causes a lot of extra memory consumption. But not all sparse matrices will have a devastating ...

3

Matlab internally uses compressed sparse column (CSC) format for sparse matrices. The design and implementation of Matlab's sparse matrices are described in this document. As a consequence of using CSC format, indexing into sparse matrices can be an expensive operation. This is discussed in the help pages on sparse matrices.

3

There are many sparse matrices in Matrix Market A visual repository of test data for use in comparative studies of algorithms for numerical linear algebra, featuring nearly 500 sparse matrices from a variety of applications, as well as matrix generation tools and services. Use Matrix Market in conjunction with MM_TO_MSM: Matrix Market File to ...

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