A message from our CEO about the future of Stack Overflow and Stack Exchange. Read now.
19

See https://math.stackexchange.com/questions/861674/decompose-a-2d-arbitrary-transform-into-only-scaling-and-rotation (sorry, I would have put that in a comment but I've registered just to post this so I can't post comments yet). But since I'm writing it as an answer, I'll also write the method: $$E=\frac{m_{00}+m_{11}}{2}; F=\frac{m_{00}-m_{11}}{2}; G=\...


9

@Pedro Gimeno "I doubt it can be any more robust than that." Challenge accepted. I noticed the usual approach is to use trig functions like atan2. Intuitively, there shouldn't be a need to use trig functions. Indeed, all the results end up as sines and cosines of arctans--which can be simplified to algebraic functions. It took quite a while, but I managed ...


8

The GSL has a 2-by-2 SVD solver underlying the QR decomposition part of the main SVD algorithm for gsl_linalg_SV_decomp. See the svdstep.c file and look for the svd2 function. The function has a few special cases, isn't exactly trivial, and looks to be doing several things to be numerically careful (e.g., using hypot to avoid overflows).


8

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{bmatrix} R_1 & R_{12}\\ 0 & R_{22} \end{bmatrix}, $$ with $R_1$ triangular of size $k\times k$, and $R_{22}$ typically not triangular (since we stopped ...


8

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 reversed. The first $P$ swaps the rows, and the last $P$ swaps the columns. We have the Cholesky decomposition $$PAP=LL^T$$ which implies $$A = (PLP)(PLP)^T,$$ ...


7

When we say "numerically robust" we usually mean an algorithm in which we do things like pivoting to avoid error propagation. However, for a 2x2 matrix, you can write the result down in terms of explicit formulas -- i.e., write down formulas for the SVD elements that state the result only in terms of the inputs, rather than in terms of intermediate values ...


7

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 is to use a randomized algorithm to find low-rank approximations. These can, at least in theory, be significantly faster than working on the entire matrix ...


6

It may help to define $N$, the number of discretization points along a 1D edge, and relate it to $n$, the number of unknowns in the system. In 2D on a square grid of points, $n = O(N^2)$. Nested dissection efficiently reduces the sparse system you usually get from discretizations by eliminating levels of "interior" points. The result is a more dense system ...


6

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 latter statement gives a recipe for computing $A$ given $B$. Given $B$, compute its eigendecomposition $B=V\Lambda V^T$, in which $V$ is the orthonormal set of of ...


6

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 done only once for all RHS), and the solution find (for each RHS) is much cheaper. Direct solvers do not suffer from poor conditioning of the matrix or, in ...


6

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 integer. For example, if you were using a higher-order finite difference discretization that used more points to calculate a derivative, you'd get higher-...


6

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 commonly, they are backward stable, and promise that the actual floating-point output will be $(Q+\delta Q)(D'+\delta D')(Q+\delta Q)^\top$, for some small ...


6

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 guess.. $\begin{bmatrix} \mathbf A_{11} & \mathbf a_{21}^T \\ \mathbf a_{21} & \alpha_{22} \end{bmatrix} = \begin{bmatrix} \mathbf L_{11}^T & \...


5

This code is based on Blinn's paper, Ellis paper, SVD lecture, and additional calculations. An algorithm is suitable for regular and singular real matrices. All previous versions works 100% as well as this one. #include <stdio.h> #include <math.h> void svd22(const double a[4], double u[4], double s[2], double v[4]) { s[0] = (sqrt(pow(a[0] -...


5

Meijerink and van der Vorst showed that incomplete Cholesky dose not break down for $M$-matrices. As for finite element mass matrices, you will have to specify a basis to have any hope of making that claim. Given any basis $\Phi = [\phi_0 | \phi_1 | \phi_2 | \phi_3]$ such that the mass matrix $A = \Phi^T \Phi$ is SPD, we can construct a new basis $\hat \Phi ...


5

Your matrix does not have an inverse, which is why dgesv returns an error.


5

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$ and $\sigma_2$ as follows: $A = USV$, which can be expanded like: $ \begin{bmatrix} a & b \\ c & d \end{bmatrix} = \begin{bmatrix} c_1 & s_1 \\ -...


5

Yes, it is possible. By construction, the $L$ matrix is lower triangular with all its diagonal entries equal to one. Therefore, $\det(L) = \det(L^T) = 1$ and, consequently, $\det(L D L^T) = \det(L) \, \det(D) \, \det(L^T) = \det(D)$. From this, you can conclude that $\log(\det(A)) = \log(\det(D)) = \log(\prod_{i = 1}^n D_{ii}) = \sum_{i = 1}^n \log(D_{ii})$...


5

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}\right] $$ Since the matrices are small and very general (do not feature any known structure, zeroes, relative scales of the elements), I think it would be ...


5

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 analysis phase to recognize consecutive columns that share the same nonzero structure, and reorder/aggregate them into one "supervariable" that will be ...


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


5

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 finding a Cholesky factorization somehow, you should look at modified Cholesky factorization algorithms that perturb the covariance as little as possible to make ...


5

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 enough that they can be replaced by other software, such as ATLAS, MKL, and so on. Another way of looking at this is that any software that does linear ...


4

The structure of $C_i$ gives a hope that $C_i$ can be described as a low-rank update. If that is the case (size of the square block is much smaller than the size of the matrix $B$), one can apply classical methods for low-rank updates based on Sherman-Morrisson-Woodbury identity. While updating the LU factorization is not the easiest task (the required ...


4

Consider what happens if you apply modified / classical Gram-Schmidt to the matrix $$ \begin{pmatrix} 1 & 1 & 1 \\ 1 & 1 & 1 \\ 0 & \varepsilon & 0 \\ 0 & 0 & \varepsilon \\ \end{pmatrix} $$ where $\varepsilon \ll 1$. The two algorithms agree on the first and second column and produce a $Q$ factor whose entries are given by ...


4

Since a SPD matrix is invertible, we can make the Cholesky decomposition $A^{-1} = PP^T$. Since $A$ is non-singular, so is $P$, and the inverse of a triangular matrix is triangular, so writing $L = P^{-1}$ we have $A^{-1}$ = $L^{-1}L^{-T}$. Inverting both sides gives $A = L^TL$.


3

It's unlikely. Demmel and co-authors have a paper called "Fast Linear Algebra is Stable" that shows that numerically stable algorithms for eigendecompositions (up to finite precision) cost at least as much as matrix multiplication, which is $O(n^{\omega})$, where $2 \leq \omega \leq 2.376$ (or so). Methods for transforming an LDL or Cholesky decomposition ...


3

The MATLAB \ operator uses UMFPACK when the input matrix is sparse, square, and unsymmetric. UMFPACK is not a parallel sparse solver. However, the UMFPACK factorization step calls BLAS level-3 routines in its computational kernel. Most modern implementations of the BLAS (e.g. Intel MKL in MATLAB) support multicore architectures so UMFPACK can obtain fairly ...


3

I'm going through some of my old StackExchange posts and came across this one. As it turns out, the answer led to a section in a published paper! As detailed in that paper (and in its notation), if you wish to minimize $\mathrm{KL}(B|AA^\top)$, the following iterations will decrease the objective---and in practice work well--- $$ \begin{array}{rl} U_k&\...


3

I have used the description at http://www.lucidarme.me/?p=4624 to create this C++ code. The Matrices are those of the Eigen library, but you can easily create your own data structure from this example: $A=U\Sigma V^T$ #include <cmath> #include <Eigen/Core> using namespace Eigen; Matrix2d A; // ... fill A double a = A(0,0); double b = A(0,1); ...


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