Factor once, solve many

Cholesky costs n³/3 once; each later solve is two triangular substitutions at O(n²). So m right-hand sides cost O(n³ + mn²), not O(mn³).

np.linalg.inv(A) @ B is both slower and less accurate than factoring, and is never the right call.

Measured on the real digit kernel matrix with 200 right-hand sides:

• factor once with cho_factor/cho_solve — 1×
• explicit inverse — about 5× slower, and about 7× less accurate on the residual
• one np.linalg.solve per column — about 340× slower

The fastest option is also the most accurate one. That does not happen often, and it is worth remembering when it does.