3

我在 Math Stackexchange 中问过这个问题，但似乎没有得到足够的关注，所以我在这里问。https://math.stackexchange.com/questions/1729946/why-do-we-say-svd-can-handle-singular-matrx-when-doing-least-square-comparison?noredirect=1#comment3530971_1729946

我从一些教程中了解到，在解决最小二乘问题时，SVD 应该比 QR 分解更稳定，并且能够处理奇异矩阵。但是下面我用matlab写的例子似乎支持了相反的结论。我对 SVD 没有深入的了解，所以如果你能在 Math StackExchange 的旧帖子中查看我的问题并向我解释，我将不胜感激。

我使用具有大条件数（e+13）的矩阵。结果显示 SVD 得到的误差（0.8）比 QR（e-27）大得多

% we do a linear regression between Y and X
data= [
47.667483331 -122.1070832;
47.667483331001 -122.1070832
];
X = data(:,1);
Y = data(:,2);

X_1 =  [ones(length(X),1),X];

%%
%SVD method
[U,D,V] = svd(X_1,'econ');
beta_svd = V*diag(1./diag(D))*U'*Y;


%% QR method(here one can also use "\" operator, which will get the same result as I tested. I just wrote down backward substitution to educate myself)
[Q,R] = qr(X_1)
%now do backward substitution
[nr nc] = size(R)
beta_qr=[]
Y_1 = Q'*Y
for i = nc:-1:1
    s = Y_1(i)
    for j = m:-1:i+1
        s = s - R(i,j)*beta_qr(j)
    end
    beta_qr(i) = s/R(i,i)
end

svd_error = 0;
qr_error = 0;
for i=1:length(X)
   svd_error = svd_error + (Y(i) - beta_svd(1) - beta_svd(2) * X(i))^2;
   qr_error = qr_error + (Y(i) - beta_qr(1) - beta_qr(2) * X(i))^2;
end
4

3 回答 3

6

您基于 SVD 的方法与MATLAB中的pinv函数基本相同（请参阅Pseudo-inverse 和 SVD）。但是（出于数字原因）您缺少的是使用容差值，使得任何小于该容差的奇异值都被视为零。

如果您参考edit pinv.m，您会看到类似以下内容（我不会在此处发布确切的代码，因为该文件的版权归 MathWorks 所有）：

[U,S,V] = svd(A,'econ');
s = diag(S);
tol = max(size(A)) * eps(norm(s,inf));
% .. use above tolerance to truncate singular values
invS = diag(1./s);
out = V*invS*U';

实际上pinv有第二种语法，pinv(A,tol)如果默认值不合适，您可以在其中显式指定容差值...

因此，在求解 形式的最小二乘问题时minimize norm(A*x-b)，您应该了解pinv和mldivide解具有不同的性质：

  • x = pinv(A)*b的特点norm(x)是小于任何其他解决方案的范数。
  • x = A\b具有尽可能少的非零分量（即稀疏）。

使用您的示例（请注意，rcond(A)在机器 epsilon 附近非常小）：

data = [
    47.667483331    -122.1070832;
    47.667483331001 -122.1070832
];
A = [ones(size(data,1),1), data(:,1)];
b = data(:,2);

让我们比较一下两种解决方案：

x1 = A\b;
x2 = pinv(A)*b;

首先，您可以看到如何mldivide返回x1具有一个零分量的解决方案（这显然是一个有效的解决方案，因为您可以通过乘以零来解决这两个方程，如 中所示b + a*0 = b）：

>> sol = [x1 x2]
sol =
 -122.1071   -0.0537
         0   -2.5605

接下来，您将看到如何pinv返回x2具有较小范数的解决方案：

>> nrm = [norm(x1) norm(x2)]
nrm =
  122.1071    2.5611

这是两种解决方案的错误，可以接受的非常小：

>> err = [norm(A*x1-b) norm(A*x2-b)]
err =
   1.0e-11 *
         0    0.1819

请注意，使用mldivide、、linsolve或qr将给出几乎相同的结果：

>> x3 = linsolve(A,b)
Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND =  2.159326e-16. 
x3 =
 -122.1071
         0

>> [Q,R] = qr(A); x4 = R\(Q'*b)
x4 =
 -122.1071
         0
于 2016-04-11T18:38:32.200 回答
2

SVD 可以处理秩不足。对角矩阵 D 在您的代码中有一个接近零的元素，您需要对 SVD 使用伪逆，即将 1./diag(D) 的第二个元素设置为 0，而不是巨大的值 (10^14)。您应该会发现 SVD 和 QR 在您的示例中具有同样好的准确性。有关详细信息，请参阅此文档http://www.cs.princeton.edu/courses/archive/fall11/cos323/notes/cos323_f11_lecture09_svd.pdf

于 2016-04-10T21:12:55.337 回答
2

试试这个称为块 SVD 的 SVD 版本——你只需将迭代设置为你想要的精度——通常 1 就足够了。如果你想要所有的因素（这有一个默认的 # 选择用于减少因素）然后如果我正确回忆我的 MATLAB，则将 k= 行编辑为 size(matrix)

A= randn(100,5000);
A=corr(A);
% A is your correlation matrix
tic
k = 1000; % number of factors to extract
bsize = k +50;
block = randn(size(A,2),bsize);
iter = 2; % could set via tolerance

[block,R] = qr(A*block,0);
for i=1:iter
    [block,R] = qr(A*(A'*block),0);
end
M = block'*A;
% Economy size dense SVD.
[U,S] = svd(M,0);
U = block*U(:,1:k);
S = S(1:k,1:k);
% Note SVD of a symmetric matrix is:
% A = U*S*U' since V=U in this case, S=eigenvalues, U=eigenvectors
V=real(U*sqrt(S)); %scaling matrix for simulation
toc
% reduced randomized matrix for simulation
sims = 2000;
randnums = randn(k,sims);
corrrandnums = V*randnums;
est_corr_matrix = corr(corrrandnums');
total_corrmatrix_difference =sum(sum(est_corr_matrix-A))
于 2016-04-15T01:03:48.600 回答