Code covered by the BSD License

# Least-square with 2-norm constraint

### Bruno Luong (view profile)

Minimize |A*x-b|^2 such that |x| = cte

spherelsq.m
```function x = spherelsq(A,b,xnorm)
% x = spherelsq(A,b,xnorm)
%
% Minimize J(x) = |A*x-b|^2 such that |x| = xnorm
%
% A is matrix, b is a vector having the same number of rows
% xnorm is a positive scalar. |.| denotes l^2 norm.
%
% Reference: W. Gander, G. H. Golub, and U. Von Matt, "A constrained
% eigenvalue problem", Linear Algebra Appl., 114-115 (1989), pp. 815839.
%
% Author: Bruno Luong <brunoluong@yahoo.com>
% History
%   12-May-2010: original version
%

if size(A,1) ~= size(b,1)
error('spherelsq: sizes of A and b do not match');
end

if size(b,2) ~= 1
error('spherelsq: cannot handle multiple RHSs');
end

if ~isscalar(xnorm) || xnorm<0
error('spherelsq: xnorm must be positive scalar');
end
% ignore the imaginary part
xnorm = real(xnorm);
if xnorm==0
x = zeros(size(A,2),1);
return
end

H = A'*A;
g = A'*b;
A2 = speye(size(H));
A1 = -2*H;
gn = g/xnorm;
A0 = H*H - gn*gn';
% QEP
[V lambda] = polyeig(A0,A1,A2); %#ok
% smallest real eigenvalue
lambda = lambda(imag(lambda)==0);
if isempty(lambda)
error('spherelsq: problem with floating point accuracy');
end
lambda = min(lambda);

% solution
x = (H-lambda*speye(size(H))) \ g;

end % spherelsq```