A lot of numerical solvers (such as this one), for some reason, only supports solving Ax=b where x, b are vectors.

What if I want to solve AX = B, where X and B are matrices.

Looping would have overhead, and I dont want to directly use numpy.linalg.solve because I do want to use the special feature of Ax=b solver. Is there a feature in python that speeds this up?