Мой пример показывает, что SVD менее численно устойчив, чем QR-разложение.

Я задал этот вопрос на 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

matlab least-squares svd matrix-decomposition qr-decomposition
person daydayup    schedule 06.04.2016    source источник

Ответы (3)


arrow_upward
6
arrow_downward

Ваш подход на основе SVD в основном такой же, как pinv function в MATLAB. (см. Псевдоинверсия и 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) очень мало рядом с машинным эпсилоном):

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
person Amro    schedule 11.04.2016
comment
кстати, если вы будете следовать диаграмме в разделе алгоритма mldivide, вы увидите, что в вашем случае не обязательно выполняется декомпозиция QR : mathworks.com/help/matlab/ref/mldivide.html#moreabout - person Amro; 11.04.2016

arrow_upward
2
arrow_downward

СВД может справиться с дефицитом ранга. Диагональная матрица D имеет почти нулевой элемент в вашем коде, и вам нужно использовать псевдоинверсию для SVD, т.е. установить 2-й элемент 1./diag(D) равным 0, отличному от огромного значения (10^14). Вы должны обнаружить, что SVD и QR имеют одинаковую точность в вашем примере. Для получения дополнительной информации см. этот документ http://www.cs.princeton.edu/courses/archive/fall11/cos323/notes/cos323_f11_lecture09_svd.pdf

person uPhone    schedule 10.04.2016

arrow_upward
2
arrow_downward

Попробуйте эту версию SVD, называемую блочной SVD — вы просто устанавливаете итерации, равные нужной точности — обычно 1 достаточно. Если вам нужны все факторы (по умолчанию # выбрано для уменьшения фактора), тогда отредактируйте строку k= до размера (матрицы), если я правильно вспомню свой MATLAB

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))
person Matt    schedule 15.04.2016