利用相关性进行SVD计算

首先在说SVD之前回忆一波特征值分解：

特征值分解

特征值分解简单来说就是把矩阵 A \mathbf{A} A分解成特征向量矩阵 × \times ×特征值形成的对角矩阵 × \times ×特征向量矩阵的逆的形式,即

A = V Λ V − 1 (1) \tag{1} \mathbf{A = V\Lambda V^{-1}} A=VΛV−1(1)

其中 V \mathbf{V} V 是 A \mathbf{A} A 的特征向量

上式还可以写成这样，即原矩阵 × \times × 特征向量矩阵 = 特征向量矩阵 × \times × 特征值矩阵，特征值分解只能应用于方阵

A V = V Λ (2) \tag{2} \mathbf{A V = V\Lambda } AV=VΛ(2)

SVD分解

关于一个向量 X \mathbf{X} X ,首先将其看成是列向量的矩阵,且假设 X ∈ R n × m ， （ n > > m ） \mathbf{X}\in\mathbb{R}^{n \times m}，（n>> m） X∈Rn×m，（n>>m）

X = [ ∣ ∣ ∣ x 1 x 2 ⋯ x m ∣ ∣ ∣ ] (3) \tag{3} \mathbf{X}=\left[\begin{array}{llll} \mid & \mid & & \mid\\ x_1 & x_{2} & \cdots & x_{m} \\ \mid & \mid & & \mid \end{array}\right] X=⎣⎡​∣x1​∣​∣x2​∣​⋯​∣xm​∣​⎦⎤​(3)

则 X X ⊤ \mathbf{XX^{\top}} XX⊤ 是一个 m × m m\times m m×m的矩阵

而这个矩阵就是矩阵中 x 1 , x 2 ⋯ x m x_1,x_2\cdots x_m x1​,x2​⋯xm​ 相互之间的相关性

X X ⊤ = [ x 1 ⊤ x 1 x 1 ⊤ x 2 ⋯ x 1 ⊤ x m x 2 ⊤ x 1 x 2 ⊤ x 2 ⋯ x 2 ⊤ x m ⋮ ⋮ ⋱ ⋮ x m ⊤ x 1 x m ⊤ x 2 − x m ⊤ x m ] (4) \tag{4} \mathbf{XX^{\top}}=\left[\begin{array}{cccc} x_{1}^{\top} x_{1} & x_{1}^{\top} x_{2} & \cdots & x_{1}^{\top} x_{m} \\ x_{2}^{\top} x_{1} & x_{2}^{\top} x_{2} & \cdots & x_{2}^{\top} x_{m} \\ \vdots & \vdots & \ddots & \vdots \\ x_{m}^{\top} x_{1} & x_{m}^{\top} x_{2} & - & x_{m}^{\top} x_{m} \end{array}\right] XX⊤=⎣⎢⎢⎢⎡​x1⊤​x1​x2⊤​x1​⋮xm⊤​x1​​x1⊤​x2​x2⊤​x2​⋮xm⊤​x2​​⋯⋯⋱−​x1⊤​xm​x2⊤​xm​⋮xm⊤​xm​​⎦⎥⎥⎥⎤​(4)

解释完 X X ⊤ \mathbf{XX^{\top}} XX⊤的含义就可以看SVD分解了

首先SVD的分解形式为：

X = U Σ V ⊤ (5) \tag{5} \mathbf{{X}={U}{\Sigma} V^{\top}} X=UΣV⊤(5)

其中 U \mathbf{U} U 称之为 X \mathbf{X} X的左奇异向量矩阵， V \mathbf{V} V称之为 X \mathbf{X} X的奇异向量矩阵， U \mathbf{U} U 和 V \mathbf{V} V 都是共轭转置矩阵，即其本身乘以其共轭转置得到单位矩阵

U U ⋆ = I V V ⋆ = I (6) \tag{6} \mathbf{U U^\star = I} \\ \mathbf{V V^\star = I} \\ UU⋆=IVV⋆=I(6)

( U ⋆ = U ⊤ 、 V ⋆ = V ⊤ ) \mathbf{(U ^\star =U^\top、V^{\star} = V^\top)} (U⋆=U⊤、V⋆=V⊤) 如果 U , V \mathbf{U,V} U,V 不是复数矩阵。

Σ \mathbf{\Sigma} Σ 称之为 X \mathbf{X} X的奇异值矩阵，其中 Σ \mathbf{\Sigma} Σ 是对角矩阵，且每个奇异值的重要性按行顺序排列。

SVD分解可以是非方阵，也就是说任何矩阵都有其奇异值矩阵、左/右奇异向量矩阵，这和特征值分解是不同的。

SVD分解和特征值分解的关系

对 X \mathbf{X} X和 X ⊤ \mathbf{X^{\top}} X⊤进行简化SVD分解可以分别得到

注：所谓简化SVD分解是

matlab

的一个函数

svd（X,'economy'）

，其得到的是近似最优，以 $ \mathbf{U}$ 的维度确定 $ \Sigma$ 的维度，下面将矩阵加帽以表示

X ˉ = U ^ Σ ^ V ⊤ X ˉ ⊤ = V Σ ^ U ^ ⊤ (7) \tag{7} \begin{array}{l} \mathbf{\bar{X}=\hat{U} \hat{\Sigma} V^{\top}}\\ \mathbf{\bar{X}^{\top}=V \hat{\Sigma} \hat{U}^{\top}} \end{array} Xˉ=U^Σ^V⊤Xˉ⊤=VΣ^U^⊤​(7)

则 X ⊤ X \mathbf{X^{\top}X} X⊤X 可以写成：

X ⊤ X = V Σ ^ U ^ ⊤ U ^ Σ ^ V ⊤ = V ^ Σ ^ 2 V ^ ⊤ ⟹ X ⊤ X V ^ = V ^ Σ ^ 2 (8) \tag{8} \mathbf{X^{\top}X}=\mathbf{V \hat{\Sigma} \hat{U}^{\top}\hat{U} \hat{\Sigma} V^{\top} =\hat{V} \hat{\Sigma}^{2}\hat{V}^{\top}}\\\\ \Longrightarrow\qquad \mathbf{X^{\top}X\hat{V}}=\mathbf{ \hat{V}\hat{\Sigma}^{2}} X⊤X=VΣ^U^⊤U^Σ^V⊤=V^Σ^2V^⊤⟹X⊤XV^=V^Σ^2(8)

则 X X ⊤ \mathbf{XX^{\top}} XX⊤ 可以写成：

X X ⊤ = U ^ Σ ^ V ⊤ V Σ ^ U ^ ⊤ = U ^ Σ ^ 2 U ^ ⊤ ⟹ X X ⊤ U ^ = U ^ Σ ^ 2 (9) \tag{9} \mathbf{XX^{\top}}=\mathbf{\hat{U} \hat{\Sigma} V^{\top}V \hat{\Sigma} \hat{U}^{\top} =\hat{U} \hat{\Sigma}^{2}\hat{U}^{\top}}\\\\ \Longrightarrow\qquad \mathbf{XX^{\top}\hat{U}}=\mathbf{\hat{U} \hat{\Sigma}^{2}} XX⊤=U^Σ^V⊤VΣ^U^⊤=U^Σ^2U^⊤⟹XX⊤U^=U^Σ^2(9)

上面两个式子，结合公式 ( 2 ) (2) (2) 可以推理得到的是 V ^ \mathbf{\hat{V}} V^ 和 U ^ \mathbf{\hat{U}} U^ 都是上述对应矩阵 X ⊤ X \mathbf{X^{\top}X} X⊤X 和 X X ⊤ \mathbf{XX^{\top}} XX⊤ 的特征向量矩阵，而 Σ 2 \mathbf{\Sigma^2} Σ2是特征值。也就是说，一个矩阵和其自身转置相乘得到其协方差阵，然后进行协方差特征值分解， X ⊤ X \mathbf{X^{\top}X} X⊤X 得到的是 X \mathbf{X} X的右特征值 V \mathbf{V} V， X X ⊤ \mathbf{XX^{\top}} XX⊤得到的是左特征值 U \mathbf{U} U，这样便实现了使用特征值分解得到奇异分解的结果了。

但是需要注意的是在公式 ( 9 ) (9) (9) 中得出来的的矩阵 X X ⊤ \mathbf{XX^{\top}} XX⊤ 一般来说都非常巨大，是 n × n n \times n n×n维矩阵，所以应该尽量避免（一般情况下数据集行大于列 n > > m n>>m n>>m，尤其是使用了Matlab命令

svd(A,'economy')

的时候）。

使用snapshot的方法进行左singular vector的计算：

当使用 ( 8 ) (8) (8) 的时候，我们可以获取 V ^ 、 Σ ^ \mathbf{ \hat{V}、\hat{\Sigma}} V^、Σ^ 的值，这个时候我们可以求出 U ^ \mathbf{\hat{U}} U^的值为：

U ^ = X ^ V Σ ^ − 1 (10) \tag{10} \mathbf{\hat{U}=\hat{X} V\hat{\Sigma}^{-1}} U^=X^VΣ^−1(10)

这个方法就是snpshot，因为矩阵 V \mathbf{V} V就相当于 X \mathbf{X} X在时间序列上的变换，每一列（ V ⊤ \mathbf{V^{\top}} V⊤的每一行）相当于一个时间点所有的数据的融合。

Underdetermined system 和overdetermined system

两种系统，用

A x = b (11) \tag{11} Ax = b Ax=b(11)

表示的话，大概是这个样子：

对于第一种不确定系统，其系统矩阵的形状如下：

这种情况下，因为 x x x的数量众多， b b b无法提供足够的数量去填充 x x x，因此系统有无数个解，通俗用一个公式解释就是 3 x 1 + 4 x 2 = 0 3x_ 1+4x_ 2=0 3x1​+4x2​=0 ,对于 x 1 、 x 2 x_ 1、x_ 2 x1​、x2​有无穷多个解满足此公式，因为 b = 0 , b b=0,b b=0,b的数量只有一个而变量有两个。

另一种情况就是 over determined system，这种情况就是变量过少，结果过多，导致没有一个准确的解。

基于SVD的PCA 分析

PCA 算法

话不多说，先上如何计算PCA。

假设 X \mathbf{X} X是一个矮胖矩阵：

1. 首先对一个矩阵 X \mathbf{X} X所有行进行的平均运算，得到一个很胖的单行向量， 1 ∗ m 1*m 1∗m 维度的矢量 x ˉ \bar{\mathbf{x}} xˉ 。

然后进行构建矩阵，用一个单位列向量构建：

X ‾ = [ 1 ⋮ 1 ] x ‾ (12) \tag{12} \overline{\mathbf{X}}=\left[\begin{array}{c} 1 \\ \vdots \\ 1 \end{array}\right] \overline{\mathbf{x}} X=⎣⎢⎡​1⋮1​⎦⎥⎤​x(12)

2. 然后用 X X X减去 X ˉ \bar{X} Xˉ得到去掉均值的矩阵 B \mathbf{B} B ：

B = X − X ‾ (13) \tag{13} \mathbf{B}=\mathbf{X}-\overline{\mathbf{X}} B=X−X(13)

1. 计算 B \mathbf{B} B的协方差矩阵 C \mathbf{C} C：

   C = 1 n − 1 B ∗ B (14) \tag{14} \mathbf{C}=\frac{1}{n-1} \mathbf{B}^{*} \mathbf{B} C=n−11​B∗B(14)

2. 计算协方差矩阵 C \mathbf{C} C的特征值和特征向量：

v 1 ⊤ B ⊤ B v 1 (15) \tag{15} \mathbf{v_{1}^{\top} B^{\top} B v_{1}} v1⊤​B⊤Bv1​(15)

其中 v 1 \mathbf{v_{1}} v1​是 C \mathbf{C} C的第一个特征向量。 V \mathbf{V} V是 C \mathbf{C} C 的特征向量矩阵

1. 可以通过特征值分解（eigen-decomposition）进行解析：

   C V = V D (16) \tag{16} \mathbf{C V} =\mathbf{V D} CV=VD(16)

   C = V D V − 1 (17) \tag{17} \mathbf{C } =\mathbf{V D V^{-1}} C=VDV−1(17)

   其中 D \mathbf{D} D是特征值， V \mathbf{V} V是特征向量

2. 最后使用 B \mathbf{B} B和 V \mathbf{V} V 的乘积求出principal components

   T = B V (18) \tag{18} \mathbf{T}=\mathbf{B V} T=BV(18)

   其中 V \mathbf{V} V 称之为loadings

   在此结果上，因为

   B = U Σ V ⊤ (19) \tag{19} \mathbf{B=U \Sigma V^{\top}} B=UΣV⊤(19)

   所以的出来principal components可以是：

   T = U Σ (20) \tag{20} \mathbf{T}=\mathbf{U \Sigma} T=UΣ(20)

   其中 Σ \mathbf{\Sigma} Σ是 B \mathbf{B} B的奇异值， U \mathbf{U} U是左奇异向量。

知识痛点补充：

1. 为什么上述的协方差计算使用 1 n − 1 \frac{1}{n-1} n−11​而不是 1 n \frac{1}{n} n1​作为系数？

   首先明确协方差和方差的区别，数据的方差反映的是一个数据集自身的变化程度，也就是这个数据集的离散度，是对一个随机变量的测量值。

   协方差则是描述对于不同的，相关的随机变量之间的变化关系，比如一个人的身高和体重两个变量的关系是正相关的，那么这两个变量的协方差就描述了身高和体重的变化之间的相互影响，就是多少身高带来多少体重、或者是多少体重的变化带来多少身高的变化。

   因为在我们使用的任何数据集中，基本上都是更大一部分数据集的子集。在这部分数据集中，我们求得的协方差并不是能反映整体数据集的自然误差状态。换句话说，如果我们用 1 n \frac{1}{n} n1​ 做系数，就相当于我们用了大数据集中的一部分数据的协方差去反映整个大数据集的测量误差，所以用 1 n − 1 \frac{1}{n-1} n−11​ 而非 1 n \frac{1}{n} n1​ 来弥补这种不合理性。比如我们在一个班里抽取学生的身高体重样本来反映整个班级的身高体重水平，班里有100个学生，我们现在取了20个学生（因为我们没有时间把100个学生的身高体重都测量一遍），使用他们的身高体重求协方差，然后用这个协方差去衡量整个班级的身高体重水平，这个时候我们的数据集只是这个大数据集的一小部分，所以我们的衡量是有偏差的，因此我们用 1 n − 1 \frac{1}{n-1} n−11​ 而不是 1 n \frac{1}{n} n1​ 来抵消这种偏差。

PCA的用处

Eigenface

PCA有很多作用，比如在人脸识别领域，拾取人脸得关键特征。最著名的一个例子是所谓的eigenface。

顾名思义，eigenface就是脸的

eigenvalue

，即脸的特征（其实任意图片都有特征，但是使用一类图片比如人脸图片较为好比较）。

使用PCA提取人脸图片的特征，并且将其重新组成图片，过程如下。

1. 获取 n n n张同一个人不同的人脸照片，比如192x168像素的照片，用灰度表示（比较简单，这样每一个pixel就只有一个值而不是3个值了）。
2. 将每一张照片变成列向量，每一列代表一个照片，构成一个矩阵 X X X，即矩阵 X X X的维度是 32256 × n 32256 \times n 32256×n。
3. 计算矩阵 X X X的每一行的平均值，得到一个 32256 × 1 32256\times 1 32256×1 的列向量，然后将这个列向量扩充成矩阵，变成 32256 × n 32256\times n 32256×n的矩阵，每一列的值都相同，这个矩阵称之为均值矩阵 X ˉ \bar{X} Xˉ 。
4. 通过公式 ( 13 ) (13) (13)算出去平均值后的矩阵 B B B
5. 将 B B B进行SVD分解得到 U , S , V U,S,V U,S,V 其中，矩阵 U U U就是eigenface矩阵，每一列就是每张照片的eigenface,每一张都记录了不同方面的信息

以上五步就可以的出来eigenface，下面发一个代码，以及本人靓照：

使用eigenfaces可以拟合出来原图，选取秩数的不同会得到不同的还原度，而且即便用于测试的图片非常不一样，下图我用了一张类似的照片做测试。

注意上图最后一张，也许是使用了

matlab

的

svd(,'econ')

命令导致的，最后一张的图像总会得到类似于马赛克的效果。

图片原图：

我用于测试的图片：

代码：

clear all, close all, clc

file_path = 'D:\OneDrive\西安交大\实验\MATLAB 练习\Data-driven Science\Practice\svd\myface\eyes\';
h = 192; % 高
w = 144; % 宽
n = length(dir(strcat(file_path,'*.jpg'))); %图片数量
m = h*w; % 单列维度
I = zeros(m,n);
for i=1:n-1
    img = strcat(num2str(i),'.jpg');
    gray_img = rgb2gray(imread(strcat(file_path,img)));
    I(:,i) = reshape(gray_img,m,1);    %将所有的像素变成列向量
%     imagesc(reshape(I(:,i),h,w)); %展示图片
%     pause(0.1)
end
%% 搞出来平均脸然后去平均后进行SVD
avg_face = mean(I,2);
% imagesc(reshape(avg_face,h,w)), colormap gray; %平均脸
%画出所有的eigenface，越往后占的重要程度越低，到最后一张变成一坨糊糊
B = I-avg_face*ones(1,size(I,2));
[U,S,V] = svd(B,'econ');
for i = 1:n
    imagesc(reshape(U(:,i),h,w));
    pause(0.1);
end

%% 画出不同r的eigenface，用最后一张照片做测试
test_img = reshape(rgb2gray(im2double(imread(strcat(file_path,strcat(num2str(n),'.jpg'))))),m,1); 
truncated = test_img - avg_face;

% for r = 1: 1 :n
%     new_img = avg_face + (U(:,1:r)*(U(:,1:r)'*truncated));
%     imagesc(reshape(new_img,h,w)),colormap gray;
%     pause(0.1)
% end

%% 将所有图画到一张上去，一定注意不要太多张图
f=figure;
row = 3; %3 行排列
col = floor(n/row); 
zoom = 2;% 缩放最后一张图的大小
set(gcf,'Position',[500 500 col*w*zoom row*h*zoom])
movegui(f,'center');

for r = 1:n
    new_img = avg_face + (U(:,1:r)*(U(:,1:r)'*truncated));
    subplot(row,col,r),imagesc(reshape(new_img,h,w));colormap gray, axis off;
    s = sprintf('%s%d','秩数为r =',r);
    title(s);
end
sgtitle('累计r的效果，累计度越高越清晰，因为叠加的特征多了');

eigenfaces = figure();
set(gcf,'Position',[200 500 col*w*zoom row*h*zoom])
for i=1:n
    subplot(row,col,i),imagesc(reshape(U(:,i),h,w));colormap gray, axis off;
    s = sprintf('%s%d%s','第',i,'张eigenface');
    title(s);
end
sgtitle('eigenfaces,每一个都记录了不同方面的信息');
           

帅的人都很细心地发现了上面的图片都只动了我的慧眼，其他部分保持了不变，这是因为SVD对于图片处理有天然的缺陷。即SVD 对图像进行分析的时候受到图像的排放方式影响很大。

图片来源于: <Data-Driven Science and Engineering Machine Learning Dynamics Systems and Control> \text{<Data-Driven Science and Engineering Machine Learning Dynamics Systems and Control>} <Data-Driven Science and Engineering Machine Learning Dynamics Systems and Control>

*SVD受到图像内主体的排列方位影响很大，因为这直接导致了整个图像矩阵奇异值改变。*而奇异值分解是对图像进行操作的主要办法，通过获取不同程度的秩的图像，可以得到不同还原度的原图像。

如何选取最优秩r作为分界线实现图片压缩

对于使用SVD方法对图像进行压缩，其实实质上是我们不需要将所有的奇异值都给弄出来用作压缩图片使用，因为奇异值和奇异向量是按照其重要程度排序的，而且其奇异值矩阵中，前面一小部分奇异值的重要性是是超过后面大部分奇异值重要性之和的。因此，在压缩图像的时候，我们只需要前面 r r r 个奇异值和奇异向量就可以实现图像的压缩。

硬分界线r

我们设法搞一个分界线，这个分界线以使用奇异值矩阵秩的数量 r r r为边界，**使其能够保留90%以上信息的前提下实现使用最小秩的数量，**也就是使 r r r最小。这个分界线 r r r就是下图中如何划分 U ^ r e m \hat{\mathbf{U}}_ {rem} U^rem​和 U ~ \tilde{\mathbf{U}} U~的关键。

有一种方法是最优硬分界线，选取秩r大于该分界线的部分最为Truncated SVD，往往最有效率。

那么最优硬分界线划分假设被分解的矩阵 X \mathbf{X} X可以被分成两个部分：

X = X true + γ X noise  (21) \tag{21} \mathbf{X}=\mathbf{X}_ {\text{true}} + \gamma \mathbf{X}_ {\text {noise }} X=Xtrue​+γXnoise ​(21)

其中 X noise \mathbf{X}_ \text{noise} Xnoise​ 是高斯噪音， γ \gamma γ是噪音的magnitude。

1. 如果 X ∈ R n × n \mathbf{X} \in \mathbb{R}^{n\times n} X∈Rn×n,即 X \mathbf{X} X是个方阵，而且$\gamma $已知则硬分界线为：

   τ = ( 4 / 3 ) n γ (22) \tag{22} \tau=(4 / \sqrt{3}) \sqrt{n} \gamma τ=(4/3

   ​)n

   ​γ(22)

2. 如果 X ∈ R n × m \mathbf{X} \in \mathbb{R}^{n\times m} X∈Rn×m,且 m < < n m<<n m<<n 即 X \mathbf{X} X是个非常细长的矩阵，则 4 / 3 4/\sqrt{3} 4/3

   ​替换成 β = m / n \beta=m / n β=m/n ：

   τ = λ ( β ) n γ (23) \tag{23} \tau=\lambda(\beta) \sqrt{n} \gamma τ=λ(β)n

   ​γ(23)

   λ ( β ) = ( 2 ( β + 1 ) + 8 β ( β + 1 ) + ( β 2 + 14 β + 1 ) 1 / 2 ) 1 / 2 (24) \tag{24} \lambda(\beta)=\left(2(\beta+1)+\frac{8 \beta}{(\beta+1)+\left(\beta^{2}+14 \beta+1\right)^{1 / 2}}\right)^{1 / 2} λ(β)=(2(β+1)+(β+1)+(β2+14β+1)1/28β​)1/2(24)

3. 如果 X ∈ R n × m \mathbf{X} \in \mathbb{R}^{n\times m} X∈Rn×m,而且$\gamma $未知，则最优硬分界线为：

   τ = ω ( β ) σ m e d (25) \tag{25} \tau=\omega(\beta) \sigma_{\mathrm{med}} τ=ω(β)σmed​(25)

   ω ( β ) = λ ( β ) / μ β (26) \tag{26} \omega(\beta)=\lambda(\beta) / \mu_{\beta} ω(β)=λ(β)/μβ​(26)

   其中 σ m e d \sigma_{\mathrm{med}} σmed​是奇异值的中位数, μ β \mu_{\beta} μβ​ 由下式解出：

∫ ( 1 − β ) 2 μ β [ ( ( 1 + β ) 2 − t ) ( t − ( 1 − β ) 2 ) ] 1 / 2 2 π t d t = 1 2 (27) \tag{27} \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} ∫(1−β)2μβ​​2πt[((1+β

​)2−t)(t−(1−β

​)2)]1/2​dt=21​(27)

例子：

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)