You can use the regress function to do this.
Here is an example:
% Generate an example
n = 60;
theta = rand(n);
% Create regressors
[M,N] = meshgrid(1:n,1:n);
X = [M(:), N(:)];
% Regress
B=regress(theta(:), X);
% Compare the results
theta_hat = reshape(X*B,n,n);
plot3(M,N,theta,'o');
hold on;
surf(M,N,theta_hat);
Notice that the regression is done on theta(:)
which is a (3600,1) vector containing the values of theta(k1,k2) uses the corresponding coordinates in X which is (3600,2). The first column of X is k1, the second is k2.
The result of calling regress gives you B=[a;b]
the coefficients that best fit the data in theta.
One final note is that the least squares could be solved directly using
B=inv(X'*X)*X'*theta(:)
which should give the same result, but regress
is the preferred MATLAB method.