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 ...
- 18.5k
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-...
- 11.1k
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:
...
- 451
21
votes
Accepted
In FEM, why is the stiffness matrix positive definite?
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 ...
- 12.2k
21
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_+=\...
- 11.1k
19
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 ...
- 18.5k
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 ...
- 2,391
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 ...
- 11.1k
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
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) .$$...
- 12.2k
14
votes
14
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 ...
- 54.1k
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 ...
- 666
12
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 (@...
- 1,330
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, ...
- 11.1k
11
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 ...
- 2,305
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.
...
- 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 ...
- 2,219
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 ...
- 1,391
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 ...
- 2,993
9
votes
Does the "cofactor technique" for inverting a matrix have any practical significance?
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 ...
- 171
9
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$ ...
- 91
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 ...
- 54.1k
9
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 ...
- 757
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=...
- 254
8
votes
Accepted
Finding the matrix inverse given a solver for the matrix equation $Ax=b$
Your two ideas make it much too complicated. If $X$ is the inverse of $A$,
$$ AX=I, $$
and $x_i$ is the $i$-th column of $X$ and $e_i$ is the $i$-th column of the identity matrix $I$ ($e_i$ is a ...
- 11.4k
8
votes
Accepted
Stabilizing a 3x3 real symmetric matrix eigenvalue calculation
This is trying to compute the eigenvalues by computing the roots of the characteristic polynomial. In this case, the characteristic polynomial is $p(t) = t^3-2t^2x$, $x=1.25\times 10^6$, and zero is a ...
- 11.4k
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 ...
- 382
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 ...
- 54.1k
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{...
- 11.1k
Only top scored, non community-wiki answers of a minimum length are eligible