截取关键成分

最优硬分界线

在SVD分解中，选取合适的秩r作为分界线至关重要。这个分界线就是下图中如何划分$\hat{\mathbf{U}}_ {rem}$和$\tilde{\mathbf{U}}$的关键。有一种方法是最优硬分界线，选取秩r大于该分界线的部分最为Truncated SVD，往往最有效率。

那么最优硬分界线划分假设被分解的矩阵$\mathbf{X}$可以被分成两个部分：
$$
\begin{equation}
\mathbf{X}=\mathbf{X}{\text{true}} + \gamma \mathbf{X} {\text {noise }}
\end{equation}
$$
其中$\mathbf{X}_ \text{noise}$ 是高斯噪音，$\gamma$是噪音的magnitude。

如果$\mathbf{X} \in \mathbb{R}^{n\times n}$,即$\mathbf{X}$是个方阵，而且$\gamma $已知则硬分界线为：
$$
\begin{equation}
\tau=(4 / \sqrt{3}) \sqrt{n} \gamma
\end{equation}
$$
如果$\mathbf{X} \in \mathbb{R}^{n\times m}$,且$m<<n$ 即$\mathbf{X}$是个非常细长的矩阵，则$\frac{4}{\sqrt{3}}$替换成$\beta=m / n$ ：
$$
\begin{equation}
\tau=\lambda(\beta) \sqrt{n} \gamma
\end{equation}
$$

$$
\begin{equation}
\lambda(\beta)=\left(2(\beta+1)+\frac{8 \beta}{(\beta+1)+\left(\beta^{2}+14 \beta+1\right)^{1 / 2}}\right)^{1 / 2}
\end{equation}
$$

如果$\mathbf{X} \in \mathbb{R}^{n\times m}$,而且$\gamma $未知，则最优硬分界线为：
$$
\begin{equation}
\tau=\omega(\beta) \sigma_{\mathrm{med}}
\end{equation}
$$

$$
\begin{equation}
\omega(\beta)=\lambda(\beta) / \mu_{\beta}
\end{equation}
$$

其中$\sigma_{\mathrm{med}}$是奇异值的中位数,$\mu_{\beta}$ 由下式解出：

$$
\begin{equation}
\int_{(1-\beta)^{2}}^{\mu_{\beta}} \frac{\left[ \left((1+\sqrt{\beta})^{2}-t\right)\left(t-(1-\sqrt{\beta})^{2}\right)\right]^{1 / 2}}{2 \pi t} d t=\frac{1}{2}
\end{equation}
$$
书中例子：

clear all, close all, clc
%% 造出来一个真实X
t = (-3:.01:3)';
Utrue = [cos(17*t).*exp(-t.^2) sin(11*t)];
Strue = [2 0; 0 .5];
Vtrue = [sin(5*t).*exp(-t.^2) cos(13*t)];
X = Utrue*Strue*Vtrue';
figure, imshow(X);
%% 加噪音
sigma = 1;
Xnoisy = X+sigma*randn(size(X));
figure, imshow(Xnoisy);
%% 使用硬分界线产生一个clean matrix
[U,S,V] = svd(Xnoisy);
N = size(Xnoisy,1);
cutoff = (4/sqrt(3))*sqrt(N)*sigma; % Hard threshold
r = max(find(diag(S)>cutoff)); % Keep modes w/ sig > cutoff
Xclean = U(:,1:r)*S(1:r,1:r)*V(:,1:r)';
figure, imshow(Xclean)
%% 找到90%累计能量
cdS = cumsum(diag(S))./sum(diag(S));  % Cumulative energy
r90 = min(find(cdS>0.90));  % Find r to capture 90% energy
X90 = U(:,1:r90)*S(1:r90,1:r90)*V(:,1:r90)';
figure, imshow(X90)

%% plot singular values 画图
semilogy(diag(S),'-ok','LineWidth',1.5), hold on, grid on
semilogy(diag(S(1:r,1:r)),'or','LineWidth',1.5)
plot([-20 N+20],[cutoff cutoff],'r--','LineWidth',2)
axis([-10 610 .003 300])
rectangle('Position',[-5,20,100,200],'LineWidth',2,'LineStyle','--')
     
figure
semilogy(diag(S),'-ok','LineWidth',1.5)
hold on, grid on
semilogy(diag(S(1:r,1:r)),'or','LineWidth',1.5)
plot([-20 N+20],[cutoff cutoff],'r--','LineWidth',2)
axis([-5 100 20 200])
     
figure
plot(cdS,'-ok','LineWidth',1.5)
hold on, grid on
plot(cdS(1:r90),'ob','LineWidth',1.5)
plot(cdS(1:r),'or','LineWidth',1.5)
set(gca,'XTick',[0 300 r90 600],'YTick',[0 .5 0.9 1.0])
xlim([-10 610])
plot([r90 r90 -10],[0 0.9 0.9],'b--','LineWidth',1.5)