Home > atmlab > math > lse.m

lse

PURPOSE ^

solves A*x = b for X in a least squares sense, given C*x = d

SYNOPSIS ^

function x = lse(A,b,C,d,solverflag,weights)

DESCRIPTION ^

 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

DOWNLOAD ^

lse.m

SOURCE CODE ^

0001 function x = lse(A,b,C,d,solverflag,weights)
0002 % solves A*x = b for X in a least squares sense, given C*x = d
0003 % usage: x = lse(A,b,C,d)
0004 % usage: x = lse(A,b,C,d,solverflag)
0005 % usage: x = lse(A,b,C,d,solverflag,weights)
0006 %
0007 % Minimizes norm(A*x - b),
0008 % subject to C*x = d
0009 %
0010 % If b has multiple columns, then so will x.
0011 %
0012 % arguments: (input)
0013 %  A - nxp array, for the least squares problem
0014 %      A may be a sparse matrix, it need not be
0015 %      of full rank.
0016 %
0017 %  b - nx1 vector (or nxq array) of right hand
0018 %      side(s) for the least squares problem
0019 %
0020 %  C - mxp array for the constraint system. C
0021 %      must be a full matrix (not sparse), but
0022 %      it may be rank deficient.
0023 %
0024 %      C (and d) may be empty, in which case no
0025 %      constraints are applied
0026 %
0027 %  d - mx1 vector - right hand side for the
0028 %      constraint system.
0029 %
0030 %  solverflag - (OPTIONAL) - character flag -
0031 %      Specifies the basic style of solver used
0032 %      for the least squares problem. Use of the
0033 %      pinv solution will produce a minimum norm
0034 %      solution when A is itself singular.
0035 %
0036 %      solverflag may be any of
0037 %
0038 %      {'\', 'pinv', 'backslash'}
0039 %
0040 %      Capitalization is ignored, and any
0041 %      shortening of the string is allowed,
0042 %      as far as {'\', 'p', 'b'}
0043 %
0044 %      DEFAULT: '\'
0045 %
0046 %  weights - nx1 vector of weights for the
0047 %      regression problem. All weights must be
0048 %      non-negative. A weight of k is equivalent
0049 %      to having replicated that data point by
0050 %      a factor of k times.
0051 %
0052 %
0053 % arguments: (output)
0054 %  x - px1 vector (or pxq array) of solutions to
0055 %      the least squares problem A*x = b, subject
0056 %      to the linear equality constraint system
0057 %      C*x = d
0058 %
0059 %
0060 % Example usage:
0061 %  A = rand(10,3);
0062 %  b = rand(10,2); % two right hand sides
0063 %  C = [1 1 1;1 -1 0.5];
0064 %  d = [1;0.5];
0065 %
0066 %  X = lse(A,b,C,d)
0067 %  X =
0068 %      0.71593      0.55371
0069 %      0.23864      0.18457
0070 %     0.045427      0.26172
0071 %
0072 %  As a test, we should recover the constraint
0073 %  right hand side d for each solution X.
0074 %
0075 %  C*X
0076 %  ans =
0077 %            1            1
0078 %          0.5          0.5
0079 %
0080 %
0081 % Example usage:
0082 %  A = rand(10,3);
0083 %  b = rand(10,1);
0084 %
0085 % with a rank deficient constraint system
0086 %  C = [1 1 1;1 1 1];
0087 %  d = [1;1];
0088 %
0089 %  X = lse(A,b,C,d)
0090 %  X =
0091 %       0.5107
0092 %      0.57451
0093 %    -0.085212
0094 %
0095 %  C*X
0096 %  ans =
0097 %            1
0098 %            1
0099 %
0100 %
0101 % Example usage: (where both A and C are rank deficient)
0102 %  A = rand(10,2);
0103 %  A = [A,A];
0104 %  b = rand(10,1);
0105 %
0106 %  C = ones(2,4);
0107 %  d = [1;1];
0108 %
0109 % The \ solution will see the singularity in A
0110 %  X = lse(A,b,C,d,'\')
0111 %  Warning: Rank deficient, rank = 1,  tol =    3.1821e-15.
0112 %  > In lse at 205
0113 %  X =
0114 %      0.17097
0115 %      0.82903
0116 %            0
0117 %            0
0118 %
0119 % The pinv version will survive the singulatity
0120 %  Xp = lse(A,b,C,d,'pinv')
0121 %  Xp =
0122 %     0.085486
0123 %      0.41451
0124 %     0.085486
0125 %      0.41451
0126 %
0127 % Use of pinv will produce the minimum norm solution
0128 %  norm(X)
0129 %  ans =
0130 %      0.84647
0131 %
0132 %  norm(Xp)
0133 %  ans =
0134 %      0.59855
0135 %
0136 %
0137 % Methodology:
0138 %  Both alternative methods offered in lse are effectively subspace
0139 %  methods. That is, the equality constraints are used to reduce the
0140 %  problem to a lower dimensional problem. The '\' method uses a pivoted
0141 %  QR factorization to choose a set of variables to be eliminated. This
0142 %  chooses the best subset of variables for elimination, avoiding small
0143 %  "pivots" where possible, as well as resolving the case where the
0144 %  supplied constraints are rank deficient. Note that when the constraint
0145 %  system is rank deficient, this method will result in one or more of
0146 %  the unknowns to be set to zero. An at length description of the QR
0147 %  based method for linear least squares subject to linear equality
0148 %  constraints is found in:
0149 %
0150 %  http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=8553&objectType=FILE
0151 %
0152 %  The 'pinv' method uses a related approach, reducing the problem to
0153 %  a least squares solution in a subspace. This method is chosen to be
0154 %  consistent with pinv as used for unconstrained least squares problems.
0155 %  The reduction to a subspace does not require the selection of specific
0156 %  variables to be eliminated. Instead the reduction is a projection as
0157 %  defined by a singular value decomposition. When the constraint system
0158 %  is rank deficient, the svd allows for a minimum norm solution, much
0159 %  as is done with pinv. This may be preferable for some users.
0160 %
0161 %  Both methods allow the application of weights if supplied. As well,
0162 %  problems with multiple right hand sides (b) are solved in a fully
0163 %  vectorized fashion.
0164 %
0165 %
0166 % See also: slash, lsqlin, lsequal
0167 %
0168 %
0169 % Author: John D'Errico
0170 % E-mail: woodchips@rochester.rr.com
0171 % Release: 3.0
0172 % Release date: 1/31/07
0173 
0174 % check sizes
0175 [n,p] = size(A);
0176 [r,nrhs] = size(b);
0177 [m,ccols] = size(C);
0178 if n~=r
0179   error 'A and b are incompatible in size (wrong number of rows)'
0180 elseif ~isempty(C) && (p~=ccols)
0181   error 'A and C must have the same number of columns'
0182 elseif ~isempty(C) && issparse(C)
0183   error 'C may not be a sparse matrix'
0184 elseif ~isempty(C) && (m~=size(d,1))
0185   error 'C and d are incompatible in size (wrong number of rows)'
0186 elseif ~isempty(C) && (size(d,2)~=1)
0187   error 'd must have only one column'
0188 elseif isempty(C) && ~isempty(d)
0189   error 'C and d are inconsistent with each other (one was empty)'
0190 elseif ~isempty(C) && isempty(d)
0191   error 'C and d are inconsistent with each other (one was empty)'
0192 end
0193 
0194 % default solver is '\'
0195 if (nargin<5) || isempty(solverflag)
0196   solverflag = '\';
0197 elseif ~ischar(solverflag)
0198   error 'If supplied, solverflag must be character'
0199 else
0200   % solverflag was supplied. Make sure it is legal.
0201   valid = {'\', 'backslash', 'pinv'};
0202   ind = strmatch(solverflag,valid);
0203   if (length(ind)==1)
0204     solverflag = valid{ind};
0205   else
0206     error(['Invalid solverflag: ',solverflag])
0207   end
0208 end
0209 
0210 % default for weights = []
0211 if (nargin<6) || isempty(weights)
0212   weights = [];
0213 else
0214   weights = weights(:);
0215   if (length(weights)~=n) || any(weights<0)
0216     error 'weights should be empty or a non-negative vector of length n'
0217   elseif all(weights==0)
0218     error 'At least some of the weights must be non-zero'
0219   else
0220     % weights was supplied. scale it to have mean value = 1
0221     weights = weights./mean(weights);
0222     % also sqrt the weights for application as an
0223     % effective replication factor. remember that
0224     % least squares will minimize the sum of squares.
0225     weights = sqrt(weights);
0226   end
0227 end
0228 
0229 % tolerance used on the solve
0230 Ctol = 1.e-13;
0231 
0232 if (nargin<3) || isempty(C)
0233   % solve with \ or pinv as desired.
0234   switch solverflag
0235     case {'\' 'backslash'}
0236       % solve with or without weights
0237       if isempty(weights)
0238         x = A\b;
0239       else
0240         x = (repmat(weights,1,size(A,2)).*A)\ ...
0241           (repmat(weights,1,nrhs).*b);
0242       end
0243       
0244     case 'pinv'
0245       % solve with or without weights
0246       if isempty(weights)
0247         ptol = Ctol*norm(A,1);
0248         x = pinv(A,ptol)*b;
0249       else
0250         Aw = repmat(weights,1,size(A,2)).*A;
0251         ptol = Ctol*norm(Aw,1);
0252         x = pinv(Aw,ptol)*(repmat(weights,1,nrhs).*b);
0253       end
0254   end
0255   
0256   % no Constraints, so we are done here.
0257   return
0258 end
0259 
0260 % Which solver do we use?
0261 switch solverflag
0262   case {'\' 'backslash'}
0263     % allow a rank deficient equality constraint matrix
0264     % column pivoted qr to eliminate variables
0265     [Q,R,E]=qr(C,0);
0266     
0267     % get the numerical rank of R (and therefore C)
0268     if m == 1
0269       rdiag = R(1,1);
0270     else
0271       rdiag = abs(diag(R));
0272     end
0273     crank = sum((rdiag(1)*Ctol) <= rdiag);
0274     if crank >= p
0275       error 'Overly constrained problem.'
0276     end
0277     
0278     % check for consistency in the constraints in
0279     % the event of rank deficiency in the constraint
0280     % system
0281     if crank < m
0282       k = Q(:,(crank+1):end)'*d;
0283       if any(k > (Ctol*norm(d)));
0284         error 'The constraint system is deficient and numerically inconsistent'
0285       end
0286     end
0287     
0288     % only need the first crank columns of Q
0289     qpd = Q(:,1:crank)'*d;
0290     
0291     % which columns of A (variables) will we eliminate?
0292     j_subs = E(1:crank);
0293     % those that remain will be estimated
0294     j_est = E((crank+1):p);
0295     
0296     r1 = R(1:crank,1:crank);
0297     r2 = R(1:crank,(crank+1):p);
0298     
0299     A1 = A(:,j_subs);
0300     qpd = qpd(1:crank,:);
0301     
0302     % note that \ is still ok here, even if pinv
0303     % is used for the main regression.
0304     bmod = b-A1*(r1\repmat(qpd,1,nrhs));
0305     Amod = A(:,j_est)-A1*(r1\r2);
0306     
0307     % now solve the reduced problem, with or without weights
0308     if isempty(weights)
0309       x2 = Amod\bmod;
0310     else
0311       x2 = (repmat(weights,1,size(Amod,2)).*Amod)\ ...
0312         (repmat(weights,1,nrhs).*bmod);
0313     end
0314     
0315     % recover eliminated unknowns
0316     x1 = r1\(repmat(qpd,1,nrhs)-r2*x2);
0317     
0318     % stuff all estimated parameters into final vector
0319     x = zeros(p,nrhs);
0320     x(j_est,:) = x2;
0321     x(j_subs,:) = x1;
0322     
0323   case 'pinv'
0324     % allow a rank deficient equality constraint matrix
0325     Ctol = 1e-13;
0326     % use svd to deal with the variables
0327     [U,S,V] = svd(C,0);
0328     
0329     % get the numerical rank of S (and therefore C)
0330     if m == 1
0331       sdiag = S(1,1);
0332     else
0333       sdiag = diag(S);
0334     end
0335     crank = sum((sdiag(1)*Ctol) <= sdiag);
0336     if crank >= p
0337       error 'Overly constrained problem.'
0338     end
0339     
0340     % check for consistency in the constraints in
0341     % the event of rank deficiency in the constraint
0342     % system
0343     if crank < m
0344       k = U(:,(crank+1):end)'*d;
0345       if any(k > (Ctol*norm(d)));
0346         error 'The constraint system is deficient and numerically inconsistent'
0347       end
0348     end
0349     
0350     % only need the first crank columns of U, and the
0351     % effectively non-zero diagonal elements of S.
0352     sinv = diag(S);
0353     sinv = diag(1./sinv(1:crank));
0354     % we will use a transformation
0355     %  Z = V'*X = inv(S)*U'*d
0356     Z = sinv*U(:,1:crank)'*d;
0357     
0358     % Rather than explicitly dropping columns of A, we will
0359     % work in a reduced space as defined by the svd.
0360     Atrans = A*V;
0361     
0362     % thus, solve (A*V)*Z = b, subject to the constraints Z = supd
0363     % use pinv for the solution here.
0364     ptol = Ctol*norm(Atrans(:,(crank+1):end),1);
0365     if isempty(weights)
0366       Zsolve = pinv(Atrans(:,(crank+1):end),ptol)* ...
0367         (b - repmat(Atrans(:,1:crank)*Z(1:crank),1,nrhs));
0368     else
0369       w = spdiags(weights,0,n,n);
0370       Zsolve = pinv(w*Atrans(:,(crank+1):end),ptol)* ...
0371         (w*(b - repmat(Atrans(:,1:crank)*Z(1:crank),1,nrhs)));
0372     end
0373     
0374     % put it back together in the transformed state
0375     Z = [repmat(Z(1:crank),1,nrhs);Zsolve];
0376     
0377     % untransform back to the original variables
0378     x = V*Z;
0379     
0380 end

Generated on Mon 15-Sep-2014 13:31:28 by m2html © 2005