25 votes
Accepted

Complexity of matrix inversion in numpy

(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 ...
Christian Clason's user avatar
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 ...
Brian Borchers's user avatar
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 ...
Wolfgang Bangerth's user avatar
13 votes
Accepted

Why does sparse linear algebra have a low arithmetic intensity?

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 ...
rchilton1980's user avatar
  • 4,862
10 votes
Accepted

Is there an iterative solver for dense matrices with possible zero diagonal entries?

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

Rule of thumb for sparse vs dense matrix storage

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 ...
Henrik Schumacher's user avatar
7 votes

Complexity of matrix inversion in numpy

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 ...
origimbo's user avatar
  • 2,239
6 votes
Accepted

Matrix exponential of a Hamiltonian matrix

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

inertia count sparse matrix with dense low-rank perturbation

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} - \...
rchilton1980's user avatar
  • 4,862
5 votes

Why does sparse linear algebra have a low arithmetic intensity?

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 ...
Reid.Atcheson's user avatar
4 votes

Is there an iterative solver for dense matrices with possible zero diagonal entries?

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 ...
niran90's user avatar
  • 233
3 votes

Fast matrix multiplication with matrix elements computed on-the-fly (without forming the matrix)

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

Solve linear system for only part of the solution vector

The only way that I'm aware of to take advantage of only needing a partial solution is in the triangular solves. But those are usually a small part of the total time in dense LU (unless you have many ...
Neil Lindquist's user avatar
3 votes

Smart way to multiply 3 matrices

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

Compute all eigenvectors and eigenvalues of small symmetric matrices

You might want to try the SelfAdjointEigenSolver class https://eigen.tuxfamily.org/dox/classEigen_1_1SelfAdjointEigenSolver.html in the Eigen C++ class library http://eigen.tuxfamily.org/index.php?...
Bill Greene's user avatar
  • 5,984
2 votes
Accepted

inertia count sparse matrix with dense low-rank perturbation

This answer is based on rchilton1980 answer and the comments that have followed. Let $K = L L^T$ be the Cholesky factorization of $K$. Let $Y = L^{-1} Z$ , computed by forward substitution. Let $\...
Olivier's user avatar
  • 81
2 votes

Tools to compare two matrices with same dimensions

Welcome to Scicomp! You might be interested in the field of (multimodal) medical image registration. In medical contexts one often wants to register a CT image with an MR or PET-Scan. The amplitudes ...
MPIchael's user avatar
  • 2,842
2 votes
Accepted

Solving a linear system whose coefficient matrix is dense but symmetric

It can be lowered; it is called packed storage and Lapack has some functions to deal with it, e.g., ?PPSVX, ?SPSVX. As this ...
Federico Poloni's user avatar
2 votes
Accepted

Trace of inverse from LU decomposition

Base case: Computing the factorization requires $\frac 2 3n^3$ operations and the inverse requires $\frac 4 3n^3$ operations. Computing the trace adds $O(n)$. Let's round that to $2n^3$. LU case: ...
Charlie S's user avatar
  • 661
1 vote

Library to solve dense linear system with GMRES

PETSC might be a good option. Not super user friendly, but its good with lots of options.
EMP's user avatar
  • 2,079
1 vote

Matrix exponential of a Hamiltonian matrix

You might have an option to use hierarchical matrices ($\mathcal H$-matrices) and the corresponding functionality of the libraries that support them. Effectively, if every matrix $A$, $G$, and $Q$ ...
Anton Menshov's user avatar
  • 8,652
1 vote

Optimal algorithm choice for mixed diagonal/dense problem

This is what I wrote in the comments, formulated as an answer. If $\hat{\beta}_j=0$, then $\hat{\alpha}_j\neq 0$ otherwise the matrix wouldn't be invertible. So you can swap column $j$ with column $n+...
Federico Poloni's user avatar
1 vote

Comparison between two matrices

There are multiple ways to judge the information lost. The easiest way would be to calculate Frobenius norm of the, so called, sparsified matrix $\hat A$ and an original dense matrix. $||A||_F=\sqrt{...
Anton Menshov's user avatar
  • 8,652
1 vote
Accepted

Discretized matrix from the integral kernel function

(1) It appears that $(i_1,i_2),(j_1,j_2)$ are coordinate pairs, such that $(j_1,j_2)$ are the integer coordinates on the source grid and $(i_1,i_2)$ are the integer coordinates on the observer grid. ...
coolguy1000000's user avatar
1 vote

Preconditioning of two step iteration for dense matrices

I would try to circumvent the problem entirely by the substitution $x=By$. Then your equation reduces to $$ \left[ (\gamma^+ + \gamma^-) AB + \frac{1}{2} (\gamma^+-\gamma^-) I \right] y = b.$$ If $AB$ ...
Carl Christian's user avatar
1 vote

Computing sparse matrix products into a dense result

Since the question was asked a long time ago, I would give a detailed overview answer. The first question is how the assembled matrix $A$ will be used later, as the construction of the matrix itself ...
Anton Menshov's user avatar
  • 8,652

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