olsmatrix

 function f = olsmatrix(X,mode)

 <X> is samples x parameters
 <mode> (optional) is
   0 means normal operation
   1 means use inv instead of \ and omit unit-length normalization.
     the point of this mode is to reduce memory usage (i think).
   default: 0.

 what we want to do is to perform OLS regression using <X>
 and obtain the parameter estimates.  this is accomplished
 by inv(X'*X)*X'*y = f*y where y is the data (samples x cases).

 what this function does is to return <f> which has dimensions
 parameters x samples. 

 to ensure well-conditioning, we unit-length normalize each column
 of <X> before doing the inv.m operation.  also, we actually use
 left-slash (\), which apparently is equivalent to inv.m but faster
 and more accurate (see code for details).  if you pass <mode> as 1,
 we omit the normalization and use the inv method instead of the \ method.

 also, in the case that one or more regressors in <X> are all zero, then
 without special handling, then this case will result in warnings and
 NaN results.  what we do is to explicitly ensure that all-zero regressors
 are ignored and that there are all-zero rows for these regressors 
 in <f>.  this makes it such that the weights estimated for these 
 regressors are simply zero.

 history:
 2011/06/27 - explicitly ignore 0
 2011/03/29 - mode 1 now omits unit length normalization.
 2011/03/28 - add <mode> input.
 2010/08/05: now, we do not try to detect weird cases with low variance.
             instead, we just blindly unit-length normalize each column.

0001 function f = olsmatrix(X,mode)
0002 
0003 % function f = olsmatrix(X,mode)
0004 %
0005 % <X> is samples x parameters
0006 % <mode> (optional) is
0007 %   0 means normal operation
0008 %   1 means use inv instead of \ and omit unit-length normalization.
0009 %     the point of this mode is to reduce memory usage (i think).
0010 %   default: 0.
0011 %
0012 % what we want to do is to perform OLS regression using <X>
0013 % and obtain the parameter estimates.  this is accomplished
0014 % by inv(X'*X)*X'*y = f*y where y is the data (samples x cases).
0015 %
0016 % what this function does is to return <f> which has dimensions
0017 % parameters x samples.
0018 %
0019 % to ensure well-conditioning, we unit-length normalize each column
0020 % of <X> before doing the inv.m operation.  also, we actually use
0021 % left-slash (\), which apparently is equivalent to inv.m but faster
0022 % and more accurate (see code for details).  if you pass <mode> as 1,
0023 % we omit the normalization and use the inv method instead of the \ method.
0024 %
0025 % also, in the case that one or more regressors in <X> are all zero, then
0026 % without special handling, then this case will result in warnings and
0027 % NaN results.  what we do is to explicitly ensure that all-zero regressors
0028 % are ignored and that there are all-zero rows for these regressors
0029 % in <f>.  this makes it such that the weights estimated for these
0030 % regressors are simply zero.
0031 %
0032 % history:
0033 % 2011/06/27 - explicitly ignore 0
0034 % 2011/03/29 - mode 1 now omits unit length normalization.
0035 % 2011/03/28 - add <mode> input.
0036 % 2010/08/05: now, we do not try to detect weird cases with low variance.
0037 %             instead, we just blindly unit-length normalize each column.
0038 
0039 % input
0040 if ~exist('mode','var') || isempty(mode)
0041   mode = 0;
0042 end
0043 
0044 % deal with degenerate regressors
0045 good = ~all(X==0,1);
0046 
0047 % initialize result
0048 f = zeros(size(X,2),size(X,1));
0049 
0050 % do it
0051 switch mode
0052 case 0
0053   [X,len] = unitlength(X(:,good),1,[],0);
0054   temp = diag(1./len)*((X'*X)\X');  % hm, suggested by M-Lint
0055 case 1
0056   temp = inv(X(:,good)'*X(:,good))*X(:,good)';
0057 end
0058 
0059 % return
0060 if any(good)
0061   f(good,:) = temp;
0062 end