www.gusucode.com > mbcmodels 工具箱 matlab 源码程序 > mbcmodels/xreglsq.m
function [B,Xinv,Ri]= xreglsq(X,y,lam,doPrecond) % XREGLSQ linear least squares for large sparse problems % % [B,Xinv,Ri]= xreglsq(X,y,lam,doPrecond); % Note xreglsq will add a ridge matrix so an answer is always returned. % % Inputs % X Regression Matrix % y Output Data % lam Ridge matrix (can be 0) % doPreCond precondition regression matrix with column norms % Outputs % B lsq coefficients % Xinv pseudo inverse % Ri inverse of R from qr % Copyright 2000-2008 The MathWorks, Inc. and Ford Global Technologies, Inc. if nargin<3 lam=0; doRidge = false; else doRidge = ~all(lam(:)==0 ); end if nargin<4 doPrecond=1; end N= size(X,2); if isscalar(lam) L= lam*speye(N,N); elseif any(size(lam)==1) % diagonal matrix L= spdiags(lam,0,N,N); else L= lam; end if doPrecond % precondition X matrix [Xs,sd]= xregprecond(X); else Xs= X; sd=1; end if ~doRidge if issparse(Xs) % use q-less decomposition for sparse [qy,R] = qr(Xs,y,0); else % full qr decomposition [Q,R]= qr(Xs,0); qy= Q'*y; end rd= abs(diag(R)); tol= size(Xs,1)*eps*max(rd); if any(rd<tol) ; doRidge = true; % use ridge regression with a small lam to make sure problem is well-posed. if issparse(Xs) lam= sqrt(eps)*svds(Xs,1); else lam= sqrt(eps)*norm(Xs); end lam = max(lam,eps); L= lam*speye(N); end end if doRidge % augment qr for ridge XL= [Xs;L]; yL= [y; zeros(size(L,1),1)]; if issparse(XL) % use q-less decomposition for sparse [qy,R] = qr(XL,yL,0); else % full qr decomposition [Q,R]= qr(XL,0); qy= Q'*yL; end % rescale coefficients % B= sd*(R\(R'\(XL'*yL))); B= sd*(R\qy); else % q-less decomposition x = R\qy; if issparse(Xs) % iterative refinement r = y - Xs*x; e = R\(R'\(Xs'*r)); % scale solution x = (x + e); end B= sd*x; end if nargout>1 % pseudo inverse Xinv= sd*((Xs/R)/R')'; end if nargout>2 % inverse of R Ri= sd*inv(R); end