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.

Filter by
Sorted by
Tagged with
168 votes
8 answers
135k 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 ...
user avatar
57 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 (...
user avatar
  • 3,085
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 ...
user avatar
  • 8,452
40 votes
4 answers
33k 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) ...
user avatar
  • 3,334
37 votes
3 answers
17k 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 ...
user avatar
  • 607
36 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,...
user avatar
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, ...
user avatar
33 votes
6 answers
30k 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 ...
user avatar
  • 501
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 ...
user avatar
30 votes
5 answers
27k 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 ...
user avatar
  • 498
29 votes
11 answers
9k 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 ...
user avatar
  • 956
29 votes
2 answers
13k 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-...
user avatar
  • 3,334
29 votes
4 answers
18k 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 ...
user avatar
28 votes
2 answers
14k 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 ...
user avatar
  • 25.3k
27 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 (...
user avatar
  • 11.8k
25 votes
3 answers
887 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$) ...
user avatar
  • 1,300
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 ...
user avatar
  • 1,842
24 votes
2 answers
10k 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 ...
user avatar
  • 1,145
23 votes
3 answers
2k 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 ...
user avatar
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 ...
user avatar
20 votes
3 answers
3k 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-...
user avatar
19 votes
2 answers
10k 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 ...
user avatar
  • 1,101
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 ...
user avatar
  • 11.8k
18 votes
5 answers
767 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 ...
user avatar
17 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 𝑘 ...
user avatar
17 votes
2 answers
4k 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$...
user avatar
  • 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 ...
user avatar
  • 11.8k
17 votes
3 answers
9k 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 (...
user avatar
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 ...
user avatar
  • 11.8k
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 ...
user avatar
  • 373
16 votes
3 answers
1k 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, ...
user avatar
  • 6,861
16 votes
4 answers
3k 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 ...
user avatar
  • 4,410
16 votes
2 answers
887 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 ...
user avatar
16 votes
1 answer
4k 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.
user avatar
  • 25.3k
15 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?
user avatar
15 votes
2 answers
518 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, ...
user avatar
  • 3,013
15 votes
2 answers
13k 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 ...
user avatar
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 ...
user avatar
  • 3,048
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?
user avatar
  • 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 ...
user avatar
  • 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$ ...
user avatar
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 ...
user avatar
  • 3,334
13 votes
3 answers
9k 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 ...
user avatar
13 votes
5 answers
538 views

Rapidly determining whether or not a dense matrix is of low rank

In a software project that I'm working on, certain computations are vastly easier for dense low-rank matrices. Some problem instances involve dense low-rank matrices, but they're given to me in full, ...
user avatar
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 ...
user avatar
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 ...
user avatar
  • 785
13 votes
3 answers
394 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, ...
user avatar
  • 1,645
13 votes
2 answers
5k 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 ...
user avatar
  • 233
13 votes
1 answer
326 views

Specialized methods for complex symmetric tridiagonal generalized eigenvalue problems

I have to solve generalized eigenvalue problems $Ax = \lambda Bx$ where $A$ and $B$ are both tridiagonal, $B$ is symmetric positive definite and real, but $A$ is only complex symmetric (not definite ...
user avatar
  • 4,410
13 votes
0 answers
637 views

Fast Eigenvalue and SVD Solver for Structured Matrices

I am looking for a fast Eigenvalue and SVD solver for small dense structured matrices (Hankel and Toeplitz). I have searched for efficient implementations in libraries like MKL but I am not able to ...
user avatar

1
2 3 4 5
22