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