Questions tagged [linear-algebra]
Questions on the algorithmic/computational aspects of linear algebra, including the solution of linear systems, least squares problems, eigenproblems, and other such matters.
1,124
questions
170
votes
8
answers
139k
views
Recommendations for a usable, fast C++ matrix library?
Does anyone have recommendations on a usable, fast C++ matrix library?
What I mean by usable is the following:
Matrix objects have an intuitive interface (ex.: I can use rows and columns while ...
- 30.3k
58
votes
4
answers
9k
views
What guidelines should I follow when choosing a sparse linear system solver?
Sparse linear systems turn up with increasing frequency in applications. One has a lot of routines to choose from for solving these systems. At the highest level, there is a watershed between direct (...
- 3,145
40
votes
23
answers
10k
views
Good examples of "two is easy, three is hard" in computational sciences
I recently encountered a formulation of the meta-phenomenon: "two is easy, three is hard" (phrased this way by Federico Poloni), which can be described, as follows:
When a certain problem is ...
- 8,602
40
votes
4
answers
35k
views
How does the MATLAB backslash operator solve $Ax=b$ for square matrices?
I was comparing a few of my codes to "stock" MATLAB codes. I am surprised at the results.
I ran a sample code (Sparse Matrix)
...
- 3,344
40
votes
3
answers
18k
views
How to choose a method for solving linear equations
To my knowledge, there are 4 ways to solving a system of linear equations (correct me if there are more):
If the system matrix is a full-rank square matrix, you can use Cramer’s Rule;
Compute the ...
- 647
38
votes
10
answers
9k
views
stupid + stupid = brilliant in scientific computing
I'm interested in examples of very effective methods in scientific computing that are the sum or naive combination of very ineffective or bad ones.
- 10.1k
37
votes
2
answers
12k
views
why is A*v+B*v faster than (A+B)*v?
$A$ and $B$ are $n \times n$ matrices and $v$ is a vector with $n$ elements. $Av$ has $\approx 2n^2$ flops and $A+B$ has $n^2$ flops. Following this logic, $(A+B)v$ should be faster than $Av+Bv$.
Yet,...
- 525
34
votes
6
answers
33k
views
What is the fastest way to calculate the largest eigenvalue of a general matrix?
EDIT: I am testing if any eigenvalues have a magnitude of one or greater.
I need to find the largest absolute eigenvalue of a large sparse, non-symmetric matrix.
I have been using R's ...
- 511
33
votes
5
answers
15k
views
Performance differences between ATLAS and MKL?
ATLAS is a free BLAS/LAPACK replacement that tunes itself to the machine when compiled. MKL is the commercial library shipped by Intel. Are these two libraries comparable when it comes to performance, ...
- 1,589
31
votes
3
answers
4k
views
What is the relationship of BLAS, LAPACK, and other linear algebra libraries?
I have been looking into C++ linear algebra libraries for a project I've been working on. Something that I still don't have any grasp on is the connection of BLAS and LAPACK to other linear algebra ...
- 469
31
votes
5
answers
29k
views
Permute a matrix in-place in numpy
I want to modify a dense square transition matrix in-place by changing the order of several of its rows and columns, using python's numpy library. Mathematically this corresponds to pre-multiplying ...
- 508
30
votes
11
answers
10k
views
Robust algorithm for $2 \times 2$ SVD
What is a simple algorithm for computing the SVD of $2 \times 2$ matrices?
Ideally, I'd like a numerically robust algorithm, but I'll like to see both simple and not-so-simple implementations. C code ...
- 966
30
votes
2
answers
15k
views
Why is my iterative linear solver not converging?
What can go wrong when using preconditoned Krylov methods from KSP (PETSc's linear solver package) to solve a sparse linear system such as those obtained by discretizing and linearizing partial ...
- 25.6k
29
votes
2
answers
15k
views
Does a tiny determinant imply ill-conditioning of a matrix?
If I have a square invertible matrix and I take its determinant, and I find that $\det(A) \approx 0$, does this imply that the matrix is poorly conditioned?
Is the converse also true? Does an ill-...
- 3,344
29
votes
4
answers
19k
views
Dealing with the inverse of a positive definite symmetric (covariance) matrix?
In statistics and its various applications, we often calculate the covariance matrix, which is positive definite (in the cases considered) and symmetric, for various uses. Sometimes, we need the ...
28
votes
3
answers
4k
views
What is the principle behind the convergence of Krylov subspace methods for solving linear systems of equations?
As I understand it, there are two major categories of iterative methods for solving linear systems of equations:
Stationary Methods (Jacobi, Gauss-Seidel, SOR, Multigrid)
Krylov Subspace methods (...
- 12k
25
votes
2
answers
11k
views
Libraries for solving sparse linear systems
There are a number of different libraries out there that solve a sparse linear system of equations, however I'm finding it difficult to figure out what the differences are.
As far as I can tell there ...
- 1,155
25
votes
3
answers
917
views
Solving $(G^TA^{-1}G)x = b$ without inverting $A$
I have matrices $A$ and $G$. $A$ is sparse and is $n\times n$ with $n$ very large (can be on the order of several million.) $G$ is an $n\times m$ tall matrix with $m$ rather small ($1 \lt m \lt 1000$) ...
- 1,320
24
votes
4
answers
5k
views
When do orthogonal transformations outperform Gaussian elimination?
As we know, orthogonal transformations methods (Givens rotations and Housholder reflections) for systems of linear equations are more expensive than Gaussian elimination, but theoretically have nicer ...
- 1,832
23
votes
3
answers
3k
views
Can diagonal plus fixed symmetric linear systems be solved in quadratic time after precomputation?
Is there an $O(n^3+n^2 k)$ method to solve $k$ linear systems of the form $(D_i + A) x_i = b_i$ where $A$ is a fixed SPD matrix and $D_i$ are positive diagonal matrices?
For example, if each $D_i$ is ...
- 3,909
21
votes
3
answers
4k
views
What guidelines should I use when searching for good preconditioning methods for a specific problem?
For the solution of large linear systems $Ax=b$ using iterative methods, it is often of interest to introduce preconditioning, e.g. solve instead $M^{-1}(Ax=b)$, where $M$ is here used for left-...
- 3,196
21
votes
1
answer
3k
views
Diagonal update of a symmetric positive definite matrix
$A$ is a $n \times n$ symmetric positive definite (SPD) sparse matrix. $G$ is a sparse diagonal matrix. $n$ is large ($n$ >10000) and the number of nonzeros in the $G$ is usually 100 ~ 1000.
$A$ has ...
luogyong
20
votes
3
answers
11k
views
Null-space of a rectangular dense matrix
Given a dense matrix $$A \in R^{m \times n}, m >> n; max(m) \approx 100000 $$ what is the best way to find its null-space basis within some tolerance $\epsilon$?
Based on that basis can I then ...
- 1,111
19
votes
3
answers
3k
views
Why do we usually not want the eigenvalues of non-symmetric matrices?
I came across this line in a class note I am reading where it discusses finding eigenvalues of matrices.
In reality we don't go all the way with Arnoldi. We stop at a decent value of 𝑘. Then the 𝑘 ...
- 383
18
votes
4
answers
4k
views
Why can't Householder reflections diagonalize a matrix?
When computing the QR factorization in practice, one uses Householder reflections to zero out the lower portion of a matrix. I know that for computing eigenvalues of symmetric matrices, the best you ...
- 4,480
18
votes
6
answers
2k
views
How to reorder variables to produce a banded matrix of minimum bandwidth?
I'm trying to solve a 2D Poisson equation by finite differences. In the process, I obtain a sparse matrix with only $5$ variables in each equation. For example, if the variables were $U$, then the ...
- 12k
18
votes
5
answers
789
views
20% performance penalty for a nice software design
I'm writing a small library for sparse matrix computations as a way to teach myself to make the best use of object-oriented programming. I've worked really hard on having a nice object model, where ...
- 10.1k
18
votes
3
answers
10k
views
Efficient computation of the matrix square root inverse
A common problem in statistics is computing the square root inverse of a symmetric positive definite matrix. What would be the most efficient way of computing this?
I came across some literature (...
- 345
17
votes
2
answers
5k
views
Stopping criteria for iterative linear solvers applied to nearly singular systems
Consider $Ax=b$ with $A$ nearly singular which means there is an eigenvalue $\lambda_0$ of $A$ that is very small. The usual stop criterion of an iterative method is based on the residual $r_n:=b-Ax_n$...
- 1,319
17
votes
1
answer
2k
views
Can a Krylov subspace method be used as a smoother for multigrid?
As far as I am aware, multigrid solvers use iterative smoothers such as Jacobi, Gauss-Seidel, and SOR to dampen the error at various frequencies. Could a Krylov subspace method (like conjugate ...
- 12k
16
votes
4
answers
1k
views
Can the solution of a linear system of equations be approximated for only the first few variables?
I have a linear system of equations of size mxm, where m is large. However, the variables that I'm interested in are just the first n variables (n is small compared to m). Is there a way I can ...
- 12k
16
votes
3
answers
2k
views
multigrid method to solve PDE
I need simple explanation of the Multigrid Method or some literature about this.
I am familiar with iterational methods including BiCGStab,CG,GS,Jacobi and preconditioning, but I am a beginner with ...
- 373
16
votes
3
answers
2k
views
Why would a computational scientist need to implement their own version of std::complex?
Many of the better-known C++ libraries in computational science such as Eigen, Trilinos, and deal.II use the standard C++ template header library object, ...
- 6,931
16
votes
2
answers
986
views
Preconditioning a Krylov method with another Krylov method
In methods like gmres or bicgstab it could be attractive to use another Krylov method as a preconditioner. After all they are easy to implement in a matrix-free way and in a parallel environment. For ...
16
votes
2
answers
2k
views
Estimation of condition numbers for very large matrices
Which approaches are used in practice for estimating the condition number of large sparse matrices?
- 3,196
16
votes
3
answers
16k
views
What is the fastest way to compute all eigenvalues of a very big and sparse adjacency matrix in python?
I'm trying to figure out if there is a faster way to compute all the eigenvalues and eigenvectors of a very big and sparse adjacency matrix than using scipy.sparse.linalg.eigsh As far as I know, this ...
- 263
16
votes
1
answer
5k
views
How can I estimate the condition number of a large sparse matrix using PETSc?
I have a PETSc Mat and would like to estimate its condition number.
- 25.6k
15
votes
2
answers
573
views
Is there any way to do "double preconditioning"
Question:
Suppose that you have two different (factored) preconditioners for a symmetric positive definite matrix $A$:
$$A \approx B^TB$$
and
$$A \approx C^TC,$$
where the inverses of the factors $B, ...
- 3,143
15
votes
2
answers
7k
views
Compute all eigenvalues of a very big and very sparse adjacency matrix
I have two graphs with nearly n~100000 nodes each. In both graphs, each node is connected to exactly 3 other nodes so the adjacency matrix is symmetric and very sparse.
The hard part is I need all ...
- 253
14
votes
3
answers
10k
views
Understanding how Numpy does SVD
I have been using different methods to calculate both the rank of a matrix and the solution of a matrix system of equations. I came across the function linalg.svd. Comparing this to my own effort of ...
- 529
14
votes
6
answers
1k
views
Approximate spectrum of a large matrix
I want to compute the spectrum (all the eigenvalues) of a large sparse matrix (hundreds of thousands of rows). This is hard.
I am willing to settle for an approximation. Are there approximation ...
- 3,068
14
votes
3
answers
2k
views
Why is my parallel solver slower than my sequential solver?
I was playing around with PETSc and noticed that when I run my program with more than one process via MPI it seems to run even slower! How can I check to see what is going on?
- 1,370
14
votes
3
answers
2k
views
SVD for finding the largest eigenvalue of a 50x50 matrix -- am I wasting significant amounts of time?
I've got a program that computes the largest eigenvalue of many real symmetric 50x50 matrices by performing singular-value decompositions on all of them. The SVD is a bottleneck in the program.
Are ...
- 243
14
votes
1
answer
1k
views
How is Krylov-accelerated Multigrid (using MG as a preconditioner) motivated?
Multigrid (MG) may be used to solve a linear system $Ax=b$ by constructing an initial guess $x_0$ and repeating the following for $i=0,1..$ until convergence:
Compute the residual $r_i = b-Ax_i$
...
- 849
14
votes
2
answers
2k
views
How useful is PETSc for Dense Matrices?
Wherever I have seen, PETSc tutorial/documents etc. say that it is useful for linear algebra and usually specifies that sparse systems will benefit. What about dense matrices? I am concerned about ...
- 3,344
13
votes
2
answers
2k
views
Does the "cofactor technique" for inverting a matrix have any practical significance?
The title is the question. This technique involves using the "matrix of cofactors", or "adjugate matrix", and gives explicit formulae for the components of the inverse of a square matrix. It is not ...
- 515
13
votes
3
answers
1k
views
Blaze linear algebra library?
The paper "Expression Templates Revisited: A Performance Analysis of Current Methodologies" in SIAM Journal of Scientific Computing references the "Blaze" linear algebra library. I haven't heard of it ...
- 805
13
votes
3
answers
4k
views
What are the fastest available implementations of BLAS/LAPACK or other linear algebra routines on GPU systems?
nVidia, for example, has CUBLAS, which promises 7-14x speedup. Naively, this is nowhere near the theoretical throughput of any of nVidia's GPU cards. What are the challenges in speeding up linear ...
- 231
13
votes
3
answers
411
views
In what application cases are additive preconditioning schemes superior to multiplicative ones?
In both domain decomposition (DD) and multigrid (MG) methods, one may compose the application of the block updates or coarse corrections as either additive or multiplicative. For pointwise solvers, ...
- 1,665
13
votes
1
answer
2k
views
Purely rotational least squares match
Could anyone recommend a method for the following least-squares problem:
find $R \in \mathbb{R}^{3 \times 3}$ that minimizes: $\sum\limits_{i=0}^N (Rx_i - b_i)^2 \rightarrow \min$, where $R$ is a ...
