-
Notifications
You must be signed in to change notification settings - Fork 308
/
subspaceSVD.m
33 lines (31 loc) · 946 Bytes
/
subspaceSVD.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
function [U1, D1, V1, r] = subspaceSVD(S)
% Returns quantities satisfying :math:`S = U1 D1 V1^T`,
% where `U1` and `V1` have `r` orthonormal columns,
% and `D1` is `r x r` diagonal and has numerical rank `r`.
%
% The matrix `S` may be diagonal.
% Matrices `U1` and `V1` may be used by `projectSVD.m`
% to project given vectors onto certain subspaces.
%
% USAGE:
%
% [U1, D1, V1, r] = subspaceSVD(S)
%
% INPUT:
% S: matrix
%
% OUTPUTS:
% U1, D1, V1: matrices
% r: numerical rank of `D1`
%
% .. Author: - Michael Saunders, 29 Jul 2009 First version of subspaceSVD.m
% written as alternative to Ronan's subspaceProjector.m.
[U1,D1,V1] = svd(full(S),'econ');
d = diag(D1);
dmax = max(d);
tol = 1e-12;
pos = find(d > dmax*tol); % pos = 1:r = indices of positive d
r = length(pos);
U1 = U1( : ,1:r);
D1 = D1(1:r,1:r);
V1 = V1( : ,1:r);