奇異值分解（Singular Value Decomposition, SVD）是矩陣的一種分解方法，與特徵值分解不同，它可以應用於長方矩陣中，並將其分解成三個特殊矩陣的乘積。此外SVD分解還具有許多優良的性質，在影像壓縮，資料降維和推薦系統等演算法中都具有廣泛的應用。
奇異值分解的引入
我們考慮二維的情形，考慮一組二維空間上的單位正交向量 \(\boldsymbol{v}_1 , \boldsymbol{v}_2\) ，設任意一個變換矩陣 \(M \in \mathbb R ^ {m \times 2}\) ，對其作變換得到另外一組正交向量 \(M\boldsymbol{v}_1, M\boldsymbol{v}_2\) ，容易知道變換後的正交向量仍然是該二維平面上的一組基底，可以對這組基底進行單位化得到 \(\boldsymbol{u}_1, \boldsymbol{u}_2\)，即單位化前後的向量存在伸縮關係：
此外，對於該二維平面上的任意一個向量 \(\boldsymbol{x} \in \mathbb R^2\) ，可以在基底 \((\boldsymbol{v}_1, \boldsymbol{v}_2)\) 下線性表示為
對該向量作上述線性變換，可以得到
若定義 \(\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}\) ，則對於任一變換矩陣 \(M \in \mathbb R ^ {2 \times 2}\) ，都可以分解成如下兩組單位正交基底與對角矩陣的乘積的形式
一般的，對於任意變換矩陣 \(A \in \mathbb R ^{m \times n}\) ，都存在單位正交陣 \(\boldsymbol{U} \in \mathbb R ^ {m \times m}， \boldsymbol{V} \in \mathbb R ^ {n \times n}\)，類對角矩陣 \(\mathbf \Sigma \in \mathbb R ^{m \times n}\) ，滿足
上式即為矩陣的奇異值分解，其中類對角矩陣 \(\mathbf \Sigma\) 中的非零對角元素稱為矩陣 \(A\) 的奇異值，矩陣 \(U,V\) 分別叫做左奇異矩陣和右奇異矩陣（酉矩陣）
奇異值分解的計算
下面考慮如何獲取一個長方矩陣的奇異值分解，並給出一個具體的奇異值分解算例。
對於矩陣 \(A \in \mathbb R ^ {m \times n}\) ，我們考察其方陣形式 \(A ^ \text T A\) 和 \(AA ^ \text T\) 。由於 \(A^\text T A\) 是 \(n \times n\) 的實對陳方陣，因此可以對其進行特徵值分解。
同時，矩陣 \(A\) 存在奇異值分解，即
對比上面兩式的結果
可以發現右奇異矩陣即為 \(A^ \text T A\) 的特徵向量矩陣，而所謂的奇異值則為 \(A^\text T A\) 的特徵值的算術平方根。
同理可以再考察 \(AA^\text T\) ，即
可以發現左奇異矩陣即為 \(AA^\text T\) 的特徵向量矩陣。因此，要求解矩陣 \(A\) 的奇異值分解，只需分別計算 \(A^\text T A\) 和 \(AA^\text T\) 的特徵值分解即可，下面透過一個具體的算例進一步說明。
（算例）計算矩陣 \(A =
\begin{pmatrix}
0 & 1 \\
1 & 1 \\
1 & 0 \\
\end{pmatrix}\)的 SVD 分解。
解：考察 \(A^\text T A\) 和 \(AA^\text T\) ，
作特徵值分解並單位化特徵向量
據此可以得到奇異值 \(\sigma_1 = \sqrt 3, \sigma_2 = 1\) ，則SVD分解為
需要注意的一點是，在求解奇異值或者特徵值矩陣時，為了方便後續處理，我們一般將特徵值或奇異值按大小降序排列，且特徵向量一般都做歸一化處理。
奇異值分解的幾何意義
單位正交矩陣可以看作空間中的旋轉矩陣，而對角矩陣則表示沿座標軸方向上的伸縮變換，因此對於一個矩陣，也可以看作是一種線性變換，所謂的奇異值分解就是將這一變換分解成兩次旋轉變換和一次沿座標軸的拉伸變換的過程，更直觀的變換過程如下圖所示：
奇異值分解的應用
奇異值分解在平面擬合，影像壓縮和降噪，主成分分析和推薦系統等演算法中都有廣泛的應用，下面以直線擬合和影像壓縮為例說明。
直線擬合問題
設直線的方程為 \(ax+by+c=0\) ，採集到一組點集 \((x_i, y_i)\) ，需要根據這組點集擬合直線方程。要求解這一問題，可以構建誤差函式
當所有點的誤差最小時，即認為擬合效果最好。這本質上是一個線性最小二乘問題，其數學形式可以表示為
構造矩陣 \(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}\) ，則可以將上述求解過程寫成矩陣形式
注意到 \(V^\text T \boldsymbol{x}\) 其實是向量 \(\boldsymbol{x}\) 在基底 \(V\) 下的一組座標，可以記作 \(\bm \alpha\) ，即
則上式可以進一步寫成
在考慮歸一化的擬合解前提下，即 \(\bm{\left| x \right|} = 1\) 時，當且僅當座標 \(\bm \alpha\) 取 \(\begin{pmatrix} 0 & 0 & \cdots & 1 \end{pmatrix}^\text T\) 時取等號，此時求解以下方程
容易得到，此時\(\boldsymbol x = \boldsymbol{v}_N\)。也就是說，要求解誤差的最小值，只需要進行奇異值分解，取右奇異矩陣的最後一維作為解 \(\bm{x}\) 即可，此時對應的擬合誤差即為 \(\sigma_N^2\) ，也就是在最小二乘意義下的最小誤差。
下面透過一個直線擬合的程式設計例項來進一步說明上述結論。
（程式設計例項）考慮直線 \(x + 2y + 5 = 0\) ，給定一系列帶有噪聲的點集，根據這組點集擬合直線，並與真值對比，得到擬合引數的均方誤差。
採用C++語言編寫程式，呼叫Eigen線性代數庫進行SVD分解，程式碼如下：
Eigen::Vector3d abc = Eigen::Vector3d(1, 2, 5).normalized(); // 直線引數真值
const double sigma = 0.6; // 噪聲方差
const size_t count = 15; // 取樣點的個數
// 利用 OpenCV 的隨機數種子生成取樣點
cv::RNG rng;
Eigen::MatrixXd A(count, 3);
for(size_t i=0; i<count; ++i) {
A.row(i) << i, -(abc[0]/abc[1])*i - (abc[2]/abc[1]) + rng.gaussian(sigma), 1;
}
// SVD求解, 並取出 V 矩陣的最後一維作為擬合解 x
Eigen::VectorXd x = Eigen::JacobiSVD<Eigen::MatrixXd>(A, Eigen::ComputeThinV).matrixV().col(2);
// 輸出求解結果
cout << "==> true value: abc = " << abc.transpose() << endl; // 列印真值
cout << "==> svd solving result: x = " << x.transpose() << endl; // 列印擬合引數
cout << " ==> error RMS = " << sqrt( (A*x).squaredNorm() / double(count) ) << endl; // 列印均方誤差
程式執行結果截圖為：
注意SVD分解時考慮的是歸一化的向量，因此得到的也是歸一化後的擬合結果，與實際直線引數之間差了一個比例係數。
可以看到在噪聲方差為0.6，取樣點數為15個點時，得到的擬合結果的均方誤差為0.2，擬合結果基本可以接受，在噪聲方差更小的情況下，可以得到更可靠的擬合結果。
影像壓縮問題
考慮一個 \(m \times n\) 的矩陣 \(A\) ，其奇異值分解式為
我們發現求得的奇異值衰減速度較快，前面幾項的奇異值較大，而越往後的奇異值則越小，也就是說我們可以透過只保留前面幾項奇異值來對 \(A\) 進行近似，這種做法可以壓縮儲存空間，在影像處理領域有著廣泛的應用。具體做法是我們可以取出前 \(k\) 個奇異值，得到近似的 \(\tilde A\) ，其中 \(\frac{k}{\min\{m,n\}}\) 稱為矩陣的壓縮比。
（程式設計例項）給定一張圖片（如下），考慮採用不同的壓縮比對其進行壓縮，並對比分析壓縮效果。
採用C++語言編寫程式，呼叫OpenCV影像處理庫來進行影像矩陣的SVD分解，採用不同的壓縮比（0.01, 0.05, 0.1, 0.2, 0.3, 0.5）重複執行程式，得到輸出結果。程式程式碼如下：
// 讀入影像
cv::Mat img_original = cv::imread("../image.png", cv::IMREAD_GRAYSCALE);
cv::Mat img = img_original.clone();
img.convertTo(img, CV_64FC1, 1/255.0);
// 對影像進行SVD分解
cv::Mat U, Vt, W;
cv::SVD::compute(img, W, U, Vt);
// 影像壓縮操作
int rank = min(img.rows, img.cols);
cv::Mat W_hat = cv::Mat(rank, rank, CV_64FC1, cv::Scalar(0)); // 取前k維的奇異值矩陣
double compression_ratio = 0.3; // 設定壓縮比
int k = rank * compression_ratio;
for(size_t i=0; i<k; ++i) {
W_hat.at<double>(i, i) = W.at<double>(i, 0);
}
cv::Mat img_compression = U * W_hat * Vt; // 計算壓縮後的影像
// 壓縮影像輸出
img_compression.convertTo(img_compression, CV_8UC1, 255.0);
cv::imwrite("./ratio=0_3.png", img_compression);
輸出壓縮結果如下所示
- 壓縮比 = 0.01
- 壓縮比 = 0.05
- 壓縮比 = 0.1
- 壓縮比 = 0.2
- 壓縮比 = 0.3
- 壓縮比 = 0.5
可以看到測試影像的基本輪廓資訊在壓縮比為0.05的時候就已經大致凸顯，而在壓縮比為0.2時就已經基本保留影像的細節資訊，壓縮效率較好。
參考資料
[1] http://www.ams.org/publicoutreach/feature-column/fcarc-svd
[2] https://blog.csdn.net/lomodays207/article/details/88687126
[3] https://www.cnblogs.com/pinard/p/6251584.html
[4] https://blog.csdn.net/u012198575/article/details/99548136