inv(A)@b vs solve: measured accuracy and speed, NumPy/OpenBLAS
n=500 with random orthogonal U,V and geometric singular values, cond 1e2 to 1e14. np.linalg.solve gives forward error 2e-14/3e-11/1e-7/9e-4 and inv(A)@b gives 1e-13/5e-10/4e-6/3e-2, so inv is 6-33x worse. Solve's scaled residual stays around 1e-16 (backward stable), while inv's grows with the condition number (1.7e-14 up to 2.6e-4). Speed at n=2000 on 2 OpenBLAS threads: solve 82 ms against 232 ms for inv@b with 1-50 right-hand sides. The one case where inv wins is applying the same A to many single vectors one at a time: inv@b took 0.48 ms per vector, but scipy.linalg.lu_solve took 1.0 ms even with check_finite=False, because the triangular solves are slower than a matvec. Batch the right-hand sides into one matrix instead: 0.086 ms per vector after a 67 ms lu_factor.
Author self-check with 5 seeds per condition number and min-of-5 timings (NumPy 2.5.4, OpenBLAS). Accuracy holds: the inv/solve forward-error ratio was 4.7-6.8x at cond 1e2, 12-20x at 1e6, 19-28x at 1e10 and 22-33x at 1e14. solve's scaled residual stays at 1.5e-16 to 1.9e-15, inv's median is 1.7e-14 at cond 1e2 and 2e-4 at 1e14. One nuance on the per-vector case: it is a threading effect. With OPENBLAS_NUM_THREADS=1, inv@v took 0.96 ms against 1.04 ms for lu_solve one at a time, so a tie. With 2 threads inv@v was 0.48 ms and lu_solve still 1.03 ms, because the matvec parallelises and the triangular solves don't. Batched lu_solve is 0.15 ms (1 thread) or 0.08 ms (2 threads) per vector, and that remains the answer. n=2000 solve vs inv@b: 98 vs 403 ms on 1 thread, 65 vs 246 ms on 2.