# Tag Info

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

27

Stability under perturbations Let $E$ be a perturbation such that $\|E\| \leq \varepsilon$. If $A$ is symmetric, then the eigenvalues of $A+E$ are at a distance $\varepsilon$ from those of $A$. (Bauer-Fike Theorem.) If $A$ is non-symmetric, then the eigenvalues of $A+E$ could be just anywhere in the complex plane. No bound can be formulated a priori. (...

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

19

See https://math.stackexchange.com/questions/861674/decompose-a-2d-arbitrary-transform-into-only-scaling-and-rotation (sorry, I would have put that in a comment but I've registered just to post this so I can't post comments yet). But since I'm writing it as an answer, I'll also write the method: $$E=\frac{m_{00}+m_{11}}{2}; F=\frac{m_{00}-m_{11}}{2}; G=\... 19 The property follows from the property of the corresponding (weak form of the) partial differential equation; this is one of the advantages of finite element methods compared to, e.g., finite difference methods. To see that, first recall that the finite element method starts from the weak form of the Poisson equation (I'm assuming Dirichlet boundary ... 18 For what it is worth, for random sparse matrices of size 10,000 by 10,000 vs. dense matrices of the same size, on my Xeon workstation using MATLAB and Intel MKL as the BLAS, the sparse matrix-vector multiply was faster for densities of 15% or less. At 67% (as proposed by another answer), the dense matrix-vector multiplication was about three time faster. 14 In fundamental C++, I find the problem here is that C++ will allocate a new object of cx_mat to store evolutionMatrix*stateMatrix, and then copy the new object to stateMatrix with operator=(). I think you're right that it's creating temporaries, which is too slow, but I think the reason for why it's doing that is wrong. Armadillo, like any good C++ linear ... 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 For problems I am interested in, the matrix dimension is 30 or less. As WolfgangBangerth notes, unless you have a large number of these matrices (millions, billions), performance of matrix inversion typically isn't an issue. Given a positive definite symmetric matrix, what is the fastest algorithm for computing the inverse matrix and its determinant? If ... 13 If the solution of Ax=b is unstable, the matrix is very ill-conditioned (i.e., has a very large condition number), and (paraphrasing Lanczos) no amount of mathematical trickery can make it stable. The best you can hope for is to solve a different problem that is a) stable and b) gives you a solution that is sufficiently close; this is called regularization.... 13 Here is a quick example which is very practical related to memory usage in PDEs. When one discretizes the Laplace operator, \Delta u, for example, in the Heat Equation$$ u_t = \Delta u + f(t,u) .$$To solve it numerically, one ends up with sparse matrices A, and a method of lines discretization then solve$$ u_t = Au + f(t,u) $$The canonical 1D ... 12 The LU factors of a sparse matrix are at least somewhat sparse. The Q matrix in QR can also somewhat preserve sparsity, and is typically used when the matrix is very long and skinny. The SVD of a sparse matrix will almost always have fully dense U and V factors, so it destroys any reason to perform the computations treating the matrix sparsely. 12 Consolidating the comments: No, you are very unlikely to beat a typical BLAS library such as Intel's MKL, AMD's Math Core Library, or OpenBLAS.1 These not only use vectorization, but also (at least for the major functions) use kernels that are hand-written in architecture-specific assembly language in order to optimally exploit available vector extensions (... 12 Normally there are some principal reasons to prefer solve a linear system respect to use the inverse. Briefly: problem with the conditional number (@GoHokies comment) problem in the sparse case (@ChrisRackauckas answer) efficiency (@Kirill comment) Anyway, as @ChristianClason remarked in comments, can be some cases where the use of the inverse is a good ... 12 All matrix operations are memory bound (and not compute bound) on today's processors. So basically, you have to ask which format stores fewer bytes. This is easy to compute: For a full matrix, you store 8 bytes (one double) per entry For a sparse matrix, you store 12 bytes per entry (one double for the value, and one integer for the column index of the ... 12 Quick sketch of an answer for the Frobenius norm: Replace X with the closest symmetric matrix, Y=\frac12 (X+X^\top). Take an eigendecomposition Y=QDQ^\top, and form the diagonal matrix D_+=\max(D,0) (elementwise maximum). The closest symmetric positive semidefinite matrix to X is Z=QD_+Q^\top. The closest positive definite matrix to X does not ... 12 If you can compute products with A and A^T, as you specify in a comment, you can run the classical sparse SVD algorithms such as scipy.sparse.linalg.svds, Matlab's svds, or Julia's Arpack.svds, which are based on Lanczos bidiagonalization. They are designed to compute singular values, and are likely to be more robust than a minimization routine coded by ... 11 The code that you've posted uses the eigenvalue decomposition of the symmetric matrix to compute A^{-1/2}. The statement d=(d+abs(d))/2 effectively takes any negative entry in d and sets it to 0, while leaving non-negative entries alone. That is, any negative eigenvalue of A is treated as though it was 0. In theory, the eigenvalues of A should ... 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 It is true that compilers are getting better and better at auto-vectorization, and for basic coefficient-wise operations like 2*A-4*B a library like Eigen cannot do much better than recent compilers. However, for slightly more complicated expressions like matrix products, reductions, transposition, powers, etc. the compiler cannot do much. On the other hand, ... 11 First: there is a must read on this topic Moler, Cleve, and Charles Van Loan. "Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later." SIAM review 45.1 (2003): 3-49. (in case you wonder, the original paper is Moler, Cleve, and Charles Van Loan. "Nineteen dubious ways to compute the exponential of a matrix." SIAM review ... 11 Since$$ A = B(I-B)^{-1} = (I-B)^{-1}(I-B)B(I-B)^{-1} = (I-B)^{-1}B(I-B)(I-B)^{-1} =(I-B)^{-1}B $$So you want to solve$$ (I-B)A=B $$You seem to need only the first three columns of A. Solve the matrix problems$$ (I-B)a_i = b_i, \qquad i=0,1,2 $$where b_0,b_1,b_2 are first three columns of B. Then a_0,a_1,a_2 are the first three columns of A. 10 @Pedro Gimeno "I doubt it can be any more robust than that." Challenge accepted. I noticed the usual approach is to use trig functions like atan2. Intuitively, there shouldn't be a need to use trig functions. Indeed, all the results end up as sines and cosines of arctans--which can be simplified to algebraic functions. It took quite a while, but I managed ... 10 I think I'd rather see it done this way. Since A is invertible, then \min_i\sigma_i>0,$$ \min_i\sigma_i = \inf_{x\neq 0}\frac{\|Ax\|_2}{\|x\|_2} \quad\Longleftrightarrow \frac{1}{\min_i\sigma_i} = \sup_{x\neq 0}\frac{\|x\|_2}{\|Ax\|_2}. $$Then we have$$ \frac{1}{\min_i\sigma_i} = \sup_{x\neq 0}\frac{\|x\|_2}{\|Ax\|_2} = \sup_{A^{-1}y\neq 0}\frac{\|...

10

This can be interpreted as summing over an index of a tensor when the vector $x$ is reshaped into a box of numbers instead of a list. In particular, if $X$ is the $d\text{-by-}d$ folded version of $x$, then the operation you are doing is, \begin{align} Dx &= \mathrm{vec}\left((I \otimes \mathbf{1})\mathrm{vec}(X)\right) \\ &= \mathrm{vec}(\mathbf{1}^...

10

The question is ill-posed -- you don't say what you mean by "desirable" or "superior". If a phenomenon I try to model involves only real numbers, then clearly using real matrices is the better way. On the other hand, if I try to model something in the Fourier domain, then complex numbers appear and I would try to use a matrix with complex elements, rather ...

9

I'm going against the crowd - the adjugate matrix is in fact very useful for some specialty applications with small dimensionality (like four or less), in particular when you need the inverse of a matrix but don't care about scale. Two examples include computation of an inverse homography and Rayleigh quotient iteration for very small problems (which in ...

9

It's an important concept when developing numerical methods for ordinary differential equations, at least as long as they are linear. Consider, for simplicity, a homogenous system of the form $\dot x(t)=Ax(t), x(0)=x_0$, then the solution is $x(t)=e^{At}x_0$. This is easily extended to the case $\dot x(t)=Ax(t)+f(t)$. The point now is not that one may ...

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

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