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