Skip to main content
68 votes
Accepted

why is A*v+B*v faster than (A+B)*v?

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 ...
Brian Borchers's user avatar
32 votes
Accepted

Why do we usually not want the eigenvalues of non-symmetric matrices?

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-...
Federico Poloni's user avatar
22 votes
Accepted

How to find the nearest/a near positive definite from a given matrix?

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_+=\...
Federico Poloni's user avatar
22 votes

why is A*v+B*v faster than (A+B)*v?

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: ...
Rainer P.'s user avatar
  • 501
20 votes

Rule of thumb for sparse vs dense matrix storage

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 ...
Brian Borchers's user avatar
16 votes
Accepted

Rule of thumb for sparse vs dense matrix storage

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 ...
Wolfgang Bangerth's user avatar
16 votes
Accepted

Why is it that SVD routines on nearly square matrices run significantly faster than if the matrix was highly non square?

You are asking for a full (dense) SVD, which also needs to generate the unitary components of $U$ and $V$ which correspond with the null space of your input. for the $1000 \times 800$ case, your input ...
helloworld922's user avatar
15 votes

Practical example of why it is not good to invert a matrix

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) .$$...
Chris Rackauckas's user avatar
15 votes

Matrix multiplication accuracy Matlab vs Python

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 ...
Federico Poloni's user avatar
14 votes
Accepted

Beating typical BLAS libraries matrix multiplication performance

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

Matrix multiplication accuracy Matlab vs Python

Here is R1, as computed in MATLAB: ...
Mark L. Stone's user avatar
13 votes
Accepted

Practical example of why it is not good to invert a matrix

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 (@...
Mauro Vanzetto's user avatar
12 votes
Accepted

Why do libraries need hand-vectorized code instead of compiler auto vectorization

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 ...
ggael's user avatar
  • 689
12 votes
Accepted

Compute all eigenvalues of a very big and very sparse adjacency matrix

You can use the shift-invert spectral transform [1] and compute the spectrum band by band. The technique is also explained in my article [2]. Besides the implementation in [1], an implementation is ...
BrunoLevy's user avatar
  • 2,315
12 votes

Checking singularity of a matrix

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, ...
Federico Poloni's user avatar
11 votes
Accepted

Approximating the exponent of a matrix $\exp(A)$ using Taylor series

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. ...
Dirk's user avatar
  • 1,738
11 votes
Accepted

Method to check for positive definite matrices

The standard way of approaching this is really to attempt a Cholesky factorization and check to see if such a factorization exists. This is both fast and realiable. Here it is in MATLAB notation: A ...
Tolga Birdal's user avatar
  • 2,249
11 votes

Computing the Cholesky decomposition based of the QR decomposition

No. In general, the QR decomposition has no relation to the Cholesky decomposition of a symmetric positive definite matrix. Moreover, the QR decomposition is substantially more expensive to compute ...
Carl Christian's user avatar
11 votes
Accepted

Do I really need to invert this matrix

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 ...
cfdlab's user avatar
  • 3,028
11 votes

What algorithm(s) do numpy and scipy use to calculate matrix inverses?

Documentation to numpy.linalg.inv and scipy.linalg.inv does not mention the algorithm used. Judging from the source, ...
Vladimir Lysikov's user avatar
11 votes

Optimized Lanczos method for finding eigenvalues of $A \otimes B$

There isn't much complicated behind this idea; it's just that since Lanczos is a black-box method you can use any method of your choice to compute the products $v\mapsto (A\otimes B)v$ needed in the ...
Federico Poloni's user avatar
10 votes

Robust algorithm for $2 \times 2$ SVD

I needed an algorithm that has little branching (hopefully CMOVs) no trigonometric function calls high numerical accuracy even with 32 bit floats We want to calculate $c_1, s_1, c_2, s_2, \sigma_1$ ...
petiaccja's user avatar
  • 101
10 votes

Method to check for positive definite matrices

Mostly, I'm leaving this answer here as a cautionary tale to not use a Choleski factorization. Most of the time, this is a fine answer. However, very specifically, it can and will fail, so if this ...
wyer33's user avatar
  • 767
9 votes

In constructing matrices to model physical phenomena, are real matrices superior to complex matrices, in terms of computational cost?

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 ...
Wolfgang Bangerth's user avatar
9 votes
Accepted

Convexity of Sum of $k$-smallest Eigenvalue

Given $A \in {\bf S}^n$ (a positive definite matrix) with eigenvalues $\lambda_1 \leq \lambda_2 \leq \ldots \leq \lambda_n $, then: $\displaystyle f_k(A)=\sum_{i=1}^{k} \lambda_i$ is concave. Why? $$...
GoHokies's user avatar
  • 2,216
9 votes
Accepted

Efficiently computing $e^{tX}$ for many different values of $t$

An anti-Hermitian matrix is diagonalizable, with orthogonal eigenvectors (ref). Hence you can write $X = PDP^{-1}$, where $D$ is a diagonal matrix. Therefore the exponential can be calculated as $e^X=...
rpm2718's user avatar
  • 254
8 votes

Compute $x = B^{-1}(2A+I)(C^{-1}+A)b$ without calculating matrix inverses

As was mentioned in the comment, calculating $x=M^{-1}y$ is equivalent to solving $Mx=y$. Here is the full solution: First, you can reformulate the equation to: $Bx=(2A+I)(C^{-1}+A)b$, and by ...
Gil's user avatar
  • 392
8 votes
Accepted

Block-matrix SVD and rank bounds

Let us begin with the exact singular value decompositions $A=U_{A}S_{A}V_{A}^{T}$, $B=U_{B}S_{B}V_{B}^{T}$, $C=U_{C}S_{C}V_{C}^{T}$, $D=U_{D}S_{D}V_{D}^{T}$. Then $$ M=\underbrace{\begin{bmatrix}U_{A} ...
Richard Zhang's user avatar
8 votes

Rapidly determining whether or not a dense matrix is of low rank

The problem, of course, is that computing the true rank (e.g., via a QR decomposition) is not really any cheaper than computing a low-rank representation of the matrix. The best you can probably do ...
Wolfgang Bangerth's user avatar
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

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