numpy.linalg.pinv精度不足，如何准确恢复目标矩阵m2？

Great question! Let's break down why you're seeing poor precision in the first two columns of your recovered m2, then walk through practical, actionable solutions tailored to your matrix structure.

Root Cause of the Precision Problem #

Your matrix m1 is a low-rank singular matrix (rank = 2, to be exact):

• The first two rows are random sequences of 2s and 3s, which are almost certainly linearly independent.
• The remaining 126 rows are incremented row-wise/column-wise, making them fully linearly dependent on each other (and not adding new information beyond the first two rows).

When you use np.linalg.pinv(m1) on the full 128x128 matrix, those 126 redundant rows introduce numerical noise that gets amplified during the pseudo-inverse calculation. This noise directly impacts the precision of the first two columns of m2, since those columns are tied to the only meaningful information in m1 (its first two rows).

Practical Solutions #

1. Extract the Linear Independent Rows of m1 (Core Fix) #

Since the first two rows of m1 contain all the unique information, we can discard the redundant rows and solve the problem with a smaller, well-conditioned submatrix:

# Extract the only meaningful rows from m1 and corresponding rows from m3
m1_core = m1[:2, :]
m3_core = m3[:2, :]

# Recover m2 using the pseudo-inverse of the core matrix
m2_recovered = np.linalg.pinv(m1_core) @ m3_core

# Alternatively, use least squares (equivalent result, sometimes more stable)
# m2_recovered, _, _, _ = np.linalg.lstsq(m1_core, m3_core, rcond=None)

This eliminates the noise from redundant rows entirely, and you'll immediately see much better precision in the first two columns of m2.

2. Refine with m2's Known Structure #

You know two key details about m2—use them to polish the recovered result:

• First two columns: Values are random 2s and 3s. Round or clip the recovered values to match this:

  for col in range(2):
      # Round to nearest integer (2 or 3)
      m2_recovered[:, col] = np.round(m2_recovered[:, col])
      # Ensure no out-of-bounds values (edge case handling)
      m2_recovered[:, col] = np.clip(m2_recovered[:, col], 2, 3)
  

• Remaining columns: Values are row-wise incrementing. Fit a linear trend to each column to fix minor numerical deviations:

  for col in range(2, 128):
      row_indices = np.arange(128)
      # Fit a line to the column values
      slope, intercept = np.polyfit(row_indices, m2_recovered[:, col], 1)
      # Replace the column with the clean linear sequence
      m2_recovered[:, col] = slope * row_indices + intercept
  

3. Boost Numerical Precision (If Needed) #

If double-precision (float64) still isn't enough, try using higher-precision floats (note: this requires CPU support for float128):

# Convert matrices to float128 for higher-precision calculation
m1_core = m1[:2, :].astype(np.float128)
m3_core = m3[:2, :].astype(np.float128)

# Compute pseudo-inverse and convert back to float64 for regular use
m2_recovered = np.linalg.pinv(m1_core) @ m3_core
m2_recovered = m2_recovered.astype(np.float64)

Why Your Previous Attempts Didn't Work #

• Adjusting rcond: The redundant rows of m1 have near-zero singular values, so changing rcond won't preserve them—they'll still be filtered out, and the noise remains.
• mpmath: This library lacks native pseudo-inverse support, and integrating it with numpy's arrays introduces more complexity than it's worth. The structural fixes above are far more efficient.

内容的提问来源于stack exchange，提问作者Blackyu Sylvean