Singular value decomposition (SVD) of matrices and its applications

Singular Value Decomposition (SVD) is a decomposition method for matrices, which, unlike eigenvalue decomposition, can be applied to rectangular matrices and decompose them into the product of three special matrices. In addition, SVD decomposition has many excellent properties and has a wide range of applications in algorithms such as image compression, data dimensionality reduction and recommender systems.

Introduction of singular value decomposition

We consider the two-dimensional case by considering a set of unit orthogonal vectors on a two-dimensional space\(\boldsymbol{v}_1 , \boldsymbol{v}_2\) Let any one of the transformation matrices\(M \in \mathbb R ^ {m \times 2}\) which is transformed to get another set of orthogonal vectors\(M\boldsymbol{v}_1, M\boldsymbol{v}_2\) It is easy to know that the transformed orthogonal vectors are still a set of bases on this 2D plane, which can be unitized to get\(\boldsymbol{u}_1, \boldsymbol{u}_2\), i.e., the vectors before and after unitization have a scaling relationship:

\[\begin{cases} \begin{aligned} M\boldsymbol{v}_1 &= \sigma_1 \boldsymbol{u}_1 \\ M\boldsymbol{v}_2 &= \sigma_2 \boldsymbol{u}_2 \\ \end{aligned} \end{cases} \]

Moreover, for any vector in this two-dimensional plane\(\boldsymbol{x} \in \mathbb R^2\) The basement can be used in the\((\boldsymbol{v}_1, \boldsymbol{v}_2)\) The lower linear representation is

\[\boldsymbol{x} = (\boldsymbol{v}_1 \cdot \boldsymbol{x}) \boldsymbol{v}_1 + (\boldsymbol{v}_2 \cdot \boldsymbol{x}) \boldsymbol{v}_2 \]

Applying the above linear transformation to this vector yields

\[\begin{aligned} M\boldsymbol{x} &= (\boldsymbol{v}_1 \cdot \boldsymbol{x}) M \boldsymbol{v}_1 + (\boldsymbol{v}_2 \cdot \boldsymbol{x}) M \boldsymbol{v}_2 \\ &= (\boldsymbol{v}_1 \cdot \boldsymbol{x}) \sigma_1 \boldsymbol{u}_1 + (\boldsymbol{v}_2 \cdot \boldsymbol{x}) \sigma_2 \boldsymbol{u}_2 \\ &= \boldsymbol{v}_1 ^ \text T \boldsymbol{x} \ \sigma_1 \boldsymbol{u}_1 + \boldsymbol{v}_2 ^ \text T \boldsymbol{x}\ \sigma_2 \boldsymbol{u}_2 \\ &= \boldsymbol{u}_1 \sigma_1 \boldsymbol{v}_1 ^ \text T \boldsymbol{x} + \boldsymbol{u}_2 \sigma_2 \boldsymbol{v}_2 ^ \text T \boldsymbol{x} \\ &= \begin{pmatrix} \boldsymbol{u}_1 & \boldsymbol{u}_2 \end{pmatrix} \begin{pmatrix} \sigma_1 & \ \\ \ & \sigma_2 \\ \end{pmatrix} \begin{pmatrix} \boldsymbol{v}_1 ^ \text T \\ \boldsymbol{v}_2 ^ \text T \\ \end{pmatrix} \boldsymbol{x} \end{aligned} \]

If the definition of\(\boldsymbol{U} = \begin{pmatrix} \boldsymbol{u}_1 & \boldsymbol{u}_2 \end{pmatrix}, \mathbf \Sigma = \begin{pmatrix} \sigma_1 & 0 \\ 0 & \sigma_2 \\ \end{pmatrix}, \boldsymbol{V} ^ \text T = \begin{pmatrix} \boldsymbol{v}_1 ^ \text T \\ \boldsymbol{v}_2 ^ \text T \\ \end{pmatrix}\) , then for any transformation matrix\(M \in \mathbb R ^ {2 \times 2}\) , both of which can be decomposed into the form of the product of the following two sets of unit orthogonal bases and diagonal matrices

\[M=\boldsymbol{U} \mathbf \Sigma \boldsymbol{V} ^ \text T \]

In general, for any transformation matrix\(A \in \mathbb R ^{m \times n}\) , all of which exist as unit orthogonal arrays\(\boldsymbol{U} \in \mathbb R ^ {m \times m}， \boldsymbol{V} \in \mathbb R ^ {n \times n}\)Class Diagonal Matrix\(\mathbf \Sigma \in \mathbb R ^{m \times n}\) fulfillment

\[A = \boldsymbol{U} \mathbf \Sigma \boldsymbol{V} ^ \text T \]

The above equation is the singular value decomposition of the matrix, where the class diagonal matrix\(\mathbf \Sigma\) The non-zero diagonal elements in are called matrices\(A\) The singular values of the matrix\(U,V\) are called the left and right singular matrices (You matrices) respectively

Calculation of singular value decomposition

Consider below how to obtain the singular value decomposition of a rectangular matrix and give a specific singular value decomposition arithmetic example.

For the matrix\(A \in \mathbb R ^ {m \times n}\) We examine the form of their squares\(A ^ \text T A\) cap (a poem)\(AA ^ \text T\) . As a result of\(A^\text T A\) be\(n \times n\) of real pairs of Chan square matrices, so that they can be eigenvalued decomposed.

\[A ^ \text T A = Q \Lambda Q ^ \text T \]

At the same time, the matrix\(A\) There exists a singular value decomposition, i.e.

\[\begin{aligned} A ^ \text T A &= (\boldsymbol{U} \mathbf \Sigma \boldsymbol{V} ^ \text T) ^ \text T \boldsymbol{U} \mathbf \Sigma \boldsymbol{V} ^ \text T \\ &= V \Sigma ^ \text T U ^ \text T U \Sigma V ^ \text T \\ &= V \Sigma ^ \text T \Sigma V ^ \text T \\ &= V \Sigma ^ 2 V ^ \text T \end{aligned} \]

Compare the results of the two equations above

\[V = Q \\ \Sigma ^ 2 = \Lambda \]

It can be found that the right singular matrix i.e.\(A^ \text T A\) of the eigenvector matrix, and the so-called singular values of the\(A^\text T A\) The arithmetic square root of the eigenvalue of the
The same reasoning can be re-examined\(AA^\text T\) namely

\[\begin{aligned} AA^\text T &= (U\Sigma V^\text T)(U\Sigma V^\text T)^\text T \\ &= U\Sigma V^\text T V \Sigma ^\text T U ^ \text T \\ &= U \Sigma \Sigma ^ \text T U ^ \text T \end{aligned} \]

It can be found that the left singular matrix i.e.\(AA^\text T\) of the eigenvector matrix. Therefore, it is required to solve the matrix\(A\) The singular value decomposition is simply computed separately for the\(A^\text T A\) cap (a poem)\(AA^\text T\) The eigenvalue decomposition is sufficient, and is further illustrated by a specific arithmetic example below.
(Example) Calculate the matrix\(A = \begin{pmatrix} 0 & 1 \\ 1 & 1 \\ 1 & 0 \\ \end{pmatrix}\)The SVD decomposition.
Solution: Examination\(A^\text T A\) cap (a poem)\(AA^\text T\) ，

\[A^\text T A= \begin{pmatrix} 0 & 1 & 1 \\ 1 & 1 & 0 \\ \end{pmatrix} \begin{pmatrix} 0 & 1 \\ 1 & 1 \\ 1 & 0 \\ \end{pmatrix} = \begin{pmatrix} 2 & 1 \\ 1 & 2 \end{pmatrix} \\ \]

\[A A^\text T= \begin{pmatrix} 0 & 1 \\ 1 & 1 \\ 1 & 0 \\ \end{pmatrix} \begin{pmatrix} 0 & 1 & 1 \\ 1 & 1 & 0 \\ \end{pmatrix} = \begin{pmatrix} 1 & 1 & 0 \\ 1 & 2 & 1 \\ 0 & 1 & 1 \\ \end{pmatrix} \\ \]

Make eigenvalue decomposition and unitize eigenvectors

\[A^\text T A= \begin{pmatrix} \frac{1}{\sqrt 2} & \frac{1}{\sqrt 2} \\ \frac{-1}{\sqrt 2} & \frac{1}{\sqrt 2} \\ \end{pmatrix} \begin{pmatrix} 3 & 0 \\ 0 & 1 \\ \end{pmatrix} \begin{pmatrix} \frac{1}{\sqrt 2} & \frac{-1}{\sqrt 2} \\ \frac{1}{\sqrt 2} & \frac{1}{\sqrt 2} \\ \end{pmatrix} \\ \]

\[A A^ \text T= \begin{pmatrix} \frac{1}{\sqrt 6} & \frac{-1}{\sqrt 2} & \frac{1}{\sqrt 3} \\ \frac{2}{\sqrt 6} & 0 & \frac{-1}{\sqrt 3} \\ \frac{1}{\sqrt 6} & \frac{1}{\sqrt 2} & \frac{1}{\sqrt 3} \\ \end{pmatrix} \begin{pmatrix} 3 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 0 \\ \end{pmatrix} \begin{pmatrix} \frac{1}{\sqrt 6} & \frac{2}{\sqrt 6} & \frac{1}{\sqrt 6} \\ \frac{-1}{\sqrt 2} & 0 & \frac{1}{\sqrt 2} \\ \frac{1}{\sqrt 3} & \frac{-1}{\sqrt 3} & \frac{1}{\sqrt 3} \\ \end{pmatrix} \]

Accordingly, the singular value can be obtained\(\sigma_1 = \sqrt 3, \sigma_2 = 1\) , then the SVD decomposition is

\[A = \begin{pmatrix} \frac{1}{\sqrt 6} & \frac{-1}{\sqrt 2} & \frac{1}{\sqrt 3} \\ \frac{2}{\sqrt 6} & 0 & \frac{-1}{\sqrt 3} \\ \frac{1}{\sqrt 6} & \frac{1}{\sqrt 2} & \frac{1}{\sqrt 3} \\ \end{pmatrix} \begin{pmatrix} \sqrt 3 & 0 \\ 0 & 1 \\ 0 & 0 \\ \end{pmatrix} \begin{pmatrix} \frac{1}{\sqrt 2} & \frac{-1}{\sqrt 2} \\ \frac{1}{\sqrt 2} & \frac{1}{\sqrt 2} \\ \end{pmatrix} \\ \]

One point to note is that when solving the singular value or eigenvalue matrix, in order to facilitate the subsequent processing, we generally arrange the eigenvalues or singular values in descending order of size, and the eigenvectors are generally normalized.

Geometric significance of singular value decomposition

A unit orthogonal matrix can be viewed as a rotation matrix in space, while a diagonal matrix represents a telescoping transformation in the direction of the coordinate axes, and thus for a matrix, it can also be viewed as a linear transformation thatThe so-called singular value decomposition is the process of decomposing this transformation into two rotational transformations and one stretching transformation along the coordinate axes, a more intuitive transformation process is shown below:

Applications of singular value decomposition

Singular value decomposition has a wide range of applications in algorithms such as plane fitting, image compression and noise reduction, principal component analysis and recommender systems, and is illustrated below as an example in linear fitting and image compression.

Linear fitting problem

Let the equation of the line be\(ax+by+c=0\) The collection of a set of point sets\((x_i, y_i)\) , the equation of the line needs to be fitted to this set of point sets. To solve this problem, the error function can be constructed

\[e_i = ax_i + by_i +c \]

The fit is considered best when the error is minimized at all points. This is essentially a linear least squares problem and its mathematical form can be expressed as

\[\begin{aligned} (a,b)^* &= \arg \min \sum_{i=0}^{N}{e_i}^2 \\ &= \arg \min \sum_{i=0}^{N}{(ax_i + by_i + c)^ 2} \end{aligned} \]

construction matrix\(A = \begin{pmatrix} x_0 & y_0 & 1 \\ x_1 & y_1 & 1 \\ \vdots & \vdots & \vdots \\ x_N & y_N & 1 \\ \end{pmatrix}, \boldsymbol{x} = \begin{pmatrix} a \\ b \\ c \\ \end{pmatrix}\) The above solution process can then be written in the matrix form

\[\begin{aligned} (a,b)^* &= \arg \min (A\boldsymbol{x})^\text T(A\boldsymbol{x}) \\ &= \arg \min \boldsymbol{x}^\text T (A^\text T A)\boldsymbol{x} \\ &= \arg \min \boldsymbol{x} ^ \text T V \Sigma^2 V^\text T \boldsymbol{x} \end{aligned} \]

perceive\(V^\text T \boldsymbol{x}\) Actually, it's a vector.\(\boldsymbol{x}\) basally\(V\) The set of coordinates under can be denoted as\(\bm \alpha\) namely

\[V^\text T \boldsymbol{x} = \begin{pmatrix} \boldsymbol{v}_1^\text T \\ \boldsymbol{v}_2^\text T \\ \vdots \\ \boldsymbol{v}_N^\text T \\ \end{pmatrix} \boldsymbol{x} = \begin{pmatrix} \alpha_1 \\ \alpha_2 \\ \vdots \\ \alpha_N \\ \end{pmatrix} = \bm{\alpha} \]

Then the above equation can be further written as

\[\begin{aligned} \boldsymbol{x}^\text T V \Sigma^2 V^\text T \boldsymbol{x} &= \bm \alpha^\text T \Sigma^2 \bm \alpha \\ &= \sigma_1^2 \alpha_1^2 + \sigma_2^2 \alpha_2^2 + \cdots + \sigma_N^2 \alpha_N^2 \\ &\ge \sigma_N^2 \end{aligned} \]

Under the premise of considering the normalized fitted solution, i.e.\(\bm{\left| x \right|} = 1\) when and only when the coordinates\(\bm \alpha\) get\(\begin{pmatrix} 0 & 0 & \cdots & 1 \end{pmatrix}^\text T\) The equal sign is taken when solving the following equation

\[\begin{pmatrix} \boldsymbol{v}_1^\text T \\ \boldsymbol{v}_2^\text T \\ \vdots \\ \boldsymbol{v}_N^\text T \\ \end{pmatrix} \boldsymbol x = \begin{pmatrix} 0 \\ 0 \\ \vdots \\ 1 \\ \end{pmatrix} \]

Easy to get, at this point\(\boldsymbol x = \boldsymbol{v}_N\)That is, the solution error is minimized. That is, to require the minimum of the solution error, it is only necessary to perform a singular value decomposition, taking the last dimension of the right singular matrix as the solution\(\bm{x}\) is sufficient, at which point the corresponding fitting error is\(\sigma_N^2\) , which is also the minimum error in the least squares sense.
The above conclusions are further illustrated below by a programming example of straight line fitting.
(Programming Example) Consider a straight line\(x + 2y + 5 = 0\) , given a series of point sets with noise, a straight line is fitted to this set of point sets and compared to the true value to obtain the mean square error of the fitted parameter.
A program is written in C++ to call the Eigen linear algebra library to perform the SVD decomposition, and the code is as follows:

Eigen::Vector3d abc = Eigen::Vector3d(1, 2, 5).normalized(); // straight line parameter truth value
const double sigma = 0.6; // noise variance
const size_t count = 15; // number of points to sample
// Use OpenCV's random number seed to generate the sample points.
cv::RNG rng.
Eigen::MatrixXd A(count, 3); for(size_t i=0.5); // Sample points using OpenCV's random number seed.
for(size_t i=0; i<count; ++i) {
    (i) << i, -(abc[0]/abc[1])*i - (abc[2]/abc[1]) + (sigma), 1;
}
// SVD solve, and take the last dimension of the V matrix as the fitted solution x
Eigen::VectorXd x = Eigen::JacobiSVD<Eigen::MatrixXd>(A, Eigen::ComputeThinV).matrixV().col(2); // Output the result.
// Output the result of the solution
cout << "==> true value: abc = " << () << endl; // print the true value
cout << "==> svd solving result: x = " << () << endl; // print fit parameters
cout << " ==> error RMS = " << sqrt( (A*x).squaredNorm() / double(count) ) << endl; // print mean square error

Screenshots of the program run results are:

Note that the SVD decomposition takes into account normalized vectors, so that what is obtained is also a normalized fit, differing by a scale factor from the actual linear parameters.
It can be seen that at a noise variance of 0.6 and a sampling point count of 15 points, the mean square error of the fitting results obtained is 0.2, which is basically acceptable, and more reliable fitting results can be obtained with a smaller noise variance.

Image compression issues

Consider a\(m \times n\) matrices\(A\) and its singular value decomposition is given by

\[A=U_{m \times m}\Sigma_{m\times n} V_{n \times n}^\text T \]

We find that the obtained singular values decay faster, with larger singular values for the first few terms and smaller singular values the farther back we go, which means that we can make an analysis by keeping only the first few singular values for the\(A\) To carry out the approximation, this practice can compress the storage space, and has a wide range of applications in the field of image processing. The specific practice is that we can take out the front\(k\) singular values to obtain an approximation of the\(\tilde A\) which\(\frac{k}{\min\{m,n\}}\) is called the compression ratio of the matrix.

\[A \approx \tilde A = U_{m \times k} \Sigma_{k \times k} V_{k \times n}^\text T \]

(Programming Example) Given an image (below), consider compressing it using different compression ratios and analyze the compression results comparatively.

The program is written in C++ language and OpenCV image processing library is called to perform SVD decomposition of image matrix and the program is executed repeatedly using different compression ratios (0.01, 0.05, 0.1, 0.2, 0.3, 0.5) to get the output results. The program code is given below:

// Read image
cv::Mat img_original = cv::imread("../", cv::IMREAD_GRAYSCALE);
cv::Mat img = img_original.clone();
(img, CV_64FC1, 1/255.0);
// PerformSVDdisassembly
cv::Mat U, Vt, W;
cv::SVD::compute(img, W, U, Vt);
// Image compression operation
int rank = min(, );
cv::Mat W_hat = cv::Mat(rank, rank, CV_64FC1, cv::Scalar(0)); // choose an antecedentksingular value matrices of dimension
double compression_ratio = 0.3; // Setting the compression ratio
int k = rank * compression_ratio;
for(size_t i=0; i<k; ++i) {
    W_hat.at<double>(i, i) = <double>(i, 0);
}
cv::Mat img_compression = U * W_hat * Vt; // Calculate the compressed image
// Compressed image output
img_compression.convertTo(img_compression, CV_8UC1, 255.0);
cv::imwrite("./ratio=0_3.png", img_compression);

The output compression result is shown below

Compression ratio = 0.01
Compression ratio = 0.05
Compression ratio = 0.1
Compression ratio = 0.2
Compression ratio = 0.3
Compression ratio = 0.5

It can be seen that the basic outline information of the test image has been roughly highlighted when the compression ratio is 0.05, while the detailed information of the image has been basically retained when the compression ratio is 0.2, which is a better compression efficiency.

bibliography

[1] /publicoutreach/feature-column/fcarc-svd
[2] /lomodays207/article/details/88687126
[3] /pinard/p/
[4] /u012198575/article/details/99548136