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


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.


14

There are two choices which, if you use reasonable data structures, can solve the problem with $n>10^6$ on a laptop and $n \approx 10^{12}$ on a supercomputer. Note that for efficiency, you should use multigrid to solve with $\Delta$. The cost in either case will be a small factor more expensive than just solving with $\Delta$. The two approaches are ...


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


13

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


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

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


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


7

Calculate the SVD in place of the spectral decomposition. The results are the same in exact arithmetic, as your matrix is symmetric positve definite, but in finite precision arithmetic, you'll get the small eigenvalues with much more accuracy. Edit: See Demmel & Kahan, Accurate Singular Values of Bidiagonal Matrices, SIAM J. Sci. Stat. Comput. 11 (1990)...


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

The computation of integer-valued matrix determinants has been a subject of considerable research. Using exact arithmetic the Smith normal form can be computed, and from this diagonal form the determinant is easily found. Saunders and Wan (2004), Smith Normal Form of Dense Integer Matrices, Fast Algorithms into Practice, say "Over the past thirty years, ...


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

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


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

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

Wikipedia actually has a nice overview of approaches: http://en.wikipedia.org/wiki/Determinant#Calculation As a general remark, people do not compute determinants of large matrices (large would here be of size >10 or >50) because this is numerically difficult and, likely, not very stable anyway. If you need to do it, I would see if maybe there are ...


4

From geometrical tools (www.geometrictools.com) http://www.geometrictools.com/Documentation/EigenSymmetric3x3.pdf I hope this helps


4

A Google search returns many relevant links for this problem since it is a common problems in mechanics. A simple algorithm for symmetric and real 3x3 matrices is presented on Wikipedia: http://en.wikipedia.org/wiki/Eigenvalue_algorithm#3.C3.973_matrices


4

[Moving the comment thread down here] A direct solve will be pretty expensive, likely on the order of an hour for the sort of machine you describe. If you want invest time in something faster using more specialized methods, I suggest using FMM and a block Krylov method. FMM applies the operator in $\mathcal O(N)$ memory and time. FMM has a big constant, but ...


4

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


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


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