svds
奇异值和向量的子集
语法
说明
使用一个或多个名称-值对组参数指定其他选项。例如,s
= svds(A
,k
,sigma
,Name,Value
)svds(A,k,sigma,'Tolerance',1e-3)
将调整算法的收敛容差。
示例
最大奇异值
矩阵 A = delsq(numgrid('C',15))
是一个对称正定矩阵,奇异值合理分布在区间 (0 8) 中。计算六个最大的奇异值。
A = delsq(numgrid('C',15));
s = svds(A)
s = 6×1
7.8666
7.7324
7.6531
7.5213
7.4480
7.3517
指定第二个输入,以计算特定数量的最大奇异值。
s = svds(A,3)
s = 3×1
7.8666
7.7324
7.6531
最小奇异值
矩阵 A = delsq(numgrid('C',15))
是一个对称正定矩阵,奇异值合理分布在区间 (0 8) 中。计算五个最小的奇异值。
A = delsq(numgrid('C',15)); s = svds(A,5,'smallest')
s = 5×1
0.5520
0.4787
0.3469
0.2676
0.1334
最小非零奇异值
创建一个 100×100 的稀疏 Neumann 矩阵。
C = gallery('neumann',100);
计算十个最小的奇异值。
ss = svds(C,10,'smallest')
ss = 10×1
0.9828
0.9049
0.5625
0.5625
0.4541
0.4506
0.2256
0.1139
0.1139
0
计算十个最小的非零奇异值。由于该矩阵有一个奇异值等于零,因此 'smallestnz'
选项将忽略它。
snz = svds(C,10,'smallestnz')
snz = 10×1
0.9828
0.9828
0.9049
0.5625
0.5625
0.4541
0.4506
0.2256
0.1139
0.1139
使用函数句柄的最大奇异值
创建两个矩阵,表示稀疏矩阵右上角和左下角的非零块。
n = 500; B = rand(500); C = rand(500);
将 Afun
保存到您的当前目录中,以便其可用于 svds
。
function y = Afun(x,tflag,B,C,n) if strcmp(tflag,'notransp') y = [B*x(n+1:end); C*x(1:n)]; else y = [C'*x(n+1:end); B'*x(1:n)]; end
函数 Afun
使用 B
和 C
来计算 A*x
或 A'*x
(具体取决于指定的标志),而不用真正建立整个稀疏矩阵 A = [zeros(n) B; C zeros(n)]
。这种方法可以在计算 A*x
和 A'*x
时充分利用矩阵的稀疏模式来节省内存。
使用 Afun
计算 A
的 10 个最大的奇异值。将 B
、C
和 n
作为额外输入传递给 Afun
。
s = svds(@(x,tflag) Afun(x,tflag,B,C,n),[1000 1000],10)
s = 250.3248 249.9914 12.7627 12.7232 12.6988 12.6608 12.6166 12.5643 12.5419 12.4512
直接计算 A
的 10 个最大的奇异值以比较结果。
A = [zeros(n) B; C zeros(n)]; s = svds(A,10)
s = 250.3248 249.9914 12.7627 12.7232 12.6988 12.6608 12.6166 12.5643 12.5419 12.4512
稀疏矩阵的奇异值分解
west0479
是一个 479×479 的实数值稀疏矩阵。该矩阵有几个较大的奇异值和许多较小的奇异值。
加载 west0479
并将其存储为 A
。
load west0479
A = west0479;
计算 A
的奇异值分解,返回六个最大的奇异值和对应的奇异向量。指定第四个输出参数,以检查奇异值的收敛。
[U,S,V,cflag] = svds(A); cflag
cflag = 0
cflag
表示所有奇异值已收敛。奇异值位于输出矩阵 S
的对角线上。
s = diag(S)
s = 6×1
105 ×
3.1895
3.1725
3.1695
3.1685
3.1669
0.3038
通过计算 A
的完整奇异值分解来检查结果。将 A
转换成满矩阵并使用 svd
。
[U1,S1,V1] = svd(full(A));
使用对数刻度绘制由 svd
和 svds
计算的 A
的六个最大奇异值。
s2 = diag(S1); semilogy(s2(1:6),'r.') hold on semilogy(s,'ro','MarkerSize',10) title('Singular Values of west0479') legend('svd','svds')
解决收敛问题
创建一个稀疏对角矩阵并计算六个最大的奇异值。
A = diag(sparse([1e4*ones(1, 8) 1e4:-1:1])); s = svds(A)
Warning: Only 2 of the 6 requested singular values converged. Singular values that did not converge are NaN.
s = 6×1
104 ×
1.0000
0.9999
NaN
NaN
NaN
NaN
由于已执行完最大迭代次数,但仍然无法满足容差要求,所以 svds
算法生成一条警告消息。
解决收敛问题的最有效方法是使用更大的 'SubspaceDimension'
值,以增加计算中使用的 Krylov 子空间的最大大小。此操作可通过传入值为 60
的名称-值对组 'SubspaceDimension'
来完成。
s = svds(A,6,'largest','SubspaceDimension',60)
s = 6×1
104 ×
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
使用 QR 分解计算接近奇异矩阵的 SVD
计算接近奇异矩阵的 10 个最小的奇异值。
rng default format shortg B = spdiags([repelem([1; 1e-7], [198, 2]) ones(200, 1)], [0 1], 200, 200); s1 = svds(B,10,'smallest')
Warning: Large residual norm detected. This is likely due to bad condition of the input matrix (condition number 1.0008e+16).
s1 = 10×1
7.0945
7.0945
7.0945
7.0945
7.0945
7.0945
7.0945
7.0945
0.25927
7.0888e-16
警告表明 svds
无法计算正确的奇异值。svds
失败的原因在于最小奇异值和次小奇异值之间存在间隔。svds(...,'smallest')
需要求 B
的逆矩阵,这会产生较大的数值误差。
要进行比较,请使用 svd
计算精确的奇异值。
s = svd(full(B)); s = s(end-9:end)
s = 10×1
0.14196
0.12621
0.11045
0.094686
0.078914
0.063137
0.047356
0.031572
0.015787
7.0888e-16
要使用 svds
再现此计算,请执行 B
的 QR 分解。三角矩阵 R
的奇异值与 B
的相同。
[Q,R,p] = qr(B,0);
绘制 R
的每行的范数。
rownormR = sqrt(diag(R*R')); semilogy(rownormR) hold on; semilogy(size(R, 1), rownormR(end), 'ro')
R
中的最后一项接近于零,这导致解不稳定。
通过将 R
的最后一行精确设置为零,可以防止此项破坏解的正确部分。
R(end,:) = 0;
使用 svds
求出 R
的 10 个最小的奇异值。结果与通过 svd
获得的值类似。
sr = svds(R,10,'smallest')
sr = 10×1
0.14196
0.12621
0.11045
0.094686
0.078914
0.063137
0.047356
0.031572
0.015787
0
要使用此方法计算 B
的奇异向量,请使用 Q
和置换向量 p
转换左奇异向量和右奇异向量。
[U,S,V] = svds(R,20,'s');
U = Q*U;
V(p,:) = V;
输入参数
A
— 输入矩阵
矩阵
输入矩阵。A
通常是(但不总是)一个大型稀疏矩阵。
数据类型: double
复数支持: 是
k
— 要计算的奇异值个数
标量
要计算的奇异值个数,指定为正整数标量。如果满足下列任一条件,svds
返回的奇异值个数将少于请求的数量:
k
大于min(size(A))
sigma = 'smallestnz'
且k
大于A
的非零奇异值个数
如果 k
太大,svds
会将其替换为 k
的最大有效值。
示例: svds(A,2)
返回 A
的两个最大的奇异值。
sigma
— 奇异值的类型
'largest'
(默认) | 'smallest'
| 'smallestnz'
| 标量
奇异值的类型,指定为下列值之一。
选项 | 描述 |
---|---|
| 最大奇异值 |
| 最小奇异值 |
| 最小非零奇异值 |
标量 | 最接近标量的奇异值 |
示例: svds(A,k,'smallest')
计算 k
个最小奇异值。
示例: svds(A,k,100)
计算与 100
最接近的 k
个奇异值。
数据类型: double
| char
| string
opts
— Options 结构体
结构体
Options 结构体,指定为包含下表中的一个或多个字段的结构体。
注意
不建议使用 options 结构体指定选项,而应使用名称-值对组。
选项字段 | 描述 | 名称-值对组 |
---|---|---|
tol | 收敛容差 | 'Tolerance' |
maxit | 最大迭代次数 | 'MaxIterations' |
p | Krylov 子空间的最大大小 | 'SubspaceDimension' |
u0 | 左初始起始向量 | 'LeftStartVector' |
v0 | 右初始起始向量 | 'RightStartVector' |
disp | 诊断信息的显示级别 | 'Display' |
fail | 输出中未收敛的奇异值的处理方式 | 'FailureTreatment' |
注意
如果使用数值标量偏移值 sigma
,svds
将忽略选项 p
。
示例: opts.tol = 1e-6, opts.maxit = 500
将创建一个结构体,其中包含字段 tol
和 maxit
的设置值。
数据类型: struct
Afun
— 矩阵函数
函数句柄
矩阵函数,指定为函数句柄。函数 Afun
必须满足下列条件:
Afun(x,'notransp')
接受向量x
并返回乘积A*x
。Afun(x,'transp')
接受向量x
并返回乘积A'*x
。
注意
只有在 sigma = 'largest'
(默认值)这种情况下,才使用函数句柄。
示例: svds(Afun,[1000 1200])
n
— Afun
使用的矩阵的大小
二元素向量
Afun
使用的 A
矩阵的大小,指定为二元素向量 [m n]
。
名称-值参数
将可选的参数对组指定为 Name1=Value1,...,NameN=ValueN
,其中 Name
是参数名称,Value
是对应的值。名称-值参数必须出现在其他参数之后,但参数对组的顺序无关紧要。
在 R2021a 之前,使用逗号分隔每个名称和值,并用引号将 Name
引起来。
示例: s = svds(A,k,sigma,'Tolerance',1e-10,'MaxIterations',100)
放宽收敛容差并使用较少的迭代。
Tolerance
— 收敛容差
1e-14
(默认) | 非负实数标量
收敛容差,以逗号分隔的对组形式指定,该对组由 'Tolerance'
和一个非负实数数值标量组成。
示例: s = svds(A,k,sigma,'Tolerance',1e-3)
MaxIterations
— 算法的最大迭代次数
300
(默认) | 正整数
算法的最大迭代次数,以逗号分隔的对组形式指定,该对组由 'MaxIterations'
和一个正整数组成。
示例: s = svds(A,k,sigma,'MaxIterations',350)
SubspaceDimension
— Krylov 子空间的最大大小
max(3*k,15)
(默认) | 非负整数
Krylov 子空间的最大大小,以逗号分隔的对组形式指定,该对组由 'SubspaceDimension'
和一个非负整数组成。'SubspaceDimension'
值必须大于或等于 k + 2
,其中 k
是奇异值的数量。
对于 svds
不能收敛的问题,增大 'SubspaceDimension'
的值可以改善收敛行为。
对于数值类型的 sigma
,此选项将被忽略。
示例: s = svds(A,k,sigma,'SubspaceDimension',25)
LeftStartVector
— 左初始起始向量
向量
左初始起始向量,以逗号分隔的对组形式指定,该对组由 'LeftStartVector'
和一个数值向量组成。
您可以指定 'LeftStartVector'
或 'RightStartVector'
,但不能同时指定两者。如果未指定任何选项,则对于 m
×n
矩阵 A
来说,默认值为:
m < n
- 左初始起始向量设置为randn(m,1)
m >= n
- 右初始起始向量设置为randn(n,1)
指定不同随机起始向量的主要原因是为了控制用于生成向量的随机数流。
注意
svds
使用专用的随机数流以可再现的方式选择起始向量。更改随机数种子不会影响 randn
的这种用法。
示例: s = svds(A,k,sigma,'LeftStartVector',randn(m,1))
使用从全局随机数流中获取值的随机起始向量。
数据类型: double
RightStartVector
— 右初始起始向量
向量
右初始起始向量,以逗号分隔的对组形式指定,该对组由 'RightStartVector'
和一个数值向量组成。
您可以指定 'LeftStartVector'
或 'RightStartVector'
,但不能同时指定两者。如果未指定任何选项,则对于 m
×n
矩阵 A
来说,默认值为:
m < n
- 左初始起始向量设置为randn(m,1)
m >= n
- 右初始起始向量设置为randn(n,1)
指定不同随机起始向量的主要原因是为了控制用于生成向量的随机数流。
注意
svds
使用专用的随机数流以可再现的方式选择起始向量。更改随机数种子不会影响 randn
的这种用法。
示例: s = svds(A,k,sigma,'RightStartVector',randn(n,1))
使用从全局随机数流中获取值的随机起始向量。
数据类型: double
FailureTreatment
— 未收敛的奇异值的处理方式
'replacenan'
| 'keep'
| 'drop'
未收敛的奇异值的处理方式,以逗号分隔的对组形式指定,该对组由 'FailureTreatment'
和下列选项之一组成:'replacenan'
、'keep'
或 'drop'
。
'FailureTreatment'
的值决定未收敛的奇异值在输出中如何显示。
选项 | 对输出的影响 |
---|---|
| 从输出中删除未收敛的奇异值,这会导致 |
| 将未收敛的奇异值替换为 |
| 在输出中包括未收敛的奇异值。 |
示例: s = svds(A,k,sigma,'FailureTreatment','drop')
将从输出中删除未收敛的奇异值。
数据类型: char
| string
Display
— 开启或关闭诊断信息的显示
false
(默认) | true
| 0
| 1
开启或关闭诊断信息的显示,指定为 false
、true
、0
或 1
。值 false
或 0
将关闭显示,值 true
或 1
将开启显示。
输出参数
s
— 奇异值
列向量
奇异值,以列向量形式返回。奇异值是以降序顺序列出的非负实数。
U
— 左奇异向量
矩阵
左奇异向量,以矩阵的列形式返回。如果 A
是一个 m
×n
矩阵,而您请求 k
个奇异值,则 U
将是一个包含正交列的 m
×k
矩阵。
不同的计算机、MATLAB® 版本或参数(例如起始向量和子空间维度)可能生成不同的奇异向量,它们在数值上依然精确。U
和 V
中的相应列可以翻转其符号,因为这不会影响表达式 A = U*S*V'
的值。
S
— 奇异值
对角矩阵
奇异值,以对角矩阵形式返回。S
的对角元素是非负奇异值。如果 A
是一个 m
×n
矩阵,而您请求 k
个奇异值,则 S
将是一个 k
×k
矩阵。
V
— 右奇异向量
矩阵
右奇异向量,以矩阵的列形式返回。如果 A
是一个 m
×n
矩阵,而您请求 k
个奇异值,则 V
将是一个包含正交列的 n
×k
矩阵。
不同的计算机、MATLAB 版本或参数(例如起始向量和子空间维度)可能生成不同的奇异向量,它们在数值上依然精确。U
和 V
中的相应列可以翻转其符号,因为这不会影响表达式 A = U*S*V'
的值。
flag
— 收敛标志
标量
收敛标志,以标量形式返回。值为 0
表示所有奇异值已收敛。否则,表示并非所有奇异值均已收敛。
使用此收敛标志输出可抑制收敛失败时的警告。
提示
如果您事先不知道用
svds
指定什么秩但知道 SVD 的逼近应满足什么容差,则svdsketch
很有用。svds
使用专用的随机数流生成默认起始向量,以确保在不同运行之间的可再现性。调用svds
之前使用rng
设置随机数生成器状态不会影响输出。要求出小型稠密矩阵的几个奇异值,使用
svds
并不是最有效的方式。对于这些问题,使用svd(full(A))
可能会更快。例如,求 500×500 矩阵中的三个奇异值相对容易,使用svd
即可轻松完成。对于给定矩阵,如果
svds
无法收敛,可以通过增大'SubspaceDimension'
值来增加 Krylov 子空间的大小。作为备用方案,调整最大迭代次数 ('MaxIterations'
) 和收敛容差 ('Tolerance'
) 也有助于改善收敛行为。增大
k
有时可以提高性能,特别是当矩阵包含重复的奇异值时。
参考
[1] Baglama, J. and L. Reichel, “Augmented Implicitly Restarted Lanczos Bidiagonalization Methods.” SIAM Journal on Scientific Computing. Vol. 27, 2005, pp. 19–42.
[2] Larsen, R. M. “Lanczos Bidiagonalization with partial reorthogonalization.” Dept. of Computer Science, Aarhus University. DAIMI PB-357, 1998.
扩展功能
基于线程的环境
使用 MATLAB® backgroundPool
在后台运行代码或使用 Parallel Computing Toolbox™ ThreadPool
加快代码运行速度。
此函数完全支持基于线程的环境。有关详细信息,请参阅Run MATLAB Functions in Thread-Based Environment。
GPU 数组
通过使用 Parallel Computing Toolbox™ 在图形处理单元 (GPU) 上运行来加快代码执行。
用法说明和限制:
如果提供
sigma
参数,则值必须为'largest'
。
有关详细信息,请参阅Run MATLAB Functions on a GPU (Parallel Computing Toolbox)。
分布式数组
使用 Parallel Computing Toolbox™ 在集群的组合内存中对大型数组进行分区。
用法说明和限制:
如果提供
sigma
参数,则值必须为'largest'
。
有关详细信息,请参阅Run MATLAB Functions with Distributed Arrays (Parallel Computing Toolbox)。
版本历史记录
在 R2006a 之前推出R2016a: 可再现性变化
可再现性
现在,连续多次调用
svds
会产生相同的结果。要更改此行为,请执行以下操作:在 R2017a 或更早版本中,将 options 结构体的
u0
或v0
字段设置为随机向量。在 R2017b 或更高版本中,优先选择将
'LeftStartVector'
或'RightStartVector'
设置为随机向量。
MATLAB 命令
您点击的链接对应于以下 MATLAB 命令:
请在 MATLAB 命令行窗口中直接输入以执行命令。Web 浏览器不支持 MATLAB 命令。
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list:
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)