Hot answers tagged

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


13

The singular value decomposition for a symmetric matrix $A=A^{T}$ is one and the same as its canonical eigendecomposition (i.e. with an orthonormal matrix-of-eigenvectors), while the same thing for a nonsymmetric matrix $M=U \Sigma V^T$ is just the canonical eigenvalue decomposition for the symmetric matrix $$ H=\begin{bmatrix}0 & M\\ M^{T} & 0 \end{...


12

The LU factors of a sparse matrix are at least somewhat sparse. The $Q$ matrix in QR can also somewhat preserve sparsity, and is typically used when the matrix is very long and skinny. The SVD of a sparse matrix will almost always have fully dense $U$ and $V$ factors, so it destroys any reason to perform the computations treating the matrix sparsely.


10

The problem is called Wahba's problem, one algorithm for it is called Kabsch algorithm, and the later more popular is called Davenport q method. It's apparently used and studied in aeronautics to determine a craft orientation. There are lots of reviews about the methods. Beware that the best fit may include reflection. Kabsch method computes a 3x3 ...


10

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


9

Algorithms for the SVD, as more or less every classical linear algebra algorithm based on orthogonal transformations, are normwise backward stable, i.e., it should be guaranteed that $\frac{\|USV^* - A \|}{\|A\|} = O(u)$, where the norms are Euclidean norms, $u$ is the machine precision, and "$O(u)$" means that the first-order term in $u$ is bounded by a ...


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

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


8

Yes. Christopher Baker has implemented his incremental SVD method in a MATLAB package called IncPACK (archived on GitHub, within the imtsl project). It implements methods that are described in his master's thesis. A brief discussion of why Brand's algorithm tends to accumulate error can be found in a 2012 paper by Baker, et al. A related method by Chahlaoui, ...


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

Think of your matrix as block-diagonal with blocks $$ C_\ell = \begin{pmatrix} 0 & & \\ & D_\ell & \\ & & 0\end{pmatrix}. $$ Then it is clear that $$ C_\ell^{-1} = \begin{pmatrix} 0^{-1} & & \\ & D_\ell^{-1} & \\ & & 0^{-1}\end{pmatrix} $$ where of course $0^{-1}$ is not well defined. But that doesn't ...


6

There are a couple of options available if you want an approximate rank-k factorizations. Strongly rank-revealing QR factorizations Interpolative decomposition (ID) and other randomized techniques. Generally speaking, they provide a factorization of the form \begin{equation}\| A - MN^T\| \leq \text{factor}\times \sigma_{k+1}(A) := \epsilon \end{equation}...


6

The Pauli matrices and the identity matrix form an orthogonal basis of the space of $2\times 2$ matrices, so finding the expansion coefficients amounts to just a projection onto this basis (i.e., you need to form the inner product of your matrix onto each element of this basis). But, for $2\times 2$ matrices, the answer is simple enough to write things down ...


6

No, there is nothing, as far as I know, unless you know approximately the location of these eigenvalues. As for methods that can compute a subset of the spectrum of a matrix, I know only of methods that can produce: eigenvalues from "the extremes of the spectrum", e.g., the ones with the largest absolute value, or with the most negative real/imaginary part. ...


6

Q1: No. Here's a counter-example: >> A = eye(4)*1e-300 A = 1.0e-300 * 1.0000 0 0 0 0 1.0000 0 0 0 0 1.0000 0 0 0 0 1.0000 >> rank(A) ans = 4 >> >> rank([A, ones(4, 1)]) ans = 1 >> so, if the added column ...


6

If you are prepared to go digging around in the fortran code: The SVD algorithm consists of a few parts: Bidiagonalization (usually using Householder reflectors) Use QR shifts to reduce the bidiagonal to diagonal (usually using Givens rotations) Change the sign of any negative diagonal elements and reorder the singular values so the are decreasing. ...


5

Yes. Take a look at the 2002 incremental SVD paper by Matthew Brand. In it, he discusses how to re-weight the SVD to dilute the influence of old data. He also has a 2006 paper on rank-1 updates of SVDs, and I believe downdates are also discussed. The best review-type paper on incremental SVDs I've seen was by CG Baker, Van Dooren, and Gallivan in 2012, which ...


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

Let us begin with the exact singular value decompositions $A=U_{A}S_{A}V_{A}^{T}$, $B=U_{B}S_{B}V_{B}^{T}$, $C=U_{C}S_{C}V_{C}^{T}$, $D=U_{D}S_{D}V_{D}^{T}$. Then $$ M=\underbrace{\begin{bmatrix}U_{A} & 0 & U_{C} & 0\\ 0 & U_{B} & 0 & U_{D} \end{bmatrix}}_{X}\underbrace{\begin{bmatrix}S_{A} & 0 & 0 & 0\\ 0 & S_{B} &...


5

Pseudoinverse can be computed using the SVD $A = USV^\top$ by: $$ A^+ = V\Sigma^+ U^\top $$ where $\Sigma^+$ is formed from $\Sigma$ by taking the reciprocal of all the non-zero elements. With that in mind, you could use MATLAB's svds function as follows: [U,S,V] = svds(A,k); Ainv = V*diag(1./diag(S))*U'; Here k refers to the rank and svds computes only a ...


5

For part 1: To my knowledge, the answer is no. For part 2: This question had several good answers, all of which were negative (there isn't really a faster way). I don't believe there is meaningful new literature on the topic.


4

LAPACK has an implementation of the svd of a 2x2 triangular matrix. It appears to be very robust. The routine is XLASV2. To apply to a regular 2x2 matrix, you can simply apply a single givens rotation from the left/right.


4

One method to compute the svd of a matrix X is to first factor X=QR using the QR decomposition (for stability, use pivoting, so this is [Q,R,E] = qr(X,0) in Matlab), and then compute the svd of R. If the matrix is very rectangular in either, then the most expensive computation is the QR factorization. Thus if you increment your matrix X with another row or ...


4

Here is a method that can handle column additions: http://pcc.byu.edu/resources.html. I updated it to handle row additions: function [Up1,Sp,Vp1] = addblock_svd_update2( Uarg, Sarg, Varg, Aarg, force_orth ) U = Varg; V = Uarg; S = Sarg; A = Aarg'; current_rank = size( U, 2 ); m = U' * A; p = A - U*m; P = orth( p ); P = [ P zeros(size(P,...


4

As I already mentioned in the comments here is a possible answer which is backed by some experiments from AlexE in a further comment. The SVD in MATLAB uses the DGESVD from LAPACK, which is based on the ideas of Gene Golub. The main thing is that is is implemented on matrices in Fortran, i.e. columwise storage. In this way processing values in the same ...


4

There are matrix free algorithms (algorithms that use only matrix-vector multiplications rather than working directly with the entries of the matrix) that can compute approximate values of a few singular values/vectors of a large matrix. Since you want a low rank approximation to this matrix, you could use such an algorithm to find the $k$ largest singular ...


4

Pretty sure that what you want is Non-negative least squares. Plenty of implementations already exist, so you shouldn't have to write any code to get going. https://en.wikipedia.org/wiki/Non-negative_least_squares


4

Nemtsov, Averbuchm, and Schclar's "Matrix compression using the Nyström method" (2016) seems relevant: The Nyström method is routinely used for out-of-sample extension of kernel matrices. We describe how this method can be applied to find the singular value decomposition (SVD) of general matrices and the eigenvalue decomposition (EVD) of square matrices. [...


4

Yes. You can run rank-revealing QR on your matrix $A$, which will stop at step $k$ (hence effectively terminating in $O(mnk)$) and produce $A = QRP$, where $R$ has nonzeros only in its first $k$ rows, and $Q,P$ are orthogonal. You can now compute and SVD of $R$, and use it to piece back the factors with a few matrix products with cost $O(\max(m,n)k^2)$.


3

Just collecting the answers contained in the comments: The left and right singular vectors of $A$ are exactly the left and right kernel vectors of $A$. Since the assumption is that the kernels have dimension $1$ and $A$ is well conditioned outside the kernel, one can compute the kernel vectors using the usual, pivoted, matrix factorizations. For the LU ...


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