Skip to main content
11 votes
Accepted

Computational method to compute both the (log) determinant and inverse of a matrix

The LU decomposition will give you what you want with only $\tfrac{2}{3}n^3 + \mathcal{O}(n^2)$ FLOPs. The linear system is solved by solving two triangular systems. The determinant is the product ...
Neil Lindquist's user avatar
9 votes

Robust algorithm for $2 \times 2$ SVD

I needed an algorithm that has little branching (hopefully CMOVs) no trigonometric function calls high numerical accuracy even with 32 bit floats We want to calculate $c_1, s_1, c_2, s_2, \sigma_1$ ...
petiaccja's user avatar
8 votes
Accepted

Solving $A=B+AB$ without matrix inverse

If $A$ and $B$ are real symmetric, then $A=B+AB$ if and only if the product $AB$ is also real symmetric. In turn, $AB=BA$ holds if any only if $A$ and $B$ share a common eigendecomposition. This ...
Richard Zhang's user avatar
8 votes

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

The problem, of course, is that computing the true rank (e.g., via a QR decomposition) is not really any cheaper than computing a low-rank representation of the matrix. The best you can probably do ...
Wolfgang Bangerth's user avatar
8 votes
Accepted

Cholesky for ill-conditioned/singular covariance matrices

If your covariance matrix is singular, then you really should consider why the matrix is singular and come up with a higher-level approach that avoids the singularity. However, if you insist on ...
Brian Borchers's user avatar
8 votes
Accepted

Does a symmetric positive definite matrix also have $\mathbf{A} = \mathbf{L}^T\mathbf{L}$ (where $\mathbf{L}$ is a lower triangular matrix)?

Let $P$ be the anti-diagonal permutation matrix, $$P = \begin{bmatrix} & & & 1 \\ & & 1 \\ & 1 \\ 1 \end{bmatrix}$$ so that $PAP$ is the version of $A$ with rows and columns ...
Nick Alger's user avatar
  • 3,143
8 votes
Accepted

Bareiss algorithm vs. LU-decomposition

Bareiss' algorithm is a better choice if you have to compute exactly determinants of integer matrices, as noted in the comments. When it comes to real (floating-point) or complex matrices, the main ...
Federico Poloni's user avatar
7 votes

Supernodal vs multifrontal factorizations

The scopes of these terms are (in my opinion) somewhat different. I'd use "supernodal" as a broad term, to describe any sparse direct solver that applies sufficient intelligence during the symbolic ...
rchilton1980's user avatar
  • 4,862
7 votes
Accepted

Term for the typical "linear in the larger dimension, quadratic in the smaller" cost for linear algebra

The book "Introduction to Applied Linear Algebra" by Boyd and Vandenberghe has an appendix about complexity of basic operations in linear algebra and they call this case big-times-small-...
Dirk's user avatar
  • 1,738
7 votes

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

There is a neat trick I have recently learned from this paper. You start doing rank-revealing QR, and stop after the first $k$ Householder reflections, when you have a matrix of the form $$ \begin{...
Federico Poloni's user avatar
7 votes
Accepted

Numerically find the nearest positive semi definite matrix to a symmetric matrix

After you compute $Q$ and $D$, form $D'=\max(D,0)$, and compute $A'=QD'Q^\top$, the algorithms involved in multiplying those matrices do not promise that $A'$ will be exactly $QD'Q^\top$. Most ...
Kirill's user avatar
  • 11.4k
7 votes

Givens rotation vs 2x2 Householder reflection

This is an old question, but I want to add some notes: The following c code applies a rotation: ...
Thijs Steel's user avatar
  • 1,693
7 votes
Accepted

Do most statistical packages and libraries in high-level programming languages rely on LAPACK for their matrix inversion operations?

tl;dr Yes. But your question doesn't make it clear that you understand what LAPACK is about. LAPACK is both a software as well as an interface. That is, the operations that LAPACK defines are standard ...
Richard's user avatar
  • 3,961
7 votes
Accepted

Algorithm to factorize matrix whose many rows are already of upper triangular form?

I believe you can accomplish what you want efficiently using the recursive LU algorithm. In brief, recursive LU on a $M \times N$ matrix $A$ proceeds by partitioning the matrix into 4 blocks: \begin{...
vibe's user avatar
  • 1,058
7 votes
Accepted

Functions from Scipy, Blas, or Lapack that compute only upper triangular matrix

I think you are overestimating the overhead of computing L. There are zero extra operations needed; the only additional cost is writing to RAM some numbers that you ...
Federico Poloni's user avatar
7 votes
Accepted

Accurate Way to Calculate Matrix Powers and Matrix Exponential for Sparse Positive Semidefinite Matrices

This is a slightly modified version of my response on math.stackexchange. One standard approach to computing matrix functions times a vector $f(M)x$ or quadratic forms $x^Tf(M)x$ when $M$ is symmetric ...
overfull hbox's user avatar
7 votes

QR decomposition with only diagonal elements changing

Variants of this question get asked regularly on [scicomp.se], but unfortunately the answer is always no: one cannot update easily any matrix factorizations after a full-rank diagonal update. Some ...
Federico Poloni's user avatar
6 votes
Accepted

Fastest way to solve a sparse unsymmetric system many times

Generally speaking, for many right-hand side (RHS) problems, a direct solver is a more feasible solution for several reasons: Major computations are performed during the factorization step (which is ...
Anton Menshov's user avatar
  • 8,672
6 votes
Accepted

Ways to solve $Ax=b$ for a sparse (banded) $A$ with updates

If the only non-zero entries of $A_{ij}$ have $j$ in $\{i - 1, i, i + 1\}$, then $A$ is a banded matrix with bandwidth 1. More generally, you can talk about matrices of bandwidth $k$ where $k$ is any ...
Daniel Shapero's user avatar
6 votes

Does a symmetric positive definite matrix also have $\mathbf{A} = \mathbf{L}^T\mathbf{L}$ (where $\mathbf{L}$ is a lower triangular matrix)?

Yes, for an SPD matrix $\mathbf A$ there are a variety of Cholesky-like decompositions, you can derive the $\mathbf A = \mathbf L^T \mathbf L$ variant by first writing down an educated/structured ...
rchilton1980's user avatar
  • 4,862
6 votes
Accepted

Block-matrix: optimal fill-in reduction for LU factorization

Unfortunately, optimal reordering is NP-complete, so all the algorithms in this space are kinda heuristic. At first glance my instinct would be to try some variant of minimum degree, mainly because ...
rchilton1980's user avatar
  • 4,862
6 votes
Accepted

Numerical Linear Algebra: When to use Direct methods versus iterative methods to solve a linear system - for PDEs in particular

It's a complicated question, which is why I've recorded a whole bunch of video lectures on the topic :-) Take a look at lectures 34 and following here: https://www.math.colostate.edu/~bangerth/videos....
Wolfgang Bangerth's user avatar
6 votes
Accepted

Algorithm for solving systems which are nearly symmetric/adjoint?

Let's express the matrix $A \in \mathbb{R}^{n \times n}$ with which we want to solve linear systems as $$ A = S + U V $$ where $S$ is a symmetric matrix, $U \in \mathbb{R}^{n \times r}$, and $V \in \...
Steven Roberts's user avatar
6 votes

Is there a way to generate a matrix-free decomposition for a matrix-free operator?

This is not efficient, primarily because (i) the factors of sparse matrices are, in general, not sparse themselves, and (ii) when computing matrix factorizations, you don't just need $A_{ij}$ but you ...
Wolfgang Bangerth's user avatar
6 votes

How to efficient solve $e^{-tA} x =b$, where A is a very sparse matrix

$\newcommand{\bm}[1]{\boldsymbol{#1}}$Continuous Formulation With the new edits in your post I believe this is a diffusion problem with a time-dependent inhomogeneity on a domain with reflecting ...
lightxbulb's user avatar
  • 2,122
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

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

Another approach, which might be of interest to you is randomized sampling. This is of particular interest if you can quickly compute matrix-vector products $x\rightarrow Ax$ and $x\rightarrow A^* x$. ...
bobo's user avatar
  • 91
5 votes

Rank-one Update to a Rank Revealing QR (RRQR) Factorization?

I don't know of a usable result for the pivoted QR decomposition, but a related class of matrix decompositions are the $UTV$ decompositions with $$ A = U T V $$ with $U$ and $V$ being orthogonal ...
Florian Goth's user avatar
5 votes

The fast, and The Backward-Stable (left) $3\times 3$ matrix inverse

I will try to give my thought on the first question regarding fast $3\times 3$ inverse. Consider $$ A=\left[ \begin{array}{ccc} a & d & g\\ b & e & h\\ c & f & i \end{array}\...
Anton Menshov's user avatar
  • 8,672
5 votes

Factorize laplacian in terms of first derivative matrix

It is sufficient if you consider a $D$ that uses forward or backward differences with reflecting boundaries: \begin{equation} D_f = \frac{1}{h}\begin{bmatrix} -1 & 1 & & \\ & \ddots &...
lightxbulb's user avatar
  • 2,122

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