matlab - 我的例子表明 SVD 在数值上不如 QR 分解稳定

Question

我在 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 得到比 QR(e-27) 大得多的误差 (0.8)

% 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

Answer 1

您基于 SVD 的方法与 pinv function 基本相同在 MATLAB 中(参见 Pseudo-inverse and 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