wolffd@0: function m = sqdist(p, q, A)
wolffd@0: % SQDIST      Squared Euclidean or Mahalanobis distance.
wolffd@0: % SQDIST(p,q)   returns m(i,j) = (p(:,i) - q(:,j))'*(p(:,i) - q(:,j)).
wolffd@0: % SQDIST(p,q,A) returns m(i,j) = (p(:,i) - q(:,j))'*A*(p(:,i) - q(:,j)).
wolffd@0: 
wolffd@0: %  From Tom Minka's lightspeed toolbox
wolffd@0: 
wolffd@0: [d, pn] = size(p);
wolffd@0: [d, qn] = size(q);
wolffd@0: 
wolffd@0: if nargin == 2
wolffd@0:   
wolffd@0:   pmag = sum(p .* p, 1);
wolffd@0:   qmag = sum(q .* q, 1);
wolffd@0:   m = repmat(qmag, pn, 1) + repmat(pmag', 1, qn) - 2*p'*q;
wolffd@0:   %m = ones(pn,1)*qmag + pmag'*ones(1,qn) - 2*p'*q;
wolffd@0:   
wolffd@0: else
wolffd@0: 
wolffd@0:   if isempty(A) | isempty(p)
wolffd@0:     error('sqdist: empty matrices');
wolffd@0:   end
wolffd@0:   Ap = A*p;
wolffd@0:   Aq = A*q;
wolffd@0:   pmag = sum(p .* Ap, 1);
wolffd@0:   qmag = sum(q .* Aq, 1);
wolffd@0:   m = repmat(qmag, pn, 1) + repmat(pmag', 1, qn) - 2*p'*Aq;
wolffd@0:   
wolffd@0: end