R、C、Python中Moore-Penrose伪逆结果差异原因及准确性问询

Let's break down why you're seeing these discrepancies, even though all three implementations are rooted in SVD, and how to evaluate which result is more numerically accurate.

Key Reasons for Numerical Differences #

While the Moore-Penrose pseudoinverse is theoretically unique, practical implementations can diverge dramatically—especially for highly ill-conditioned matrices like your correlation matrix (all values are near 1, giving it an extremely large condition number). Here's the breakdown:

1. Singular Value Truncation (Thresholding) #

All robust pseudoinverse implementations truncate tiny singular values (to avoid dividing by near-zero numbers, which causes catastrophic numerical instability). But each library uses a slightly different threshold calculation:

• R's corpcor::pseudoinverse: Uses a default threshold of tol = max(dim(X)) * max(singular_values) * .Machine$double.eps (where .Machine$double.eps is ~2.2e-16 for double-precision floats).
• NumPy's numpy.linalg.pinv: Defaults to tol = max(M, N) * eps * max(s) (similar in theory, but may use subtle differences in eps definition or implementation).
• Your C implementation: If you're not applying the same threshold logic as R/NumPy, or using a different formula for tol, this will drastically alter the pseudoinverse result. For ill-conditioned matrices, even tiny changes in which singular values you keep/discard lead to huge swings in the final matrix.

2. Underlying SVD Implementation Differences #

Different libraries rely on distinct BLAS/LAPACK backends for SVD computation:

• R typically uses the reference LAPACK implementation (or a tuned version like OpenBLAS) for corpcor.
• NumPy may use OpenBLAS, MKL, or another optimized BLAS/LAPACK library, which can produce slightly different singular values/vectors due to algorithmic optimizations (e.g., pivot strategies, floating-point accumulation order).
• If your C implementation uses a custom SVD routine or a different LAPACK version, the raw SVD outputs (U, S, V) will have subtle differences that get amplified when computing the pseudoinverse (since we take 1/S for non-truncated values).

3. Floating-Point Precision & Computation Order #

Floating-point arithmetic has inherent rounding errors. For ill-conditioned matrices, these errors get magnified exponentially. Small differences in the order of operations (e.g., how matrix multiplications are accumulated) or precision handling between libraries can lead to noticeable differences in the final pseudoinverse.

How to Assess Which Result Is More Accurate #

Don't judge accuracy by raw numerical differences—instead, verify which implementation satisfies the Moore-Penrose conditions (the formal definition of a pseudoinverse) with the smallest residual error:
For a matrix ( A ) and its pseudoinverse ( A^+ ), the four conditions are:

1. ( A A^+ A = A )
2. ( A^+ A A^+ = A^+ )
3. ( (A A⁺⁾T = A A^+ )
4. ( (A^+ A)^T = A^+ A )

Calculate the Frobenius norm of the residuals for each condition (e.g., norm(A %*% A_plus %*% A - A, "F") in R). The implementation with the smallest residuals is the most accurate.

Additionally, check if ( A^+ A ) and ( A A^+ ) are orthogonal projection matrices—their eigenvalues should be either 0 or 1 (up to acceptable numerical noise).

Fixing the C Implementation to Match R #

To align your C result with corpcor::pseudoinverse:

1. Replicate the threshold logic: Implement the same tol calculation as R:

   // Example: Compute tol like corpcor
   int rows = 6, cols = 6;
   double max_s = /* maximum singular value from your SVD */;
   double eps = 2.220446049250313e-16; // .Machine$double.eps in R
   double tol = (double)fmax(rows, cols) * max_s * eps;
   // Then, set any singular value s < tol to 0 (or skip dividing by them)
   

2. Use the same BLAS/LAPACK backend: Link your C code against the same LAPACK library that R uses (e.g., reference LAPACK or OpenBLAS) to ensure consistent SVD outputs.

内容的提问来源于stack exchange，提问作者David Morocho