64
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
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
The column major layout is the scheme used by Fortran and that's why it's used in LAPACK and other libraries.
In general it is much more efficient in terms of memory bandwidth usage and cache performance to access the elements of an array in the order in which they're laid out in memory. Depending on how your matrices are stored, you'll want to pick ...
18
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=\...
17
The first thing is to recognize that you can do this using BLAS. If you data matrix is $X = [x_1 x_2 x_3 ...] \in \mathbb{R}^{m\times n}$ (each $x$ is a column vector corresponding to one measurement; rows are trials), then you can write the covariance as:
$$ C_{ij} = E[x_i,x_j] - E[x_i] E[x_j] = \frac{1}{n} \sum_k x_{ik} x_{jk} - \frac{1}{n^2} \left(\sum_k ...
17
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 ...
16
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.
15
The Cholesky factorization $C=R^TR$ leads to a Cholesky-like factorization of the inverse $C^{-1}=SS^T$ with the upper triangular matrix $S=R^{-1}$.
In practice, is best to keep the inverse factored. If $R$ is sparse then it is usually even better to keep $S$ implicit, as matrix-vector products $y=C^{-1}x$ can be computed by solving the two triangular ...
14
Matrix Market is a terrible format for reading in parallel, therefore it is better to preprocess to a better parallel format. Your matrix size is extremely small so performance is not an issue, but the easiest and most general thing is to use Python or Matlab/Octave to write the Matrix Market file in PETSc binary format, which can be read efficiently in ...
14
The most obvious thing you can do is to precompute
[L,U] = lu(A) ~ O(n^3)
Then you just compute
x = U \ (L \ b) ~ O(2 n^2)
This would reduce the cost enormously and make it faster. Accuracy would be the same.
14
Apart from the (now classical) Golub-Reinsch paper Brian notes in his answer (I have linked to the Handbook version of the paper), as well as the (also now classical) predecessor paper of Golub-Kahan, there have been a number of important developments in computing the SVD since then. First, I have to summarize how the usual method works.
The idea in ...
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
I actually wrote the original code in Matlab for A*B, both A and B sparse. Pre-allocation of space for the result was indeed the interesting part. We observed what Godric points out -- that knowing the number of nonzeros in AB is as costly as computing AB.
We did the initial implementaion of sparse Matlab around 1990, before the Edith Cohen paper that ...
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....
12
You're right -- it has absolutely no practical relevance for computing. Even if computing the determinant was an $O(n)$ operation, the complexity of the method would be at least $O(n^3)$ and, consequently, of the same complexity as Gaussian elimination. In practice, computing the determinant of a matrix is actually of exponential complexity, making this ...
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
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
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
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
In exact arithmetic you shouldn't need to reorthogonalize regularly, but practically you do. Your u1 and u2 are close to (but not exactly) the true eigenvectors, so your initial deflation almost (but not entirely) removed the true eigenvectors from u3. The tiny components you left behind will be amplified by repeated multiplication by A, you will need to ...
11
In vacuum without considering any existing software, there's no reason to prefer column major over row major from the code point of view. However, most mathematical literature is written in a way that groups vectors into a matrix by storing them as columns instead of rows. For example when you write the full eigenvalue equation $AX=X\Lambda$, the $X$ matrix ...
11
This is typically done using the Golub-Reinsch algorithm, and no, it doesn't involve computing eigenvalues and eigenvectors of $AA^{T}$.
See
G. H. Golub and C. Reinsch. Singular Value Decomposition and Least Squares Solutions. Numerische Mathematik 14:403-420, 1970.
This material is discussed in many textbooks on numerical linear algebra.
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 ...
Only top voted, non community-wiki answers of a minimum length are eligible
