Some least squares Matlab code

Related to a climate statistics post here.

function [tau,mu]=roman3(x,h,w);
% Finds tau and mu to minimize the sum of squares
% \sum_n \sum_m h(m,n)*(x(m,n)-tau(n)-mu(m))^2,
% subject to the constraint \sum_m w(m)*mu(m)=0. Here x and h are MxN,
% and w is Mx1. Normally, one would have non-negative h, positive w, and
% \sum_m w(m)=1. Upon return, tau will be 1xN, mu will be Mx1.
% The method uses a Lagrange multiplier (lambda) and Gaussian elimination.
% The matrix h must not contain an all-zero row or column.

[M,N]=size(x); % x and h have M rows, N columns
A=zeros(M+N+1,M+N+1); % A will contain the linear system matrix
A(1:N,1:N+M)=[diag(sum(h,1)) h']; % dL/dtau, N equations
A(N+1:N+M,:)=[h diag(sum(h,2)) -w/2]; % dL/dmu, M equations
A(end,N+1:N+M)=w'; % dL/dlambda, one equation (the constraint)
b=[sum(h.*x,1)'; sum(h.*x,2); 0]; % b is the "right hand side"
y=A\b; % solve Ay=b for y, using Gaussian elim.
tau=y(1:N)';
mu=y(N+1:N+M); % note the multiplier lambda is always zero
return

Ann says: Man, you beat me to it. That's just what I was thinking!