Solve 3D least squares in numpy/scipy
In fact the answer was simple, I just needed to create bigger matrices Y and X by horizontally stacking the Y_k (to create Y) and the  X_k (to create X). Then I can just solve a regular 2d least squares problem: minimize norm(Y - A.dot(X))
For some integer K around 100, I have 2 * K (n, n) arrays: X_1, ..., X_K and Y_1, ..., Y_K.
I would like to perform K least squares simultaneously, i.e. find the n by n matrix A minimizing the sum of squares over k: \sum_k norm(Y_k - A.dot(X_k), ord='fro') ** 2 (A must not depend on k).
I am looking for an easy way to do this with numpy or scipy.
I know the function I want to minimize is a quadratic form in A so I could do it by hand, but I'm looking for an off-the-shelf way of doing it. Is there one?