function x = lse(A,b,C,d,solverflag,weights)
% solves A*x = b for X in a least squares sense, given C*x = d
% usage: x = lse(A,b,C,d)
% usage: x = lse(A,b,C,d,solverflag)
% usage: x = lse(A,b,C,d,solverflag,weights)
%
% Minimizes norm(A*x - b),
% subject to C*x = d
%
% If b has multiple columns, then so will x.
%
% arguments: (input)
% A - nxp array, for the least squares problem
% A may be a sparse matrix, it need not be
% of full rank.
%
% b - nx1 vector (or nxq array) of right hand
% side(s) for the least squares problem
%
% C - mxp array for the constraint system. C
% must be a full matrix (not sparse), but
% it may be rank deficient.
%
% C (and d) may be empty, in which case no
% constraints are applied
%
% d - mx1 vector - right hand side for the
% constraint system.
%
% solverflag - (OPTIONAL) - character flag -
% Specifies the basic style of solver used
% for the least squares problem. Use of the
% pinv solution will produce a minimum norm
% solution when A is itself singular.
%
% solverflag may be any of
%
% {'\', 'pinv', 'backslash'}
%
% Capitalization is ignored, and any
% shortening of the string is allowed,
% as far as {'\', 'p', 'b'}
%
% DEFAULT: '\'
%
% weights - nx1 vector of weights for the
% regression problem. All weights must be
% non-negative. A weight of k is equivalent
% to having replicated that data point by
% a factor of k times.
%
%
% arguments: (output)
% x - px1 vector (or pxq array) of solutions to
% the least squares problem A*x = b, subject
% to the linear equality constraint system
% C*x = d
%
%
% Example usage:
% A = rand(10,3);
% b = rand(10,2); % two right hand sides
% C = [1 1 1;1 -1 0.5];
% d = [1;0.5];
%
% X = lse(A,b,C,d)
% X =
% 0.71593 0.55371
% 0.23864 0.18457
% 0.045427 0.26172
%
% As a test, we should recover the constraint
% right hand side d for each solution X.
%
% C*X
% ans =
% 1 1
% 0.5 0.5
%
%
% Example usage:
% A = rand(10,3);
% b = rand(10,1);
%
% with a rank deficient constraint system
% C = [1 1 1;1 1 1];
% d = [1;1];
%
% X = lse(A,b,C,d)
% X =
% 0.5107
% 0.57451
% -0.085212
%
% C*X
% ans =
% 1
% 1
%
%
% Example usage: (where both A and C are rank deficient)
% A = rand(10,2);
% A = [A,A];
% b = rand(10,1);
%
% C = ones(2,4);
% d = [1;1];
%
% The \ solution will see the singularity in A
% X = lse(A,b,C,d,'\')
% Warning: Rank deficient, rank = 1, tol = 3.1821e-15.
% > In lse at 205
% X =
% 0.17097
% 0.82903
% 0
% 0
%
% The pinv version will survive the singulatity
% Xp = lse(A,b,C,d,'pinv')
% Xp =
% 0.085486
% 0.41451
% 0.085486
% 0.41451
%
% Use of pinv will produce the minimum norm solution
% norm(X)
% ans =
% 0.84647
%
% norm(Xp)
% ans =
% 0.59855
%
%
% Methodology:
% Both alternative methods offered in lse are effectively subspace
% methods. That is, the equality constraints are used to reduce the
% problem to a lower dimensional problem. The '\' method uses a pivoted
% QR factorization to choose a set of variables to be eliminated. This
% chooses the best subset of variables for elimination, avoiding small
% "pivots" where possible, as well as resolving the case where the
% supplied constraints are rank deficient. Note that when the constraint
% system is rank deficient, this method will result in one or more of
% the unknowns to be set to zero. An at length description of the QR
% based method for linear least squares subject to linear equality
% constraints is found in:
%
% http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=8553&objectType=FILE
%
% The 'pinv' method uses a related approach, reducing the problem to
% a least squares solution in a subspace. This method is chosen to be
% consistent with pinv as used for unconstrained least squares problems.
% The reduction to a subspace does not require the selection of specific
% variables to be eliminated. Instead the reduction is a projection as
% defined by a singular value decomposition. When the constraint system
% is rank deficient, the svd allows for a minimum norm solution, much
% as is done with pinv. This may be preferable for some users.
%
% Both methods allow the application of weights if supplied. As well,
% problems with multiple right hand sides (b) are solved in a fully
% vectorized fashion.
%
%
% See also: slash, lsqlin, lsequal
%
%
% Author: John D'Errico
% E-mail: woodchips@rochester.rr.com
% Release: 3.0
% Release date: 1/31/07
% check sizes
[n,p] = size(A);
[r,nrhs] = size(b);
[m,ccols] = size(C);
if n~=r
error 'A and b are incompatible in size (wrong number of rows)'
elseif ~isempty(C) && (p~=ccols)
error 'A and C must have the same number of columns'
elseif ~isempty(C) && issparse(C)
error 'C may not be a sparse matrix'
elseif ~isempty(C) && (m~=size(d,1))
error 'C and d are incompatible in size (wrong number of rows)'
elseif ~isempty(C) && (size(d,2)~=1)
error 'd must have only one column'
elseif isempty(C) && ~isempty(d)
error 'C and d are inconsistent with each other (one was empty)'
elseif ~isempty(C) && isempty(d)
error 'C and d are inconsistent with each other (one was empty)'
end
% default solver is '\'
if (nargin<5) || isempty(solverflag)
solverflag = '\';
elseif ~ischar(solverflag)
error 'If supplied, solverflag must be character'
else
% solverflag was supplied. Make sure it is legal.
valid = {'\', 'backslash', 'pinv'};
ind = strmatch(solverflag,valid);
if (length(ind)==1)
solverflag = valid{ind};
else
error(['Invalid solverflag: ',solverflag])
end
end
% default for weights = []
if (nargin<6) || isempty(weights)
weights = [];
else
weights = weights(:);
if (length(weights)~=n) || any(weights<0)
error 'weights should be empty or a non-negative vector of length n'
elseif all(weights==0)
error 'At least some of the weights must be non-zero'
else
% weights was supplied. scale it to have mean value = 1
weights = weights./mean(weights);
% also sqrt the weights for application as an
% effective replication factor. remember that
% least squares will minimize the sum of squares.
weights = sqrt(weights);
end
end
% tolerance used on the solve
Ctol = 1.e-13;
if (nargin<3) || isempty(C)
% solve with \ or pinv as desired.
switch solverflag
case {'\' 'backslash'}
% solve with or without weights
if isempty(weights)
x = A\b;
else
x = (repmat(weights,1,size(A,2)).*A)\ ...
(repmat(weights,1,nrhs).*b);
end
case 'pinv'
% solve with or without weights
if isempty(weights)
ptol = Ctol*norm(A,1);
x = pinv(A,ptol)*b;
else
Aw = repmat(weights,1,size(A,2)).*A;
ptol = Ctol*norm(Aw,1);
x = pinv(Aw,ptol)*(repmat(weights,1,nrhs).*b);
end
end
% no Constraints, so we are done here.
return
end
% Which solver do we use?
switch solverflag
case {'\' 'backslash'}
% allow a rank deficient equality constraint matrix
% column pivoted qr to eliminate variables
[Q,R,E]=qr(C,0);
% get the numerical rank of R (and therefore C)
if m == 1
rdiag = R(1,1);
else
rdiag = abs(diag(R));
end
crank = sum((rdiag(1)*Ctol) <= rdiag);
if crank >= p
error 'Overly constrained problem.'
end
% check for consistency in the constraints in
% the event of rank deficiency in the constraint
% system
if crank < m
k = Q(:,(crank+1):end)'*d;
if any(k > (Ctol*norm(d)));
error 'The constraint system is deficient and numerically inconsistent'
end
end
% only need the first crank columns of Q
qpd = Q(:,1:crank)'*d;
% which columns of A (variables) will we eliminate?
j_subs = E(1:crank);
% those that remain will be estimated
j_est = E((crank+1):p);
r1 = R(1:crank,1:crank);
r2 = R(1:crank,(crank+1):p);
A1 = A(:,j_subs);
qpd = qpd(1:crank,:);
% note that \ is still ok here, even if pinv
% is used for the main regression.
bmod = b-A1*(r1\repmat(qpd,1,nrhs));
Amod = A(:,j_est)-A1*(r1\r2);
% now solve the reduced problem, with or without weights
if isempty(weights)
x2 = Amod\bmod;
else
x2 = (repmat(weights,1,size(Amod,2)).*Amod)\ ...
(repmat(weights,1,nrhs).*bmod);
end
% recover eliminated unknowns
x1 = r1\(repmat(qpd,1,nrhs)-r2*x2);
% stuff all estimated parameters into final vector
x = zeros(p,nrhs);
x(j_est,:) = x2;
x(j_subs,:) = x1;
case 'pinv'
% allow a rank deficient equality constraint matrix
Ctol = 1e-13;
% use svd to deal with the variables
[U,S,V] = svd(C,0);
% get the numerical rank of S (and therefore C)
if m == 1
sdiag = S(1,1);
else
sdiag = diag(S);
end
crank = sum((sdiag(1)*Ctol) <= sdiag);
if crank >= p
error 'Overly constrained problem.'
end
% check for consistency in the constraints in
% the event of rank deficiency in the constraint
% system
if crank < m
k = U(:,(crank+1):end)'*d;
if any(k > (Ctol*norm(d)));
error 'The constraint system is deficient and numerically inconsistent'
end
end
% only need the first crank columns of U, and the
% effectively non-zero diagonal elements of S.
sinv = diag(S);
sinv = diag(1./sinv(1:crank));
% we will use a transformation
% Z = V'*X = inv(S)*U'*d
Z = sinv*U(:,1:crank)'*d;
% Rather than explicitly dropping columns of A, we will
% work in a reduced space as defined by the svd.
Atrans = A*V;
% thus, solve (A*V)*Z = b, subject to the constraints Z = supd
% use pinv for the solution here.
ptol = Ctol*norm(Atrans(:,(crank+1):end),1);
if isempty(weights)
Zsolve = pinv(Atrans(:,(crank+1):end),ptol)* ...
(b - repmat(Atrans(:,1:crank)*Z(1:crank),1,nrhs));
else
w = spdiags(weights,0,n,n);
Zsolve = pinv(w*Atrans(:,(crank+1):end),ptol)* ...
(w*(b - repmat(Atrans(:,1:crank)*Z(1:crank),1,nrhs)));
end
% put it back together in the transformed state
Z = [repmat(Z(1:crank),1,nrhs);Zsolve];
% untransform back to the original variables
x = V*Z;
end