In my opinion you are better off creating a pseudo inverse from the singular value decomposition. Moore-Penrose pseudo inverse will work, but I think it gives some weird results sometimes. The least squares can be unstable, especially since the rank of your example matrix is not full.
Also; don't calculate the inverse for every iteration!
Here is an example, where rank-deficiency of A is taken care of:
A = [2,1,5;-2,-1,-5];
% pseudo inverse of A
[m,n] = size(A);
% get SVD
[U,S,V] = svd(A);
% check rank
r = rank(S);
SR = S(1:r,1:r);
% make complete if rank is not full
SRc = [SR^-1 zeros(r,m-r);zeros(n-r,r) zeros(n-r,m-r)];
% create inverse
A_inv = V*SRc*U.';
X=[];
for i = 1:1600
% this actually takes most of the time
k31 = rand(1000, 1);
k32 = rand(1000, 1);
B = [k31';k32'];
% X is overwritten on every loop...
X = [X A_inv*B];
end
N_xy = 1000;
x = rand(N_xy,1);
y = rand(N_xy,1);
g1 = X(1,1:N_xy);
figure(1), clf
scatter(x,y,5,g1)
I didn't plot all 1600000 points since that can't be what you want