function [A,b,Error]=LinearNet(X,Y);
% function [A, b, Error]=LinearNet(X,Y)
%  Solves AX+b=Y for each column of X and Y- So X, Y are dim x numpts
%  If X is n x p, Y is m x p, then A is m x n and b is n x 1.

[n,p]=size(X);
[m,k]=size(Y);
if k~=p
    error('X and Y should be dim x number of points in LinearNet (X and Y have the same number of columns)');
end

Xhat=[X;ones(1,p)];
C=Y*pinv(Xhat);  
A=C(:,1:n);
b=C(:,n+1);
Error=norm(C*Xhat-Y);