Inverse a non-square matrix with high precision

Question

I am trying to solve a problem of A*b = c, where A and c are:

A = np.block([[X_matrix - rho_0, np.zeros((3, 3)), np.zeros((3, 3))],
              [np.zeros((3, 3)), X_matrix - rho_1, np.zeros((3, 3))],
              [np.zeros((3, 3)), np.zeros((3, 3)), X_matrix - rho_2],
              [I, I, I]])

c = np.block([[np.zeros((3, 3))],
              [np.zeros((3, 3))],
              [np.zeros((3, 3))],
              [I]])

Since A is a non-square matrix here, I can only use

A_inv = np.linalg.pinv(A)
# Solve for b
b = np.matmul(A_inv, c)

However, this method is not very precise as print(np.linalg.norm(np.matmul(A, b) - c)) will give 0.212871998042824. May I ask whether there is a more precise method to obtain b?

I have tried np.linalg.lstsq but it seems to give the same result as np.linalg.pinv(A)

Also, it seems that A**-1 from mpmath will arrive at the error: ValueError: only powers of square matrices are defined.