Listing applications of the SVD

Question

The SVD (singular value decomposition) is taught in many linear algebra courses. It's taken for granted that it's important. I have helped teach a linear algebra course before, and I feel like I need to provide a better motivation for the SVD. This has motivated me to ask what the applications of the SVD might be.

I have considered the following applications, but I am sceptical of some of the ones I list:

• Low-rank approximation in image compression. Aren't images usually compressed using an algorithm like JPEG which is based on the discrete cosine transform (DCT) and not on low-rank approximation? Basically, is this something that's used in competitive image compression algorithms?

• Low-rank approximation in recommender systems. The problem to solve is that you have a matrix $M$ where the columns are users and the rows are products. The entries of the matrix tell you what rating each user gave to each product. Notice that for many use-cases, some of the entries of the matrix will be missing. Picture those missing entries as taking on the special value "?" (that is, a question mark). I am almost certain that you can't use a default value like $0$. Supposedly, given a low-rank approximation $M \approx U\Sigma V^T$, the columns of $V$ represent products, and columns which have a high cosine-similarity represent similar products. "Cosine-similarity" is simply the cosine of the angle between between vectors. The problem is that you can't take the SVD of a matrix with missing values. This point is well known among practitioners who have developed approaches to low-rank approximation of matrices with missing entries that don't directly involve using the SVD (and I'm not sure how they even could use the SVD).

• Finding the nearest orthogonal matrix to a given matrix in machine vision. The objective is, given a matrix $M$, find the orthogonal matrix $\Omega$ for which $\|M - \Omega\|$ is minimised. For some matrix norms, this problem can be solved optimally (in the sense of finding the optimal $\Omega$) using SVD. This is potentially useful in machine vision for turning a matrix which ought to be orthogonal, but which isn't because of various sources of error, into one that is orthogonal. This is helpful for reconstructing the position of a camera in space based on photographs taken with that camera.

• Principal component analysis (PCA). The objective is to reduce the dimensionality of a dataset in order to plot the data or use it in machine learning algorithms that might not work well with high-dimensional data. The idea here is that you imagine that your data is drawn from a multivariate normal distribution. What PCA does is infer the underlying multivariate normal distribution and approximate it using a lower-dimensional multivariate normal. Note that if a dimension gets scaled in the sample set, then the output of PCA can be changed. In many cases in which PCA is used (most cases at the moment?), the dataset is not actually drawn from a normal distribution; and I don't know how these deviations from normality affect the suitability of PCA.

• Linear regression. My understanding is that SVD can be used to find Moore-Penrose pseudoinverses, and Moore-Penrose pseudoinverses can in turn be used to fit linear regression models. It's worth pointing out that linear regression is usually done alongside regularisation in order to prevent overfitting. One form of regularised linear regression is called ridge regression. There is indeed a way to fit a Ridge regression model using SVD. Some properties of Ridge regression can be studied this way.

I'd love to see applications in pure mathematics or physics as well. The above summarises my investigations into SVD thus far.

Answer 1

I teach a course on Applied Linear Algebra, intended for Engineers, where the final project is always to give a presentation on an application of linear algebra in the student's field of study. Since I myself am not an engineer, I often learn new things here, but there is the risk that my level of understanding is no better than an undergraduate engineer.

With this caveat, one of the common topics my students talk about is image identication. You have, say, $M$ reference photos of people's faces, each of them $K \times K$ grey-scale pixels and positioned similarly in the camera frame. You get a new photo and you want to identify which face it is most like.

The naive thing to due is to treat each of the $M$ images as a vector with $K^2$ entries, and find the reference photo in $K^2$ dimensional space which is closest to your photo. The trouble is that you will be thrown off by random noise introduced by the photography process, and you will get a face which is closest to your photo in the directions of random variation, and not in the directions that faces naturally vary.

A better idea is to use SVD on the $M \times K^2$ matrix whose columns are your reference photos, to find a lower dimensional plane $V$ in $K^2$-dimensional space, such that all of your reference photos are close to $V$. Then project all the reference photos, and the input photo, onto $V$, before finding the closest one. This way, the output photo will be close to the input photo in the ways in which faces vary, while the random noise will be projected away.

This is called the eigenfaces method. It can of course apply to images other than faces, such as optical character recognition.

Answer 2

It is a fun little exercise to show using SVD that a linear transformation $A\colon\mathbb{R}^n\to\mathbb{R}^n$ maps ellipses to ellipses (not necessarily centered at the origin).

Given an ellipse $E$, find a linear transformation $B$ such that $B(E)$ is a sphere, say, by rotating $E$ around the origin such that the axes of the ellipse become parallel to the coordinate axes, and then by scaling along them. Then $A=U\Sigma V$ certainly maps the sphere $B(E)$ to an ellipse, because $VB(E)$ is a sphere, so $\Sigma VB(E)$ is an ellipse, as well as $AB(E)=U\Sigma VB(E)$. Now apply this to a linear map $AB^{-1}$.

Answer 3

In quantum physics, one often studies the entanglement between to parts of the system, in terms of the entanglement entropy, which can be expressed in terms of the Schmidt coefficients occurring in the Schmidt decomposition (more or less an SVD) of a pure state. In particular the entanglement entropy $S$ is given by $$ S = - \sum_i |\alpha_i|^2 \log (|\alpha_i|^2) \ , $$ for a pure state $| \psi \rangle$ divided in two parts $A$ and $B$ as $$ | \psi \rangle = \sum_{i} \alpha_i | u_i \rangle_A \otimes | v_i \rangle_B \ , $$ where $|u_i \rangle_A$ and $| v_i \rangle_B$ are orthonormal bases for subsystem $A$ and $B$ respectively.

This answer is heavily based on this wikipedia page.

Answer 4

A huge amount of the work in analyzing measurements of earth's atmosphere in order to assess climate variability consists of clever applications of the SVD. See for example this paper by Coughlin and Tung or these notes from a course taught by Dennis Hartmann (which was my intro to the subject).

Answer 5

In Positron Emission Imaging, scans take long time (>10 mins) and obviously this is affected by patient motion, one of the main factors being respiratory motion. Patients breathe, so they "blur" the image.

A good way to fix this is to "bin" your data in different parts of the respiratory phase, i.e. reconstruct 10 images instead of 1, of the full respiratory cycle. But people breathe non-uniformly in amplitude and period, so we need to find a way to obtain a "respiratory surrogate signal", a plot that tells us for each data point where it was acquired in the respiratory phase.

Turns out that if you bin the data in e.g. 0.5 second pieces and apply Principal Component Analysis (PCA, "same thing" as SVD) to the data, and find the principal component that produces a change in value with a period of a common respiratory signal (3~6s or so) when applied to the data, you can very accurately obtain the respiratory phase.

This (or a variation of this) is nowadays clinically used in PET scanners all across the globe to correct for respiratory motion.

K. Thielemans, S. Rathore, F. Engbrant and P. Razifar, "Device-less gating for PET/CT using PCA," 2011 IEEE Nuclear Science Symposium Conference Record, 2011, pp. 3904-3910, doi: 10.1109/NSSMIC.2011.6153742.

Answer 6

I will assume that using just the singular values (rather then SVD as a whole) counts.

Identify the space of symmetric positive-definite $n\times n$ real matrices with $\operatorname{GL}(n,\mathbb{R})/\operatorname{O}(n)$ by mapping $A\mapsto AA^\top$. Then the heat kernel $K(X,Y)$ with respect to the Riemannian structure carried over from the homogeneous space is in fact a function of the singular values of $XY^{-1}$. The explicit formulas are only available in dimensions $2$ and (less explicit) $3$, see this paper.

Answer 7

In control theory, specifially in $H_2/H_\infty$ or $\mu$ synthesis and analysis, pretty much everything relies on the SVD. For further info, see e.g.

SVD controllers for $H_2$, $H_\infty$ and $\mu$ optimal control

Answer 8

$\newcommand{\RR}{\mathbb{R}}$Filtering an image $u\in\RR^{n\times m}$ by some filter $h\in \RR^{2r+1\times 2s+1}$ means computing $$ \sum_{k=-r}^r\sum_{l=-s}^s u_{i+k,j+l}h_{k,l} $$ at every pixel $(i,j)$. This is basically a convolution (up the reflection of $h$). To compute one pixel of the filtered image, one needs $O((2r+1)(2s+1))$ operations. If the filter is separable, i.e. it is of the form $h=ab^T$ with $a\in\RR^{2r+1}$ and $b\in\RR^{2s+1}$, the filtering can be done by filtering with $a$ first and then with $b$ and the amount of operations is $O((2r+1)+(2s+1))$ which may be much smaller. Of course, a separable filter is just a rank-one matrix, so some software uses the SVD of the filter to check, how good the filter $h$ can be approximated by a rank-one matrix.

One can get a little more: Using the SVD one can write $$ h = \sigma_1 a_1b_i^T + \cdots + \sigma_k a_k b_k^T $$ and represent $h$ as the sum of $k$ separable filters and, depending on $k$, using this representation may still save a few computations (and, of course one may just omit very small singular valued $\sigma_j$ to same more).

Answer 9

The SVD is used to analyze linear regularization methods for linear inverse problems.

Here is very short introduction: A linear inverse problem is the challenge to find a good approximation of a linear operator equation $Ax=b$ when $b$ is only given with an error, i.e. you only get $b^\delta$ with $\|b-b^\delta\|\leq\delta$ for some known $\delta>0$. Here $A$ is a linear and continuous operator between two Hilbert spaces $X$ and $Y$. If $A$ is compact, it has a singular value decomposition $$ Ax = \sum_i \sigma_i \langle u_i,x\rangle v_i $$ where $(u_i)$ is an orthonormal basis of $R(A^T)$, $(v_i)$ is an orthonormal basis of $R(A)$, and $\sigma_i>0$.

In theory, you can use the Moore-Penrose pseudo-inverse to obtain a least-squares minimum-norm solution of $Ax=b^\delta$, i.e. $x = A^\dagger b^\delta$ and this can be computed with the SVD as $$ A^\dagger b = \sum_i \sigma_i^{-1} \langle b,v_i\rangle u_i $$ (that's already a nice application, btw). Alas, this sum may not converge. The reason is, that the pseudo-inverse is usually not a bounded operator (it's only defined on $R(A)\oplus R(A)^T$ and in the case of compact operators, $R(A)$ either finite dimensional or not closed). The SVD characterizes the range of $A$: $y\in R(A)$ holds exactly if the Picard condition is satisfied $$ \sum_i \frac{|\langle y,v_i\rangle|^2}{\sigma_i^2}<\infty. $$

To stabilize the (pseudo-)inversion one can do many things, and the linear methods are usually written as operators $R_\alpha:Y\to X$ defined by $$ R_\alpha y = \sum_i F_\alpha(\sigma_i^2)\sigma_i\langle y,v_i\rangle u_i $$ for $\alpha>0$ (the regularization parameter) and a function $F_\alpha$. It can be shown that if $F_\alpha$ fulfills

1. $F_\alpha(\lambda) \to 1/\lambda$ for $\lambda>0$
2. $\lambda |F_\alpha(\lambda)| \leq C$

then $R_\alpha$ is indeed a regularization of the pseudo-inverse (you can find more in the lecture notes by Christian Clason in Section 5, for example). So, the SVD allows for a quite general theory which the allows to analyze widely used methods such as Tikhonov regularization, the Landweber method, or the truncated singular valued decomposition as special cases of the general theory.

Answer 10

In statistics, the SVD can be used to assess the conditioning of a design matrix, and thus the stability of the parameter coefficient estimates. One key advantage of this approach over the more common Variance Inflation Factors (VIF) is that it also allows detecting collinearity with the intercept column. One disadvantage is that it is much harder to explain to a nontechnical audience.

More information can be found here:

• Belsley, David A.; Kuh, E.; Welsch, R. E. Regression Diagnostics: Identifying Influential Data and Sources of Collinearity (Wiley, 1980)
• Belsley, David A. Conditioning Diagnostics: Collinearity and Weak Data in Regression (Wiley, 1991)
• Belsley, David A. (1991). A Guide to using the collinearity diagnostics. Computer Science in Economics and Management 4, 33-50

Answer 11

See my other answer for the caveat that I am writing about things I have learned from student presentations, so they may be flawed.

Suppose we have a block of raw material which has some shape $X \subset \mathbb{R}^3$. We want to mold it into a new shape $Y \subset \mathbb{R}^3$, but, if the material is stretched too much, it will break. How can we model this?

A possible way of molding $X$ into $Y$ is a smooth bijection $f : X \to Y$. At any point, this map has a Jacobian matrix $Df$. If $Df$ is orthogonal, then we are just rotating the material, so there is no problem at all. How can we numerically impose that $Df$ is nearly orthogonal?

We choose some parameter $\epsilon$, and insist that all singular values of $Df$ are in the range $(1-\epsilon, 1+\epsilon)$. Then we have a problem in PDE's: To map $f : X \to Y$ while imposing that the singular values of $Df$ stay in a bounded range.

Answer 12

Google Scholar turns up a lot of applications in cryptography. One reference that is very accessible to students (but not particularly deep) is "Singular Value Analysis of Cryptograms", Cleve Moler and Donald Morrison, MAA Monthly, Feb. 1983 vol 90., no. 2, pp. 78-87. It's available on JSTOR. http://www.jstor.com/stable/2975804

Answer 13

For an image compression example which can be used directly in teaching, see my answer at What is the intuition behind SVD?.

I hope this example also serves to give some interesting insights into the svd!

Answer 14

SVD is often used to perform tensor decompositions in Tucker and Tensor Train formats. HOSVD (Higher-order SVD) is an algorithm that approximates a given tensor (a multidimensional array of real or complex numbers) with a smaller tensor and a bunch of matrices - this is called the Tucker format. Sometimes HOSVD is used to initialize such a decomposition which is then improved using ALS (Alternating Least Squares) - the resulting algorithm that combines these two parts is called HOOI (Higher Order Orthogonal Iteration). All these three algorithms as well as Tucker format are described in Tensor Decompositions and Applications by T. G. Kolda, B. W. Bader (http://www.kolda.net/publication/koba09/). Also, these algorithms are implemented in http://tensorly.org/stable/modules/generated/tensorly.decomposition.Tucker.html. Some applications of these decompositions are given in the aforementioned paper. Another application I would like to mention is making neural networks work faster by approximating their tensors in Tucker format using these decompositions, see Automated Multi-Stage Compression of Neural Networks by Gusak et al. (https://www.semanticscholar.org/paper/Automated-Multi-Stage-Compression-of-Neural-Gusak-Kholyavchenko/b8bbda9bc5e0861a64a54057af7f6a88b49498c7).

To approximate a tensor using Tensor Train format, one typically uses a bunch of SVDs and nothing else. The format and the algorithm are described in Tensor-train decomposition by Oseledets (http://pitt.edu/~sjh95/related_papers/tensor_train_decomposition.pdf). I don't know much about its applications. I know that quantum physicists call this format Matrix Product State when it's used to represent quantum states and Matrix Product Operator when it's used to represent a Hamiltonian, more can be read about this in Hand-waving and Interpretive Dance: An Introductory Course on Tensor Networks by Jacob C. Bridgeman, Christopher T. Chubb (https://arxiv.org/abs/1603.03039).

Answer 15

It was shown that many quantum algorithms (e.g. linear systems, Hamiltonian simulation, quantum phase estimation) can be thought of essentially as applications of the Quantum Singular Value Transform (QSVT). Even more recently, it was proposed that QSVT could lead to a Grand Unification of Quantum Algorithms. Here is a recent talk given on the latter paper.

Answer 16

A speech recognition practitioner here. SVD is used in certain algorithms, although not very ubiquitous.

I am reminded of the fast speaker adaptation algorithm fMLLR[1], dated, but involving a clever and beautiful shuffling of matrices. SVD may not be the very core of it, but does make an appearance; also, the algorithm as such is indeed inspired by the PCA.

First, instead of a model space transformation, an affine transformation in the feature space is learned. The feature space in speech recognition is typically $\mathbb{R}^D, D=13\ldots 80$, with the frequency around 100 vectors/s. An affine transformation matrix dimensions are $D\times(D+1)$ by combining $\mathbf{W} = [\mathbf{A};\mathbf{b}]$ of the affine $\mathbf{A}\,\mathbf{x}+\mathbf{b}$, with an extra $1$ appended to the feature vector $\mathbf{x}$ to reduce the transform to a single matrix operation.

The transformation is a weighted sum of $N$ basis matrices $\mathbf{W} = \mathbf{d}\,\mathbf{W}_{(N)}$, with the weights $d_i$ adjusted at the decode time, which constitutes the adaptation process proper. The basis size $N$ is all-important: a larger value means more values of the $d_i$ to learn at decode time, and better but slower-converging adaptation. The key insight is to vary the $N$ with the utterance length.

This is done by precomputing a sequence of $D\,(D+1)$ basis matrices $\mathbf{W}_i$, ordered such that for any basis size $N \leq D\,(D+1)$, the prefix of $N$ matrices in such a sequence form a basis maximizing the likelihood of training-time data (in other words, a prefix of any length $N$ gives the "best" $N$ matrices for the task). This sequence is computed at the training time, and is convenient for fast decode-time access. The problem is analytically intractable, and computationally barely so, but under certain assumptions and a second-order Taylor approximation all $\mathbf{W}_i$ can be computed at once relatively quickly with what I think of as a "stacking trick".

All rows of each (yet unknown) $\mathbf{W}_i$ are stacked as the $i$-th column of a new matrix $\mathbf{X}$, which is then computed to minimize a cost function under the constraint that columns of the corresponding matrix $\mathbf{\hat X}$ are orthonormal vectors, where $\mathbf{\hat X}$ is composed of similarly stacked matrices $\mathbf{\hat W}_i$, which in turn are obtained from $\mathbf{W}_i$ each transformed by the action of a certain matrix. The orthonormality constraint translates into a representation of the objective function in the form $\mathrm{tr}(\mathbf{X}^{\mathrm{T}} \mathbf{\hat M}\,\mathbf{X})$, where the hat in $\mathbf{\hat X}$ stands the same matrix action. (I'm skipping on the details and the training of $\mathbf{M}$ because I forgot how it's done it's not related to SVD). The $\mathbf{X}$ is then obtained by eigendecomposition of $\mathbf{\hat M}$, which in this case, thanks to the symmetries of the model and the "hat transform", is the same as SVD and is in practice done as such. Columns of the $\mathbf{X}$ thus are the eigenvectors of $\mathbf{\hat M}$ in the falling eigenvalue order. When "unstacked" by the reversal of stacking, the resulting $\mathbf{W}_i$ naturally form the sequence with the desired prefix basis property!

[1] Povey D., Yao K. A basis representation of constrained MLLR transforms for robust adaptation. Comp. Speech and Lang., 26:35–51, 2012.

Answer 17

This is a very important result in wireless communication. I don't remember the exact details, but it went something like this. Information transmission in a Multi-Input Multi-Output (MIMO) wireless channel can be represented by the equation \begin{align} \mathbf{y}=\mathbf{H}\mathbf{x} + \mathbf{n} \end{align}where $\mathbf{H}$ represents the wireless channel, $\mathbf{x}$ is the input information vector from the transmitter, $\mathbf{n}$ is the random noise in the channel and $\mathbf{y}$ is the vector observed by the receiver. Then, it can be shown that the strategy that achieves the maximum capacity that the channel can support is to transmit the vector $\mathbf{x}=\mathbf{V}\mathbf{u}$ where $\mathbf{V}$ is the right singular matrix of $\mathbf{H}$ and $\mathbf{u}$ is the original information vector that needs to be transmitted. Also, the receiver will recover his information by doing the transformation $$\hat{\mathbf{y}}=\mathbf{U}^H\mathbf{y}=\Sigma\mathbf{u} + \mathbf{n}$$Note that $\mathbf{\Sigma}$ is the diagonal singular value matrix. The original proof can be found in this paper