R ginv và Matlab pinv tạo ra các kết quả khác nhau

Question

ginv() chức năng từ MASS gói trong R tạo ra các giá trị hoàn toàn khác so với hàm MATLAB pinv(). Cả hai đều yêu cầu sản xuất nghịch đảo tổng quát Moore-Penrose của một ma trận.

Tôi đã cố gắng đặt cùng một dung sai cho việc triển khai R nhưng sự khác biệt vẫn tồn tại.

• MATLAB tol mặc định: max(size(A)) * norm(A) * eps(class(A))
• R mặc định tol: sqrt(.Machine$double.eps)

sinh sản:

R:

library(MASS) 
A <- matrix(c(47,94032,149, 94032, 217179406,313679,149,313679,499),3,3) 
ginv(A) 

kết quả đầu ra:

012.

   [,1]   [,2]   [,3] 
[1,] 1.675667e-03 -8.735203e-06 5.545605e-03 
[2,] -8.735203e-06 5.014084e-08 -2.890907e-05 
[3,] 5.545605e-03 -2.890907e-05 1.835313e-02 

svd(A)

kết quả đầu ra:

$d 
[1] 2.171799e+08 4.992800e+01 2.302544e+00 

$u 
       [,1]   [,2]   [,3] 
[1,] -0.0004329688 0.289245088 -9.572550e-01 
[2,] -0.9999988632 -0.001507826 -3.304234e-06 
[3,] -0.0014443299 0.957253888 2.892454e-01 

$v 
       [,1]   [,2]   [,3] 
[1,] -0.0004329688 0.289245088 -9.572550e-01 
[2,] -0.9999988632 -0.001507826 -3.304234e-06 
[3,] -0.0014443299 0.957253888 2.892454e-01 

MATLAB:

A = [47 94032 149; 94032 217179406 313679; 149 313679 499] 
pinv(A) 

kết quả đầu ra:

ans = 

    0.3996 -0.0000 -0.1147 
    -0.0000 0.0000 -0.0000 
    -0.1147 -0.0000 0.0547 

svd:

[U, S, V] = svd(A) 

U = 

    -0.0004 0.2892 -0.9573 
    -1.0000 -0.0015 -0.0000 
    -0.0014 0.9573 0.2892 


S = 

    1.0e+008 * 

    2.1718   0   0 
     0 0.0000   0 
     0   0 0.0000 


V = 

    -0.0004 0.2892 -0.9573 
    -1.0000 -0.0015 -0.0000 
    -0.0014 0.9573 0.2892 

Answer 1

Chạy debugonce(MASS::ginv), chúng ta thấy rằng sự khác biệt nằm ở chỗ những gì được thực hiện với sự phân hủy giá trị duy nhất.

Cụ thể, R kiểm tra như sau:

Xsvd <- svd(A) 
Positive <- Xsvd$d > max(tol * Xsvd$d[1L], 0) 
Positive 
# [1] TRUE TRUE FALSE 

Nếu các yếu tố thứ ba là đúng (mà chúng tôi có thể buộc bằng cách thiết lập tol = 0, theo đề nghị của @nicola), MASS::ginv sẽ trở lại:

Xsvd$v %*% (1/Xsvd$d * t(Xsvd$u)) 
#    [,1]   [,2]   [,3] 
# [1,] 3.996430e-01 -7.361507e-06 -1.147047e-01 
# [2,] -7.361507e-06 5.014558e-08 -2.932415e-05 
# [3,] -1.147047e-01 -2.932415e-05 5.468812e-02 

(ví dụ, giống như MATLAB).

Thay vào đó, nó sẽ trả về:

Xsvd$v[, Positive, drop = FALSE] %*% ((1/Xsvd$d[Positive]) * 
    t(Xsvd$u[, Positive, drop = FALSE])) 
#    [,1]   [,2]   [,3] 
# [1,] 1.675667e-03 -8.735203e-06 5.545605e-03 
# [2,] -8.735203e-06 5.014084e-08 -2.890907e-05 
# [3,] 5.545605e-03 -2.890907e-05 1.835313e-02 

Nhờ @FaridCher để chỉ ra mã nguồn của pinv.

Tôi không chắc mình hiểu 100% mã MATLAB, nhưng tôi nghĩ nó có sự khác biệt về cách sử dụng tol.Các MATLAB thư để Positive trong R là:

`r = sum(s>tol)` 

đâu tol là những gì được cung cấp bởi người sử dụng; nếu không được cung cấp, chúng tôi nhận được:

m = 0; 
% I don't get the point of this for loop -- why not just `m = max(size(A))`? 
for i = 1:n 
    m = max(m,length(A(:,i))); 
end 
% contrast with simply `tol * Xsvd$d[1L]` in R 
% (note: i believe the elements of d are sorted largest to smallest) 
tol = m*eps(max(s));