# Tag Info

22

(This is getting too long for comments...) I'll assume you actually need to compute an inverse in your algorithm.1 First, it is important to note that these alternative algorithms are not actually claimed to be faster, just that they have better asymptotic complexity (meaning the required number of elementary operations grows more slowly). In fact, in ...

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

In general, there is no shortcut other than completely re-factoring the matrix. There have been a few similar questions on this SE that cover the topic in more depth than I can: Can diagonal plus fixed symmetric linear systems be solved in quadratic time after precomputation? LU Decom of PSD Matrix + Diagonal Matrix Perturbation of Cholesky decomposition ...

13

BLAS1-operations, BLAS2-operations, and sparse-operations share the same curse of low arithmetic intensity, that they perform $O(1)$ flops for each memory read (contrast this to a BLAS3-operation like gemm, which performs $O(N^3)$ flops over $O(N^2)$ reads and only becomes more and more arithmetic-intensive/compute-bound for large $N$). However, sparse ...

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

11

Use a Cholesky decomposition or an LDL decomposition instead of LU. Judging from answers to: Diagonal update of a symmetric positive definite matrix Solving a system with a small rank diagonal update Can diagonal plus fixed symmetric linear systems be solved in quadratic time after precomputation? there's no good way to update any of those decompositions. ...

10

Iterative Krylov-subspace solvers generally only require matrix-vector products and don't care whether or where there are zeros in the matrix. In your case, unless you have other information about the matrix (e.g., symmetry), you could for example use GMRES. What you probably had in mind is the question of preconditioning, and that you can't use things such ...

8

@BillGreene points to the "return value optimization" as a way around the fundamental problem, but this actually only helps for one half of it. Assume you have code of this form: struct ExpensiveObject { ExpensiveObject(); ~ExpensiveObject(); }; ExpensiveObject operator+ (ExpensiveObject &obj1, ExpensiveObject &obj2) { ...

7

One option here would be to form the normal equations $A^{T}Ax=A^{T}b$ and solve them by Cholesky factorization of the resulting $n$ by $n$ matrix. This squares the condition number of the problem which could potentially be a significant problem. Forming $B=A^{T}A$ doesn't require more than $O(n^2)$ memory, assuming that you can access the rows of $A$ ...

6

You should probably note that, buried deep inside the numpy source code (see https://github.com/numpy/numpy/blob/master/numpy/linalg/umath_linalg.c.src) the inv routine attempts to call the dgetrf function from your system LAPACK package, which then performs an LU decomposition of your original matrix. This is morally equivalent to Gaussian elimination, but ...

6

In double precision, a system of 20,000 equations in 20,000 variables will require 3.2 gigabytes of RAM to store the matrix. This isn't "big" by contemporary standards, and is reasonably easy to solve in practice. To make this run fast, you'll want to use of a well tuned linear algebra library. The standard for this kind of work would be to use LAPACK ...

6

Very quick answer... The exponential of a Hamiltonian matrix is symplectic, a property that you probably wish to preserve, otherwise you would simply use a non-structure-preserving method. Indeed, there is no real speed advantage in using structured method, just structure preservation. A possible way to solve your problem is the following. First find a ...

6

Even if a matrix is very sparse, its matrix product with itself can be dense. Take for example a diagonal matrix and fill its first row and column with nonzero entries; its product with itself will be completely dense. Such a matrix can arise, for examle, as graph Laplacian of a graph in which there is a vertex that is connected to all other vertices. In ...

5

Stackoverflow (https://stackoverflow.com/) is probably a better discussion forum for this question. However, here is a short answer. I doubt that the C++ compiler is generating code for this expression like you describe above. All modern C++ compilers implement an optimization called "return value optimization" (http://en.wikipedia.org/wiki/...

5

I'd like to refute your premise that this is “very, very inefficient”. Performance-wise, it's actually pretty irrelevant, because matrix multiplication is1 $\mathcal{O}(n^{2.8})$, whereas copying one of those matrices is only $\mathcal{O}(n^2)$. So for $n$ as big as you have, the multiplication will probably completely outweigh the copying time. ...

5

Modern C++ has a solution for the problem by using "move constructors" and "rvalue references". A "move constructor" is a constructor for a class, for example a matrix class, which takes another instance of the same class and moves the data from the other instance into the new instance, leaving the original instance empty. Typically, a matrix object will ...

5

For dense distributed memory linear algebra, you can't do better than Elemental these days. Although, it does not currently handle sparse matrices, I believe. Dense-sparse operations will probably require you to do some custom coding.

5

This really depends on the operations you are including in your question. If you took the sparse equivalent of any level 1 BLAS or level 2 BLAS algorithm, then yes they are memory bound (not compute bound), but that is also true of level 1 BLAS and level 2 BLAS for dense matrices. To answer this we need to understand what makes an algorithm compute bound. ...

5

Upon further examination, I do think the Woodbury identity can be used to solve this problem. With it we can write: $\left( \mathbf K - \sigma \mathbf Z \mathbf Z^T \right)^{-1} = \mathbf K^{-1} - \mathbf K^{-1} \mathbf Z \left(-\frac{1}{\sigma}\mathbf I + \mathbf Z^{T} \mathbf K^{-1} \mathbf Z\right)^{-1}\mathbf Z^T \mathbf K^{-1}$ Since $\mathbf K$ is ...

4

SLEPc is a library designed for the solution of large sparse eigenvalue problems, while your matrix is dense. The algorithms used by SLEPc are based on matrix-vector multiplication, and the performance of these algorithms is primarily limited by the available memory bandwidth rather than the number of processor cores. You might be better off using a ...

4

In addition to the QR algorithm, the divide and conquer method is also worth mentioning. It is applicable to symmetric tridiagonal matrices, but any matrix can be reduced to such a form via the Lanczos method*. It hinges on the observation that a tridiagonal matrix is, up to a rank 1 perturbation, a block diagonal matrix. One can then find the ...

4

In my experience the answer to this question is not clear-cut; it is dependent on the span of the eigenvalues and relative matrix structure itself. That said, your current approach evokes an implictly shifted QR solver that is essentially the standard for this exact type of problem. With some experimentation, however, you may squeak out a modest performance ...

3

There's no element wise multiplication operation in the BLAS library. Your best approach is probably to just implement the operation yourself using (e.g.) OpenMP threading. Before you do this, you should consider Amdahl's law and whether speeding up this bit of your code is really going to help- chances are that these elementwise multiplications are not ...

3

I have been using pari/gp (confuration displayed in output below) to calculate determinants of large, dense matrices. For example, I use a quad core i7, 16GB RAM, 32GB swap space on internal ssd, linux laptop, to calculate a matrix 1700x1700 to a precision of a little less than 3000 digits. That takes about 17 hours, using all eight threads, 100% of the time,...

3

Have you considered working with the Cholesky (or low-rank) factor of $\rho(t_0)$ rather than with the matrix itself? This might reduce the number of products that you need to make, and it has the additional benefit of preserving positive semidefiniteness across your computations. It is already cheaper to use Cholesky factors for this computation alone, if ...

3

If the matrix is small enough to fit into memory, then there is of course no cost associated with actually forming the elements: You will have to compute the elements at least once anyway to perform the matrix-vector product, so you might as well store the values you computed in memory and then run an optimized matrix-vector product function over it. But ...

3

Conclusion: I have tried CG, BiCG, and GMRES with no pre-conditioner, and these are all at least twice as slow as LU decomposition for the 121x121 matrix. Further, even for convergence tolerances near machine-precision, I do not get the accuracy that I get with LU. So I guess I'm sticking with LU :)

2

Here are two examples: http://math.nist.gov/MatrixMarket/formats.html http://tedlab.mit.edu/~dr/SVDLIBC/SVD_F_DB.html I don't claim either is a good format...

2

A differential operator will always result in rows with sum zero. Furthermore, if you want to preserve the locality of the differential operator, the entries with opposite signs will be neighbors (possibly separated by zeros). Smoothly varying entries come up if you discretize an integral operator with smooth kernel by a collocation or quadrature method. I ...

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