R ginv e Matlab pinv producono risultati diversi

Question

ginv() la funzione dal pacchetto MASS in R produce valori completamente diversi rispetto alla funzione MATLAB pinv(). Entrambi dichiarano di produrre un'inversione generalizzata di Moore-Penrose di una matrice.

Ho provato a impostare la stessa tolleranza per l'implementazione R ma la differenza persiste.

• MATLAB predefinito tol: max(size(A)) * norm(A) * eps(class(A))
• R predefinito tol: sqrt(.Machine$double.eps)

Riproduzione:

R:

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

uscite:

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)

uscite:

$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) 

uscite:

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

corso debugonce(MASS::ginv), vediamo che la differenza sta in ciò che si fa con la decomposizione in valori singolari.

In particolare, R verifica il seguente:

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

Se il terzo elemento fosse vero (che possiamo costringere impostando tol = 0, come suggerito da @nicola), MASS::ginv restituisce:

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 

(cioè uguale a MATLAB).

contrario, restituisce:

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 

Grazie a @FaridCher per indicare il codice sorgente di pinv.

Non sono sicuro di aver compreso al 100% il codice MATLAB, ma penso che si tratti di una differenza nel modo in cui viene utilizzato tol.La corrispondenza MATLAB per Positive in R è:

`r = sum(s>tol)` 

Dove tol è ciò che è fornita dall'utente; se nessuno viene fornito, otteniamo:

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));