So I have this line of code that involves a matrix inversion

X = A @ B @ np.linalg.pinv(S)

A is an n by n matrix, B is an n by m matrix and S is an m by m matrix. m is smaller than n but usually not orders of magnitude smaller. Usually m is about half of n. S is a symmetrical matrix. How do I make this line of code run faster in python?

I know one should never invert a matrix, but all the examples I have seen are about solving a linear system Ax = b, not in my case when I am solving XS = C or SX’ = C’ (if I just let C = AB, and $S$ is symmetrical). I know it’s just conceptually similar to solving SX_i = C_i where X_i is the i-th column of X. But it seems difficult to glue them together without a for-loop, which tends to be slow in python?

Or are there other tricks I missed? Note I am more looking for a practical python solution. Or is it worthwhile at all given m is half as big as n?