function g = lsq (f, delta)

% LSQ    lsq smoothes the data in the vector f by means of a
%        Least Squares smoothing filter.
%        The data in the vector f is assumed to be sampled equidistantly.
%        lsq computes a smooth function g such that the mean
%        value of  (g(i) - f(i))^2  is less than or equal to delta^2.
%
%        W. Gander,
%        On the Linear Least Squares Problem with a Quadratic Constraint,
%        Habilitationsschrift ETH Zuerich, STAN-CS-78-697,
%        Stanford University, 1978, pp. 99-111.

n = length (f);
alpha = sqrt (n) * delta;
e = ones (n, 1);
A = spdiags ([e -2*e e], 0:2, n-2, n);

lambda = 0;
while 1,
   % z = [A'; sqrt (lambda) * eye (n-2)] \ [f; 0];
   [Qbar, R] = spqr (n, sqrt (lambda));
   z = QTx (Qbar, [f; zeros(n-2, 1)]);
   z = R \ z (1:n-2);

   F = norm (A' * z)^2;
   if (F <= alpha^2), break; end;

   if (lambda > 0),
      % zprime = [A'; sqrt (lambda) * eye (n-2)] \ [0; -z / sqrt (lambda)];
      zp = QTx (Qbar, [zeros(n, 1); -z/sqrt(lambda)]);
      zp = R \ zp (1:n-2);
      Fp = -2 * z' * (z + lambda * zp);
   else
      Fp = -2 * z' * z;
   end;

   % Newton step
   lambdaold = lambda;
   lambda = lambda - 2 * (F / Fp) * (sqrt (F) / alpha - 1);
   if (lambda <= lambdaold), break; end;
end;
g = f - A' * z;

% check Lagrange equations
% norm ((A*A' + lambda * speye(n-2)) * z - A*f) / norm (A*f)
% (norm (A'*z) - alpha) / alpha
