Hi everyone, I want to vectorize multiple matrix-vector products and avoid a for loop, a little bit like np.linalg.solve does, for instance : /nDOF = 100// //M = 4// // //mat = np.random.rand(nDOF, M, M)// //u = np.random.rand(nDOF, M)/ /vecMatInv = np.linalg.solve(mat, u)  # => shape (nDOF, M) / /for i in range(nDOF):// //        assert np.allclose(vecMatInv[i], np.linalg.solve(mat[i], u[i]))  # => OK / This is quite straightforward and also works if u is only a flat vector of size M (which is great). But this is not that easy for matrix vector multiplication, and the solution I found yet is to use np.matmul with some manipulations over the array dimensions : vecMatMul = np.squeeze(np.matmul(mat, u[..., None]), axis=-1) /# => shape (nDOF, M)/ for i in range(nDOF):         assert np.allclose(vecMatMul[i], mat[i] @ u[i])  # => OK That way we get the same behavior as for np.linalg.solve. But this seems not very natural, either for np.matmul or np.dot when using arrays of dimension > 2. Since np.linalg.solve does this vectorization naturally, I wonder if there is a way to get the same behavior with matrix-vector multiplication already in numpy, and why np.matmul does not behave like np.linalg.solve does ? Best, Thibaut -- *Dr. Thibaut LUNET* Chair Computational Mathematics, /Institute of Mathematics (E-10),/ /Office 3.040, tel: +49 40 42878 3428/ *Hamburg University of Technology*