# Tag Info

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.

11

The standard tolerance for forming a pseudoinverse is to only invert singular values that are at least $\max(m,n) \epsilon \|A\|_2$, where $A$ is $m \times n$, $\epsilon$ is the machine precision, and $\|A\|_2$ coincides with the largest singular value of $A$. With that said, as J.M. mentioned, it is much more stable to avoid forming $A^H A$: First, we ...

11

Augh!! No, no, a thousand times, no! The reason people use SVD is precisely to avoid having to form the cross-product matrix $\mathbf A^\top\mathbf A$, since the formation of this matrix is a nice recipe for forming ill-conditioned linear systems! The decomposition is meant to be applied directly to $\mathbf A$. (See also some of my previous answers.) I ...

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

9

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

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

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

"Randomized algorithms" have recently become quite popular for partial svds. A header only implementation can be downloaded here: http://code.google.com/p/redsvd/ A review of the current methods can be found here: http://arxiv.org/abs/0909.4061 For full svds I am not sure if you can do better than Householder.

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

In theory, you can "square up" $A$ as Jan illustrates in his answer, but you really don't want to do that when actually computing the SVD, as it comes at a significant cost of accuracy. A better approach that leverages standard Lanczos tridiagonalization is to perform the tridiagonalization of the matrix $\left[\begin{smallmatrix} 0 & A \\ A^T & 0 \... 6 If you only want a few singular values/vectors, ARPACK should do the trick. The SVD docs aren't great, and this distribution is more up to date. EDIT: If you want to do this in python, SciPy has a wrapper. Since your matrix is dense, you could try the block sparse row (BSR) format. 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. ... 5 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 $$\| A - MN^T\| \leq \text{factor}\times \sigma_{k+1}(A) := \epsilon$$... 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

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

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

4

First, never calculate an explicit inverse, allthough sometimes you cannot avoid it. Most of the time however, the solution can be restated as multiplication of a vector with an inverse matrix. In that case several methods exist, which don't exlpicitely calculate the inverse. That said, common regularization techniques are truncated singular value ...

4

I'm not an expert in this field, but being your $A$ and $B$ sparse, matlab and LAPACK are not a good choice. For sparse matrices a quick literature search confirmed that algorithms for the computation of a few extremal generalized singular values have been described, as one might expect. (Some random results from google scholar: http://www.win.tue.nl/~...

4

I read here that research was heading towards a inverse-iteration method that guarantees orthogonality with $\mathcal{O}(N)$ complexity. I'd be interested in hearing about that or other advances. (I wanted to just make a few comments since I don't have the time to write out details, but it got rather big for the comment box.) That I believe would be the ...

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

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

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