function X = MatLS(A1,A2,B)
% A1 and A2 are each m-by-n with m > n = rank(A1) = rank(A2)
% B is m-by-m
% X minimizes norm(A1*X*A2' - B,'fro')

[U1,S1,V1] = svd(A1,0); % -1 if svd(A1), svd(A2)
[U2,S2,V2] = svd(A2,0);
% Since
% norm(A1*X*A2' - B,'fro') = norm(U1*S1*V1)*X*(U2*S2*V2')' - B,'fro')
% = norm(S1(V1*X*V2')S2 - U1'*B*U2,'fro')
% = norm(S1*Y*S2 - C,'fro')
% where
C = U1'*B*U2;
% and Y = inv(S1)*C*inv(S2), i.e.,
Y = C./(diag(S1)*diag(S2)');
X = V1*Y*V2';