why use svd() to invert a matrix?

Question

In MATLAB, i compared elapsed time to invert a Hermitian matrix using inverse(), svd(), and chol(). svd() took the longest. So is there any reason to prefer svd() to the other two methods?

Answer 1

The two methods differ, above all, by their applicability to matrix classes.

col (cholesky) decomposes Hermitian, positive-definite rectangular matrices into the product of a lower triangular matrix and its conjugate transpose;

svd (singular value decomposition) factorizes any m×n matrix into the form UΣV^*, where U and V are square real or compex unitary matrices, m×m and n×n, respectively, and Σ is an m×n rectangular diagonal matrix with non-negative real numbers on the diagonal.

The method inv internally performs an LU decomposition of the input matrix (or an LDL decomposition if the input matrix is Hermitian), but outputs only the inverse of square matrix only.

Both SVD and Cholesky can be used for computing pseudoinverse of a matrix, provided the matrix satisfies requirement for the method used. The pseudoinverse operation is used to solve linear least squares problems and the other signal processing, image processing, and big data problems.

UPDATE on OP's comment

The matrix can be both Hermitian and not a positive/negative (semi)definite, in which case it is called an Hermitian indefinite matrix. The (generalized) Cholesky decomposes Hermitian positive semidefinite matrices. You can also decompose a negative (semi)definite matrix, say, A_negdef, just call chol()) on -A_negdef, but you cannot compute an inverse of an indefinite matrix with Cholesky, because of unavoidable square roots of negative numbers.

Rank-deficient matrices (having zero singular values) are not invertible, as is the case for most general-case matrices. That is, there exists no matrix A^-1 such that A^-1A = AA^-1 = I, for rank-deficient matrices. Still, many problems which you solve through matrix inversion can be solved for indefinite (and consequently non-invertible) matrices with a generalization of matrix inversion, pseudoinverses of matrices.

I won't bore you with the description of Moore-Penrose pseudoinverse concept, on which subject many articles are readily available online and in textbooks: in short, the matrix A⁺ is a pseudoinverse of matrix A, if AA⁺A = A (consequently, A⁺AA⁺ = A⁺). If A has linearly independent columns, the pseudoinverse is a left inverse, because in this case A⁺A = I, you can check it with a sample low-dimension matrix constructed for this purpose. If A has linearly independent rows, the pseudoinverse is a right inverse, because in this case AA⁺ = I. In general case, A may have no left/right inverses, but still have pseudoinverse A⁺ satisfying the equation AA⁺A = A.

Back to quadratic programming problems: Below, all matrices and vectors are real-valued, so instead of Hermitian conjugation H we have a matrix transpose T here. For an invertible (in a "classic" sense) matrix A, a quadratic form $$ f(x) = {1\over2}x^{\textrm{T}}Ax-x^{\textrm{T}}b $$ has a minimum value at x_opt = A^-1b: $$ f(x_{opt}) = -{1\over2}b^{\textrm{T}}A^{\textrm{-1}}b $$

To compute an inverse of A, you use function chol() here.

You can solve a minimization problem for a quadratic form with a non-invertible matrix A, provided A is positive semidefinite, even if A has no inverse in this case. For the singular matrix A, the minimum value is $$ f(x_{opt}) = -{1\over2}b^{\textrm{T}}A^{\textrm{+}}b $$ which is reached for all x given by $$ x = A^{\textrm{+}}b + U^{\textrm{T}}{\left(\begin{align}{0\\z}\end{align}\right)}, \hspace{1cm} z \in \mathbf{R}^{\textrm{n-r}} $$ , n is the dimension of x, r is the rank of A.

Matrix U in this solution is the factor of singular value decomposition of A = UΣV^T, which you compute with the function svd(). So, maybe not all "the intermediate products, whether they be triangular matrices, eigenvectors, diagonal matrices, etc. are irrelevant when the only things that matter are the computation time(FLOPS)"?

In these examples, all matrices and vectors are real-valued, and I am not sure this example can be readily applied to your beamforming problem. Still, SVD is certainly useful in signal processing, e.g., see Image Reconstruction using SVD.

For the sake of completeness: 1) there exists a generalization of Cholesky decomposition useful for computation of a positive semidefinite matrix, and the pseudoinverse can be computed with the help of Greville's algorithm. But this algorithm is implemented through recursion on matrix columns, and the Cholesky speed can be offset by this recursion; 2) the "stability issues" are addressed by the LDL decomposition.

Any way, these algebraic technicalities are not "too much information" for those who work in signal processing: if, even without a preliminary analysis on positive-definiteness of matrices, you start compute using a chol() function and it turns out that you are getting weird results, you now know what might be the culprit.